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