| blob | bdaac8c06d0c66c3bce54af20c8a604d23ac9c31 |
1 SUBROUTINE DSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
2 !
3 ! -- LAPACK routine (version 3.1) --
4 ! Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 ! November 2006
6 !
7 ! .. Scalar Arguments ..
8 CHARACTER COMPZ
9 INTEGER INFO, LDZ, N
10 ! ..
11 ! .. Array Arguments ..
12 DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * )
13 ! ..
14 !
15 ! Purpose
16 ! =======
17 !
18 ! DSTEQR computes all eigenvalues and, optionally, eigenvectors of a
19 ! symmetric tridiagonal matrix using the implicit QL or QR method.
20 ! The eigenvectors of a full or band symmetric matrix can also be found
21 ! if DSYTRD or DSPTRD or DSBTRD has been used to reduce this matrix to
22 ! tridiagonal form.
23 !
24 ! Arguments
25 ! =========
26 !
27 ! COMPZ (input) CHARACTER*1
28 ! = 'N': Compute eigenvalues only.
29 ! = 'V': Compute eigenvalues and eigenvectors of the original
30 ! symmetric matrix. On entry, Z must contain the
31 ! orthogonal matrix used to reduce the original matrix
32 ! to tridiagonal form.
33 ! = 'I': Compute eigenvalues and eigenvectors of the
34 ! tridiagonal matrix. Z is initialized to the identity
35 ! matrix.
36 !
37 ! N (input) INTEGER
38 ! The order of the matrix. N >= 0.
39 !
40 ! D (input/output) DOUBLE PRECISION array, dimension (N)
41 ! On entry, the diagonal elements of the tridiagonal matrix.
42 ! On exit, if INFO = 0, the eigenvalues in ascending order.
43 !
44 ! E (input/output) DOUBLE PRECISION array, dimension (N-1)
45 ! On entry, the (n-1) subdiagonal elements of the tridiagonal
46 ! matrix.
47 ! On exit, E has been destroyed.
48 !
49 ! Z (input/output) DOUBLE PRECISION array, dimension (LDZ, N)
50 ! On entry, if COMPZ = 'V', then Z contains the orthogonal
51 ! matrix used in the reduction to tridiagonal form.
52 ! On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
53 ! orthonormal eigenvectors of the original symmetric matrix,
54 ! and if COMPZ = 'I', Z contains the orthonormal eigenvectors
55 ! of the symmetric tridiagonal matrix.
56 ! If COMPZ = 'N', then Z is not referenced.
57 !
58 ! LDZ (input) INTEGER
59 ! The leading dimension of the array Z. LDZ >= 1, and if
60 ! eigenvectors are desired, then LDZ >= max(1,N).
61 !
62 ! WORK (workspace) DOUBLE PRECISION array, dimension (max(1,2*N-2))
63 ! If COMPZ = 'N', then WORK is not referenced.
64 !
65 ! INFO (output) INTEGER
66 ! = 0: successful exit
67 ! < 0: if INFO = -i, the i-th argument had an illegal value
68 ! > 0: the algorithm has failed to find all the eigenvalues in
69 ! a total of 30*N iterations; if INFO = i, then i
70 ! elements of E have not converged to zero; on exit, D
71 ! and E contain the elements of a symmetric tridiagonal
72 ! matrix which is orthogonally similar to the original
73 ! matrix.
74 !
75 ! =====================================================================
76 !
77 ! .. Parameters ..
78 DOUBLE PRECISION ZERO, ONE, TWO, THREE
79 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, &
80 THREE = 3.0D0 )
81 INTEGER MAXIT
82 PARAMETER ( MAXIT = 30 )
83 ! ..
84 ! .. Local Scalars ..
85 INTEGER I, ICOMPZ, II, ISCALE, J, JTOT, K, L, L1, LEND, &
86 LENDM1, LENDP1, LENDSV, LM1, LSV, M, MM, MM1, &
87 NM1, NMAXIT
88 DOUBLE PRECISION ANORM, B, C, EPS, EPS2, F, G, P, R, RT1, RT2, &
89 S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST
90 ! ..
91 ! .. External Functions ..
92 ! LOGICAL LSAME
93 ! DOUBLE PRECISION DLAMCH, DLANST, DLAPY2
94 ! EXTERNAL LSAME, DLAMCH, DLANST, DLAPY2
95 ! ..
96 ! .. External Subroutines ..
97 ! EXTERNAL DLAE2, DLAEV2, DLARTG, DLASCL, DLASET, DLASR, &
98 ! DLASRT, DSWAP, XERBLA
99 ! ..
100 ! .. Intrinsic Functions ..
101 INTRINSIC ABS, MAX, SIGN, SQRT
102 ! ..
103 ! .. Executable Statements ..
104 !
105 ! Test the input parameters.
106 !
107 INFO = 0
108 !
109 IF( LSAME( COMPZ, 'N' ) ) THEN
110 ICOMPZ = 0
111 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
112 ICOMPZ = 1
113 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
114 ! LOGICAL LSAME
115 ICOMPZ = 2
116 ELSE
117 ICOMPZ = -1
118 END IF
119 IF( ICOMPZ.LT.0 ) THEN
120 INFO = -1
121 ELSE IF( N.LT.0 ) THEN
122 INFO = -2
123 ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, &
124 N ) ) ) THEN
125 INFO = -6
126 END IF
127 IF( INFO.NE.0 ) THEN
128 CALL XERBLA( 'DSTEQR', -INFO )
129 RETURN
130 END IF
131 !
132 ! Quick return if possible
133 !
134 IF( N.EQ.0 ) &
135 RETURN
136 !
137 IF( N.EQ.1 ) THEN
138 IF( ICOMPZ.EQ.2 ) &
139 Z( 1, 1 ) = ONE
140 RETURN
141 END IF
142 !
143 ! Determine the unit roundoff and over/underflow thresholds.
144 !
145 EPS = DLAMCH( 'E' )
146 EPS2 = EPS**2
147 SAFMIN = DLAMCH( 'S' )
148 SAFMAX = ONE / SAFMIN
149 SSFMAX = SQRT( SAFMAX ) / THREE
150 SSFMIN = SQRT( SAFMIN ) / EPS2
151 !
152 ! Compute the eigenvalues and eigenvectors of the tridiagonal
153 ! matrix.
154 !
155 IF( ICOMPZ.EQ.2 ) &
156 CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
157 !
158 NMAXIT = N*MAXIT
159 JTOT = 0
160 !
161 ! Determine where the matrix splits and choose QL or QR iteration
162 ! for each block, according to whether top or bottom diagonal
163 ! element is smaller.
164 !
165 L1 = 1
166 NM1 = N - 1
167 !
168 10 CONTINUE
169 IF( L1.GT.N ) &
170 GO TO 160
171 IF( L1.GT.1 ) &
172 E( L1-1 ) = ZERO
173 IF( L1.LE.NM1 ) THEN
174 DO 20 M = L1, NM1
175 TST = ABS( E( M ) )
176 IF( TST.EQ.ZERO ) &
177 GO TO 30
178 IF( TST.LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+ &
179 1 ) ) ) )*EPS ) THEN
180 E( M ) = ZERO
181 GO TO 30
182 END IF
183 20 CONTINUE
184 END IF
185 M = N
186 !
187 30 CONTINUE
188 L = L1
189 LSV = L
190 LEND = M
191 LENDSV = LEND
192 L1 = M + 1
193 IF( LEND.EQ.L ) &
194 GO TO 10
195 !
196 ! Scale submatrix in rows and columns L to LEND
197 !
198 ANORM = DLANST( 'I', LEND-L+1, D( L ), E( L ) )
199 ISCALE = 0
200 IF( ANORM.EQ.ZERO ) &
201 GO TO 10
202 IF( ANORM.GT.SSFMAX ) THEN
203 ISCALE = 1
204 CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N, &
205 INFO )
206 CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N, &
207 INFO )
208 ELSE IF( ANORM.LT.SSFMIN ) THEN
209 ISCALE = 2
210 CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N, &
211 INFO )
212 CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N, &
213 INFO )
214 END IF
215 !
216 ! Choose between QL and QR iteration
217 !
218 IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
219 LEND = LSV
220 L = LENDSV
221 END IF
222 !
223 IF( LEND.GT.L ) THEN
224 !
225 ! QL Iteration
226 !
227 ! Look for small subdiagonal element.
228 !
229 40 CONTINUE
230 IF( L.NE.LEND ) THEN
231 LENDM1 = LEND - 1
232 DO 50 M = L, LENDM1
233 TST = ABS( E( M ) )**2
234 IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M+1 ) )+ &
235 SAFMIN )GO TO 60
236 50 CONTINUE
237 END IF
238 !
239 M = LEND
240 !
241 60 CONTINUE
242 IF( M.LT.LEND ) &
243 E( M ) = ZERO
244 P = D( L )
245 IF( M.EQ.L ) &
246 GO TO 80
247 !
248 ! If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
249 ! to compute its eigensystem.
250 !
251 IF( M.EQ.L+1 ) THEN
252 IF( ICOMPZ.GT.0 ) THEN
253 CALL DLAEV2( D( L ), E( L ), D( L+1 ), RT1, RT2, C, S )
254 WORK( L ) = C
255 WORK( N-1+L ) = S
256 CALL DLASR( 'R', 'V', 'B', N, 2, WORK( L ), &
257 WORK( N-1+L ), Z( 1, L ), LDZ )
258 ELSE
259 CALL DLAE2( D( L ), E( L ), D( L+1 ), RT1, RT2 )
260 END IF
261 D( L ) = RT1
262 D( L+1 ) = RT2
263 E( L ) = ZERO
264 L = L + 2
265 IF( L.LE.LEND ) &
266 GO TO 40
267 GO TO 140
268 END IF
269 !
270 IF( JTOT.EQ.NMAXIT ) &
271 GO TO 140
272 JTOT = JTOT + 1
273 !
274 ! Form shift.
275 !
276 G = ( D( L+1 )-P ) / ( TWO*E( L ) )
277 R = DLAPY2( G, ONE )
278 G = D( M ) - P + ( E( L ) / ( G+SIGN( R, G ) ) )
279 !
280 S = ONE
281 C = ONE
282 P = ZERO
283 !
284 ! Inner loop
285 !
286 MM1 = M - 1
287 DO 70 I = MM1, L, -1
288 F = S*E( I )
289 B = C*E( I )
290 CALL DLARTG( G, F, C, S, R )
291 IF( I.NE.M-1 ) &
292 E( I+1 ) = R
293 G = D( I+1 ) - P
294 R = ( D( I )-G )*S + TWO*C*B
295 P = S*R
296 D( I+1 ) = G + P
297 G = C*R - B
298 !
299 ! If eigenvectors are desired, then save rotations.
300 !
301 IF( ICOMPZ.GT.0 ) THEN
302 WORK( I ) = C
303 WORK( N-1+I ) = -S
304 END IF
305 !
306 70 CONTINUE
307 !
308 ! If eigenvectors are desired, then apply saved rotations.
309 !
310 IF( ICOMPZ.GT.0 ) THEN
311 MM = M - L + 1
312 CALL DLASR( 'R', 'V', 'B', N, MM, WORK( L ), WORK( N-1+L ), &
313 Z( 1, L ), LDZ )
314 END IF
315 !
316 D( L ) = D( L ) - P
317 E( L ) = G
318 GO TO 40
319 !
320 ! Eigenvalue found.
321 !
322 80 CONTINUE
323 D( L ) = P
324 !
325 L = L + 1
326 IF( L.LE.LEND ) &
327 GO TO 40
328 GO TO 140
329 !
330 ELSE
331 !
332 ! QR Iteration
333 !
334 ! Look for small superdiagonal element.
335 !
336 90 CONTINUE
337 IF( L.NE.LEND ) THEN
338 LENDP1 = LEND + 1
339 DO 100 M = L, LENDP1, -1
340 TST = ABS( E( M-1 ) )**2
341 IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M-1 ) )+ &
342 SAFMIN )GO TO 110
343 100 CONTINUE
344 END IF
345 !
346 M = LEND
347 !
348 110 CONTINUE
349 IF( M.GT.LEND ) &
350 E( M-1 ) = ZERO
351 P = D( L )
352 IF( M.EQ.L ) &
353 GO TO 130
354 !
355 ! If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
356 ! to compute its eigensystem.
357 !
358 IF( M.EQ.L-1 ) THEN
359 IF( ICOMPZ.GT.0 ) THEN
360 CALL DLAEV2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2, C, S )
361 WORK( M ) = C
362 WORK( N-1+M ) = S
363 CALL DLASR( 'R', 'V', 'F', N, 2, WORK( M ), &
364 WORK( N-1+M ), Z( 1, L-1 ), LDZ )
365 ELSE
366 CALL DLAE2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2 )
367 END IF
368 D( L-1 ) = RT1
369 D( L ) = RT2
370 E( L-1 ) = ZERO
371 L = L - 2
372 IF( L.GE.LEND ) &
373 GO TO 90
374 GO TO 140
375 END IF
376 !
377 IF( JTOT.EQ.NMAXIT ) &
378 GO TO 140
379 JTOT = JTOT + 1
380 !
381 ! Form shift.
382 !
383 G = ( D( L-1 )-P ) / ( TWO*E( L-1 ) )
384 R = DLAPY2( G, ONE )
385 G = D( M ) - P + ( E( L-1 ) / ( G+SIGN( R, G ) ) )
386 !
387 S = ONE
388 C = ONE
389 P = ZERO
390 !
391 ! Inner loop
392 !
393 LM1 = L - 1
394 DO 120 I = M, LM1
395 F = S*E( I )
396 B = C*E( I )
397 CALL DLARTG( G, F, C, S, R )
398 IF( I.NE.M ) &
399 E( I-1 ) = R
400 G = D( I ) - P
401 R = ( D( I+1 )-G )*S + TWO*C*B
402 P = S*R
403 D( I ) = G + P
404 G = C*R - B
405 !
406 ! If eigenvectors are desired, then save rotations.
407 !
408 IF( ICOMPZ.GT.0 ) THEN
409 WORK( I ) = C
410 WORK( N-1+I ) = S
411 END IF
412 !
413 120 CONTINUE
414 !
415 ! If eigenvectors are desired, then apply saved rotations.
416 !
417 IF( ICOMPZ.GT.0 ) THEN
418 MM = L - M + 1
419 CALL DLASR( 'R', 'V', 'F', N, MM, WORK( M ), WORK( N-1+M ), &
420 Z( 1, M ), LDZ )
421 END IF
422 !
423 D( L ) = D( L ) - P
424 E( LM1 ) = G
425 GO TO 90
426 !
427 ! Eigenvalue found.
428 !
429 130 CONTINUE
430 D( L ) = P
431 !
432 L = L - 1
433 IF( L.GE.LEND ) &
434 GO TO 90
435 GO TO 140
436 !
437 END IF
438 !
439 ! Undo scaling if necessary
440 !
441 140 CONTINUE
442 IF( ISCALE.EQ.1 ) THEN
443 CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1, &
444 D( LSV ), N, INFO )
445 CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV, 1, E( LSV ), &
446 N, INFO )
447 ELSE IF( ISCALE.EQ.2 ) THEN
448 CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1, &
449 D( LSV ), N, INFO )
450 CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV, 1, E( LSV ), &
451 N, INFO )
452 END IF
453 !
454 ! Check for no convergence to an eigenvalue after a total
455 ! of N*MAXIT iterations.
456 !
457 IF( JTOT.LT.NMAXIT ) &
458 GO TO 10
459 DO 150 I = 1, N - 1
460 IF( E( I ).NE.ZERO ) &
461 INFO = INFO + 1
462 150 CONTINUE
463 GO TO 190
464 !
465 ! Order eigenvalues and eigenvectors.
466 !
467 160 CONTINUE
468 IF( ICOMPZ.EQ.0 ) THEN
469 !
470 ! Use Quick Sort
471 !
472 CALL DLASRT( 'I', N, D, INFO )
473 !
474 ELSE
475 !
476 ! Use Selection Sort to minimize swaps of eigenvectors
477 !
478 DO 180 II = 2, N
479 I = II - 1
480 K = I
481 P = D( I )
482 DO 170 J = II, N
483 IF( D( J ).LT.P ) THEN
484 K = J
485 P = D( J )
486 END IF
487 170 CONTINUE
488 IF( K.NE.I ) THEN
489 D( K ) = D( I )
490 D( I ) = P
491 CALL DSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
492 END IF
493 180 CONTINUE
494 END IF
495 !
496 190 CONTINUE
497 RETURN
498 !
499 ! End of DSTEQR
500 !
501 END SUBROUTINE DSTEQR