Annotation of rpl/lapack/lapack/zlahqr.f, revision 1.1
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>