|
2329
|
1 SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, |
|
|
2 $ ILOZ, IHIZ, Z, LDZ, INFO ) |
|
|
3 * |
|
7034
|
4 * -- LAPACK auxiliary routine (version 3.1) -- |
|
|
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. |
|
|
6 * November 2006 |
|
2329
|
7 * |
|
|
8 * .. Scalar Arguments .. |
|
7034
|
9 INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N |
|
2329
|
10 LOGICAL WANTT, WANTZ |
|
|
11 * .. |
|
|
12 * .. Array Arguments .. |
|
|
13 DOUBLE PRECISION H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * ) |
|
|
14 * .. |
|
|
15 * |
|
7034
|
16 * Purpose |
|
|
17 * ======= |
|
2329
|
18 * |
|
7034
|
19 * DLAHQR is an auxiliary routine called by DHSEQR to update the |
|
|
20 * eigenvalues and Schur decomposition already computed by DHSEQR, by |
|
|
21 * dealing with the Hessenberg submatrix in rows and columns ILO to |
|
|
22 * IHI. |
|
2329
|
23 * |
|
7034
|
24 * Arguments |
|
|
25 * ========= |
|
2329
|
26 * |
|
7034
|
27 * WANTT (input) LOGICAL |
|
2329
|
28 * = .TRUE. : the full Schur form T is required; |
|
|
29 * = .FALSE.: only eigenvalues are required. |
|
|
30 * |
|
7034
|
31 * WANTZ (input) LOGICAL |
|
2329
|
32 * = .TRUE. : the matrix of Schur vectors Z is required; |
|
|
33 * = .FALSE.: Schur vectors are not required. |
|
|
34 * |
|
7034
|
35 * N (input) INTEGER |
|
2329
|
36 * The order of the matrix H. N >= 0. |
|
|
37 * |
|
7034
|
38 * ILO (input) INTEGER |
|
|
39 * IHI (input) INTEGER |
|
2329
|
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 * |
|
7034
|
47 * H (input/output) DOUBLE PRECISION array, dimension (LDH,N) |
|
2329
|
48 * On entry, the upper Hessenberg matrix H. |
|
7034
|
49 * On exit, if INFO is zero and if WANTT is .TRUE., H is upper |
|
|
50 * quasi-triangular in rows and columns ILO:IHI, with any |
|
|
51 * 2-by-2 diagonal blocks in standard form. If INFO is zero |
|
|
52 * and WANTT is .FALSE., the contents of H are unspecified on |
|
|
53 * exit. The output state of H if INFO is nonzero is given |
|
|
54 * below under the description of INFO. |
|
2329
|
55 * |
|
7034
|
56 * LDH (input) INTEGER |
|
2329
|
57 * The leading dimension of the array H. LDH >= max(1,N). |
|
|
58 * |
|
7034
|
59 * WR (output) DOUBLE PRECISION array, dimension (N) |
|
|
60 * WI (output) DOUBLE PRECISION array, dimension (N) |
|
2329
|
61 * The real and imaginary parts, respectively, of the computed |
|
|
62 * eigenvalues ILO to IHI are stored in the corresponding |
|
|
63 * elements of WR and WI. If two eigenvalues are computed as a |
|
|
64 * complex conjugate pair, they are stored in consecutive |
|
|
65 * elements of WR and WI, say the i-th and (i+1)th, with |
|
|
66 * WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the |
|
|
67 * eigenvalues are stored in the same order as on the diagonal |
|
|
68 * of the Schur form returned in H, with WR(i) = H(i,i), and, if |
|
|
69 * H(i:i+1,i:i+1) is a 2-by-2 diagonal block, |
|
|
70 * WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i). |
|
|
71 * |
|
7034
|
72 * ILOZ (input) INTEGER |
|
|
73 * IHIZ (input) INTEGER |
|
2329
|
74 * Specify the rows of Z to which transformations must be |
|
|
75 * applied if WANTZ is .TRUE.. |
|
|
76 * 1 <= ILOZ <= ILO; IHI <= IHIZ <= N. |
|
|
77 * |
|
7034
|
78 * Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N) |
|
2329
|
79 * If WANTZ is .TRUE., on entry Z must contain the current |
|
|
80 * matrix Z of transformations accumulated by DHSEQR, and on |
|
|
81 * exit Z has been updated; transformations are applied only to |
|
|
82 * the submatrix Z(ILOZ:IHIZ,ILO:IHI). |
|
|
83 * If WANTZ is .FALSE., Z is not referenced. |
|
|
84 * |
|
7034
|
85 * LDZ (input) INTEGER |
|
2329
|
86 * The leading dimension of the array Z. LDZ >= max(1,N). |
|
|
87 * |
|
7034
|
88 * INFO (output) INTEGER |
|
|
89 * = 0: successful exit |
|
|
90 * .GT. 0: If INFO = i, DLAHQR failed to compute all the |
|
|
91 * eigenvalues ILO to IHI in a total of 30 iterations |
|
|
92 * per eigenvalue; elements i+1:ihi of WR and WI |
|
|
93 * contain those eigenvalues which have been |
|
|
94 * successfully computed. |
|
|
95 * |
|
|
96 * If INFO .GT. 0 and WANTT is .FALSE., then on exit, |
|
|
97 * the remaining unconverged eigenvalues are the |
|
|
98 * eigenvalues of the upper Hessenberg matrix rows |
|
|
99 * and columns ILO thorugh INFO of the final, output |
|
|
100 * value of H. |
|
2329
|
101 * |
|
7034
|
102 * If INFO .GT. 0 and WANTT is .TRUE., then on exit |
|
|
103 * (*) (initial value of H)*U = U*(final value of H) |
|
|
104 * where U is an orthognal matrix. The final |
|
|
105 * value of H is upper Hessenberg and triangular in |
|
|
106 * rows and columns INFO+1 through IHI. |
|
3333
|
107 * |
|
7034
|
108 * If INFO .GT. 0 and WANTZ is .TRUE., then on exit |
|
|
109 * (final value of Z) = (initial value of Z)*U |
|
|
110 * where U is the orthogonal matrix in (*) |
|
|
111 * (regardless of the value of WANTT.) |
|
|
112 * |
|
|
113 * Further Details |
|
|
114 * =============== |
|
|
115 * |
|
|
116 * 02-96 Based on modifications by |
|
3333
|
117 * David Day, Sandia National Laboratory, USA |
|
|
118 * |
|
7034
|
119 * 12-04 Further modifications by |
|
|
120 * Ralph Byers, University of Kansas, USA |
|
|
121 * |
|
|
122 * This is a modified version of DLAHQR from LAPACK version 3.0. |
|
|
123 * It is (1) more robust against overflow and underflow and |
|
|
124 * (2) adopts the more conservative Ahues & Tisseur stopping |
|
|
125 * criterion (LAWN 122, 1997). |
|
|
126 * |
|
|
127 * ========================================================= |
|
2329
|
128 * |
|
|
129 * .. Parameters .. |
|
7034
|
130 INTEGER ITMAX |
|
|
131 PARAMETER ( ITMAX = 30 ) |
|
|
132 DOUBLE PRECISION ZERO, ONE, TWO |
|
|
133 PARAMETER ( ZERO = 0.0d0, ONE = 1.0d0, TWO = 2.0d0 ) |
|
2329
|
134 DOUBLE PRECISION DAT1, DAT2 |
|
7034
|
135 PARAMETER ( DAT1 = 3.0d0 / 4.0d0, DAT2 = -0.4375d0 ) |
|
2329
|
136 * .. |
|
|
137 * .. Local Scalars .. |
|
7034
|
138 DOUBLE PRECISION AA, AB, BA, BB, CS, DET, H11, H12, H21, H21S, |
|
|
139 $ H22, RT1I, RT1R, RT2I, RT2R, RTDISC, S, SAFMAX, |
|
|
140 $ SAFMIN, SMLNUM, SN, SUM, T1, T2, T3, TR, TST, |
|
|
141 $ ULP, V2, V3 |
|
|
142 INTEGER I, I1, I2, ITS, J, K, L, M, NH, NR, NZ |
|
2329
|
143 * .. |
|
|
144 * .. Local Arrays .. |
|
7034
|
145 DOUBLE PRECISION V( 3 ) |
|
2329
|
146 * .. |
|
|
147 * .. External Functions .. |
|
7034
|
148 DOUBLE PRECISION DLAMCH |
|
|
149 EXTERNAL DLAMCH |
|
2329
|
150 * .. |
|
|
151 * .. External Subroutines .. |
|
7034
|
152 EXTERNAL DCOPY, DLABAD, DLANV2, DLARFG, DROT |
|
2329
|
153 * .. |
|
|
154 * .. Intrinsic Functions .. |
|
7034
|
155 INTRINSIC ABS, DBLE, MAX, MIN, SQRT |
|
2329
|
156 * .. |
|
|
157 * .. Executable Statements .. |
|
|
158 * |
|
|
159 INFO = 0 |
|
|
160 * |
|
|
161 * Quick return if possible |
|
|
162 * |
|
|
163 IF( N.EQ.0 ) |
|
|
164 $ RETURN |
|
|
165 IF( ILO.EQ.IHI ) THEN |
|
|
166 WR( ILO ) = H( ILO, ILO ) |
|
|
167 WI( ILO ) = ZERO |
|
|
168 RETURN |
|
|
169 END IF |
|
|
170 * |
|
7034
|
171 * ==== clear out the trash ==== |
|
|
172 DO 10 J = ILO, IHI - 3 |
|
|
173 H( J+2, J ) = ZERO |
|
|
174 H( J+3, J ) = ZERO |
|
|
175 10 CONTINUE |
|
|
176 IF( ILO.LE.IHI-2 ) |
|
|
177 $ H( IHI, IHI-2 ) = ZERO |
|
|
178 * |
|
2329
|
179 NH = IHI - ILO + 1 |
|
|
180 NZ = IHIZ - ILOZ + 1 |
|
|
181 * |
|
|
182 * Set machine-dependent constants for the stopping criterion. |
|
|
183 * |
|
7034
|
184 SAFMIN = DLAMCH( 'SAFE MINIMUM' ) |
|
|
185 SAFMAX = ONE / SAFMIN |
|
|
186 CALL DLABAD( SAFMIN, SAFMAX ) |
|
|
187 ULP = DLAMCH( 'PRECISION' ) |
|
|
188 SMLNUM = SAFMIN*( DBLE( NH ) / ULP ) |
|
2329
|
189 * |
|
|
190 * I1 and I2 are the indices of the first row and last column of H |
|
|
191 * to which transformations must be applied. If eigenvalues only are |
|
|
192 * being computed, I1 and I2 are set inside the main loop. |
|
|
193 * |
|
|
194 IF( WANTT ) THEN |
|
|
195 I1 = 1 |
|
|
196 I2 = N |
|
|
197 END IF |
|
|
198 * |
|
|
199 * The main loop begins here. I is the loop index and decreases from |
|
|
200 * IHI to ILO in steps of 1 or 2. Each iteration of the loop works |
|
|
201 * with the active submatrix in rows and columns L to I. |
|
|
202 * Eigenvalues I+1 to IHI have already converged. Either L = ILO or |
|
|
203 * H(L,L-1) is negligible so that the matrix splits. |
|
|
204 * |
|
|
205 I = IHI |
|
7034
|
206 20 CONTINUE |
|
2329
|
207 L = ILO |
|
|
208 IF( I.LT.ILO ) |
|
7034
|
209 $ GO TO 160 |
|
2329
|
210 * |
|
|
211 * Perform QR iterations on rows and columns ILO to I until a |
|
|
212 * submatrix of order 1 or 2 splits off at the bottom because a |
|
|
213 * subdiagonal element has become negligible. |
|
|
214 * |
|
7034
|
215 DO 140 ITS = 0, ITMAX |
|
2329
|
216 * |
|
|
217 * Look for a single small subdiagonal element. |
|
|
218 * |
|
7034
|
219 DO 30 K = I, L + 1, -1 |
|
|
220 IF( ABS( H( K, K-1 ) ).LE.SMLNUM ) |
|
|
221 $ GO TO 40 |
|
|
222 TST = ABS( H( K-1, K-1 ) ) + ABS( H( K, K ) ) |
|
|
223 IF( TST.EQ.ZERO ) THEN |
|
|
224 IF( K-2.GE.ILO ) |
|
|
225 $ TST = TST + ABS( H( K-1, K-2 ) ) |
|
|
226 IF( K+1.LE.IHI ) |
|
|
227 $ TST = TST + ABS( H( K+1, K ) ) |
|
|
228 END IF |
|
|
229 * ==== The following is a conservative small subdiagonal |
|
|
230 * . deflation criterion due to Ahues & Tisseur (LAWN 122, |
|
|
231 * . 1997). It has better mathematical foundation and |
|
|
232 * . improves accuracy in some cases. ==== |
|
|
233 IF( ABS( H( K, K-1 ) ).LE.ULP*TST ) THEN |
|
|
234 AB = MAX( ABS( H( K, K-1 ) ), ABS( H( K-1, K ) ) ) |
|
|
235 BA = MIN( ABS( H( K, K-1 ) ), ABS( H( K-1, K ) ) ) |
|
|
236 AA = MAX( ABS( H( K, K ) ), |
|
|
237 $ ABS( H( K-1, K-1 )-H( K, K ) ) ) |
|
|
238 BB = MIN( ABS( H( K, K ) ), |
|
|
239 $ ABS( H( K-1, K-1 )-H( K, K ) ) ) |
|
|
240 S = AA + AB |
|
|
241 IF( BA*( AB / S ).LE.MAX( SMLNUM, |
|
|
242 $ ULP*( BB*( AA / S ) ) ) )GO TO 40 |
|
|
243 END IF |
|
2329
|
244 30 CONTINUE |
|
7034
|
245 40 CONTINUE |
|
2329
|
246 L = K |
|
|
247 IF( L.GT.ILO ) THEN |
|
|
248 * |
|
|
249 * H(L,L-1) is negligible |
|
|
250 * |
|
|
251 H( L, L-1 ) = ZERO |
|
|
252 END IF |
|
|
253 * |
|
|
254 * Exit from loop if a submatrix of order 1 or 2 has split off. |
|
|
255 * |
|
|
256 IF( L.GE.I-1 ) |
|
7034
|
257 $ GO TO 150 |
|
2329
|
258 * |
|
|
259 * Now the active submatrix is in rows and columns L to I. If |
|
|
260 * eigenvalues only are being computed, only the active submatrix |
|
|
261 * need be transformed. |
|
|
262 * |
|
|
263 IF( .NOT.WANTT ) THEN |
|
|
264 I1 = L |
|
|
265 I2 = I |
|
|
266 END IF |
|
|
267 * |
|
|
268 IF( ITS.EQ.10 .OR. ITS.EQ.20 ) THEN |
|
|
269 * |
|
|
270 * Exceptional shift. |
|
|
271 * |
|
7034
|
272 H11 = DAT1*S + H( I, I ) |
|
|
273 H12 = DAT2*S |
|
|
274 H21 = S |
|
|
275 H22 = H11 |
|
2329
|
276 ELSE |
|
|
277 * |
|
3333
|
278 * Prepare to use Francis' double shift |
|
|
279 * (i.e. 2nd degree generalized Rayleigh quotient) |
|
2329
|
280 * |
|
7034
|
281 H11 = H( I-1, I-1 ) |
|
|
282 H21 = H( I, I-1 ) |
|
|
283 H12 = H( I-1, I ) |
|
|
284 H22 = H( I, I ) |
|
|
285 END IF |
|
|
286 S = ABS( H11 ) + ABS( H12 ) + ABS( H21 ) + ABS( H22 ) |
|
|
287 IF( S.EQ.ZERO ) THEN |
|
|
288 RT1R = ZERO |
|
|
289 RT1I = ZERO |
|
|
290 RT2R = ZERO |
|
|
291 RT2I = ZERO |
|
|
292 ELSE |
|
|
293 H11 = H11 / S |
|
|
294 H21 = H21 / S |
|
|
295 H12 = H12 / S |
|
|
296 H22 = H22 / S |
|
|
297 TR = ( H11+H22 ) / TWO |
|
|
298 DET = ( H11-TR )*( H22-TR ) - H12*H21 |
|
|
299 RTDISC = SQRT( ABS( DET ) ) |
|
|
300 IF( DET.GE.ZERO ) THEN |
|
3333
|
301 * |
|
7034
|
302 * ==== complex conjugate shifts ==== |
|
|
303 * |
|
|
304 RT1R = TR*S |
|
|
305 RT2R = RT1R |
|
|
306 RT1I = RTDISC*S |
|
|
307 RT2I = -RT1I |
|
|
308 ELSE |
|
|
309 * |
|
|
310 * ==== real shifts (use only one of them) ==== |
|
|
311 * |
|
|
312 RT1R = TR + RTDISC |
|
|
313 RT2R = TR - RTDISC |
|
|
314 IF( ABS( RT1R-H22 ).LE.ABS( RT2R-H22 ) ) THEN |
|
|
315 RT1R = RT1R*S |
|
|
316 RT2R = RT1R |
|
3333
|
317 ELSE |
|
7034
|
318 RT2R = RT2R*S |
|
|
319 RT1R = RT2R |
|
3333
|
320 END IF |
|
7034
|
321 RT1I = ZERO |
|
|
322 RT2I = ZERO |
|
3333
|
323 END IF |
|
2329
|
324 END IF |
|
|
325 * |
|
|
326 * Look for two consecutive small subdiagonal elements. |
|
|
327 * |
|
7034
|
328 DO 50 M = I - 2, L, -1 |
|
2329
|
329 * Determine the effect of starting the double-shift QR |
|
|
330 * iteration at row M, and see if this would make H(M,M-1) |
|
7034
|
331 * negligible. (The following uses scaling to avoid |
|
|
332 * overflows and most underflows.) |
|
2329
|
333 * |
|
7034
|
334 H21S = H( M+1, M ) |
|
|
335 S = ABS( H( M, M )-RT2R ) + ABS( RT2I ) + ABS( H21S ) |
|
|
336 H21S = H( M+1, M ) / S |
|
|
337 V( 1 ) = H21S*H( M, M+1 ) + ( H( M, M )-RT1R )* |
|
|
338 $ ( ( H( M, M )-RT2R ) / S ) - RT1I*( RT2I / S ) |
|
|
339 V( 2 ) = H21S*( H( M, M )+H( M+1, M+1 )-RT1R-RT2R ) |
|
|
340 V( 3 ) = H21S*H( M+2, M+1 ) |
|
|
341 S = ABS( V( 1 ) ) + ABS( V( 2 ) ) + ABS( V( 3 ) ) |
|
|
342 V( 1 ) = V( 1 ) / S |
|
|
343 V( 2 ) = V( 2 ) / S |
|
|
344 V( 3 ) = V( 3 ) / S |
|
2329
|
345 IF( M.EQ.L ) |
|
7034
|
346 $ GO TO 60 |
|
|
347 IF( ABS( H( M, M-1 ) )*( ABS( V( 2 ) )+ABS( V( 3 ) ) ).LE. |
|
|
348 $ ULP*ABS( V( 1 ) )*( ABS( H( M-1, M-1 ) )+ABS( H( M, |
|
|
349 $ M ) )+ABS( H( M+1, M+1 ) ) ) )GO TO 60 |
|
2329
|
350 50 CONTINUE |
|
7034
|
351 60 CONTINUE |
|
2329
|
352 * |
|
|
353 * Double-shift QR step |
|
|
354 * |
|
7034
|
355 DO 130 K = M, I - 1 |
|
2329
|
356 * |
|
|
357 * The first iteration of this loop determines a reflection G |
|
|
358 * from the vector V and applies it from left and right to H, |
|
|
359 * thus creating a nonzero bulge below the subdiagonal. |
|
|
360 * |
|
|
361 * Each subsequent iteration determines a reflection G to |
|
|
362 * restore the Hessenberg form in the (K-1)th column, and thus |
|
|
363 * chases the bulge one step toward the bottom of the active |
|
|
364 * submatrix. NR is the order of G. |
|
|
365 * |
|
|
366 NR = MIN( 3, I-K+1 ) |
|
|
367 IF( K.GT.M ) |
|
|
368 $ CALL DCOPY( NR, H( K, K-1 ), 1, V, 1 ) |
|
|
369 CALL DLARFG( NR, V( 1 ), V( 2 ), 1, T1 ) |
|
|
370 IF( K.GT.M ) THEN |
|
|
371 H( K, K-1 ) = V( 1 ) |
|
|
372 H( K+1, K-1 ) = ZERO |
|
|
373 IF( K.LT.I-1 ) |
|
|
374 $ H( K+2, K-1 ) = ZERO |
|
|
375 ELSE IF( M.GT.L ) THEN |
|
|
376 H( K, K-1 ) = -H( K, K-1 ) |
|
|
377 END IF |
|
|
378 V2 = V( 2 ) |
|
|
379 T2 = T1*V2 |
|
|
380 IF( NR.EQ.3 ) THEN |
|
|
381 V3 = V( 3 ) |
|
|
382 T3 = T1*V3 |
|
|
383 * |
|
|
384 * Apply G from the left to transform the rows of the matrix |
|
|
385 * in columns K to I2. |
|
|
386 * |
|
7034
|
387 DO 70 J = K, I2 |
|
2329
|
388 SUM = H( K, J ) + V2*H( K+1, J ) + V3*H( K+2, J ) |
|
|
389 H( K, J ) = H( K, J ) - SUM*T1 |
|
|
390 H( K+1, J ) = H( K+1, J ) - SUM*T2 |
|
|
391 H( K+2, J ) = H( K+2, J ) - SUM*T3 |
|
7034
|
392 70 CONTINUE |
|
2329
|
393 * |
|
|
394 * Apply G from the right to transform the columns of the |
|
|
395 * matrix in rows I1 to min(K+3,I). |
|
|
396 * |
|
7034
|
397 DO 80 J = I1, MIN( K+3, I ) |
|
2329
|
398 SUM = H( J, K ) + V2*H( J, K+1 ) + V3*H( J, K+2 ) |
|
|
399 H( J, K ) = H( J, K ) - SUM*T1 |
|
|
400 H( J, K+1 ) = H( J, K+1 ) - SUM*T2 |
|
|
401 H( J, K+2 ) = H( J, K+2 ) - SUM*T3 |
|
7034
|
402 80 CONTINUE |
|
2329
|
403 * |
|
|
404 IF( WANTZ ) THEN |
|
|
405 * |
|
|
406 * Accumulate transformations in the matrix Z |
|
|
407 * |
|
7034
|
408 DO 90 J = ILOZ, IHIZ |
|
2329
|
409 SUM = Z( J, K ) + V2*Z( J, K+1 ) + V3*Z( J, K+2 ) |
|
|
410 Z( J, K ) = Z( J, K ) - SUM*T1 |
|
|
411 Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2 |
|
|
412 Z( J, K+2 ) = Z( J, K+2 ) - SUM*T3 |
|
7034
|
413 90 CONTINUE |
|
2329
|
414 END IF |
|
|
415 ELSE IF( NR.EQ.2 ) THEN |
|
|
416 * |
|
|
417 * Apply G from the left to transform the rows of the matrix |
|
|
418 * in columns K to I2. |
|
|
419 * |
|
7034
|
420 DO 100 J = K, I2 |
|
2329
|
421 SUM = H( K, J ) + V2*H( K+1, J ) |
|
|
422 H( K, J ) = H( K, J ) - SUM*T1 |
|
|
423 H( K+1, J ) = H( K+1, J ) - SUM*T2 |
|
7034
|
424 100 CONTINUE |
|
2329
|
425 * |
|
|
426 * Apply G from the right to transform the columns of the |
|
|
427 * matrix in rows I1 to min(K+3,I). |
|
|
428 * |
|
7034
|
429 DO 110 J = I1, I |
|
2329
|
430 SUM = H( J, K ) + V2*H( J, K+1 ) |
|
|
431 H( J, K ) = H( J, K ) - SUM*T1 |
|
|
432 H( J, K+1 ) = H( J, K+1 ) - SUM*T2 |
|
7034
|
433 110 CONTINUE |
|
2329
|
434 * |
|
|
435 IF( WANTZ ) THEN |
|
|
436 * |
|
|
437 * Accumulate transformations in the matrix Z |
|
|
438 * |
|
7034
|
439 DO 120 J = ILOZ, IHIZ |
|
2329
|
440 SUM = Z( J, K ) + V2*Z( J, K+1 ) |
|
|
441 Z( J, K ) = Z( J, K ) - SUM*T1 |
|
|
442 Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2 |
|
7034
|
443 120 CONTINUE |
|
2329
|
444 END IF |
|
|
445 END IF |
|
7034
|
446 130 CONTINUE |
|
2329
|
447 * |
|
7034
|
448 140 CONTINUE |
|
2329
|
449 * |
|
|
450 * Failure to converge in remaining number of iterations |
|
|
451 * |
|
|
452 INFO = I |
|
|
453 RETURN |
|
|
454 * |
|
7034
|
455 150 CONTINUE |
|
2329
|
456 * |
|
|
457 IF( L.EQ.I ) THEN |
|
|
458 * |
|
|
459 * H(I,I-1) is negligible: one eigenvalue has converged. |
|
|
460 * |
|
|
461 WR( I ) = H( I, I ) |
|
|
462 WI( I ) = ZERO |
|
|
463 ELSE IF( L.EQ.I-1 ) THEN |
|
|
464 * |
|
|
465 * H(I-1,I-2) is negligible: a pair of eigenvalues have converged. |
|
|
466 * |
|
|
467 * Transform the 2-by-2 submatrix to standard Schur form, |
|
|
468 * and compute and store the eigenvalues. |
|
|
469 * |
|
|
470 CALL DLANV2( H( I-1, I-1 ), H( I-1, I ), H( I, I-1 ), |
|
|
471 $ H( I, I ), WR( I-1 ), WI( I-1 ), WR( I ), WI( I ), |
|
|
472 $ CS, SN ) |
|
|
473 * |
|
|
474 IF( WANTT ) THEN |
|
|
475 * |
|
|
476 * Apply the transformation to the rest of H. |
|
|
477 * |
|
|
478 IF( I2.GT.I ) |
|
|
479 $ CALL DROT( I2-I, H( I-1, I+1 ), LDH, H( I, I+1 ), LDH, |
|
|
480 $ CS, SN ) |
|
|
481 CALL DROT( I-I1-1, H( I1, I-1 ), 1, H( I1, I ), 1, CS, SN ) |
|
|
482 END IF |
|
|
483 IF( WANTZ ) THEN |
|
|
484 * |
|
|
485 * Apply the transformation to Z. |
|
|
486 * |
|
|
487 CALL DROT( NZ, Z( ILOZ, I-1 ), 1, Z( ILOZ, I ), 1, CS, SN ) |
|
|
488 END IF |
|
|
489 END IF |
|
|
490 * |
|
7034
|
491 * return to start of the main loop with new value of I. |
|
2329
|
492 * |
|
|
493 I = L - 1 |
|
7034
|
494 GO TO 20 |
|
2329
|
495 * |
|
7034
|
496 160 CONTINUE |
|
2329
|
497 RETURN |
|
|
498 * |
|
|
499 * End of DLAHQR |
|
|
500 * |
|
|
501 END |
