![[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: *
3: * -- LAPACK routine (version 3.2) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * November 2006
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,
57: $ SIGMA, SSFMAX, SSFMIN
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
93: *
94: * Compute the eigenvalues of the tridiagonal matrix.
95: *
96: NMAXIT = N*MAXIT
97: SIGMA = ZERO
98: JTOT = 0
99: *
100: * Determine where the matrix splits and choose QL or QR iteration
101: * for each block, according to whether top or bottom diagonal
102: * element is smaller.
103: *
104: L1 = 1
105: *
106: 10 CONTINUE
107: IF( L1.GT.N )
108: $ GO TO 170
109: IF( L1.GT.1 )
110: $ E( L1-1 ) = ZERO
111: DO 20 M = L1, N - 1
112: IF( ABS( E( M ) ).LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+
113: $ 1 ) ) ) )*EPS ) THEN
114: E( M ) = ZERO
115: GO TO 30
116: END IF
117: 20 CONTINUE
118: M = N
119: *
120: 30 CONTINUE
121: L = L1
122: LSV = L
123: LEND = M
124: LENDSV = LEND
125: L1 = M + 1
126: IF( LEND.EQ.L )
127: $ GO TO 10
128: *
129: * Scale submatrix in rows and columns L to LEND
130: *
131: ANORM = DLANST( 'I', LEND-L+1, D( L ), E( L ) )
132: ISCALE = 0
133: IF( ANORM.GT.SSFMAX ) THEN
134: ISCALE = 1
135: CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,
136: $ INFO )
137: CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
138: $ INFO )
139: ELSE IF( ANORM.LT.SSFMIN ) THEN
140: ISCALE = 2
141: CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,
142: $ INFO )
143: CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
144: $ INFO )
145: END IF
146: *
147: DO 40 I = L, LEND - 1
148: E( I ) = E( I )**2
149: 40 CONTINUE
150: *
151: * Choose between QL and QR iteration
152: *
153: IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
154: LEND = LSV
155: L = LENDSV
156: END IF
157: *
158: IF( LEND.GE.L ) THEN
159: *
160: * QL Iteration
161: *
162: * Look for small subdiagonal element.
163: *
164: 50 CONTINUE
165: IF( L.NE.LEND ) THEN
166: DO 60 M = L, LEND - 1
167: IF( ABS( E( M ) ).LE.EPS2*ABS( D( M )*D( M+1 ) ) )
168: $ GO TO 70
169: 60 CONTINUE
170: END IF
171: M = LEND
172: *
173: 70 CONTINUE
174: IF( M.LT.LEND )
175: $ E( M ) = ZERO
176: P = D( L )
177: IF( M.EQ.L )
178: $ GO TO 90
179: *
180: * If remaining matrix is 2 by 2, use DLAE2 to compute its
181: * eigenvalues.
182: *
183: IF( M.EQ.L+1 ) THEN
184: RTE = SQRT( E( L ) )
185: CALL DLAE2( D( L ), RTE, D( L+1 ), RT1, RT2 )
186: D( L ) = RT1
187: D( L+1 ) = RT2
188: E( L ) = ZERO
189: L = L + 2
190: IF( L.LE.LEND )
191: $ GO TO 50
192: GO TO 150
193: END IF
194: *
195: IF( JTOT.EQ.NMAXIT )
196: $ GO TO 150
197: JTOT = JTOT + 1
198: *
199: * Form shift.
200: *
201: RTE = SQRT( E( L ) )
202: SIGMA = ( D( L+1 )-P ) / ( TWO*RTE )
203: R = DLAPY2( SIGMA, ONE )
204: SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
205: *
206: C = ONE
207: S = ZERO
208: GAMMA = D( M ) - SIGMA
209: P = GAMMA*GAMMA
210: *
211: * Inner loop
212: *
213: DO 80 I = M - 1, L, -1
214: BB = E( I )
215: R = P + BB
216: IF( I.NE.M-1 )
217: $ E( I+1 ) = S*R
218: OLDC = C
219: C = P / R
220: S = BB / R
221: OLDGAM = GAMMA
222: ALPHA = D( I )
223: GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
224: D( I+1 ) = OLDGAM + ( ALPHA-GAMMA )
225: IF( C.NE.ZERO ) THEN
226: P = ( GAMMA*GAMMA ) / C
227: ELSE
228: P = OLDC*BB
229: END IF
230: 80 CONTINUE
231: *
232: E( L ) = S*P
233: D( L ) = SIGMA + GAMMA
234: GO TO 50
235: *
236: * Eigenvalue found.
237: *
238: 90 CONTINUE
239: D( L ) = P
240: *
241: L = L + 1
242: IF( L.LE.LEND )
243: $ GO TO 50
244: GO TO 150
245: *
246: ELSE
247: *
248: * QR Iteration
249: *
250: * Look for small superdiagonal element.
251: *
252: 100 CONTINUE
253: DO 110 M = L, LEND + 1, -1
254: IF( ABS( E( M-1 ) ).LE.EPS2*ABS( D( M )*D( M-1 ) ) )
255: $ GO TO 120
256: 110 CONTINUE
257: M = LEND
258: *
259: 120 CONTINUE
260: IF( M.GT.LEND )
261: $ E( M-1 ) = ZERO
262: P = D( L )
263: IF( M.EQ.L )
264: $ GO TO 140
265: *
266: * If remaining matrix is 2 by 2, use DLAE2 to compute its
267: * eigenvalues.
268: *
269: IF( M.EQ.L-1 ) THEN
270: RTE = SQRT( E( L-1 ) )
271: CALL DLAE2( D( L ), RTE, D( L-1 ), RT1, RT2 )
272: D( L ) = RT1
273: D( L-1 ) = RT2
274: E( L-1 ) = ZERO
275: L = L - 2
276: IF( L.GE.LEND )
277: $ GO TO 100
278: GO TO 150
279: END IF
280: *
281: IF( JTOT.EQ.NMAXIT )
282: $ GO TO 150
283: JTOT = JTOT + 1
284: *
285: * Form shift.
286: *
287: RTE = SQRT( E( L-1 ) )
288: SIGMA = ( D( L-1 )-P ) / ( TWO*RTE )
289: R = DLAPY2( SIGMA, ONE )
290: SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
291: *
292: C = ONE
293: S = ZERO
294: GAMMA = D( M ) - SIGMA
295: P = GAMMA*GAMMA
296: *
297: * Inner loop
298: *
299: DO 130 I = M, L - 1
300: BB = E( I )
301: R = P + BB
302: IF( I.NE.M )
303: $ E( I-1 ) = S*R
304: OLDC = C
305: C = P / R
306: S = BB / R
307: OLDGAM = GAMMA
308: ALPHA = D( I+1 )
309: GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
310: D( I ) = OLDGAM + ( ALPHA-GAMMA )
311: IF( C.NE.ZERO ) THEN
312: P = ( GAMMA*GAMMA ) / C
313: ELSE
314: P = OLDC*BB
315: END IF
316: 130 CONTINUE
317: *
318: E( L-1 ) = S*P
319: D( L ) = SIGMA + GAMMA
320: GO TO 100
321: *
322: * Eigenvalue found.
323: *
324: 140 CONTINUE
325: D( L ) = P
326: *
327: L = L - 1
328: IF( L.GE.LEND )
329: $ GO TO 100
330: GO TO 150
331: *
332: END IF
333: *
334: * Undo scaling if necessary
335: *
336: 150 CONTINUE
337: IF( ISCALE.EQ.1 )
338: $ CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,
339: $ D( LSV ), N, INFO )
340: IF( ISCALE.EQ.2 )
341: $ CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,
342: $ D( LSV ), N, INFO )
343: *
344: * Check for no convergence to an eigenvalue after a total
345: * of N*MAXIT iterations.
346: *
347: IF( JTOT.LT.NMAXIT )
348: $ GO TO 10
349: DO 160 I = 1, N - 1
350: IF( E( I ).NE.ZERO )
351: $ INFO = INFO + 1
352: 160 CONTINUE
353: GO TO 180
354: *
355: * Sort eigenvalues in increasing order.
356: *
357: 170 CONTINUE
358: CALL DLASRT( 'I', N, D, INFO )
359: *
360: 180 CONTINUE
361: RETURN
362: *
363: * End of DSTERF
364: *
365: END