![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zsteqr.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack|
Return to zsteqr.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 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