![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsterf.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 dsterf.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE DSTERF( N, D, E, INFO )
2: *
1.8 ! bertrand 3: * -- LAPACK routine (version 3.3.1) --
1.1 bertrand 4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8 ! bertrand 6: * -- April 2011 --
1.1 bertrand 7: *
8: * .. Scalar Arguments ..
9: INTEGER INFO, N
10: * ..
11: * .. Array Arguments ..
12: DOUBLE PRECISION D( * ), E( * )
13: * ..
14: *
15: * Purpose
16: * =======
17: *
18: * DSTERF computes all eigenvalues of a symmetric tridiagonal matrix
19: * using the Pal-Walker-Kahan variant of the QL or QR algorithm.
20: *
21: * Arguments
22: * =========
23: *
24: * N (input) INTEGER
25: * The order of the matrix. N >= 0.
26: *
27: * D (input/output) DOUBLE PRECISION array, dimension (N)
28: * On entry, the n diagonal elements of the tridiagonal matrix.
29: * On exit, if INFO = 0, the eigenvalues in ascending order.
30: *
31: * E (input/output) DOUBLE PRECISION array, dimension (N-1)
32: * On entry, the (n-1) subdiagonal elements of the tridiagonal
33: * matrix.
34: * On exit, E has been destroyed.
35: *
36: * INFO (output) INTEGER
37: * = 0: successful exit
38: * < 0: if INFO = -i, the i-th argument had an illegal value
39: * > 0: the algorithm failed to find all of the eigenvalues in
40: * a total of 30*N iterations; if INFO = i, then i
41: * elements of E have not converged to zero.
42: *
43: * =====================================================================
44: *
45: * .. Parameters ..
46: DOUBLE PRECISION ZERO, ONE, TWO, THREE
47: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
48: $ THREE = 3.0D0 )
49: INTEGER MAXIT
50: PARAMETER ( MAXIT = 30 )
51: * ..
52: * .. Local Scalars ..
53: INTEGER I, ISCALE, JTOT, L, L1, LEND, LENDSV, LSV, M,
54: $ NMAXIT
55: DOUBLE PRECISION ALPHA, ANORM, BB, C, EPS, EPS2, GAMMA, OLDC,
56: $ OLDGAM, P, R, RT1, RT2, RTE, S, SAFMAX, SAFMIN,
1.8 ! bertrand 57: $ SIGMA, SSFMAX, SSFMIN, RMAX
1.1 bertrand 58: * ..
59: * .. External Functions ..
60: DOUBLE PRECISION DLAMCH, DLANST, DLAPY2
61: EXTERNAL DLAMCH, DLANST, DLAPY2
62: * ..
63: * .. External Subroutines ..
64: EXTERNAL DLAE2, DLASCL, DLASRT, XERBLA
65: * ..
66: * .. Intrinsic Functions ..
67: INTRINSIC ABS, SIGN, SQRT
68: * ..
69: * .. Executable Statements ..
70: *
71: * Test the input parameters.
72: *
73: INFO = 0
74: *
75: * Quick return if possible
76: *
77: IF( N.LT.0 ) THEN
78: INFO = -1
79: CALL XERBLA( 'DSTERF', -INFO )
80: RETURN
81: END IF
82: IF( N.LE.1 )
83: $ RETURN
84: *
85: * Determine the unit roundoff for this environment.
86: *
87: EPS = DLAMCH( 'E' )
88: EPS2 = EPS**2
89: SAFMIN = DLAMCH( 'S' )
90: SAFMAX = ONE / SAFMIN
91: SSFMAX = SQRT( SAFMAX ) / THREE
92: SSFMIN = SQRT( SAFMIN ) / EPS2
1.8 ! bertrand 93: RMAX = DLAMCH( 'O' )
1.1 bertrand 94: *
95: * Compute the eigenvalues of the tridiagonal matrix.
96: *
97: NMAXIT = N*MAXIT
98: SIGMA = ZERO
99: JTOT = 0
100: *
101: * Determine where the matrix splits and choose QL or QR iteration
102: * for each block, according to whether top or bottom diagonal
103: * element is smaller.
104: *
105: L1 = 1
106: *
107: 10 CONTINUE
108: IF( L1.GT.N )
109: $ GO TO 170
110: IF( L1.GT.1 )
111: $ E( L1-1 ) = ZERO
112: DO 20 M = L1, N - 1
113: IF( ABS( E( M ) ).LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+
114: $ 1 ) ) ) )*EPS ) THEN
115: E( M ) = ZERO
116: GO TO 30
117: END IF
118: 20 CONTINUE
119: M = N
120: *
121: 30 CONTINUE
122: L = L1
123: LSV = L
124: LEND = M
125: LENDSV = LEND
126: L1 = M + 1
127: IF( LEND.EQ.L )
128: $ GO TO 10
129: *
130: * Scale submatrix in rows and columns L to LEND
131: *
1.8 ! bertrand 132: ANORM = DLANST( 'M', LEND-L+1, D( L ), E( L ) )
1.1 bertrand 133: ISCALE = 0
1.8 ! bertrand 134: IF( ANORM.EQ.ZERO )
! 135: $ GO TO 10
! 136: IF( (ANORM.GT.SSFMAX) ) THEN
1.1 bertrand 137: ISCALE = 1
138: CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,
139: $ INFO )
140: CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
141: $ INFO )
142: ELSE IF( ANORM.LT.SSFMIN ) THEN
143: ISCALE = 2
144: CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,
145: $ INFO )
146: CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
147: $ INFO )
148: END IF
149: *
150: DO 40 I = L, LEND - 1
151: E( I ) = E( I )**2
152: 40 CONTINUE
153: *
154: * Choose between QL and QR iteration
155: *
156: IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
157: LEND = LSV
158: L = LENDSV
159: END IF
160: *
161: IF( LEND.GE.L ) THEN
162: *
163: * QL Iteration
164: *
165: * Look for small subdiagonal element.
166: *
167: 50 CONTINUE
168: IF( L.NE.LEND ) THEN
169: DO 60 M = L, LEND - 1
170: IF( ABS( E( M ) ).LE.EPS2*ABS( D( M )*D( M+1 ) ) )
171: $ GO TO 70
172: 60 CONTINUE
173: END IF
174: M = LEND
175: *
176: 70 CONTINUE
177: IF( M.LT.LEND )
178: $ E( M ) = ZERO
179: P = D( L )
180: IF( M.EQ.L )
181: $ GO TO 90
182: *
183: * If remaining matrix is 2 by 2, use DLAE2 to compute its
184: * eigenvalues.
185: *
186: IF( M.EQ.L+1 ) THEN
187: RTE = SQRT( E( L ) )
188: CALL DLAE2( D( L ), RTE, D( L+1 ), RT1, RT2 )
189: D( L ) = RT1
190: D( L+1 ) = RT2
191: E( L ) = ZERO
192: L = L + 2
193: IF( L.LE.LEND )
194: $ GO TO 50
195: GO TO 150
196: END IF
197: *
198: IF( JTOT.EQ.NMAXIT )
199: $ GO TO 150
200: JTOT = JTOT + 1
201: *
202: * Form shift.
203: *
204: RTE = SQRT( E( L ) )
205: SIGMA = ( D( L+1 )-P ) / ( TWO*RTE )
206: R = DLAPY2( SIGMA, ONE )
207: SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
208: *
209: C = ONE
210: S = ZERO
211: GAMMA = D( M ) - SIGMA
212: P = GAMMA*GAMMA
213: *
214: * Inner loop
215: *
216: DO 80 I = M - 1, L, -1
217: BB = E( I )
218: R = P + BB
219: IF( I.NE.M-1 )
220: $ E( I+1 ) = S*R
221: OLDC = C
222: C = P / R
223: S = BB / R
224: OLDGAM = GAMMA
225: ALPHA = D( I )
226: GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
227: D( I+1 ) = OLDGAM + ( ALPHA-GAMMA )
228: IF( C.NE.ZERO ) THEN
229: P = ( GAMMA*GAMMA ) / C
230: ELSE
231: P = OLDC*BB
232: END IF
233: 80 CONTINUE
234: *
235: E( L ) = S*P
236: D( L ) = SIGMA + GAMMA
237: GO TO 50
238: *
239: * Eigenvalue found.
240: *
241: 90 CONTINUE
242: D( L ) = P
243: *
244: L = L + 1
245: IF( L.LE.LEND )
246: $ GO TO 50
247: GO TO 150
248: *
249: ELSE
250: *
251: * QR Iteration
252: *
253: * Look for small superdiagonal element.
254: *
255: 100 CONTINUE
256: DO 110 M = L, LEND + 1, -1
257: IF( ABS( E( M-1 ) ).LE.EPS2*ABS( D( M )*D( M-1 ) ) )
258: $ GO TO 120
259: 110 CONTINUE
260: M = LEND
261: *
262: 120 CONTINUE
263: IF( M.GT.LEND )
264: $ E( M-1 ) = ZERO
265: P = D( L )
266: IF( M.EQ.L )
267: $ GO TO 140
268: *
269: * If remaining matrix is 2 by 2, use DLAE2 to compute its
270: * eigenvalues.
271: *
272: IF( M.EQ.L-1 ) THEN
273: RTE = SQRT( E( L-1 ) )
274: CALL DLAE2( D( L ), RTE, D( L-1 ), RT1, RT2 )
275: D( L ) = RT1
276: D( L-1 ) = RT2
277: E( L-1 ) = ZERO
278: L = L - 2
279: IF( L.GE.LEND )
280: $ GO TO 100
281: GO TO 150
282: END IF
283: *
284: IF( JTOT.EQ.NMAXIT )
285: $ GO TO 150
286: JTOT = JTOT + 1
287: *
288: * Form shift.
289: *
290: RTE = SQRT( E( L-1 ) )
291: SIGMA = ( D( L-1 )-P ) / ( TWO*RTE )
292: R = DLAPY2( SIGMA, ONE )
293: SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
294: *
295: C = ONE
296: S = ZERO
297: GAMMA = D( M ) - SIGMA
298: P = GAMMA*GAMMA
299: *
300: * Inner loop
301: *
302: DO 130 I = M, L - 1
303: BB = E( I )
304: R = P + BB
305: IF( I.NE.M )
306: $ E( I-1 ) = S*R
307: OLDC = C
308: C = P / R
309: S = BB / R
310: OLDGAM = GAMMA
311: ALPHA = D( I+1 )
312: GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
313: D( I ) = OLDGAM + ( ALPHA-GAMMA )
314: IF( C.NE.ZERO ) THEN
315: P = ( GAMMA*GAMMA ) / C
316: ELSE
317: P = OLDC*BB
318: END IF
319: 130 CONTINUE
320: *
321: E( L-1 ) = S*P
322: D( L ) = SIGMA + GAMMA
323: GO TO 100
324: *
325: * Eigenvalue found.
326: *
327: 140 CONTINUE
328: D( L ) = P
329: *
330: L = L - 1
331: IF( L.GE.LEND )
332: $ GO TO 100
333: GO TO 150
334: *
335: END IF
336: *
337: * Undo scaling if necessary
338: *
339: 150 CONTINUE
340: IF( ISCALE.EQ.1 )
341: $ CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,
342: $ D( LSV ), N, INFO )
343: IF( ISCALE.EQ.2 )
344: $ CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,
345: $ D( LSV ), N, INFO )
346: *
347: * Check for no convergence to an eigenvalue after a total
348: * of N*MAXIT iterations.
349: *
350: IF( JTOT.LT.NMAXIT )
351: $ GO TO 10
352: DO 160 I = 1, N - 1
353: IF( E( I ).NE.ZERO )
354: $ INFO = INFO + 1
355: 160 CONTINUE
356: GO TO 180
357: *
358: * Sort eigenvalues in increasing order.
359: *
360: 170 CONTINUE
361: CALL DLASRT( 'I', N, D, INFO )
362: *
363: 180 CONTINUE
364: RETURN
365: *
366: * End of DSTERF
367: *
368: END