![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsteqr.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 dsteqr.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE DSTEQR( 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( * ), 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( '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