1 SUBROUTINE ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
2 $ 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 COMPLEX*16 H( LDH, * ), W( * ), Z( LDZ, * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * ZLAHQR is an auxiliary routine called by CHSEQR to update the
20 * eigenvalues and Schur decomposition already computed by CHSEQR, 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 triangular in rows and
41 * columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1).
42 * ZLAHQR works primarily with the Hessenberg submatrix in rows
43 * and columns ILO to IHI, but applies transformations to all of
44 * H if WANTT is .TRUE..
45 * 1 <= ILO <= max(1,IHI); IHI <= N.
46 *
47 * H (input/output) COMPLEX*16 array, dimension (LDH,N)
48 * On entry, the upper Hessenberg matrix H.
49 * On exit, if INFO is zero and if WANTT is .TRUE., then H
50 * is upper triangular in rows and columns ILO:IHI. If INFO
51 * is zero and if WANTT is .FALSE., then the contents of H
52 * are unspecified on exit. The output state of H in case
53 * INF is positive is below under the description of INFO.
54 *
55 * LDH (input) INTEGER
56 * The leading dimension of the array H. LDH >= max(1,N).
57 *
58 * W (output) COMPLEX*16 array, dimension (N)
59 * The computed eigenvalues ILO to IHI are stored in the
60 * corresponding elements of W. If WANTT is .TRUE., the
61 * eigenvalues are stored in the same order as on the diagonal
62 * of the Schur form returned in H, with W(i) = H(i,i).
63 *
64 * ILOZ (input) INTEGER
65 * IHIZ (input) INTEGER
66 * Specify the rows of Z to which transformations must be
67 * applied if WANTZ is .TRUE..
68 * 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
69 *
70 * Z (input/output) COMPLEX*16 array, dimension (LDZ,N)
71 * If WANTZ is .TRUE., on entry Z must contain the current
72 * matrix Z of transformations accumulated by CHSEQR, and on
73 * exit Z has been updated; transformations are applied only to
74 * the submatrix Z(ILOZ:IHIZ,ILO:IHI).
75 * If WANTZ is .FALSE., Z is not referenced.
76 *
77 * LDZ (input) INTEGER
78 * The leading dimension of the array Z. LDZ >= max(1,N).
79 *
80 * INFO (output) INTEGER
81 * = 0: successful exit
82 * .GT. 0: if INFO = i, ZLAHQR failed to compute all the
83 * eigenvalues ILO to IHI in a total of 30 iterations
84 * per eigenvalue; elements i+1:ihi of W contain
85 * those eigenvalues which have been successfully
86 * computed.
87 *
88 * If INFO .GT. 0 and WANTT is .FALSE., then on exit,
89 * the remaining unconverged eigenvalues are the
90 * eigenvalues of the upper Hessenberg matrix
91 * rows and columns ILO thorugh INFO of the final,
92 * output value of H.
93 *
94 * If INFO .GT. 0 and WANTT is .TRUE., then on exit
95 * (*) (initial value of H)*U = U*(final value of H)
96 * where U is an orthognal matrix. The final
97 * value of H is upper Hessenberg and triangular in
98 * rows and columns INFO+1 through IHI.
99 *
100 * If INFO .GT. 0 and WANTZ is .TRUE., then on exit
101 * (final value of Z) = (initial value of Z)*U
102 * where U is the orthogonal matrix in (*)
103 * (regardless of the value of WANTT.)
104 *
105 * Further Details
106 * ===============
107 *
108 * 02-96 Based on modifications by
109 * David Day, Sandia National Laboratory, USA
110 *
111 * 12-04 Further modifications by
112 * Ralph Byers, University of Kansas, USA
113 * This is a modified version of ZLAHQR from LAPACK version 3.0.
114 * It is (1) more robust against overflow and underflow and
115 * (2) adopts the more conservative Ahues & Tisseur stopping
116 * criterion (LAWN 122, 1997).
117 *
118 * =========================================================
119 *
120 * .. Parameters ..
121 INTEGER ITMAX
122 PARAMETER ( ITMAX = 30 )
123 COMPLEX*16 ZERO, ONE
124 PARAMETER ( ZERO = ( 0.0d0, 0.0d0 ),
125 $ ONE = ( 1.0d0, 0.0d0 ) )
126 DOUBLE PRECISION RZERO, RONE, HALF
127 PARAMETER ( RZERO = 0.0d0, RONE = 1.0d0, HALF = 0.5d0 )
128 DOUBLE PRECISION DAT1
129 PARAMETER ( DAT1 = 3.0d0 / 4.0d0 )
130 * ..
131 * .. Local Scalars ..
132 COMPLEX*16 CDUM, H11, H11S, H22, SC, SUM, T, T1, TEMP, U,
133 $ V2, X, Y
134 DOUBLE PRECISION AA, AB, BA, BB, H10, H21, RTEMP, S, SAFMAX,
135 $ SAFMIN, SMLNUM, SX, T2, TST, ULP
136 INTEGER I, I1, I2, ITS, J, JHI, JLO, K, L, M, NH, NZ
137 * ..
138 * .. Local Arrays ..
139 COMPLEX*16 V( 2 )
140 * ..
141 * .. External Functions ..
142 COMPLEX*16 ZLADIV
143 DOUBLE PRECISION DLAMCH
144 EXTERNAL ZLADIV, DLAMCH
145 * ..
146 * .. External Subroutines ..
147 EXTERNAL DLABAD, ZCOPY, ZLARFG, ZSCAL
148 * ..
149 * .. Statement Functions ..
150 DOUBLE PRECISION CABS1
151 * ..
152 * .. Intrinsic Functions ..
153 INTRINSIC ABS, DBLE, DCONJG, DIMAG, MAX, MIN, SQRT
154 * ..
155 * .. Statement Function definitions ..
156 CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
157 * ..
158 * .. Executable Statements ..
159 *
160 INFO = 0
161 *
162 * Quick return if possible
163 *
164 IF( N.EQ.0 )
165 $ RETURN
166 IF( ILO.EQ.IHI ) THEN
167 W( ILO ) = H( ILO, ILO )
168 RETURN
169 END IF
170 *
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 * ==== ensure that subdiagonal entries are real ====
179 IF( WANTT ) THEN
180 JLO = 1
181 JHI = N
182 ELSE
183 JLO = ILO
184 JHI = IHI
185 END IF
186 DO 20 I = ILO + 1, IHI
187 IF( DIMAG( H( I, I-1 ) ).NE.RZERO ) THEN
188 * ==== The following redundant normalization
189 * . avoids problems with both gradual and
190 * . sudden underflow in ABS(H(I,I-1)) ====
191 SC = H( I, I-1 ) / CABS1( H( I, I-1 ) )
192 SC = DCONJG( SC ) / ABS( SC )
193 H( I, I-1 ) = ABS( H( I, I-1 ) )
194 CALL ZSCAL( JHI-I+1, SC, H( I, I ), LDH )
195 CALL ZSCAL( MIN( JHI, I+1 )-JLO+1, DCONJG( SC ),
196 $ H( JLO, I ), 1 )
197 IF( WANTZ )
198 $ CALL ZSCAL( IHIZ-ILOZ+1, DCONJG( SC ), Z( ILOZ, I ), 1 )
199 END IF
200 20 CONTINUE
201 *
202 NH = IHI - ILO + 1
203 NZ = IHIZ - ILOZ + 1
204 *
205 * Set machine-dependent constants for the stopping criterion.
206 *
207 SAFMIN = DLAMCH( 'SAFE MINIMUM' )
208 SAFMAX = RONE / SAFMIN
209 CALL DLABAD( SAFMIN, SAFMAX )
210 ULP = DLAMCH( 'PRECISION' )
211 SMLNUM = SAFMIN*( DBLE( NH ) / ULP )
212 *
213 * I1 and I2 are the indices of the first row and last column of H
214 * to which transformations must be applied. If eigenvalues only are
215 * being computed, I1 and I2 are set inside the main loop.
216 *
217 IF( WANTT ) THEN
218 I1 = 1
219 I2 = N
220 END IF
221 *
222 * The main loop begins here. I is the loop index and decreases from
223 * IHI to ILO in steps of 1. Each iteration of the loop works
224 * with the active submatrix in rows and columns L to I.
225 * Eigenvalues I+1 to IHI have already converged. Either L = ILO, or
226 * H(L,L-1) is negligible so that the matrix splits.
227 *
228 I = IHI
229 30 CONTINUE
230 IF( I.LT.ILO )
231 $ GO TO 150
232 *
233 * Perform QR iterations on rows and columns ILO to I until a
234 * submatrix of order 1 splits off at the bottom because a
235 * subdiagonal element has become negligible.
236 *
237 L = ILO
238 DO 130 ITS = 0, ITMAX
239 *
240 * Look for a single small subdiagonal element.
241 *
242 DO 40 K = I, L + 1, -1
243 IF( CABS1( H( K, K-1 ) ).LE.SMLNUM )
244 $ GO TO 50
245 TST = CABS1( H( K-1, K-1 ) ) + CABS1( H( K, K ) )
246 IF( TST.EQ.ZERO ) THEN
247 IF( K-2.GE.ILO )
248 $ TST = TST + ABS( DBLE( H( K-1, K-2 ) ) )
249 IF( K+1.LE.IHI )
250 $ TST = TST + ABS( DBLE( H( K+1, K ) ) )
251 END IF
252 * ==== The following is a conservative small subdiagonal
253 * . deflation criterion due to Ahues & Tisseur (LAWN 122,
254 * . 1997). It has better mathematical foundation and
255 * . improves accuracy in some examples. ====
256 IF( ABS( DBLE( H( K, K-1 ) ) ).LE.ULP*TST ) THEN
257 AB = MAX( CABS1( H( K, K-1 ) ), CABS1( H( K-1, K ) ) )
258 BA = MIN( CABS1( H( K, K-1 ) ), CABS1( H( K-1, K ) ) )
259 AA = MAX( CABS1( H( K, K ) ),
260 $ CABS1( H( K-1, K-1 )-H( K, K ) ) )
261 BB = MIN( CABS1( H( K, K ) ),
262 $ CABS1( H( K-1, K-1 )-H( K, K ) ) )
263 S = AA + AB
264 IF( BA*( AB / S ).LE.MAX( SMLNUM,
265 $ ULP*( BB*( AA / S ) ) ) )GO TO 50
266 END IF
267 40 CONTINUE
268 50 CONTINUE
269 L = K
270 IF( L.GT.ILO ) THEN
271 *
272 * H(L,L-1) is negligible
273 *
274 H( L, L-1 ) = ZERO
275 END IF
276 *
277 * Exit from loop if a submatrix of order 1 has split off.
278 *
279 IF( L.GE.I )
280 $ GO TO 140
281 *
282 * Now the active submatrix is in rows and columns L to I. If
283 * eigenvalues only are being computed, only the active submatrix
284 * need be transformed.
285 *
286 IF( .NOT.WANTT ) THEN
287 I1 = L
288 I2 = I
289 END IF
290 *
291 IF( ITS.EQ.10 ) THEN
292 *
293 * Exceptional shift.
294 *
295 S = DAT1*ABS( DBLE( H( L+1, L ) ) )
296 T = S + H( L, L )
297 ELSE IF( ITS.EQ.20 ) THEN
298 *
299 * Exceptional shift.
300 *
301 S = DAT1*ABS( DBLE( H( I, I-1 ) ) )
302 T = S + H( I, I )
303 ELSE
304 *
305 * Wilkinson's shift.
306 *
307 T = H( I, I )
308 U = SQRT( H( I-1, I ) )*SQRT( H( I, I-1 ) )
309 S = CABS1( U )
310 IF( S.NE.RZERO ) THEN
311 X = HALF*( H( I-1, I-1 )-T )
312 SX = CABS1( X )
313 S = MAX( S, CABS1( X ) )
314 Y = S*SQRT( ( X / S )**2+( U / S )**2 )
315 IF( SX.GT.RZERO ) THEN
316 IF( DBLE( X / SX )*DBLE( Y )+DIMAG( X / SX )*
317 $ DIMAG( Y ).LT.RZERO )Y = -Y
318 END IF
319 T = T - U*ZLADIV( U, ( X+Y ) )
320 END IF
321 END IF
322 *
323 * Look for two consecutive small subdiagonal elements.
324 *
325 DO 60 M = I - 1, L + 1, -1
326 *
327 * Determine the effect of starting the single-shift QR
328 * iteration at row M, and see if this would make H(M,M-1)
329 * negligible.
330 *
331 H11 = H( M, M )
332 H22 = H( M+1, M+1 )
333 H11S = H11 - T
334 H21 = DBLE( H( M+1, M ) )
335 S = CABS1( H11S ) + ABS( H21 )
336 H11S = H11S / S
337 H21 = H21 / S
338 V( 1 ) = H11S
339 V( 2 ) = H21
340 H10 = DBLE( H( M, M-1 ) )
341 IF( ABS( H10 )*ABS( H21 ).LE.ULP*
342 $ ( CABS1( H11S )*( CABS1( H11 )+CABS1( H22 ) ) ) )
343 $ GO TO 70
344 60 CONTINUE
345 H11 = H( L, L )
346 H22 = H( L+1, L+1 )
347 H11S = H11 - T
348 H21 = DBLE( H( L+1, L ) )
349 S = CABS1( H11S ) + ABS( H21 )
350 H11S = H11S / S
351 H21 = H21 / S
352 V( 1 ) = H11S
353 V( 2 ) = H21
354 70 CONTINUE
355 *
356 * Single-shift QR step
357 *
358 DO 120 K = M, I - 1
359 *
360 * The first iteration of this loop determines a reflection G
361 * from the vector V and applies it from left and right to H,
362 * thus creating a nonzero bulge below the subdiagonal.
363 *
364 * Each subsequent iteration determines a reflection G to
365 * restore the Hessenberg form in the (K-1)th column, and thus
366 * chases the bulge one step toward the bottom of the active
367 * submatrix.
368 *
369 * V(2) is always real before the call to ZLARFG, and hence
370 * after the call T2 ( = T1*V(2) ) is also real.
371 *
372 IF( K.GT.M )
373 $ CALL ZCOPY( 2, H( K, K-1 ), 1, V, 1 )
374 CALL ZLARFG( 2, V( 1 ), V( 2 ), 1, T1 )
375 IF( K.GT.M ) THEN
376 H( K, K-1 ) = V( 1 )
377 H( K+1, K-1 ) = ZERO
378 END IF
379 V2 = V( 2 )
380 T2 = DBLE( T1*V2 )
381 *
382 * Apply G from the left to transform the rows of the matrix
383 * in columns K to I2.
384 *
385 DO 80 J = K, I2
386 SUM = DCONJG( T1 )*H( K, J ) + T2*H( K+1, J )
387 H( K, J ) = H( K, J ) - SUM
388 H( K+1, J ) = H( K+1, J ) - SUM*V2
389 80 CONTINUE
390 *
391 * Apply G from the right to transform the columns of the
392 * matrix in rows I1 to min(K+2,I).
393 *
394 DO 90 J = I1, MIN( K+2, I )
395 SUM = T1*H( J, K ) + T2*H( J, K+1 )
396 H( J, K ) = H( J, K ) - SUM
397 H( J, K+1 ) = H( J, K+1 ) - SUM*DCONJG( V2 )
398 90 CONTINUE
399 *
400 IF( WANTZ ) THEN
401 *
402 * Accumulate transformations in the matrix Z
403 *
404 DO 100 J = ILOZ, IHIZ
405 SUM = T1*Z( J, K ) + T2*Z( J, K+1 )
406 Z( J, K ) = Z( J, K ) - SUM
407 Z( J, K+1 ) = Z( J, K+1 ) - SUM*DCONJG( V2 )
408 100 CONTINUE
409 END IF
410 *
411 IF( K.EQ.M .AND. M.GT.L ) THEN
412 *
413 * If the QR step was started at row M > L because two
414 * consecutive small subdiagonals were found, then extra
415 * scaling must be performed to ensure that H(M,M-1) remains
416 * real.
417 *
418 TEMP = ONE - T1
419 TEMP = TEMP / ABS( TEMP )
420 H( M+1, M ) = H( M+1, M )*DCONJG( TEMP )
421 IF( M+2.LE.I )
422 $ H( M+2, M+1 ) = H( M+2, M+1 )*TEMP
423 DO 110 J = M, I
424 IF( J.NE.M+1 ) THEN
425 IF( I2.GT.J )
426 $ CALL ZSCAL( I2-J, TEMP, H( J, J+1 ), LDH )
427 CALL ZSCAL( J-I1, DCONJG( TEMP ), H( I1, J ), 1 )
428 IF( WANTZ ) THEN
429 CALL ZSCAL( NZ, DCONJG( TEMP ), Z( ILOZ, J ),
430 $ 1 )
431 END IF
432 END IF
433 110 CONTINUE
434 END IF
435 120 CONTINUE
436 *
437 * Ensure that H(I,I-1) is real.
438 *
439 TEMP = H( I, I-1 )
440 IF( DIMAG( TEMP ).NE.RZERO ) THEN
441 RTEMP = ABS( TEMP )
442 H( I, I-1 ) = RTEMP
443 TEMP = TEMP / RTEMP
444 IF( I2.GT.I )
445 $ CALL ZSCAL( I2-I, DCONJG( TEMP ), H( I, I+1 ), LDH )
446 CALL ZSCAL( I-I1, TEMP, H( I1, I ), 1 )
447 IF( WANTZ ) THEN
448 CALL ZSCAL( NZ, TEMP, Z( ILOZ, I ), 1 )
449 END IF
450 END IF
451 *
452 130 CONTINUE
453 *
454 * Failure to converge in remaining number of iterations
455 *
456 INFO = I
457 RETURN
458 *
459 140 CONTINUE
460 *
461 * H(I,I-1) is negligible: one eigenvalue has converged.
462 *
463 W( I ) = H( I, I )
464 *
465 * return to start of the main loop with new value of I.
466 *
467 I = L - 1
468 GO TO 30
469 *
470 150 CONTINUE
471 RETURN
472 *
473 * End of ZLAHQR
474 *
475 END
2 $ 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 COMPLEX*16 H( LDH, * ), W( * ), Z( LDZ, * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * ZLAHQR is an auxiliary routine called by CHSEQR to update the
20 * eigenvalues and Schur decomposition already computed by CHSEQR, 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 triangular in rows and
41 * columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1).
42 * ZLAHQR works primarily with the Hessenberg submatrix in rows
43 * and columns ILO to IHI, but applies transformations to all of
44 * H if WANTT is .TRUE..
45 * 1 <= ILO <= max(1,IHI); IHI <= N.
46 *
47 * H (input/output) COMPLEX*16 array, dimension (LDH,N)
48 * On entry, the upper Hessenberg matrix H.
49 * On exit, if INFO is zero and if WANTT is .TRUE., then H
50 * is upper triangular in rows and columns ILO:IHI. If INFO
51 * is zero and if WANTT is .FALSE., then the contents of H
52 * are unspecified on exit. The output state of H in case
53 * INF is positive is below under the description of INFO.
54 *
55 * LDH (input) INTEGER
56 * The leading dimension of the array H. LDH >= max(1,N).
57 *
58 * W (output) COMPLEX*16 array, dimension (N)
59 * The computed eigenvalues ILO to IHI are stored in the
60 * corresponding elements of W. If WANTT is .TRUE., the
61 * eigenvalues are stored in the same order as on the diagonal
62 * of the Schur form returned in H, with W(i) = H(i,i).
63 *
64 * ILOZ (input) INTEGER
65 * IHIZ (input) INTEGER
66 * Specify the rows of Z to which transformations must be
67 * applied if WANTZ is .TRUE..
68 * 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
69 *
70 * Z (input/output) COMPLEX*16 array, dimension (LDZ,N)
71 * If WANTZ is .TRUE., on entry Z must contain the current
72 * matrix Z of transformations accumulated by CHSEQR, and on
73 * exit Z has been updated; transformations are applied only to
74 * the submatrix Z(ILOZ:IHIZ,ILO:IHI).
75 * If WANTZ is .FALSE., Z is not referenced.
76 *
77 * LDZ (input) INTEGER
78 * The leading dimension of the array Z. LDZ >= max(1,N).
79 *
80 * INFO (output) INTEGER
81 * = 0: successful exit
82 * .GT. 0: if INFO = i, ZLAHQR failed to compute all the
83 * eigenvalues ILO to IHI in a total of 30 iterations
84 * per eigenvalue; elements i+1:ihi of W contain
85 * those eigenvalues which have been successfully
86 * computed.
87 *
88 * If INFO .GT. 0 and WANTT is .FALSE., then on exit,
89 * the remaining unconverged eigenvalues are the
90 * eigenvalues of the upper Hessenberg matrix
91 * rows and columns ILO thorugh INFO of the final,
92 * output value of H.
93 *
94 * If INFO .GT. 0 and WANTT is .TRUE., then on exit
95 * (*) (initial value of H)*U = U*(final value of H)
96 * where U is an orthognal matrix. The final
97 * value of H is upper Hessenberg and triangular in
98 * rows and columns INFO+1 through IHI.
99 *
100 * If INFO .GT. 0 and WANTZ is .TRUE., then on exit
101 * (final value of Z) = (initial value of Z)*U
102 * where U is the orthogonal matrix in (*)
103 * (regardless of the value of WANTT.)
104 *
105 * Further Details
106 * ===============
107 *
108 * 02-96 Based on modifications by
109 * David Day, Sandia National Laboratory, USA
110 *
111 * 12-04 Further modifications by
112 * Ralph Byers, University of Kansas, USA
113 * This is a modified version of ZLAHQR from LAPACK version 3.0.
114 * It is (1) more robust against overflow and underflow and
115 * (2) adopts the more conservative Ahues & Tisseur stopping
116 * criterion (LAWN 122, 1997).
117 *
118 * =========================================================
119 *
120 * .. Parameters ..
121 INTEGER ITMAX
122 PARAMETER ( ITMAX = 30 )
123 COMPLEX*16 ZERO, ONE
124 PARAMETER ( ZERO = ( 0.0d0, 0.0d0 ),
125 $ ONE = ( 1.0d0, 0.0d0 ) )
126 DOUBLE PRECISION RZERO, RONE, HALF
127 PARAMETER ( RZERO = 0.0d0, RONE = 1.0d0, HALF = 0.5d0 )
128 DOUBLE PRECISION DAT1
129 PARAMETER ( DAT1 = 3.0d0 / 4.0d0 )
130 * ..
131 * .. Local Scalars ..
132 COMPLEX*16 CDUM, H11, H11S, H22, SC, SUM, T, T1, TEMP, U,
133 $ V2, X, Y
134 DOUBLE PRECISION AA, AB, BA, BB, H10, H21, RTEMP, S, SAFMAX,
135 $ SAFMIN, SMLNUM, SX, T2, TST, ULP
136 INTEGER I, I1, I2, ITS, J, JHI, JLO, K, L, M, NH, NZ
137 * ..
138 * .. Local Arrays ..
139 COMPLEX*16 V( 2 )
140 * ..
141 * .. External Functions ..
142 COMPLEX*16 ZLADIV
143 DOUBLE PRECISION DLAMCH
144 EXTERNAL ZLADIV, DLAMCH
145 * ..
146 * .. External Subroutines ..
147 EXTERNAL DLABAD, ZCOPY, ZLARFG, ZSCAL
148 * ..
149 * .. Statement Functions ..
150 DOUBLE PRECISION CABS1
151 * ..
152 * .. Intrinsic Functions ..
153 INTRINSIC ABS, DBLE, DCONJG, DIMAG, MAX, MIN, SQRT
154 * ..
155 * .. Statement Function definitions ..
156 CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
157 * ..
158 * .. Executable Statements ..
159 *
160 INFO = 0
161 *
162 * Quick return if possible
163 *
164 IF( N.EQ.0 )
165 $ RETURN
166 IF( ILO.EQ.IHI ) THEN
167 W( ILO ) = H( ILO, ILO )
168 RETURN
169 END IF
170 *
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 * ==== ensure that subdiagonal entries are real ====
179 IF( WANTT ) THEN
180 JLO = 1
181 JHI = N
182 ELSE
183 JLO = ILO
184 JHI = IHI
185 END IF
186 DO 20 I = ILO + 1, IHI
187 IF( DIMAG( H( I, I-1 ) ).NE.RZERO ) THEN
188 * ==== The following redundant normalization
189 * . avoids problems with both gradual and
190 * . sudden underflow in ABS(H(I,I-1)) ====
191 SC = H( I, I-1 ) / CABS1( H( I, I-1 ) )
192 SC = DCONJG( SC ) / ABS( SC )
193 H( I, I-1 ) = ABS( H( I, I-1 ) )
194 CALL ZSCAL( JHI-I+1, SC, H( I, I ), LDH )
195 CALL ZSCAL( MIN( JHI, I+1 )-JLO+1, DCONJG( SC ),
196 $ H( JLO, I ), 1 )
197 IF( WANTZ )
198 $ CALL ZSCAL( IHIZ-ILOZ+1, DCONJG( SC ), Z( ILOZ, I ), 1 )
199 END IF
200 20 CONTINUE
201 *
202 NH = IHI - ILO + 1
203 NZ = IHIZ - ILOZ + 1
204 *
205 * Set machine-dependent constants for the stopping criterion.
206 *
207 SAFMIN = DLAMCH( 'SAFE MINIMUM' )
208 SAFMAX = RONE / SAFMIN
209 CALL DLABAD( SAFMIN, SAFMAX )
210 ULP = DLAMCH( 'PRECISION' )
211 SMLNUM = SAFMIN*( DBLE( NH ) / ULP )
212 *
213 * I1 and I2 are the indices of the first row and last column of H
214 * to which transformations must be applied. If eigenvalues only are
215 * being computed, I1 and I2 are set inside the main loop.
216 *
217 IF( WANTT ) THEN
218 I1 = 1
219 I2 = N
220 END IF
221 *
222 * The main loop begins here. I is the loop index and decreases from
223 * IHI to ILO in steps of 1. Each iteration of the loop works
224 * with the active submatrix in rows and columns L to I.
225 * Eigenvalues I+1 to IHI have already converged. Either L = ILO, or
226 * H(L,L-1) is negligible so that the matrix splits.
227 *
228 I = IHI
229 30 CONTINUE
230 IF( I.LT.ILO )
231 $ GO TO 150
232 *
233 * Perform QR iterations on rows and columns ILO to I until a
234 * submatrix of order 1 splits off at the bottom because a
235 * subdiagonal element has become negligible.
236 *
237 L = ILO
238 DO 130 ITS = 0, ITMAX
239 *
240 * Look for a single small subdiagonal element.
241 *
242 DO 40 K = I, L + 1, -1
243 IF( CABS1( H( K, K-1 ) ).LE.SMLNUM )
244 $ GO TO 50
245 TST = CABS1( H( K-1, K-1 ) ) + CABS1( H( K, K ) )
246 IF( TST.EQ.ZERO ) THEN
247 IF( K-2.GE.ILO )
248 $ TST = TST + ABS( DBLE( H( K-1, K-2 ) ) )
249 IF( K+1.LE.IHI )
250 $ TST = TST + ABS( DBLE( H( K+1, K ) ) )
251 END IF
252 * ==== The following is a conservative small subdiagonal
253 * . deflation criterion due to Ahues & Tisseur (LAWN 122,
254 * . 1997). It has better mathematical foundation and
255 * . improves accuracy in some examples. ====
256 IF( ABS( DBLE( H( K, K-1 ) ) ).LE.ULP*TST ) THEN
257 AB = MAX( CABS1( H( K, K-1 ) ), CABS1( H( K-1, K ) ) )
258 BA = MIN( CABS1( H( K, K-1 ) ), CABS1( H( K-1, K ) ) )
259 AA = MAX( CABS1( H( K, K ) ),
260 $ CABS1( H( K-1, K-1 )-H( K, K ) ) )
261 BB = MIN( CABS1( H( K, K ) ),
262 $ CABS1( H( K-1, K-1 )-H( K, K ) ) )
263 S = AA + AB
264 IF( BA*( AB / S ).LE.MAX( SMLNUM,
265 $ ULP*( BB*( AA / S ) ) ) )GO TO 50
266 END IF
267 40 CONTINUE
268 50 CONTINUE
269 L = K
270 IF( L.GT.ILO ) THEN
271 *
272 * H(L,L-1) is negligible
273 *
274 H( L, L-1 ) = ZERO
275 END IF
276 *
277 * Exit from loop if a submatrix of order 1 has split off.
278 *
279 IF( L.GE.I )
280 $ GO TO 140
281 *
282 * Now the active submatrix is in rows and columns L to I. If
283 * eigenvalues only are being computed, only the active submatrix
284 * need be transformed.
285 *
286 IF( .NOT.WANTT ) THEN
287 I1 = L
288 I2 = I
289 END IF
290 *
291 IF( ITS.EQ.10 ) THEN
292 *
293 * Exceptional shift.
294 *
295 S = DAT1*ABS( DBLE( H( L+1, L ) ) )
296 T = S + H( L, L )
297 ELSE IF( ITS.EQ.20 ) THEN
298 *
299 * Exceptional shift.
300 *
301 S = DAT1*ABS( DBLE( H( I, I-1 ) ) )
302 T = S + H( I, I )
303 ELSE
304 *
305 * Wilkinson's shift.
306 *
307 T = H( I, I )
308 U = SQRT( H( I-1, I ) )*SQRT( H( I, I-1 ) )
309 S = CABS1( U )
310 IF( S.NE.RZERO ) THEN
311 X = HALF*( H( I-1, I-1 )-T )
312 SX = CABS1( X )
313 S = MAX( S, CABS1( X ) )
314 Y = S*SQRT( ( X / S )**2+( U / S )**2 )
315 IF( SX.GT.RZERO ) THEN
316 IF( DBLE( X / SX )*DBLE( Y )+DIMAG( X / SX )*
317 $ DIMAG( Y ).LT.RZERO )Y = -Y
318 END IF
319 T = T - U*ZLADIV( U, ( X+Y ) )
320 END IF
321 END IF
322 *
323 * Look for two consecutive small subdiagonal elements.
324 *
325 DO 60 M = I - 1, L + 1, -1
326 *
327 * Determine the effect of starting the single-shift QR
328 * iteration at row M, and see if this would make H(M,M-1)
329 * negligible.
330 *
331 H11 = H( M, M )
332 H22 = H( M+1, M+1 )
333 H11S = H11 - T
334 H21 = DBLE( H( M+1, M ) )
335 S = CABS1( H11S ) + ABS( H21 )
336 H11S = H11S / S
337 H21 = H21 / S
338 V( 1 ) = H11S
339 V( 2 ) = H21
340 H10 = DBLE( H( M, M-1 ) )
341 IF( ABS( H10 )*ABS( H21 ).LE.ULP*
342 $ ( CABS1( H11S )*( CABS1( H11 )+CABS1( H22 ) ) ) )
343 $ GO TO 70
344 60 CONTINUE
345 H11 = H( L, L )
346 H22 = H( L+1, L+1 )
347 H11S = H11 - T
348 H21 = DBLE( H( L+1, L ) )
349 S = CABS1( H11S ) + ABS( H21 )
350 H11S = H11S / S
351 H21 = H21 / S
352 V( 1 ) = H11S
353 V( 2 ) = H21
354 70 CONTINUE
355 *
356 * Single-shift QR step
357 *
358 DO 120 K = M, I - 1
359 *
360 * The first iteration of this loop determines a reflection G
361 * from the vector V and applies it from left and right to H,
362 * thus creating a nonzero bulge below the subdiagonal.
363 *
364 * Each subsequent iteration determines a reflection G to
365 * restore the Hessenberg form in the (K-1)th column, and thus
366 * chases the bulge one step toward the bottom of the active
367 * submatrix.
368 *
369 * V(2) is always real before the call to ZLARFG, and hence
370 * after the call T2 ( = T1*V(2) ) is also real.
371 *
372 IF( K.GT.M )
373 $ CALL ZCOPY( 2, H( K, K-1 ), 1, V, 1 )
374 CALL ZLARFG( 2, V( 1 ), V( 2 ), 1, T1 )
375 IF( K.GT.M ) THEN
376 H( K, K-1 ) = V( 1 )
377 H( K+1, K-1 ) = ZERO
378 END IF
379 V2 = V( 2 )
380 T2 = DBLE( T1*V2 )
381 *
382 * Apply G from the left to transform the rows of the matrix
383 * in columns K to I2.
384 *
385 DO 80 J = K, I2
386 SUM = DCONJG( T1 )*H( K, J ) + T2*H( K+1, J )
387 H( K, J ) = H( K, J ) - SUM
388 H( K+1, J ) = H( K+1, J ) - SUM*V2
389 80 CONTINUE
390 *
391 * Apply G from the right to transform the columns of the
392 * matrix in rows I1 to min(K+2,I).
393 *
394 DO 90 J = I1, MIN( K+2, I )
395 SUM = T1*H( J, K ) + T2*H( J, K+1 )
396 H( J, K ) = H( J, K ) - SUM
397 H( J, K+1 ) = H( J, K+1 ) - SUM*DCONJG( V2 )
398 90 CONTINUE
399 *
400 IF( WANTZ ) THEN
401 *
402 * Accumulate transformations in the matrix Z
403 *
404 DO 100 J = ILOZ, IHIZ
405 SUM = T1*Z( J, K ) + T2*Z( J, K+1 )
406 Z( J, K ) = Z( J, K ) - SUM
407 Z( J, K+1 ) = Z( J, K+1 ) - SUM*DCONJG( V2 )
408 100 CONTINUE
409 END IF
410 *
411 IF( K.EQ.M .AND. M.GT.L ) THEN
412 *
413 * If the QR step was started at row M > L because two
414 * consecutive small subdiagonals were found, then extra
415 * scaling must be performed to ensure that H(M,M-1) remains
416 * real.
417 *
418 TEMP = ONE - T1
419 TEMP = TEMP / ABS( TEMP )
420 H( M+1, M ) = H( M+1, M )*DCONJG( TEMP )
421 IF( M+2.LE.I )
422 $ H( M+2, M+1 ) = H( M+2, M+1 )*TEMP
423 DO 110 J = M, I
424 IF( J.NE.M+1 ) THEN
425 IF( I2.GT.J )
426 $ CALL ZSCAL( I2-J, TEMP, H( J, J+1 ), LDH )
427 CALL ZSCAL( J-I1, DCONJG( TEMP ), H( I1, J ), 1 )
428 IF( WANTZ ) THEN
429 CALL ZSCAL( NZ, DCONJG( TEMP ), Z( ILOZ, J ),
430 $ 1 )
431 END IF
432 END IF
433 110 CONTINUE
434 END IF
435 120 CONTINUE
436 *
437 * Ensure that H(I,I-1) is real.
438 *
439 TEMP = H( I, I-1 )
440 IF( DIMAG( TEMP ).NE.RZERO ) THEN
441 RTEMP = ABS( TEMP )
442 H( I, I-1 ) = RTEMP
443 TEMP = TEMP / RTEMP
444 IF( I2.GT.I )
445 $ CALL ZSCAL( I2-I, DCONJG( TEMP ), H( I, I+1 ), LDH )
446 CALL ZSCAL( I-I1, TEMP, H( I1, I ), 1 )
447 IF( WANTZ ) THEN
448 CALL ZSCAL( NZ, TEMP, Z( ILOZ, I ), 1 )
449 END IF
450 END IF
451 *
452 130 CONTINUE
453 *
454 * Failure to converge in remaining number of iterations
455 *
456 INFO = I
457 RETURN
458 *
459 140 CONTINUE
460 *
461 * H(I,I-1) is negligible: one eigenvalue has converged.
462 *
463 W( I ) = H( I, I )
464 *
465 * return to start of the main loop with new value of I.
466 *
467 I = L - 1
468 GO TO 30
469 *
470 150 CONTINUE
471 RETURN
472 *
473 * End of ZLAHQR
474 *
475 END