Annotation of rpl/lapack/lapack/zlahqr.f, revision 1.18
1.11 bertrand 1: *> \brief \b ZLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm.
1.8 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
1.16 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.8 bertrand 7: *
8: *> \htmlonly
1.16 bertrand 9: *> Download ZLAHQR + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlahqr.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlahqr.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlahqr.f">
1.8 bertrand 15: *> [TXT]</a>
1.16 bertrand 16: *> \endhtmlonly
1.8 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
22: * IHIZ, Z, LDZ, INFO )
1.16 bertrand 23: *
1.8 bertrand 24: * .. Scalar Arguments ..
25: * INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
26: * LOGICAL WANTT, WANTZ
27: * ..
28: * .. Array Arguments ..
29: * COMPLEX*16 H( LDH, * ), W( * ), Z( LDZ, * )
30: * ..
1.16 bertrand 31: *
1.8 bertrand 32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> ZLAHQR is an auxiliary routine called by CHSEQR to update the
39: *> eigenvalues and Schur decomposition already computed by CHSEQR, by
40: *> dealing with the Hessenberg submatrix in rows and columns ILO to
41: *> IHI.
42: *> \endverbatim
43: *
44: * Arguments:
45: * ==========
46: *
47: *> \param[in] WANTT
48: *> \verbatim
49: *> WANTT is LOGICAL
50: *> = .TRUE. : the full Schur form T is required;
51: *> = .FALSE.: only eigenvalues are required.
52: *> \endverbatim
53: *>
54: *> \param[in] WANTZ
55: *> \verbatim
56: *> WANTZ is LOGICAL
57: *> = .TRUE. : the matrix of Schur vectors Z is required;
58: *> = .FALSE.: Schur vectors are not required.
59: *> \endverbatim
60: *>
61: *> \param[in] N
62: *> \verbatim
63: *> N is INTEGER
64: *> The order of the matrix H. N >= 0.
65: *> \endverbatim
66: *>
67: *> \param[in] ILO
68: *> \verbatim
69: *> ILO is INTEGER
70: *> \endverbatim
71: *>
72: *> \param[in] IHI
73: *> \verbatim
74: *> IHI is INTEGER
75: *> It is assumed that H is already upper triangular in rows and
76: *> columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1).
77: *> ZLAHQR works primarily with the Hessenberg submatrix in rows
78: *> and columns ILO to IHI, but applies transformations to all of
79: *> H if WANTT is .TRUE..
80: *> 1 <= ILO <= max(1,IHI); IHI <= N.
81: *> \endverbatim
82: *>
83: *> \param[in,out] H
84: *> \verbatim
85: *> H is COMPLEX*16 array, dimension (LDH,N)
86: *> On entry, the upper Hessenberg matrix H.
87: *> On exit, if INFO is zero and if WANTT is .TRUE., then H
88: *> is upper triangular in rows and columns ILO:IHI. If INFO
89: *> is zero and if WANTT is .FALSE., then the contents of H
90: *> are unspecified on exit. The output state of H in case
91: *> INF is positive is below under the description of INFO.
92: *> \endverbatim
93: *>
94: *> \param[in] LDH
95: *> \verbatim
96: *> LDH is INTEGER
97: *> The leading dimension of the array H. LDH >= max(1,N).
98: *> \endverbatim
99: *>
100: *> \param[out] W
101: *> \verbatim
102: *> W is COMPLEX*16 array, dimension (N)
103: *> The computed eigenvalues ILO to IHI are stored in the
104: *> corresponding elements of W. If WANTT is .TRUE., the
105: *> eigenvalues are stored in the same order as on the diagonal
106: *> of the Schur form returned in H, with W(i) = H(i,i).
107: *> \endverbatim
108: *>
109: *> \param[in] ILOZ
110: *> \verbatim
111: *> ILOZ is INTEGER
112: *> \endverbatim
113: *>
114: *> \param[in] IHIZ
115: *> \verbatim
116: *> IHIZ is INTEGER
117: *> Specify the rows of Z to which transformations must be
118: *> applied if WANTZ is .TRUE..
119: *> 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
120: *> \endverbatim
121: *>
122: *> \param[in,out] Z
123: *> \verbatim
124: *> Z is COMPLEX*16 array, dimension (LDZ,N)
125: *> If WANTZ is .TRUE., on entry Z must contain the current
126: *> matrix Z of transformations accumulated by CHSEQR, and on
127: *> exit Z has been updated; transformations are applied only to
128: *> the submatrix Z(ILOZ:IHIZ,ILO:IHI).
129: *> If WANTZ is .FALSE., Z is not referenced.
130: *> \endverbatim
131: *>
132: *> \param[in] LDZ
133: *> \verbatim
134: *> LDZ is INTEGER
135: *> The leading dimension of the array Z. LDZ >= max(1,N).
136: *> \endverbatim
137: *>
138: *> \param[out] INFO
139: *> \verbatim
140: *> INFO is INTEGER
141: *> = 0: successful exit
142: *> .GT. 0: if INFO = i, ZLAHQR failed to compute all the
143: *> eigenvalues ILO to IHI in a total of 30 iterations
144: *> per eigenvalue; elements i+1:ihi of W contain
145: *> those eigenvalues which have been successfully
146: *> computed.
147: *>
148: *> If INFO .GT. 0 and WANTT is .FALSE., then on exit,
149: *> the remaining unconverged eigenvalues are the
150: *> eigenvalues of the upper Hessenberg matrix
151: *> rows and columns ILO thorugh INFO of the final,
152: *> output value of H.
153: *>
154: *> If INFO .GT. 0 and WANTT is .TRUE., then on exit
155: *> (*) (initial value of H)*U = U*(final value of H)
156: *> where U is an orthognal matrix. The final
157: *> value of H is upper Hessenberg and triangular in
158: *> rows and columns INFO+1 through IHI.
159: *>
160: *> If INFO .GT. 0 and WANTZ is .TRUE., then on exit
161: *> (final value of Z) = (initial value of Z)*U
162: *> where U is the orthogonal matrix in (*)
163: *> (regardless of the value of WANTT.)
164: *> \endverbatim
165: *
166: * Authors:
167: * ========
168: *
1.16 bertrand 169: *> \author Univ. of Tennessee
170: *> \author Univ. of California Berkeley
171: *> \author Univ. of Colorado Denver
172: *> \author NAG Ltd.
1.8 bertrand 173: *
1.16 bertrand 174: *> \date December 2016
1.8 bertrand 175: *
176: *> \ingroup complex16OTHERauxiliary
177: *
178: *> \par Contributors:
179: * ==================
180: *>
181: *> \verbatim
182: *>
183: *> 02-96 Based on modifications by
184: *> David Day, Sandia National Laboratory, USA
185: *>
186: *> 12-04 Further modifications by
187: *> Ralph Byers, University of Kansas, USA
188: *> This is a modified version of ZLAHQR from LAPACK version 3.0.
189: *> It is (1) more robust against overflow and underflow and
190: *> (2) adopts the more conservative Ahues & Tisseur stopping
191: *> criterion (LAWN 122, 1997).
192: *> \endverbatim
193: *>
194: * =====================================================================
1.1 bertrand 195: SUBROUTINE ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
196: $ IHIZ, Z, LDZ, INFO )
197: *
1.16 bertrand 198: * -- LAPACK auxiliary routine (version 3.7.0) --
1.8 bertrand 199: * -- LAPACK is a software package provided by Univ. of Tennessee, --
200: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.16 bertrand 201: * December 2016
1.1 bertrand 202: *
203: * .. Scalar Arguments ..
204: INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
205: LOGICAL WANTT, WANTZ
206: * ..
207: * .. Array Arguments ..
208: COMPLEX*16 H( LDH, * ), W( * ), Z( LDZ, * )
209: * ..
210: *
1.8 bertrand 211: * =========================================================
1.1 bertrand 212: *
213: * .. Parameters ..
214: COMPLEX*16 ZERO, ONE
215: PARAMETER ( ZERO = ( 0.0d0, 0.0d0 ),
216: $ ONE = ( 1.0d0, 0.0d0 ) )
217: DOUBLE PRECISION RZERO, RONE, HALF
218: PARAMETER ( RZERO = 0.0d0, RONE = 1.0d0, HALF = 0.5d0 )
219: DOUBLE PRECISION DAT1
220: PARAMETER ( DAT1 = 3.0d0 / 4.0d0 )
221: * ..
222: * .. Local Scalars ..
223: COMPLEX*16 CDUM, H11, H11S, H22, SC, SUM, T, T1, TEMP, U,
224: $ V2, X, Y
225: DOUBLE PRECISION AA, AB, BA, BB, H10, H21, RTEMP, S, SAFMAX,
226: $ SAFMIN, SMLNUM, SX, T2, TST, ULP
1.14 bertrand 227: INTEGER I, I1, I2, ITS, ITMAX, J, JHI, JLO, K, L, M,
228: $ NH, NZ
1.1 bertrand 229: * ..
230: * .. Local Arrays ..
231: COMPLEX*16 V( 2 )
232: * ..
233: * .. External Functions ..
234: COMPLEX*16 ZLADIV
235: DOUBLE PRECISION DLAMCH
236: EXTERNAL ZLADIV, DLAMCH
237: * ..
238: * .. External Subroutines ..
239: EXTERNAL DLABAD, ZCOPY, ZLARFG, ZSCAL
240: * ..
241: * .. Statement Functions ..
242: DOUBLE PRECISION CABS1
243: * ..
244: * .. Intrinsic Functions ..
245: INTRINSIC ABS, DBLE, DCONJG, DIMAG, MAX, MIN, SQRT
246: * ..
247: * .. Statement Function definitions ..
248: CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
249: * ..
250: * .. Executable Statements ..
251: *
252: INFO = 0
253: *
254: * Quick return if possible
255: *
256: IF( N.EQ.0 )
257: $ RETURN
258: IF( ILO.EQ.IHI ) THEN
259: W( ILO ) = H( ILO, ILO )
260: RETURN
261: END IF
262: *
263: * ==== clear out the trash ====
264: DO 10 J = ILO, IHI - 3
265: H( J+2, J ) = ZERO
266: H( J+3, J ) = ZERO
267: 10 CONTINUE
268: IF( ILO.LE.IHI-2 )
269: $ H( IHI, IHI-2 ) = ZERO
270: * ==== ensure that subdiagonal entries are real ====
271: IF( WANTT ) THEN
272: JLO = 1
273: JHI = N
274: ELSE
275: JLO = ILO
276: JHI = IHI
277: END IF
278: DO 20 I = ILO + 1, IHI
279: IF( DIMAG( H( I, I-1 ) ).NE.RZERO ) THEN
280: * ==== The following redundant normalization
281: * . avoids problems with both gradual and
282: * . sudden underflow in ABS(H(I,I-1)) ====
283: SC = H( I, I-1 ) / CABS1( H( I, I-1 ) )
284: SC = DCONJG( SC ) / ABS( SC )
285: H( I, I-1 ) = ABS( H( I, I-1 ) )
286: CALL ZSCAL( JHI-I+1, SC, H( I, I ), LDH )
287: CALL ZSCAL( MIN( JHI, I+1 )-JLO+1, DCONJG( SC ),
288: $ H( JLO, I ), 1 )
289: IF( WANTZ )
290: $ CALL ZSCAL( IHIZ-ILOZ+1, DCONJG( SC ), Z( ILOZ, I ), 1 )
291: END IF
292: 20 CONTINUE
293: *
294: NH = IHI - ILO + 1
295: NZ = IHIZ - ILOZ + 1
296: *
297: * Set machine-dependent constants for the stopping criterion.
298: *
299: SAFMIN = DLAMCH( 'SAFE MINIMUM' )
300: SAFMAX = RONE / SAFMIN
301: CALL DLABAD( SAFMIN, SAFMAX )
302: ULP = DLAMCH( 'PRECISION' )
303: SMLNUM = SAFMIN*( DBLE( NH ) / ULP )
304: *
305: * I1 and I2 are the indices of the first row and last column of H
306: * to which transformations must be applied. If eigenvalues only are
307: * being computed, I1 and I2 are set inside the main loop.
308: *
309: IF( WANTT ) THEN
310: I1 = 1
311: I2 = N
312: END IF
313: *
1.14 bertrand 314: * ITMAX is the total number of QR iterations allowed.
315: *
1.16 bertrand 316: ITMAX = 30 * MAX( 10, NH )
1.14 bertrand 317: *
1.1 bertrand 318: * The main loop begins here. I is the loop index and decreases from
319: * IHI to ILO in steps of 1. Each iteration of the loop works
320: * with the active submatrix in rows and columns L to I.
321: * Eigenvalues I+1 to IHI have already converged. Either L = ILO, or
322: * H(L,L-1) is negligible so that the matrix splits.
323: *
324: I = IHI
325: 30 CONTINUE
326: IF( I.LT.ILO )
327: $ GO TO 150
328: *
329: * Perform QR iterations on rows and columns ILO to I until a
330: * submatrix of order 1 splits off at the bottom because a
331: * subdiagonal element has become negligible.
332: *
333: L = ILO
334: DO 130 ITS = 0, ITMAX
335: *
336: * Look for a single small subdiagonal element.
337: *
338: DO 40 K = I, L + 1, -1
339: IF( CABS1( H( K, K-1 ) ).LE.SMLNUM )
340: $ GO TO 50
341: TST = CABS1( H( K-1, K-1 ) ) + CABS1( H( K, K ) )
342: IF( TST.EQ.ZERO ) THEN
343: IF( K-2.GE.ILO )
344: $ TST = TST + ABS( DBLE( H( K-1, K-2 ) ) )
345: IF( K+1.LE.IHI )
346: $ TST = TST + ABS( DBLE( H( K+1, K ) ) )
347: END IF
348: * ==== The following is a conservative small subdiagonal
349: * . deflation criterion due to Ahues & Tisseur (LAWN 122,
350: * . 1997). It has better mathematical foundation and
351: * . improves accuracy in some examples. ====
352: IF( ABS( DBLE( H( K, K-1 ) ) ).LE.ULP*TST ) THEN
353: AB = MAX( CABS1( H( K, K-1 ) ), CABS1( H( K-1, K ) ) )
354: BA = MIN( CABS1( H( K, K-1 ) ), CABS1( H( K-1, K ) ) )
355: AA = MAX( CABS1( H( K, K ) ),
356: $ CABS1( H( K-1, K-1 )-H( K, K ) ) )
357: BB = MIN( CABS1( H( K, K ) ),
358: $ CABS1( H( K-1, K-1 )-H( K, K ) ) )
359: S = AA + AB
360: IF( BA*( AB / S ).LE.MAX( SMLNUM,
361: $ ULP*( BB*( AA / S ) ) ) )GO TO 50
362: END IF
363: 40 CONTINUE
364: 50 CONTINUE
365: L = K
366: IF( L.GT.ILO ) THEN
367: *
368: * H(L,L-1) is negligible
369: *
370: H( L, L-1 ) = ZERO
371: END IF
372: *
373: * Exit from loop if a submatrix of order 1 has split off.
374: *
375: IF( L.GE.I )
376: $ GO TO 140
377: *
378: * Now the active submatrix is in rows and columns L to I. If
379: * eigenvalues only are being computed, only the active submatrix
380: * need be transformed.
381: *
382: IF( .NOT.WANTT ) THEN
383: I1 = L
384: I2 = I
385: END IF
386: *
387: IF( ITS.EQ.10 ) THEN
388: *
389: * Exceptional shift.
390: *
391: S = DAT1*ABS( DBLE( H( L+1, L ) ) )
392: T = S + H( L, L )
393: ELSE IF( ITS.EQ.20 ) THEN
394: *
395: * Exceptional shift.
396: *
397: S = DAT1*ABS( DBLE( H( I, I-1 ) ) )
398: T = S + H( I, I )
399: ELSE
400: *
401: * Wilkinson's shift.
402: *
403: T = H( I, I )
404: U = SQRT( H( I-1, I ) )*SQRT( H( I, I-1 ) )
405: S = CABS1( U )
406: IF( S.NE.RZERO ) THEN
407: X = HALF*( H( I-1, I-1 )-T )
408: SX = CABS1( X )
409: S = MAX( S, CABS1( X ) )
410: Y = S*SQRT( ( X / S )**2+( U / S )**2 )
411: IF( SX.GT.RZERO ) THEN
412: IF( DBLE( X / SX )*DBLE( Y )+DIMAG( X / SX )*
413: $ DIMAG( Y ).LT.RZERO )Y = -Y
414: END IF
415: T = T - U*ZLADIV( U, ( X+Y ) )
416: END IF
417: END IF
418: *
419: * Look for two consecutive small subdiagonal elements.
420: *
421: DO 60 M = I - 1, L + 1, -1
422: *
423: * Determine the effect of starting the single-shift QR
424: * iteration at row M, and see if this would make H(M,M-1)
425: * negligible.
426: *
427: H11 = H( M, M )
428: H22 = H( M+1, M+1 )
429: H11S = H11 - T
430: H21 = DBLE( H( M+1, M ) )
431: S = CABS1( H11S ) + ABS( H21 )
432: H11S = H11S / S
433: H21 = H21 / S
434: V( 1 ) = H11S
435: V( 2 ) = H21
436: H10 = DBLE( H( M, M-1 ) )
437: IF( ABS( H10 )*ABS( H21 ).LE.ULP*
438: $ ( CABS1( H11S )*( CABS1( H11 )+CABS1( H22 ) ) ) )
439: $ GO TO 70
440: 60 CONTINUE
441: H11 = H( L, L )
442: H22 = H( L+1, L+1 )
443: H11S = H11 - T
444: H21 = DBLE( H( L+1, L ) )
445: S = CABS1( H11S ) + ABS( H21 )
446: H11S = H11S / S
447: H21 = H21 / S
448: V( 1 ) = H11S
449: V( 2 ) = H21
450: 70 CONTINUE
451: *
452: * Single-shift QR step
453: *
454: DO 120 K = M, I - 1
455: *
456: * The first iteration of this loop determines a reflection G
457: * from the vector V and applies it from left and right to H,
458: * thus creating a nonzero bulge below the subdiagonal.
459: *
460: * Each subsequent iteration determines a reflection G to
461: * restore the Hessenberg form in the (K-1)th column, and thus
462: * chases the bulge one step toward the bottom of the active
463: * submatrix.
464: *
465: * V(2) is always real before the call to ZLARFG, and hence
466: * after the call T2 ( = T1*V(2) ) is also real.
467: *
468: IF( K.GT.M )
469: $ CALL ZCOPY( 2, H( K, K-1 ), 1, V, 1 )
470: CALL ZLARFG( 2, V( 1 ), V( 2 ), 1, T1 )
471: IF( K.GT.M ) THEN
472: H( K, K-1 ) = V( 1 )
473: H( K+1, K-1 ) = ZERO
474: END IF
475: V2 = V( 2 )
476: T2 = DBLE( T1*V2 )
477: *
478: * Apply G from the left to transform the rows of the matrix
479: * in columns K to I2.
480: *
481: DO 80 J = K, I2
482: SUM = DCONJG( T1 )*H( K, J ) + T2*H( K+1, J )
483: H( K, J ) = H( K, J ) - SUM
484: H( K+1, J ) = H( K+1, J ) - SUM*V2
485: 80 CONTINUE
486: *
487: * Apply G from the right to transform the columns of the
488: * matrix in rows I1 to min(K+2,I).
489: *
490: DO 90 J = I1, MIN( K+2, I )
491: SUM = T1*H( J, K ) + T2*H( J, K+1 )
492: H( J, K ) = H( J, K ) - SUM
493: H( J, K+1 ) = H( J, K+1 ) - SUM*DCONJG( V2 )
494: 90 CONTINUE
495: *
496: IF( WANTZ ) THEN
497: *
498: * Accumulate transformations in the matrix Z
499: *
500: DO 100 J = ILOZ, IHIZ
501: SUM = T1*Z( J, K ) + T2*Z( J, K+1 )
502: Z( J, K ) = Z( J, K ) - SUM
503: Z( J, K+1 ) = Z( J, K+1 ) - SUM*DCONJG( V2 )
504: 100 CONTINUE
505: END IF
506: *
507: IF( K.EQ.M .AND. M.GT.L ) THEN
508: *
509: * If the QR step was started at row M > L because two
510: * consecutive small subdiagonals were found, then extra
511: * scaling must be performed to ensure that H(M,M-1) remains
512: * real.
513: *
514: TEMP = ONE - T1
515: TEMP = TEMP / ABS( TEMP )
516: H( M+1, M ) = H( M+1, M )*DCONJG( TEMP )
517: IF( M+2.LE.I )
518: $ H( M+2, M+1 ) = H( M+2, M+1 )*TEMP
519: DO 110 J = M, I
520: IF( J.NE.M+1 ) THEN
521: IF( I2.GT.J )
522: $ CALL ZSCAL( I2-J, TEMP, H( J, J+1 ), LDH )
523: CALL ZSCAL( J-I1, DCONJG( TEMP ), H( I1, J ), 1 )
524: IF( WANTZ ) THEN
525: CALL ZSCAL( NZ, DCONJG( TEMP ), Z( ILOZ, J ),
526: $ 1 )
527: END IF
528: END IF
529: 110 CONTINUE
530: END IF
531: 120 CONTINUE
532: *
533: * Ensure that H(I,I-1) is real.
534: *
535: TEMP = H( I, I-1 )
536: IF( DIMAG( TEMP ).NE.RZERO ) THEN
537: RTEMP = ABS( TEMP )
538: H( I, I-1 ) = RTEMP
539: TEMP = TEMP / RTEMP
540: IF( I2.GT.I )
541: $ CALL ZSCAL( I2-I, DCONJG( TEMP ), H( I, I+1 ), LDH )
542: CALL ZSCAL( I-I1, TEMP, H( I1, I ), 1 )
543: IF( WANTZ ) THEN
544: CALL ZSCAL( NZ, TEMP, Z( ILOZ, I ), 1 )
545: END IF
546: END IF
547: *
548: 130 CONTINUE
549: *
550: * Failure to converge in remaining number of iterations
551: *
552: INFO = I
553: RETURN
554: *
555: 140 CONTINUE
556: *
557: * H(I,I-1) is negligible: one eigenvalue has converged.
558: *
559: W( I ) = H( I, I )
560: *
561: * return to start of the main loop with new value of I.
562: *
563: I = L - 1
564: GO TO 30
565: *
566: 150 CONTINUE
567: RETURN
568: *
569: * End of ZLAHQR
570: *
571: 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