Annotation of rpl/lapack/lapack/zlahqr.f, revision 1.7
1.1 bertrand 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlahqr.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack