|
2814
|
1 SUBROUTINE DSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO ) |
|
|
2 * |
|
7034
|
3 * -- LAPACK routine (version 3.1) -- |
|
|
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. |
|
|
5 * November 2006 |
|
2814
|
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 ICOMPZ = 2 |
|
|
115 ELSE |
|
|
116 ICOMPZ = -1 |
|
|
117 END IF |
|
|
118 IF( ICOMPZ.LT.0 ) THEN |
|
|
119 INFO = -1 |
|
|
120 ELSE IF( N.LT.0 ) THEN |
|
|
121 INFO = -2 |
|
|
122 ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, |
|
|
123 $ N ) ) ) THEN |
|
|
124 INFO = -6 |
|
|
125 END IF |
|
|
126 IF( INFO.NE.0 ) THEN |
|
|
127 CALL XERBLA( 'DSTEQR', -INFO ) |
|
|
128 RETURN |
|
|
129 END IF |
|
|
130 * |
|
|
131 * Quick return if possible |
|
|
132 * |
|
|
133 IF( N.EQ.0 ) |
|
|
134 $ RETURN |
|
|
135 * |
|
|
136 IF( N.EQ.1 ) THEN |
|
|
137 IF( ICOMPZ.EQ.2 ) |
|
|
138 $ Z( 1, 1 ) = ONE |
|
|
139 RETURN |
|
|
140 END IF |
|
|
141 * |
|
|
142 * Determine the unit roundoff and over/underflow thresholds. |
|
|
143 * |
|
|
144 EPS = DLAMCH( 'E' ) |
|
|
145 EPS2 = EPS**2 |
|
|
146 SAFMIN = DLAMCH( 'S' ) |
|
|
147 SAFMAX = ONE / SAFMIN |
|
|
148 SSFMAX = SQRT( SAFMAX ) / THREE |
|
|
149 SSFMIN = SQRT( SAFMIN ) / EPS2 |
|
|
150 * |
|
|
151 * Compute the eigenvalues and eigenvectors of the tridiagonal |
|
|
152 * matrix. |
|
|
153 * |
|
|
154 IF( ICOMPZ.EQ.2 ) |
|
|
155 $ CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ ) |
|
|
156 * |
|
|
157 NMAXIT = N*MAXIT |
|
|
158 JTOT = 0 |
|
|
159 * |
|
|
160 * Determine where the matrix splits and choose QL or QR iteration |
|
|
161 * for each block, according to whether top or bottom diagonal |
|
|
162 * element is smaller. |
|
|
163 * |
|
|
164 L1 = 1 |
|
|
165 NM1 = N - 1 |
|
|
166 * |
|
|
167 10 CONTINUE |
|
|
168 IF( L1.GT.N ) |
|
|
169 $ GO TO 160 |
|
|
170 IF( L1.GT.1 ) |
|
|
171 $ E( L1-1 ) = ZERO |
|
|
172 IF( L1.LE.NM1 ) THEN |
|
|
173 DO 20 M = L1, NM1 |
|
|
174 TST = ABS( E( M ) ) |
|
|
175 IF( TST.EQ.ZERO ) |
|
|
176 $ GO TO 30 |
|
|
177 IF( TST.LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+ |
|
|
178 $ 1 ) ) ) )*EPS ) THEN |
|
|
179 E( M ) = ZERO |
|
|
180 GO TO 30 |
|
|
181 END IF |
|
|
182 20 CONTINUE |
|
|
183 END IF |
|
|
184 M = N |
|
|
185 * |
|
|
186 30 CONTINUE |
|
|
187 L = L1 |
|
|
188 LSV = L |
|
|
189 LEND = M |
|
|
190 LENDSV = LEND |
|
|
191 L1 = M + 1 |
|
|
192 IF( LEND.EQ.L ) |
|
|
193 $ GO TO 10 |
|
|
194 * |
|
|
195 * Scale submatrix in rows and columns L to LEND |
|
|
196 * |
|
|
197 ANORM = DLANST( 'I', LEND-L+1, D( L ), E( L ) ) |
|
|
198 ISCALE = 0 |
|
|
199 IF( ANORM.EQ.ZERO ) |
|
|
200 $ GO TO 10 |
|
|
201 IF( ANORM.GT.SSFMAX ) THEN |
|
|
202 ISCALE = 1 |
|
|
203 CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N, |
|
|
204 $ INFO ) |
|
|
205 CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N, |
|
|
206 $ INFO ) |
|
|
207 ELSE IF( ANORM.LT.SSFMIN ) THEN |
|
|
208 ISCALE = 2 |
|
|
209 CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N, |
|
|
210 $ INFO ) |
|
|
211 CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N, |
|
|
212 $ INFO ) |
|
|
213 END IF |
|
|
214 * |
|
|
215 * Choose between QL and QR iteration |
|
|
216 * |
|
|
217 IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN |
|
|
218 LEND = LSV |
|
|
219 L = LENDSV |
|
|
220 END IF |
|
|
221 * |
|
|
222 IF( LEND.GT.L ) THEN |
|
|
223 * |
|
|
224 * QL Iteration |
|
|
225 * |
|
|
226 * Look for small subdiagonal element. |
|
|
227 * |
|
|
228 40 CONTINUE |
|
|
229 IF( L.NE.LEND ) THEN |
|
|
230 LENDM1 = LEND - 1 |
|
|
231 DO 50 M = L, LENDM1 |
|
|
232 TST = ABS( E( M ) )**2 |
|
|
233 IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M+1 ) )+ |
|
|
234 $ SAFMIN )GO TO 60 |
|
|
235 50 CONTINUE |
|
|
236 END IF |
|
|
237 * |
|
|
238 M = LEND |
|
|
239 * |
|
|
240 60 CONTINUE |
|
|
241 IF( M.LT.LEND ) |
|
|
242 $ E( M ) = ZERO |
|
|
243 P = D( L ) |
|
|
244 IF( M.EQ.L ) |
|
|
245 $ GO TO 80 |
|
|
246 * |
|
|
247 * If remaining matrix is 2-by-2, use DLAE2 or SLAEV2 |
|
|
248 * to compute its eigensystem. |
|
|
249 * |
|
|
250 IF( M.EQ.L+1 ) THEN |
|
|
251 IF( ICOMPZ.GT.0 ) THEN |
|
|
252 CALL DLAEV2( D( L ), E( L ), D( L+1 ), RT1, RT2, C, S ) |
|
|
253 WORK( L ) = C |
|
|
254 WORK( N-1+L ) = S |
|
|
255 CALL DLASR( 'R', 'V', 'B', N, 2, WORK( L ), |
|
|
256 $ WORK( N-1+L ), Z( 1, L ), LDZ ) |
|
|
257 ELSE |
|
|
258 CALL DLAE2( D( L ), E( L ), D( L+1 ), RT1, RT2 ) |
|
|
259 END IF |
|
|
260 D( L ) = RT1 |
|
|
261 D( L+1 ) = RT2 |
|
|
262 E( L ) = ZERO |
|
|
263 L = L + 2 |
|
|
264 IF( L.LE.LEND ) |
|
|
265 $ GO TO 40 |
|
|
266 GO TO 140 |
|
|
267 END IF |
|
|
268 * |
|
|
269 IF( JTOT.EQ.NMAXIT ) |
|
|
270 $ GO TO 140 |
|
|
271 JTOT = JTOT + 1 |
|
|
272 * |
|
|
273 * Form shift. |
|
|
274 * |
|
|
275 G = ( D( L+1 )-P ) / ( TWO*E( L ) ) |
|
|
276 R = DLAPY2( G, ONE ) |
|
|
277 G = D( M ) - P + ( E( L ) / ( G+SIGN( R, G ) ) ) |
|
|
278 * |
|
|
279 S = ONE |
|
|
280 C = ONE |
|
|
281 P = ZERO |
|
|
282 * |
|
|
283 * Inner loop |
|
|
284 * |
|
|
285 MM1 = M - 1 |
|
|
286 DO 70 I = MM1, L, -1 |
|
|
287 F = S*E( I ) |
|
|
288 B = C*E( I ) |
|
|
289 CALL DLARTG( G, F, C, S, R ) |
|
|
290 IF( I.NE.M-1 ) |
|
|
291 $ E( I+1 ) = R |
|
|
292 G = D( I+1 ) - P |
|
|
293 R = ( D( I )-G )*S + TWO*C*B |
|
|
294 P = S*R |
|
|
295 D( I+1 ) = G + P |
|
|
296 G = C*R - B |
|
|
297 * |
|
|
298 * If eigenvectors are desired, then save rotations. |
|
|
299 * |
|
|
300 IF( ICOMPZ.GT.0 ) THEN |
|
|
301 WORK( I ) = C |
|
|
302 WORK( N-1+I ) = -S |
|
|
303 END IF |
|
|
304 * |
|
|
305 70 CONTINUE |
|
|
306 * |
|
|
307 * If eigenvectors are desired, then apply saved rotations. |
|
|
308 * |
|
|
309 IF( ICOMPZ.GT.0 ) THEN |
|
|
310 MM = M - L + 1 |
|
|
311 CALL DLASR( 'R', 'V', 'B', N, MM, WORK( L ), WORK( N-1+L ), |
|
|
312 $ Z( 1, L ), LDZ ) |
|
|
313 END IF |
|
|
314 * |
|
|
315 D( L ) = D( L ) - P |
|
|
316 E( L ) = G |
|
|
317 GO TO 40 |
|
|
318 * |
|
|
319 * Eigenvalue found. |
|
|
320 * |
|
|
321 80 CONTINUE |
|
|
322 D( L ) = P |
|
|
323 * |
|
|
324 L = L + 1 |
|
|
325 IF( L.LE.LEND ) |
|
|
326 $ GO TO 40 |
|
|
327 GO TO 140 |
|
|
328 * |
|
|
329 ELSE |
|
|
330 * |
|
|
331 * QR Iteration |
|
|
332 * |
|
|
333 * Look for small superdiagonal element. |
|
|
334 * |
|
|
335 90 CONTINUE |
|
|
336 IF( L.NE.LEND ) THEN |
|
|
337 LENDP1 = LEND + 1 |
|
|
338 DO 100 M = L, LENDP1, -1 |
|
|
339 TST = ABS( E( M-1 ) )**2 |
|
|
340 IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M-1 ) )+ |
|
|
341 $ SAFMIN )GO TO 110 |
|
|
342 100 CONTINUE |
|
|
343 END IF |
|
|
344 * |
|
|
345 M = LEND |
|
|
346 * |
|
|
347 110 CONTINUE |
|
|
348 IF( M.GT.LEND ) |
|
|
349 $ E( M-1 ) = ZERO |
|
|
350 P = D( L ) |
|
|
351 IF( M.EQ.L ) |
|
|
352 $ GO TO 130 |
|
|
353 * |
|
|
354 * If remaining matrix is 2-by-2, use DLAE2 or SLAEV2 |
|
|
355 * to compute its eigensystem. |
|
|
356 * |
|
|
357 IF( M.EQ.L-1 ) THEN |
|
|
358 IF( ICOMPZ.GT.0 ) THEN |
|
|
359 CALL DLAEV2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2, C, S ) |
|
|
360 WORK( M ) = C |
|
|
361 WORK( N-1+M ) = S |
|
|
362 CALL DLASR( 'R', 'V', 'F', N, 2, WORK( M ), |
|
|
363 $ WORK( N-1+M ), Z( 1, L-1 ), LDZ ) |
|
|
364 ELSE |
|
|
365 CALL DLAE2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2 ) |
|
|
366 END IF |
|
|
367 D( L-1 ) = RT1 |
|
|
368 D( L ) = RT2 |
|
|
369 E( L-1 ) = ZERO |
|
|
370 L = L - 2 |
|
|
371 IF( L.GE.LEND ) |
|
|
372 $ GO TO 90 |
|
|
373 GO TO 140 |
|
|
374 END IF |
|
|
375 * |
|
|
376 IF( JTOT.EQ.NMAXIT ) |
|
|
377 $ GO TO 140 |
|
|
378 JTOT = JTOT + 1 |
|
|
379 * |
|
|
380 * Form shift. |
|
|
381 * |
|
|
382 G = ( D( L-1 )-P ) / ( TWO*E( L-1 ) ) |
|
|
383 R = DLAPY2( G, ONE ) |
|
|
384 G = D( M ) - P + ( E( L-1 ) / ( G+SIGN( R, G ) ) ) |
|
|
385 * |
|
|
386 S = ONE |
|
|
387 C = ONE |
|
|
388 P = ZERO |
|
|
389 * |
|
|
390 * Inner loop |
|
|
391 * |
|
|
392 LM1 = L - 1 |
|
|
393 DO 120 I = M, LM1 |
|
|
394 F = S*E( I ) |
|
|
395 B = C*E( I ) |
|
|
396 CALL DLARTG( G, F, C, S, R ) |
|
|
397 IF( I.NE.M ) |
|
|
398 $ E( I-1 ) = R |
|
|
399 G = D( I ) - P |
|
|
400 R = ( D( I+1 )-G )*S + TWO*C*B |
|
|
401 P = S*R |
|
|
402 D( I ) = G + P |
|
|
403 G = C*R - B |
|
|
404 * |
|
|
405 * If eigenvectors are desired, then save rotations. |
|
|
406 * |
|
|
407 IF( ICOMPZ.GT.0 ) THEN |
|
|
408 WORK( I ) = C |
|
|
409 WORK( N-1+I ) = S |
|
|
410 END IF |
|
|
411 * |
|
|
412 120 CONTINUE |
|
|
413 * |
|
|
414 * If eigenvectors are desired, then apply saved rotations. |
|
|
415 * |
|
|
416 IF( ICOMPZ.GT.0 ) THEN |
|
|
417 MM = L - M + 1 |
|
|
418 CALL DLASR( 'R', 'V', 'F', N, MM, WORK( M ), WORK( N-1+M ), |
|
|
419 $ Z( 1, M ), LDZ ) |
|
|
420 END IF |
|
|
421 * |
|
|
422 D( L ) = D( L ) - P |
|
|
423 E( LM1 ) = G |
|
|
424 GO TO 90 |
|
|
425 * |
|
|
426 * Eigenvalue found. |
|
|
427 * |
|
|
428 130 CONTINUE |
|
|
429 D( L ) = P |
|
|
430 * |
|
|
431 L = L - 1 |
|
|
432 IF( L.GE.LEND ) |
|
|
433 $ GO TO 90 |
|
|
434 GO TO 140 |
|
|
435 * |
|
|
436 END IF |
|
|
437 * |
|
|
438 * Undo scaling if necessary |
|
|
439 * |
|
|
440 140 CONTINUE |
|
|
441 IF( ISCALE.EQ.1 ) THEN |
|
|
442 CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1, |
|
|
443 $ D( LSV ), N, INFO ) |
|
|
444 CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV, 1, E( LSV ), |
|
|
445 $ N, INFO ) |
|
|
446 ELSE IF( ISCALE.EQ.2 ) THEN |
|
|
447 CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1, |
|
|
448 $ D( LSV ), N, INFO ) |
|
|
449 CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV, 1, E( LSV ), |
|
|
450 $ N, INFO ) |
|
|
451 END IF |
|
|
452 * |
|
|
453 * Check for no convergence to an eigenvalue after a total |
|
|
454 * of N*MAXIT iterations. |
|
|
455 * |
|
|
456 IF( JTOT.LT.NMAXIT ) |
|
|
457 $ GO TO 10 |
|
|
458 DO 150 I = 1, N - 1 |
|
|
459 IF( E( I ).NE.ZERO ) |
|
|
460 $ INFO = INFO + 1 |
|
|
461 150 CONTINUE |
|
|
462 GO TO 190 |
|
|
463 * |
|
|
464 * Order eigenvalues and eigenvectors. |
|
|
465 * |
|
|
466 160 CONTINUE |
|
|
467 IF( ICOMPZ.EQ.0 ) THEN |
|
|
468 * |
|
|
469 * Use Quick Sort |
|
|
470 * |
|
|
471 CALL DLASRT( 'I', N, D, INFO ) |
|
|
472 * |
|
|
473 ELSE |
|
|
474 * |
|
|
475 * Use Selection Sort to minimize swaps of eigenvectors |
|
|
476 * |
|
|
477 DO 180 II = 2, N |
|
|
478 I = II - 1 |
|
|
479 K = I |
|
|
480 P = D( I ) |
|
|
481 DO 170 J = II, N |
|
|
482 IF( D( J ).LT.P ) THEN |
|
|
483 K = J |
|
|
484 P = D( J ) |
|
|
485 END IF |
|
|
486 170 CONTINUE |
|
|
487 IF( K.NE.I ) THEN |
|
|
488 D( K ) = D( I ) |
|
|
489 D( I ) = P |
|
|
490 CALL DSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 ) |
|
|
491 END IF |
|
|
492 180 CONTINUE |
|
|
493 END IF |
|
|
494 * |
|
|
495 190 CONTINUE |
|
|
496 RETURN |
|
|
497 * |
|
|
498 * End of DSTEQR |
|
|
499 * |
|
|
500 END |
