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