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