|
2329
|
1 SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, |
|
|
2 $ ILOZ, IHIZ, Z, LDZ, INFO ) |
|
|
3 * |
|
3333
|
4 * -- LAPACK auxiliary routine (version 3.0) -- |
|
2329
|
5 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., |
|
|
6 * Courant Institute, Argonne National Lab, and Rice University |
|
3333
|
7 * June 30, 1999 |
|
2329
|
8 * |
|
|
9 * .. Scalar Arguments .. |
|
|
10 LOGICAL WANTT, WANTZ |
|
|
11 INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N |
|
|
12 * .. |
|
|
13 * .. Array Arguments .. |
|
|
14 DOUBLE PRECISION H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * ) |
|
|
15 * .. |
|
|
16 * |
|
|
17 * Purpose |
|
|
18 * ======= |
|
|
19 * |
|
|
20 * DLAHQR is an auxiliary routine called by DHSEQR to update the |
|
|
21 * eigenvalues and Schur decomposition already computed by DHSEQR, by |
|
|
22 * dealing with the Hessenberg submatrix in rows and columns ILO to IHI. |
|
|
23 * |
|
|
24 * Arguments |
|
|
25 * ========= |
|
|
26 * |
|
|
27 * WANTT (input) LOGICAL |
|
|
28 * = .TRUE. : the full Schur form T is required; |
|
|
29 * = .FALSE.: only eigenvalues are required. |
|
|
30 * |
|
|
31 * WANTZ (input) LOGICAL |
|
|
32 * = .TRUE. : the matrix of Schur vectors Z is required; |
|
|
33 * = .FALSE.: Schur vectors are not required. |
|
|
34 * |
|
|
35 * N (input) INTEGER |
|
|
36 * The order of the matrix H. N >= 0. |
|
|
37 * |
|
|
38 * ILO (input) INTEGER |
|
|
39 * IHI (input) INTEGER |
|
|
40 * It is assumed that H is already upper quasi-triangular in |
|
|
41 * rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless |
|
|
42 * ILO = 1). DLAHQR works primarily with the Hessenberg |
|
|
43 * submatrix in rows and columns ILO to IHI, but applies |
|
|
44 * transformations to all of H if WANTT is .TRUE.. |
|
|
45 * 1 <= ILO <= max(1,IHI); IHI <= N. |
|
|
46 * |
|
|
47 * H (input/output) DOUBLE PRECISION array, dimension (LDH,N) |
|
|
48 * On entry, the upper Hessenberg matrix H. |
|
|
49 * On exit, if WANTT is .TRUE., H is upper quasi-triangular in |
|
|
50 * rows and columns ILO:IHI, with any 2-by-2 diagonal blocks in |
|
|
51 * standard form. If WANTT is .FALSE., the contents of H are |
|
|
52 * unspecified on exit. |
|
|
53 * |
|
|
54 * LDH (input) INTEGER |
|
|
55 * The leading dimension of the array H. LDH >= max(1,N). |
|
|
56 * |
|
|
57 * WR (output) DOUBLE PRECISION array, dimension (N) |
|
|
58 * WI (output) DOUBLE PRECISION array, dimension (N) |
|
|
59 * The real and imaginary parts, respectively, of the computed |
|
|
60 * eigenvalues ILO to IHI are stored in the corresponding |
|
|
61 * elements of WR and WI. If two eigenvalues are computed as a |
|
|
62 * complex conjugate pair, they are stored in consecutive |
|
|
63 * elements of WR and WI, say the i-th and (i+1)th, with |
|
|
64 * WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the |
|
|
65 * eigenvalues are stored in the same order as on the diagonal |
|
|
66 * of the Schur form returned in H, with WR(i) = H(i,i), and, if |
|
|
67 * H(i:i+1,i:i+1) is a 2-by-2 diagonal block, |
|
|
68 * WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i). |
|
|
69 * |
|
|
70 * ILOZ (input) INTEGER |
|
|
71 * IHIZ (input) INTEGER |
|
|
72 * Specify the rows of Z to which transformations must be |
|
|
73 * applied if WANTZ is .TRUE.. |
|
|
74 * 1 <= ILOZ <= ILO; IHI <= IHIZ <= N. |
|
|
75 * |
|
|
76 * Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N) |
|
|
77 * If WANTZ is .TRUE., on entry Z must contain the current |
|
|
78 * matrix Z of transformations accumulated by DHSEQR, and on |
|
|
79 * exit Z has been updated; transformations are applied only to |
|
|
80 * the submatrix Z(ILOZ:IHIZ,ILO:IHI). |
|
|
81 * If WANTZ is .FALSE., Z is not referenced. |
|
|
82 * |
|
|
83 * LDZ (input) INTEGER |
|
|
84 * The leading dimension of the array Z. LDZ >= max(1,N). |
|
|
85 * |
|
|
86 * INFO (output) INTEGER |
|
|
87 * = 0: successful exit |
|
|
88 * > 0: DLAHQR failed to compute all the eigenvalues ILO to IHI |
|
|
89 * in a total of 30*(IHI-ILO+1) iterations; if INFO = i, |
|
|
90 * elements i+1:ihi of WR and WI contain those eigenvalues |
|
|
91 * which have been successfully computed. |
|
|
92 * |
|
3333
|
93 * Further Details |
|
|
94 * =============== |
|
|
95 * |
|
|
96 * 2-96 Based on modifications by |
|
|
97 * David Day, Sandia National Laboratory, USA |
|
|
98 * |
|
2329
|
99 * ===================================================================== |
|
|
100 * |
|
|
101 * .. Parameters .. |
|
3333
|
102 DOUBLE PRECISION ZERO, ONE, HALF |
|
|
103 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, HALF = 0.5D0 ) |
|
2329
|
104 DOUBLE PRECISION DAT1, DAT2 |
|
|
105 PARAMETER ( DAT1 = 0.75D+0, DAT2 = -0.4375D+0 ) |
|
|
106 * .. |
|
|
107 * .. Local Scalars .. |
|
|
108 INTEGER I, I1, I2, ITN, ITS, J, K, L, M, NH, NR, NZ |
|
3333
|
109 DOUBLE PRECISION AVE, CS, DISC, H00, H10, H11, H12, H21, H22, |
|
|
110 $ H33, H33S, H43H34, H44, H44S, OVFL, S, SMLNUM, |
|
|
111 $ SN, SUM, T1, T2, T3, TST1, ULP, UNFL, V1, V2, |
|
|
112 $ V3 |
|
2329
|
113 * .. |
|
|
114 * .. Local Arrays .. |
|
|
115 DOUBLE PRECISION V( 3 ), WORK( 1 ) |
|
|
116 * .. |
|
|
117 * .. External Functions .. |
|
|
118 DOUBLE PRECISION DLAMCH, DLANHS |
|
|
119 EXTERNAL DLAMCH, DLANHS |
|
|
120 * .. |
|
|
121 * .. External Subroutines .. |
|
3333
|
122 EXTERNAL DCOPY, DLANV2, DLARFG, DROT |
|
2329
|
123 * .. |
|
|
124 * .. Intrinsic Functions .. |
|
3333
|
125 INTRINSIC ABS, MAX, MIN, SIGN, SQRT |
|
2329
|
126 * .. |
|
|
127 * .. Executable Statements .. |
|
|
128 * |
|
|
129 INFO = 0 |
|
|
130 * |
|
|
131 * Quick return if possible |
|
|
132 * |
|
|
133 IF( N.EQ.0 ) |
|
|
134 $ RETURN |
|
|
135 IF( ILO.EQ.IHI ) THEN |
|
|
136 WR( ILO ) = H( ILO, ILO ) |
|
|
137 WI( ILO ) = ZERO |
|
|
138 RETURN |
|
|
139 END IF |
|
|
140 * |
|
|
141 NH = IHI - ILO + 1 |
|
|
142 NZ = IHIZ - ILOZ + 1 |
|
|
143 * |
|
|
144 * Set machine-dependent constants for the stopping criterion. |
|
|
145 * If norm(H) <= sqrt(OVFL), overflow should not occur. |
|
|
146 * |
|
|
147 UNFL = DLAMCH( 'Safe minimum' ) |
|
|
148 OVFL = ONE / UNFL |
|
|
149 CALL DLABAD( UNFL, OVFL ) |
|
|
150 ULP = DLAMCH( 'Precision' ) |
|
|
151 SMLNUM = UNFL*( NH / ULP ) |
|
|
152 * |
|
|
153 * I1 and I2 are the indices of the first row and last column of H |
|
|
154 * to which transformations must be applied. If eigenvalues only are |
|
|
155 * being computed, I1 and I2 are set inside the main loop. |
|
|
156 * |
|
|
157 IF( WANTT ) THEN |
|
|
158 I1 = 1 |
|
|
159 I2 = N |
|
|
160 END IF |
|
|
161 * |
|
|
162 * ITN is the total number of QR iterations allowed. |
|
|
163 * |
|
|
164 ITN = 30*NH |
|
|
165 * |
|
|
166 * The main loop begins here. I is the loop index and decreases from |
|
|
167 * IHI to ILO in steps of 1 or 2. Each iteration of the loop works |
|
|
168 * with the active submatrix in rows and columns L to I. |
|
|
169 * Eigenvalues I+1 to IHI have already converged. Either L = ILO or |
|
|
170 * H(L,L-1) is negligible so that the matrix splits. |
|
|
171 * |
|
|
172 I = IHI |
|
|
173 10 CONTINUE |
|
|
174 L = ILO |
|
|
175 IF( I.LT.ILO ) |
|
|
176 $ GO TO 150 |
|
|
177 * |
|
|
178 * Perform QR iterations on rows and columns ILO to I until a |
|
|
179 * submatrix of order 1 or 2 splits off at the bottom because a |
|
|
180 * subdiagonal element has become negligible. |
|
|
181 * |
|
|
182 DO 130 ITS = 0, ITN |
|
|
183 * |
|
|
184 * Look for a single small subdiagonal element. |
|
|
185 * |
|
|
186 DO 20 K = I, L + 1, -1 |
|
|
187 TST1 = ABS( H( K-1, K-1 ) ) + ABS( H( K, K ) ) |
|
|
188 IF( TST1.EQ.ZERO ) |
|
|
189 $ TST1 = DLANHS( '1', I-L+1, H( L, L ), LDH, WORK ) |
|
|
190 IF( ABS( H( K, K-1 ) ).LE.MAX( ULP*TST1, SMLNUM ) ) |
|
|
191 $ GO TO 30 |
|
|
192 20 CONTINUE |
|
|
193 30 CONTINUE |
|
|
194 L = K |
|
|
195 IF( L.GT.ILO ) THEN |
|
|
196 * |
|
|
197 * H(L,L-1) is negligible |
|
|
198 * |
|
|
199 H( L, L-1 ) = ZERO |
|
|
200 END IF |
|
|
201 * |
|
|
202 * Exit from loop if a submatrix of order 1 or 2 has split off. |
|
|
203 * |
|
|
204 IF( L.GE.I-1 ) |
|
|
205 $ GO TO 140 |
|
|
206 * |
|
|
207 * Now the active submatrix is in rows and columns L to I. If |
|
|
208 * eigenvalues only are being computed, only the active submatrix |
|
|
209 * need be transformed. |
|
|
210 * |
|
|
211 IF( .NOT.WANTT ) THEN |
|
|
212 I1 = L |
|
|
213 I2 = I |
|
|
214 END IF |
|
|
215 * |
|
|
216 IF( ITS.EQ.10 .OR. ITS.EQ.20 ) THEN |
|
|
217 * |
|
|
218 * Exceptional shift. |
|
|
219 * |
|
|
220 S = ABS( H( I, I-1 ) ) + ABS( H( I-1, I-2 ) ) |
|
3333
|
221 H44 = DAT1*S + H( I, I ) |
|
2329
|
222 H33 = H44 |
|
|
223 H43H34 = DAT2*S*S |
|
|
224 ELSE |
|
|
225 * |
|
3333
|
226 * Prepare to use Francis' double shift |
|
|
227 * (i.e. 2nd degree generalized Rayleigh quotient) |
|
2329
|
228 * |
|
|
229 H44 = H( I, I ) |
|
|
230 H33 = H( I-1, I-1 ) |
|
|
231 H43H34 = H( I, I-1 )*H( I-1, I ) |
|
3333
|
232 S = H( I-1, I-2 )*H( I-1, I-2 ) |
|
|
233 DISC = ( H33-H44 )*HALF |
|
|
234 DISC = DISC*DISC + H43H34 |
|
|
235 IF( DISC.GT.ZERO ) THEN |
|
|
236 * |
|
|
237 * Real roots: use Wilkinson's shift twice |
|
|
238 * |
|
|
239 DISC = SQRT( DISC ) |
|
|
240 AVE = HALF*( H33+H44 ) |
|
|
241 IF( ABS( H33 )-ABS( H44 ).GT.ZERO ) THEN |
|
|
242 H33 = H33*H44 - H43H34 |
|
|
243 H44 = H33 / ( SIGN( DISC, AVE )+AVE ) |
|
|
244 ELSE |
|
|
245 H44 = SIGN( DISC, AVE ) + AVE |
|
|
246 END IF |
|
|
247 H33 = H44 |
|
|
248 H43H34 = ZERO |
|
|
249 END IF |
|
2329
|
250 END IF |
|
|
251 * |
|
|
252 * Look for two consecutive small subdiagonal elements. |
|
|
253 * |
|
|
254 DO 40 M = I - 2, L, -1 |
|
|
255 * Determine the effect of starting the double-shift QR |
|
|
256 * iteration at row M, and see if this would make H(M,M-1) |
|
|
257 * negligible. |
|
|
258 * |
|
|
259 H11 = H( M, M ) |
|
|
260 H22 = H( M+1, M+1 ) |
|
|
261 H21 = H( M+1, M ) |
|
|
262 H12 = H( M, M+1 ) |
|
|
263 H44S = H44 - H11 |
|
|
264 H33S = H33 - H11 |
|
|
265 V1 = ( H33S*H44S-H43H34 ) / H21 + H12 |
|
|
266 V2 = H22 - H11 - H33S - H44S |
|
|
267 V3 = H( M+2, M+1 ) |
|
|
268 S = ABS( V1 ) + ABS( V2 ) + ABS( V3 ) |
|
|
269 V1 = V1 / S |
|
|
270 V2 = V2 / S |
|
|
271 V3 = V3 / S |
|
|
272 V( 1 ) = V1 |
|
|
273 V( 2 ) = V2 |
|
|
274 V( 3 ) = V3 |
|
|
275 IF( M.EQ.L ) |
|
|
276 $ GO TO 50 |
|
|
277 H00 = H( M-1, M-1 ) |
|
|
278 H10 = H( M, M-1 ) |
|
|
279 TST1 = ABS( V1 )*( ABS( H00 )+ABS( H11 )+ABS( H22 ) ) |
|
|
280 IF( ABS( H10 )*( ABS( V2 )+ABS( V3 ) ).LE.ULP*TST1 ) |
|
|
281 $ GO TO 50 |
|
|
282 40 CONTINUE |
|
|
283 50 CONTINUE |
|
|
284 * |
|
|
285 * Double-shift QR step |
|
|
286 * |
|
|
287 DO 120 K = M, I - 1 |
|
|
288 * |
|
|
289 * The first iteration of this loop determines a reflection G |
|
|
290 * from the vector V and applies it from left and right to H, |
|
|
291 * thus creating a nonzero bulge below the subdiagonal. |
|
|
292 * |
|
|
293 * Each subsequent iteration determines a reflection G to |
|
|
294 * restore the Hessenberg form in the (K-1)th column, and thus |
|
|
295 * chases the bulge one step toward the bottom of the active |
|
|
296 * submatrix. NR is the order of G. |
|
|
297 * |
|
|
298 NR = MIN( 3, I-K+1 ) |
|
|
299 IF( K.GT.M ) |
|
|
300 $ CALL DCOPY( NR, H( K, K-1 ), 1, V, 1 ) |
|
|
301 CALL DLARFG( NR, V( 1 ), V( 2 ), 1, T1 ) |
|
|
302 IF( K.GT.M ) THEN |
|
|
303 H( K, K-1 ) = V( 1 ) |
|
|
304 H( K+1, K-1 ) = ZERO |
|
|
305 IF( K.LT.I-1 ) |
|
|
306 $ H( K+2, K-1 ) = ZERO |
|
|
307 ELSE IF( M.GT.L ) THEN |
|
|
308 H( K, K-1 ) = -H( K, K-1 ) |
|
|
309 END IF |
|
|
310 V2 = V( 2 ) |
|
|
311 T2 = T1*V2 |
|
|
312 IF( NR.EQ.3 ) THEN |
|
|
313 V3 = V( 3 ) |
|
|
314 T3 = T1*V3 |
|
|
315 * |
|
|
316 * Apply G from the left to transform the rows of the matrix |
|
|
317 * in columns K to I2. |
|
|
318 * |
|
|
319 DO 60 J = K, I2 |
|
|
320 SUM = H( K, J ) + V2*H( K+1, J ) + V3*H( K+2, J ) |
|
|
321 H( K, J ) = H( K, J ) - SUM*T1 |
|
|
322 H( K+1, J ) = H( K+1, J ) - SUM*T2 |
|
|
323 H( K+2, J ) = H( K+2, J ) - SUM*T3 |
|
|
324 60 CONTINUE |
|
|
325 * |
|
|
326 * Apply G from the right to transform the columns of the |
|
|
327 * matrix in rows I1 to min(K+3,I). |
|
|
328 * |
|
|
329 DO 70 J = I1, MIN( K+3, I ) |
|
|
330 SUM = H( J, K ) + V2*H( J, K+1 ) + V3*H( J, K+2 ) |
|
|
331 H( J, K ) = H( J, K ) - SUM*T1 |
|
|
332 H( J, K+1 ) = H( J, K+1 ) - SUM*T2 |
|
|
333 H( J, K+2 ) = H( J, K+2 ) - SUM*T3 |
|
|
334 70 CONTINUE |
|
|
335 * |
|
|
336 IF( WANTZ ) THEN |
|
|
337 * |
|
|
338 * Accumulate transformations in the matrix Z |
|
|
339 * |
|
|
340 DO 80 J = ILOZ, IHIZ |
|
|
341 SUM = Z( J, K ) + V2*Z( J, K+1 ) + V3*Z( J, K+2 ) |
|
|
342 Z( J, K ) = Z( J, K ) - SUM*T1 |
|
|
343 Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2 |
|
|
344 Z( J, K+2 ) = Z( J, K+2 ) - SUM*T3 |
|
|
345 80 CONTINUE |
|
|
346 END IF |
|
|
347 ELSE IF( NR.EQ.2 ) THEN |
|
|
348 * |
|
|
349 * Apply G from the left to transform the rows of the matrix |
|
|
350 * in columns K to I2. |
|
|
351 * |
|
|
352 DO 90 J = K, I2 |
|
|
353 SUM = H( K, J ) + V2*H( K+1, J ) |
|
|
354 H( K, J ) = H( K, J ) - SUM*T1 |
|
|
355 H( K+1, J ) = H( K+1, J ) - SUM*T2 |
|
|
356 90 CONTINUE |
|
|
357 * |
|
|
358 * Apply G from the right to transform the columns of the |
|
|
359 * matrix in rows I1 to min(K+3,I). |
|
|
360 * |
|
|
361 DO 100 J = I1, I |
|
|
362 SUM = H( J, K ) + V2*H( J, K+1 ) |
|
|
363 H( J, K ) = H( J, K ) - SUM*T1 |
|
|
364 H( J, K+1 ) = H( J, K+1 ) - SUM*T2 |
|
|
365 100 CONTINUE |
|
|
366 * |
|
|
367 IF( WANTZ ) THEN |
|
|
368 * |
|
|
369 * Accumulate transformations in the matrix Z |
|
|
370 * |
|
|
371 DO 110 J = ILOZ, IHIZ |
|
|
372 SUM = Z( J, K ) + V2*Z( J, K+1 ) |
|
|
373 Z( J, K ) = Z( J, K ) - SUM*T1 |
|
|
374 Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2 |
|
|
375 110 CONTINUE |
|
|
376 END IF |
|
|
377 END IF |
|
|
378 120 CONTINUE |
|
|
379 * |
|
|
380 130 CONTINUE |
|
|
381 * |
|
|
382 * Failure to converge in remaining number of iterations |
|
|
383 * |
|
|
384 INFO = I |
|
|
385 RETURN |
|
|
386 * |
|
|
387 140 CONTINUE |
|
|
388 * |
|
|
389 IF( L.EQ.I ) THEN |
|
|
390 * |
|
|
391 * H(I,I-1) is negligible: one eigenvalue has converged. |
|
|
392 * |
|
|
393 WR( I ) = H( I, I ) |
|
|
394 WI( I ) = ZERO |
|
|
395 ELSE IF( L.EQ.I-1 ) THEN |
|
|
396 * |
|
|
397 * H(I-1,I-2) is negligible: a pair of eigenvalues have converged. |
|
|
398 * |
|
|
399 * Transform the 2-by-2 submatrix to standard Schur form, |
|
|
400 * and compute and store the eigenvalues. |
|
|
401 * |
|
|
402 CALL DLANV2( H( I-1, I-1 ), H( I-1, I ), H( I, I-1 ), |
|
|
403 $ H( I, I ), WR( I-1 ), WI( I-1 ), WR( I ), WI( I ), |
|
|
404 $ CS, SN ) |
|
|
405 * |
|
|
406 IF( WANTT ) THEN |
|
|
407 * |
|
|
408 * Apply the transformation to the rest of H. |
|
|
409 * |
|
|
410 IF( I2.GT.I ) |
|
|
411 $ CALL DROT( I2-I, H( I-1, I+1 ), LDH, H( I, I+1 ), LDH, |
|
|
412 $ CS, SN ) |
|
|
413 CALL DROT( I-I1-1, H( I1, I-1 ), 1, H( I1, I ), 1, CS, SN ) |
|
|
414 END IF |
|
|
415 IF( WANTZ ) THEN |
|
|
416 * |
|
|
417 * Apply the transformation to Z. |
|
|
418 * |
|
|
419 CALL DROT( NZ, Z( ILOZ, I-1 ), 1, Z( ILOZ, I ), 1, CS, SN ) |
|
|
420 END IF |
|
|
421 END IF |
|
|
422 * |
|
|
423 * Decrement number of remaining iterations, and return to start of |
|
|
424 * the main loop with new value of I. |
|
|
425 * |
|
|
426 ITN = ITN - ITS |
|
|
427 I = L - 1 |
|
|
428 GO TO 10 |
|
|
429 * |
|
|
430 150 CONTINUE |
|
|
431 RETURN |
|
|
432 * |
|
|
433 * End of DLAHQR |
|
|
434 * |
|
|
435 END |
