LAPACK 3.12.1
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*> Download DSTERF + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsterf.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsterf.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsterf.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE DSTERF( N, D, E, INFO )
20*
21* .. Scalar Arguments ..
22* INTEGER INFO, N
23* ..
24* .. Array Arguments ..
25* DOUBLE PRECISION D( * ), E( * )
26* ..
27*
28*
29*> \par Purpose:
30* =============
31*>
32*> \verbatim
33*>
34*> DSTERF computes all eigenvalues of a symmetric tridiagonal matrix
35*> using the Pal-Walker-Kahan variant of the QL or QR algorithm.
36*> \endverbatim
37*
38* Arguments:
39* ==========
40*
41*> \param[in] N
42*> \verbatim
43*> N is INTEGER
44*> The order of the matrix. N >= 0.
45*> \endverbatim
46*>
47*> \param[in,out] D
48*> \verbatim
49*> D is DOUBLE PRECISION array, dimension (N)
50*> On entry, the n diagonal elements of the tridiagonal matrix.
51*> On exit, if INFO = 0, the eigenvalues in ascending order.
52*> \endverbatim
53*>
54*> \param[in,out] E
55*> \verbatim
56*> E is DOUBLE PRECISION array, dimension (N-1)
57*> On entry, the (n-1) subdiagonal elements of the tridiagonal
58*> matrix.
59*> On exit, E has been destroyed.
60*> \endverbatim
61*>
62*> \param[out] INFO
63*> \verbatim
64*> INFO is INTEGER
65*> = 0: successful exit
66*> < 0: if INFO = -i, the i-th argument had an illegal value
67*> > 0: the algorithm failed to find all of the eigenvalues in
68*> a total of 30*N iterations; if INFO = i, then i
69*> elements of E have not converged to zero.
70*> \endverbatim
71*
72* Authors:
73* ========
74*
75*> \author Univ. of Tennessee
76*> \author Univ. of California Berkeley
77*> \author Univ. of Colorado Denver
78*> \author NAG Ltd.
79*
80*> \ingroup sterf
81*
82* =====================================================================
83 SUBROUTINE dsterf( N, D, E, INFO )
84*
85* -- LAPACK computational routine --
86* -- LAPACK is a software package provided by Univ. of Tennessee, --
87* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
88*
89* .. Scalar Arguments ..
90 INTEGER INFO, N
91* ..
92* .. Array Arguments ..
93 DOUBLE PRECISION D( * ), E( * )
94* ..
95*
96* =====================================================================
97*
98* .. Parameters ..
99 DOUBLE PRECISION ZERO, ONE, TWO, THREE
100 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
101 $ three = 3.0d0 )
102 INTEGER MAXIT
103 parameter( maxit = 30 )
104* ..
105* .. Local Scalars ..
106 INTEGER I, ISCALE, JTOT, L, L1, LEND, LENDSV, LSV, M,
107 $ NMAXIT
108 DOUBLE PRECISION ALPHA, ANORM, BB, C, EPS, EPS2, GAMMA, OLDC,
109 $ OLDGAM, P, R, RT1, RT2, RTE, S, SAFMAX, SAFMIN,
110 $ SIGMA, SSFMAX, SSFMIN, RMAX
111* ..
112* .. External Functions ..
113 DOUBLE PRECISION DLAMCH, DLANST, DLAPY2
114 EXTERNAL dlamch, dlanst, dlapy2
115* ..
116* .. External Subroutines ..
117 EXTERNAL dlae2, dlascl, dlasrt, xerbla
118* ..
119* .. Intrinsic Functions ..
120 INTRINSIC abs, sign, sqrt
121* ..
122* .. Executable Statements ..
123*
124* Test the input parameters.
125*
126 info = 0
127*
128* Quick return if possible
129*
130 IF( n.LT.0 ) THEN
131 info = -1
132 CALL xerbla( 'DSTERF', -info )
133 RETURN
134 END IF
135 IF( n.LE.1 )
136 $ RETURN
137*
138* Determine the unit roundoff for this environment.
139*
140 eps = dlamch( 'E' )
141 eps2 = eps**2
142 safmin = dlamch( 'S' )
143 safmax = one / safmin
144 ssfmax = sqrt( safmax ) / three
145 ssfmin = sqrt( safmin ) / eps2
146 rmax = dlamch( 'O' )
147*
148* Compute the eigenvalues of the tridiagonal matrix.
149*
150 nmaxit = n*maxit
151 sigma = zero
152 jtot = 0
153*
154* Determine where the matrix splits and choose QL or QR iteration
155* for each block, according to whether top or bottom diagonal
156* element is smaller.
157*
158 l1 = 1
159*
160 10 CONTINUE
161 IF( l1.GT.n )
162 $ GO TO 170
163 IF( l1.GT.1 )
164 $ e( l1-1 ) = zero
165 DO 20 m = l1, n - 1
166 IF( abs( e( m ) ).LE.( sqrt( abs( d( m ) ) )*sqrt( abs( d( m+
167 $ 1 ) ) ) )*eps ) THEN
168 e( m ) = zero
169 GO TO 30
170 END IF
171 20 CONTINUE
172 m = n
173*
174 30 CONTINUE
175 l = l1
176 lsv = l
177 lend = m
178 lendsv = lend
179 l1 = m + 1
180 IF( lend.EQ.l )
181 $ GO TO 10
182*
183* Scale submatrix in rows and columns L to LEND
184*
185 anorm = dlanst( 'M', lend-l+1, d( l ), e( l ) )
186 iscale = 0
187 IF( anorm.EQ.zero )
188 $ GO TO 10
189 IF( (anorm.GT.ssfmax) ) THEN
190 iscale = 1
191 CALL dlascl( 'G', 0, 0, anorm, ssfmax, lend-l+1, 1, d( l ),
192 $ n,
193 $ info )
194 CALL dlascl( 'G', 0, 0, anorm, ssfmax, lend-l, 1, e( l ), n,
195 $ info )
196 ELSE IF( anorm.LT.ssfmin ) THEN
197 iscale = 2
198 CALL dlascl( 'G', 0, 0, anorm, ssfmin, lend-l+1, 1, d( l ),
199 $ 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:100
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:142
dlasrt
subroutine dlasrt(id, n, d, info)
DLASRT sorts numbers in increasing or decreasing order.
Definition dlasrt.f:86
dsterf
subroutine dsterf(n, d, e, info)
DSTERF
Definition dsterf.f:84