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