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: *> \date November 2011
129: *
130: *> \ingroup complex16OTHERcomputational
131: *
132: * =====================================================================
133: SUBROUTINE ZSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
134: *
135: * -- LAPACK computational routine (version 3.4.0) --
136: * -- LAPACK is a software package provided by Univ. of Tennessee, --
137: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
138: * November 2011
139: *
140: * .. Scalar Arguments ..
141: CHARACTER COMPZ
142: INTEGER INFO, LDZ, N
143: * ..
144: * .. Array Arguments ..
145: DOUBLE PRECISION D( * ), E( * ), WORK( * )
146: COMPLEX*16 Z( LDZ, * )
147: * ..
148: *
149: * =====================================================================
150: *
151: * .. Parameters ..
152: DOUBLE PRECISION ZERO, ONE, TWO, THREE
153: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
154: $ THREE = 3.0D0 )
155: COMPLEX*16 CZERO, CONE
156: PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ),
157: $ CONE = ( 1.0D0, 0.0D0 ) )
158: INTEGER MAXIT
159: PARAMETER ( MAXIT = 30 )
160: * ..
161: * .. Local Scalars ..
162: INTEGER I, ICOMPZ, II, ISCALE, J, JTOT, K, L, L1, LEND,
163: $ LENDM1, LENDP1, LENDSV, LM1, LSV, M, MM, MM1,
164: $ NM1, NMAXIT
165: DOUBLE PRECISION ANORM, B, C, EPS, EPS2, F, G, P, R, RT1, RT2,
166: $ S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST
167: * ..
168: * .. External Functions ..
169: LOGICAL LSAME
170: DOUBLE PRECISION DLAMCH, DLANST, DLAPY2
171: EXTERNAL LSAME, DLAMCH, DLANST, DLAPY2
172: * ..
173: * .. External Subroutines ..
174: EXTERNAL DLAE2, DLAEV2, DLARTG, DLASCL, DLASRT, XERBLA,
175: $ ZLASET, ZLASR, ZSWAP
176: * ..
177: * .. Intrinsic Functions ..
178: INTRINSIC ABS, MAX, SIGN, SQRT
179: * ..
180: * .. Executable Statements ..
181: *
182: * Test the input parameters.
183: *
184: INFO = 0
185: *
186: IF( LSAME( COMPZ, 'N' ) ) THEN
187: ICOMPZ = 0
188: ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
189: ICOMPZ = 1
190: ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
191: ICOMPZ = 2
192: ELSE
193: ICOMPZ = -1
194: END IF
195: IF( ICOMPZ.LT.0 ) THEN
196: INFO = -1
197: ELSE IF( N.LT.0 ) THEN
198: INFO = -2
199: ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1,
200: $ N ) ) ) THEN
201: INFO = -6
202: END IF
203: IF( INFO.NE.0 ) THEN
204: CALL XERBLA( 'ZSTEQR', -INFO )
205: RETURN
206: END IF
207: *
208: * Quick return if possible
209: *
210: IF( N.EQ.0 )
211: $ RETURN
212: *
213: IF( N.EQ.1 ) THEN
214: IF( ICOMPZ.EQ.2 )
215: $ Z( 1, 1 ) = CONE
216: RETURN
217: END IF
218: *
219: * Determine the unit roundoff and over/underflow thresholds.
220: *
221: EPS = DLAMCH( 'E' )
222: EPS2 = EPS**2
223: SAFMIN = DLAMCH( 'S' )
224: SAFMAX = ONE / SAFMIN
225: SSFMAX = SQRT( SAFMAX ) / THREE
226: SSFMIN = SQRT( SAFMIN ) / EPS2
227: *
228: * Compute the eigenvalues and eigenvectors of the tridiagonal
229: * matrix.
230: *
231: IF( ICOMPZ.EQ.2 )
232: $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDZ )
233: *
234: NMAXIT = N*MAXIT
235: JTOT = 0
236: *
237: * Determine where the matrix splits and choose QL or QR iteration
238: * for each block, according to whether top or bottom diagonal
239: * element is smaller.
240: *
241: L1 = 1
242: NM1 = N - 1
243: *
244: 10 CONTINUE
245: IF( L1.GT.N )
246: $ GO TO 160
247: IF( L1.GT.1 )
248: $ E( L1-1 ) = ZERO
249: IF( L1.LE.NM1 ) THEN
250: DO 20 M = L1, NM1
251: TST = ABS( E( M ) )
252: IF( TST.EQ.ZERO )
253: $ GO TO 30
254: IF( TST.LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+
255: $ 1 ) ) ) )*EPS ) THEN
256: E( M ) = ZERO
257: GO TO 30
258: END IF
259: 20 CONTINUE
260: END IF
261: M = N
262: *
263: 30 CONTINUE
264: L = L1
265: LSV = L
266: LEND = M
267: LENDSV = LEND
268: L1 = M + 1
269: IF( LEND.EQ.L )
270: $ GO TO 10
271: *
272: * Scale submatrix in rows and columns L to LEND
273: *
274: ANORM = DLANST( 'I', LEND-L+1, D( L ), E( L ) )
275: ISCALE = 0
276: IF( ANORM.EQ.ZERO )
277: $ GO TO 10
278: IF( ANORM.GT.SSFMAX ) THEN
279: ISCALE = 1
280: CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,
281: $ INFO )
282: CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
283: $ INFO )
284: ELSE IF( ANORM.LT.SSFMIN ) THEN
285: ISCALE = 2
286: CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,
287: $ INFO )
288: CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
289: $ INFO )
290: END IF
291: *
292: * Choose between QL and QR iteration
293: *
294: IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
295: LEND = LSV
296: L = LENDSV
297: END IF
298: *
299: IF( LEND.GT.L ) THEN
300: *
301: * QL Iteration
302: *
303: * Look for small subdiagonal element.
304: *
305: 40 CONTINUE
306: IF( L.NE.LEND ) THEN
307: LENDM1 = LEND - 1
308: DO 50 M = L, LENDM1
309: TST = ABS( E( M ) )**2
310: IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M+1 ) )+
311: $ SAFMIN )GO TO 60
312: 50 CONTINUE
313: END IF
314: *
315: M = LEND
316: *
317: 60 CONTINUE
318: IF( M.LT.LEND )
319: $ E( M ) = ZERO
320: P = D( L )
321: IF( M.EQ.L )
322: $ GO TO 80
323: *
324: * If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
325: * to compute its eigensystem.
326: *
327: IF( M.EQ.L+1 ) THEN
328: IF( ICOMPZ.GT.0 ) THEN
329: CALL DLAEV2( D( L ), E( L ), D( L+1 ), RT1, RT2, C, S )
330: WORK( L ) = C
331: WORK( N-1+L ) = S
332: CALL ZLASR( 'R', 'V', 'B', N, 2, WORK( L ),
333: $ WORK( N-1+L ), Z( 1, L ), LDZ )
334: ELSE
335: CALL DLAE2( D( L ), E( L ), D( L+1 ), RT1, RT2 )
336: END IF
337: D( L ) = RT1
338: D( L+1 ) = RT2
339: E( L ) = ZERO
340: L = L + 2
341: IF( L.LE.LEND )
342: $ GO TO 40
343: GO TO 140
344: END IF
345: *
346: IF( JTOT.EQ.NMAXIT )
347: $ GO TO 140
348: JTOT = JTOT + 1
349: *
350: * Form shift.
351: *
352: G = ( D( L+1 )-P ) / ( TWO*E( L ) )
353: R = DLAPY2( G, ONE )
354: G = D( M ) - P + ( E( L ) / ( G+SIGN( R, G ) ) )
355: *
356: S = ONE
357: C = ONE
358: P = ZERO
359: *
360: * Inner loop
361: *
362: MM1 = M - 1
363: DO 70 I = MM1, L, -1
364: F = S*E( I )
365: B = C*E( I )
366: CALL DLARTG( G, F, C, S, R )
367: IF( I.NE.M-1 )
368: $ E( I+1 ) = R
369: G = D( I+1 ) - P
370: R = ( D( I )-G )*S + TWO*C*B
371: P = S*R
372: D( I+1 ) = G + P
373: G = C*R - B
374: *
375: * If eigenvectors are desired, then save rotations.
376: *
377: IF( ICOMPZ.GT.0 ) THEN
378: WORK( I ) = C
379: WORK( N-1+I ) = -S
380: END IF
381: *
382: 70 CONTINUE
383: *
384: * If eigenvectors are desired, then apply saved rotations.
385: *
386: IF( ICOMPZ.GT.0 ) THEN
387: MM = M - L + 1
388: CALL ZLASR( 'R', 'V', 'B', N, MM, WORK( L ), WORK( N-1+L ),
389: $ Z( 1, L ), LDZ )
390: END IF
391: *
392: D( L ) = D( L ) - P
393: E( L ) = G
394: GO TO 40
395: *
396: * Eigenvalue found.
397: *
398: 80 CONTINUE
399: D( L ) = P
400: *
401: L = L + 1
402: IF( L.LE.LEND )
403: $ GO TO 40
404: GO TO 140
405: *
406: ELSE
407: *
408: * QR Iteration
409: *
410: * Look for small superdiagonal element.
411: *
412: 90 CONTINUE
413: IF( L.NE.LEND ) THEN
414: LENDP1 = LEND + 1
415: DO 100 M = L, LENDP1, -1
416: TST = ABS( E( M-1 ) )**2
417: IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M-1 ) )+
418: $ SAFMIN )GO TO 110
419: 100 CONTINUE
420: END IF
421: *
422: M = LEND
423: *
424: 110 CONTINUE
425: IF( M.GT.LEND )
426: $ E( M-1 ) = ZERO
427: P = D( L )
428: IF( M.EQ.L )
429: $ GO TO 130
430: *
431: * If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
432: * to compute its eigensystem.
433: *
434: IF( M.EQ.L-1 ) THEN
435: IF( ICOMPZ.GT.0 ) THEN
436: CALL DLAEV2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2, C, S )
437: WORK( M ) = C
438: WORK( N-1+M ) = S
439: CALL ZLASR( 'R', 'V', 'F', N, 2, WORK( M ),
440: $ WORK( N-1+M ), Z( 1, L-1 ), LDZ )
441: ELSE
442: CALL DLAE2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2 )
443: END IF
444: D( L-1 ) = RT1
445: D( L ) = RT2
446: E( L-1 ) = ZERO
447: L = L - 2
448: IF( L.GE.LEND )
449: $ GO TO 90
450: GO TO 140
451: END IF
452: *
453: IF( JTOT.EQ.NMAXIT )
454: $ GO TO 140
455: JTOT = JTOT + 1
456: *
457: * Form shift.
458: *
459: G = ( D( L-1 )-P ) / ( TWO*E( L-1 ) )
460: R = DLAPY2( G, ONE )
461: G = D( M ) - P + ( E( L-1 ) / ( G+SIGN( R, G ) ) )
462: *
463: S = ONE
464: C = ONE
465: P = ZERO
466: *
467: * Inner loop
468: *
469: LM1 = L - 1
470: DO 120 I = M, LM1
471: F = S*E( I )
472: B = C*E( I )
473: CALL DLARTG( G, F, C, S, R )
474: IF( I.NE.M )
475: $ E( I-1 ) = R
476: G = D( I ) - P
477: R = ( D( I+1 )-G )*S + TWO*C*B
478: P = S*R
479: D( I ) = G + P
480: G = C*R - B
481: *
482: * If eigenvectors are desired, then save rotations.
483: *
484: IF( ICOMPZ.GT.0 ) THEN
485: WORK( I ) = C
486: WORK( N-1+I ) = S
487: END IF
488: *
489: 120 CONTINUE
490: *
491: * If eigenvectors are desired, then apply saved rotations.
492: *
493: IF( ICOMPZ.GT.0 ) THEN
494: MM = L - M + 1
495: CALL ZLASR( 'R', 'V', 'F', N, MM, WORK( M ), WORK( N-1+M ),
496: $ Z( 1, M ), LDZ )
497: END IF
498: *
499: D( L ) = D( L ) - P
500: E( LM1 ) = G
501: GO TO 90
502: *
503: * Eigenvalue found.
504: *
505: 130 CONTINUE
506: D( L ) = P
507: *
508: L = L - 1
509: IF( L.GE.LEND )
510: $ GO TO 90
511: GO TO 140
512: *
513: END IF
514: *
515: * Undo scaling if necessary
516: *
517: 140 CONTINUE
518: IF( ISCALE.EQ.1 ) THEN
519: CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,
520: $ D( LSV ), N, INFO )
521: CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV, 1, E( LSV ),
522: $ N, INFO )
523: ELSE IF( ISCALE.EQ.2 ) THEN
524: CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,
525: $ D( LSV ), N, INFO )
526: CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV, 1, E( LSV ),
527: $ N, INFO )
528: END IF
529: *
530: * Check for no convergence to an eigenvalue after a total
531: * of N*MAXIT iterations.
532: *
533: IF( JTOT.EQ.NMAXIT ) THEN
534: DO 150 I = 1, N - 1
535: IF( E( I ).NE.ZERO )
536: $ INFO = INFO + 1
537: 150 CONTINUE
538: RETURN
539: END IF
540: GO TO 10
541: *
542: * Order eigenvalues and eigenvectors.
543: *
544: 160 CONTINUE
545: IF( ICOMPZ.EQ.0 ) THEN
546: *
547: * Use Quick Sort
548: *
549: CALL DLASRT( 'I', N, D, INFO )
550: *
551: ELSE
552: *
553: * Use Selection Sort to minimize swaps of eigenvectors
554: *
555: DO 180 II = 2, N
556: I = II - 1
557: K = I
558: P = D( I )
559: DO 170 J = II, N
560: IF( D( J ).LT.P ) THEN
561: K = J
562: P = D( J )
563: END IF
564: 170 CONTINUE
565: IF( K.NE.I ) THEN
566: D( K ) = D( I )
567: D( I ) = P
568: CALL ZSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
569: END IF
570: 180 CONTINUE
571: END IF
572: RETURN
573: *
574: * End of ZSTEQR
575: *
576: 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