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