LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dsterf.f
Go to the documentation of this file.
1*> \brief \b DSTERF
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DSTERF + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsterf.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsterf.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsterf.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DSTERF( N, D, E, INFO )
22*
23* .. Scalar Arguments ..
24* INTEGER INFO, N
25* ..
26* .. Array Arguments ..
27* DOUBLE PRECISION D( * ), E( * )
28* ..
29*
30*
31*> \par Purpose:
32* =============
33*>
34*> \verbatim
35*>
36*> DSTERF computes all eigenvalues of a symmetric tridiagonal matrix
37*> using the Pal-Walker-Kahan variant of the QL or QR algorithm.
38*> \endverbatim
39*
40* Arguments:
41* ==========
42*
43*> \param[in] N
44*> \verbatim
45*> N is INTEGER
46*> The order of the matrix. N >= 0.
47*> \endverbatim
48*>
49*> \param[in,out] D
50*> \verbatim
51*> D is DOUBLE PRECISION array, dimension (N)
52*> On entry, the n diagonal elements of the tridiagonal matrix.
53*> On exit, if INFO = 0, the eigenvalues in ascending order.
54*> \endverbatim
55*>
56*> \param[in,out] E
57*> \verbatim
58*> E is DOUBLE PRECISION array, dimension (N-1)
59*> On entry, the (n-1) subdiagonal elements of the tridiagonal
60*> matrix.
61*> On exit, E has been destroyed.
62*> \endverbatim
63*>
64*> \param[out] INFO
65*> \verbatim
66*> INFO is INTEGER
67*> = 0: successful exit
68*> < 0: if INFO = -i, the i-th argument had an illegal value
69*> > 0: the algorithm failed to find all of the eigenvalues in
70*> a total of 30*N iterations; if INFO = i, then i
71*> elements of E have not converged to zero.
72*> \endverbatim
73*
74* Authors:
75* ========
76*
77*> \author Univ. of Tennessee
78*> \author Univ. of California Berkeley
79*> \author Univ. of Colorado Denver
80*> \author NAG Ltd.
81*
82*> \ingroup sterf
83*
84* =====================================================================
85 SUBROUTINE dsterf( N, D, E, INFO )
86*
87* -- LAPACK computational routine --
88* -- LAPACK is a software package provided by Univ. of Tennessee, --
89* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
90*
91* .. Scalar Arguments ..
92 INTEGER INFO, N
93* ..
94* .. Array Arguments ..
95 DOUBLE PRECISION D( * ), E( * )
96* ..
97*
98* =====================================================================
99*
100* .. Parameters ..
101 DOUBLE PRECISION ZERO, ONE, TWO, THREE
102 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
103 $ three = 3.0d0 )
104 INTEGER MAXIT
105 parameter( maxit = 30 )
106* ..
107* .. Local Scalars ..
108 INTEGER I, ISCALE, JTOT, L, L1, LEND, LENDSV, LSV, M,
109 $ NMAXIT
110 DOUBLE PRECISION ALPHA, ANORM, BB, C, EPS, EPS2, GAMMA, OLDC,
111 $ OLDGAM, P, R, RT1, RT2, RTE, S, SAFMAX, SAFMIN,
112 $ SIGMA, SSFMAX, SSFMIN, RMAX
113* ..
114* .. External Functions ..
115 DOUBLE PRECISION DLAMCH, DLANST, DLAPY2
116 EXTERNAL dlamch, dlanst, dlapy2
117* ..
118* .. External Subroutines ..
119 EXTERNAL dlae2, dlascl, dlasrt, xerbla
120* ..
121* .. Intrinsic Functions ..
122 INTRINSIC abs, sign, sqrt
123* ..
124* .. Executable Statements ..
125*
126* Test the input parameters.
127*
128 info = 0
129*
130* Quick return if possible
131*
132 IF( n.LT.0 ) THEN
133 info = -1
134 CALL xerbla( 'DSTERF', -info )
135 RETURN
136 END IF
137 IF( n.LE.1 )
138 $ RETURN
139*
140* Determine the unit roundoff for this environment.
141*
142 eps = dlamch( 'E' )
143 eps2 = eps**2
144 safmin = dlamch( 'S' )
145 safmax = one / safmin
146 ssfmax = sqrt( safmax ) / three
147 ssfmin = sqrt( safmin ) / eps2
148 rmax = dlamch( 'O' )
149*
150* Compute the eigenvalues of the tridiagonal matrix.
151*
152 nmaxit = n*maxit
153 sigma = zero
154 jtot = 0
155*
156* Determine where the matrix splits and choose QL or QR iteration
157* for each block, according to whether top or bottom diagonal
158* element is smaller.
159*
160 l1 = 1
161*
162 10 CONTINUE
163 IF( l1.GT.n )
164 $ GO TO 170
165 IF( l1.GT.1 )
166 $ e( l1-1 ) = zero
167 DO 20 m = l1, n - 1
168 IF( abs( e( m ) ).LE.( sqrt( abs( d( m ) ) )*sqrt( abs( d( m+
169 $ 1 ) ) ) )*eps ) THEN
170 e( m ) = zero
171 GO TO 30
172 END IF
173 20 CONTINUE
174 m = n
175*
176 30 CONTINUE
177 l = l1
178 lsv = l
179 lend = m
180 lendsv = lend
181 l1 = m + 1
182 IF( lend.EQ.l )
183 $ GO TO 10
184*
185* Scale submatrix in rows and columns L to LEND
186*
187 anorm = dlanst( 'M', lend-l+1, d( l ), e( l ) )
188 iscale = 0
189 IF( anorm.EQ.zero )
190 $ GO TO 10
191 IF( (anorm.GT.ssfmax) ) THEN
192 iscale = 1
193 CALL dlascl( 'G', 0, 0, anorm, ssfmax, lend-l+1, 1, d( l ), n,
194 $ info )
195 CALL dlascl( 'G', 0, 0, anorm, ssfmax, lend-l, 1, e( l ), n,
196 $ info )
197 ELSE IF( anorm.LT.ssfmin ) THEN
198 iscale = 2
199 CALL dlascl( 'G', 0, 0, anorm, ssfmin, lend-l+1, 1, d( l ), n,
200 $ info )
201 CALL dlascl( 'G', 0, 0, anorm, ssfmin, lend-l, 1, e( l ), n,
202 $ info )
203 END IF
204*
205 DO 40 i = l, lend - 1
206 e( i ) = e( i )**2
207 40 CONTINUE
208*
209* Choose between QL and QR iteration
210*
211 IF( abs( d( lend ) ).LT.abs( d( l ) ) ) THEN
212 lend = lsv
213 l = lendsv
214 END IF
215*
216 IF( lend.GE.l ) THEN
217*
218* QL Iteration
219*
220* Look for small subdiagonal element.
221*
222 50 CONTINUE
223 IF( l.NE.lend ) THEN
224 DO 60 m = l, lend - 1
225 IF( abs( e( m ) ).LE.eps2*abs( d( m )*d( m+1 ) ) )
226 $ GO TO 70
227 60 CONTINUE
228 END IF
229 m = lend
230*
231 70 CONTINUE
232 IF( m.LT.lend )
233 $ e( m ) = zero
234 p = d( l )
235 IF( m.EQ.l )
236 $ GO TO 90
237*
238* If remaining matrix is 2 by 2, use DLAE2 to compute its
239* eigenvalues.
240*
241 IF( m.EQ.l+1 ) THEN
242 rte = sqrt( e( l ) )
243 CALL dlae2( d( l ), rte, d( l+1 ), rt1, rt2 )
244 d( l ) = rt1
245 d( l+1 ) = rt2
246 e( l ) = zero
247 l = l + 2
248 IF( l.LE.lend )
249 $ GO TO 50
250 GO TO 150
251 END IF
252*
253 IF( jtot.EQ.nmaxit )
254 $ GO TO 150
255 jtot = jtot + 1
256*
257* Form shift.
258*
259 rte = sqrt( e( l ) )
260 sigma = ( d( l+1 )-p ) / ( two*rte )
261 r = dlapy2( sigma, one )
262 sigma = p - ( rte / ( sigma+sign( r, sigma ) ) )
263*
264 c = one
265 s = zero
266 gamma = d( m ) - sigma
267 p = gamma*gamma
268*
269* Inner loop
270*
271 DO 80 i = m - 1, l, -1
272 bb = e( i )
273 r = p + bb
274 IF( i.NE.m-1 )
275 $ e( i+1 ) = s*r
276 oldc = c
277 c = p / r
278 s = bb / r
279 oldgam = gamma
280 alpha = d( i )
281 gamma = c*( alpha-sigma ) - s*oldgam
282 d( i+1 ) = oldgam + ( alpha-gamma )
283 IF( c.NE.zero ) THEN
284 p = ( gamma*gamma ) / c
285 ELSE
286 p = oldc*bb
287 END IF
288 80 CONTINUE
289*
290 e( l ) = s*p
291 d( l ) = sigma + gamma
292 GO TO 50
293*
294* Eigenvalue found.
295*
296 90 CONTINUE
297 d( l ) = p
298*
299 l = l + 1
300 IF( l.LE.lend )
301 $ GO TO 50
302 GO TO 150
303*
304 ELSE
305*
306* QR Iteration
307*
308* Look for small superdiagonal element.
309*
310 100 CONTINUE
311 DO 110 m = l, lend + 1, -1
312 IF( abs( e( m-1 ) ).LE.eps2*abs( d( m )*d( m-1 ) ) )
313 $ GO TO 120
314 110 CONTINUE
315 m = lend
316*
317 120 CONTINUE
318 IF( m.GT.lend )
319 $ e( m-1 ) = zero
320 p = d( l )
321 IF( m.EQ.l )
322 $ GO TO 140
323*
324* If remaining matrix is 2 by 2, use DLAE2 to compute its
325* eigenvalues.
326*
327 IF( m.EQ.l-1 ) THEN
328 rte = sqrt( e( l-1 ) )
329 CALL dlae2( d( l ), rte, d( l-1 ), rt1, rt2 )
330 d( l ) = rt1
331 d( l-1 ) = rt2
332 e( l-1 ) = zero
333 l = l - 2
334 IF( l.GE.lend )
335 $ GO TO 100
336 GO TO 150
337 END IF
338*
339 IF( jtot.EQ.nmaxit )
340 $ GO TO 150
341 jtot = jtot + 1
342*
343* Form shift.
344*
345 rte = sqrt( e( l-1 ) )
346 sigma = ( d( l-1 )-p ) / ( two*rte )
347 r = dlapy2( sigma, one )
348 sigma = p - ( rte / ( sigma+sign( r, sigma ) ) )
349*
350 c = one
351 s = zero
352 gamma = d( m ) - sigma
353 p = gamma*gamma
354*
355* Inner loop
356*
357 DO 130 i = m, l - 1
358 bb = e( i )
359 r = p + bb
360 IF( i.NE.m )
361 $ e( i-1 ) = s*r
362 oldc = c
363 c = p / r
364 s = bb / r
365 oldgam = gamma
366 alpha = d( i+1 )
367 gamma = c*( alpha-sigma ) - s*oldgam
368 d( i ) = oldgam + ( alpha-gamma )
369 IF( c.NE.zero ) THEN
370 p = ( gamma*gamma ) / c
371 ELSE
372 p = oldc*bb
373 END IF
374 130 CONTINUE
375*
376 e( l-1 ) = s*p
377 d( l ) = sigma + gamma
378 GO TO 100
379*
380* Eigenvalue found.
381*
382 140 CONTINUE
383 d( l ) = p
384*
385 l = l - 1
386 IF( l.GE.lend )
387 $ GO TO 100
388 GO TO 150
389*
390 END IF
391*
392* Undo scaling if necessary
393*
394 150 CONTINUE
395 IF( iscale.EQ.1 )
396 $ CALL dlascl( 'G', 0, 0, ssfmax, anorm, lendsv-lsv+1, 1,
397 $ d( lsv ), n, info )
398 IF( iscale.EQ.2 )
399 $ CALL dlascl( 'G', 0, 0, ssfmin, anorm, lendsv-lsv+1, 1,
400 $ d( lsv ), n, info )
401*
402* Check for no convergence to an eigenvalue after a total
403* of N*MAXIT iterations.
404*
405 IF( jtot.LT.nmaxit )
406 $ GO TO 10
407 DO 160 i = 1, n - 1
408 IF( e( i ).NE.zero )
409 $ info = info + 1
410 160 CONTINUE
411 GO TO 180
412*
413* Sort eigenvalues in increasing order.
414*
415 170 CONTINUE
416 CALL dlasrt( 'I', n, d, info )
417*
418 180 CONTINUE
419 RETURN
420*
421* End of DSTERF
422*
423 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dlae2
subroutine dlae2(a, b, c, rt1, rt2)
DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix.
Definition dlae2.f:102
dlascl
subroutine dlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition dlascl.f:143
dlasrt
subroutine dlasrt(id, n, d, info)
DLASRT sorts numbers in increasing or decreasing order.
Definition dlasrt.f:88
dsterf
subroutine dsterf(n, d, e, info)
DSTERF
Definition dsterf.f:86