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