1 SUBROUTINE DSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 CHARACTER COMPZ
10 INTEGER INFO, LDZ, N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DSTEQR computes all eigenvalues and, optionally, eigenvectors of a
20 * symmetric tridiagonal matrix using the implicit QL or QR method.
21 * The eigenvectors of a full or band symmetric matrix can also be found
22 * if DSYTRD or DSPTRD or DSBTRD has been used to reduce this matrix to
23 * tridiagonal form.
24 *
25 * Arguments
26 * =========
27 *
28 * COMPZ (input) CHARACTER*1
29 * = 'N': Compute eigenvalues only.
30 * = 'V': Compute eigenvalues and eigenvectors of the original
31 * symmetric matrix. On entry, Z must contain the
32 * orthogonal matrix used to reduce the original matrix
33 * to tridiagonal form.
34 * = 'I': Compute eigenvalues and eigenvectors of the
35 * tridiagonal matrix. Z is initialized to the identity
36 * matrix.
37 *
38 * N (input) INTEGER
39 * The order of the matrix. N >= 0.
40 *
41 * D (input/output) DOUBLE PRECISION array, dimension (N)
42 * On entry, the diagonal elements of the tridiagonal matrix.
43 * On exit, if INFO = 0, the eigenvalues in ascending order.
44 *
45 * E (input/output) DOUBLE PRECISION array, dimension (N-1)
46 * On entry, the (n-1) subdiagonal elements of the tridiagonal
47 * matrix.
48 * On exit, E has been destroyed.
49 *
50 * Z (input/output) DOUBLE PRECISION array, dimension (LDZ, N)
51 * On entry, if COMPZ = 'V', then Z contains the orthogonal
52 * matrix used in the reduction to tridiagonal form.
53 * On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
54 * orthonormal eigenvectors of the original symmetric matrix,
55 * and if COMPZ = 'I', Z contains the orthonormal eigenvectors
56 * of the symmetric tridiagonal matrix.
57 * If COMPZ = 'N', then Z is not referenced.
58 *
59 * LDZ (input) INTEGER
60 * The leading dimension of the array Z. LDZ >= 1, and if
61 * eigenvectors are desired, then LDZ >= max(1,N).
62 *
63 * WORK (workspace) DOUBLE PRECISION array, dimension (max(1,2*N-2))
64 * If COMPZ = 'N', then WORK is not referenced.
65 *
66 * INFO (output) INTEGER
67 * = 0: successful exit
68 * < 0: if INFO = -i, the i-th argument had an illegal value
69 * > 0: the algorithm has failed to find all the eigenvalues in
70 * a total of 30*N iterations; if INFO = i, then i
71 * elements of E have not converged to zero; on exit, D
72 * and E contain the elements of a symmetric tridiagonal
73 * matrix which is orthogonally similar to the original
74 * matrix.
75 *
76 * =====================================================================
77 *
78 * .. Parameters ..
79 DOUBLE PRECISION ZERO, ONE, TWO, THREE
80 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
81 $ THREE = 3.0D0 )
82 INTEGER MAXIT
83 PARAMETER ( MAXIT = 30 )
84 * ..
85 * .. Local Scalars ..
86 INTEGER I, ICOMPZ, II, ISCALE, J, JTOT, K, L, L1, LEND,
87 $ LENDM1, LENDP1, LENDSV, LM1, LSV, M, MM, MM1,
88 $ NM1, NMAXIT
89 DOUBLE PRECISION ANORM, B, C, EPS, EPS2, F, G, P, R, RT1, RT2,
90 $ S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST
91 * ..
92 * .. External Functions ..
93 LOGICAL LSAME
94 DOUBLE PRECISION DLAMCH, DLANST, DLAPY2
95 EXTERNAL LSAME, DLAMCH, DLANST, DLAPY2
96 * ..
97 * .. External Subroutines ..
98 EXTERNAL DLAE2, DLAEV2, DLARTG, DLASCL, DLASET, DLASR,
99 $ DLASRT, DSWAP, XERBLA
100 * ..
101 * .. Intrinsic Functions ..
102 INTRINSIC ABS, MAX, SIGN, SQRT
103 * ..
104 * .. Executable Statements ..
105 *
106 * Test the input parameters.
107 *
108 INFO = 0
109 *
110 IF( LSAME( COMPZ, 'N' ) ) THEN
111 ICOMPZ = 0
112 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
113 ICOMPZ = 1
114 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
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( 'M', 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
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 CHARACTER COMPZ
10 INTEGER INFO, LDZ, N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DSTEQR computes all eigenvalues and, optionally, eigenvectors of a
20 * symmetric tridiagonal matrix using the implicit QL or QR method.
21 * The eigenvectors of a full or band symmetric matrix can also be found
22 * if DSYTRD or DSPTRD or DSBTRD has been used to reduce this matrix to
23 * tridiagonal form.
24 *
25 * Arguments
26 * =========
27 *
28 * COMPZ (input) CHARACTER*1
29 * = 'N': Compute eigenvalues only.
30 * = 'V': Compute eigenvalues and eigenvectors of the original
31 * symmetric matrix. On entry, Z must contain the
32 * orthogonal matrix used to reduce the original matrix
33 * to tridiagonal form.
34 * = 'I': Compute eigenvalues and eigenvectors of the
35 * tridiagonal matrix. Z is initialized to the identity
36 * matrix.
37 *
38 * N (input) INTEGER
39 * The order of the matrix. N >= 0.
40 *
41 * D (input/output) DOUBLE PRECISION array, dimension (N)
42 * On entry, the diagonal elements of the tridiagonal matrix.
43 * On exit, if INFO = 0, the eigenvalues in ascending order.
44 *
45 * E (input/output) DOUBLE PRECISION array, dimension (N-1)
46 * On entry, the (n-1) subdiagonal elements of the tridiagonal
47 * matrix.
48 * On exit, E has been destroyed.
49 *
50 * Z (input/output) DOUBLE PRECISION array, dimension (LDZ, N)
51 * On entry, if COMPZ = 'V', then Z contains the orthogonal
52 * matrix used in the reduction to tridiagonal form.
53 * On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
54 * orthonormal eigenvectors of the original symmetric matrix,
55 * and if COMPZ = 'I', Z contains the orthonormal eigenvectors
56 * of the symmetric tridiagonal matrix.
57 * If COMPZ = 'N', then Z is not referenced.
58 *
59 * LDZ (input) INTEGER
60 * The leading dimension of the array Z. LDZ >= 1, and if
61 * eigenvectors are desired, then LDZ >= max(1,N).
62 *
63 * WORK (workspace) DOUBLE PRECISION array, dimension (max(1,2*N-2))
64 * If COMPZ = 'N', then WORK is not referenced.
65 *
66 * INFO (output) INTEGER
67 * = 0: successful exit
68 * < 0: if INFO = -i, the i-th argument had an illegal value
69 * > 0: the algorithm has failed to find all the eigenvalues in
70 * a total of 30*N iterations; if INFO = i, then i
71 * elements of E have not converged to zero; on exit, D
72 * and E contain the elements of a symmetric tridiagonal
73 * matrix which is orthogonally similar to the original
74 * matrix.
75 *
76 * =====================================================================
77 *
78 * .. Parameters ..
79 DOUBLE PRECISION ZERO, ONE, TWO, THREE
80 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
81 $ THREE = 3.0D0 )
82 INTEGER MAXIT
83 PARAMETER ( MAXIT = 30 )
84 * ..
85 * .. Local Scalars ..
86 INTEGER I, ICOMPZ, II, ISCALE, J, JTOT, K, L, L1, LEND,
87 $ LENDM1, LENDP1, LENDSV, LM1, LSV, M, MM, MM1,
88 $ NM1, NMAXIT
89 DOUBLE PRECISION ANORM, B, C, EPS, EPS2, F, G, P, R, RT1, RT2,
90 $ S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST
91 * ..
92 * .. External Functions ..
93 LOGICAL LSAME
94 DOUBLE PRECISION DLAMCH, DLANST, DLAPY2
95 EXTERNAL LSAME, DLAMCH, DLANST, DLAPY2
96 * ..
97 * .. External Subroutines ..
98 EXTERNAL DLAE2, DLAEV2, DLARTG, DLASCL, DLASET, DLASR,
99 $ DLASRT, DSWAP, XERBLA
100 * ..
101 * .. Intrinsic Functions ..
102 INTRINSIC ABS, MAX, SIGN, SQRT
103 * ..
104 * .. Executable Statements ..
105 *
106 * Test the input parameters.
107 *
108 INFO = 0
109 *
110 IF( LSAME( COMPZ, 'N' ) ) THEN
111 ICOMPZ = 0
112 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
113 ICOMPZ = 1
114 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
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( 'M', 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