Annotation of rpl/lapack/lapack/zlaqr5.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZLAQR5( WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, S,
! 2: $ H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, LDU, NV,
! 3: $ WV, LDWV, NH, WH, LDWH )
! 4: *
! 5: * -- LAPACK auxiliary routine (version 3.2) --
! 6: * Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
! 7: * November 2006
! 8: *
! 9: * .. Scalar Arguments ..
! 10: INTEGER IHIZ, ILOZ, KACC22, KBOT, KTOP, LDH, LDU, LDV,
! 11: $ LDWH, LDWV, LDZ, N, NH, NSHFTS, NV
! 12: LOGICAL WANTT, WANTZ
! 13: * ..
! 14: * .. Array Arguments ..
! 15: COMPLEX*16 H( LDH, * ), S( * ), U( LDU, * ), V( LDV, * ),
! 16: $ WH( LDWH, * ), WV( LDWV, * ), Z( LDZ, * )
! 17: * ..
! 18: *
! 19: * This auxiliary subroutine called by ZLAQR0 performs a
! 20: * single small-bulge multi-shift QR sweep.
! 21: *
! 22: * WANTT (input) logical scalar
! 23: * WANTT = .true. if the triangular Schur factor
! 24: * is being computed. WANTT is set to .false. otherwise.
! 25: *
! 26: * WANTZ (input) logical scalar
! 27: * WANTZ = .true. if the unitary Schur factor is being
! 28: * computed. WANTZ is set to .false. otherwise.
! 29: *
! 30: * KACC22 (input) integer with value 0, 1, or 2.
! 31: * Specifies the computation mode of far-from-diagonal
! 32: * orthogonal updates.
! 33: * = 0: ZLAQR5 does not accumulate reflections and does not
! 34: * use matrix-matrix multiply to update far-from-diagonal
! 35: * matrix entries.
! 36: * = 1: ZLAQR5 accumulates reflections and uses matrix-matrix
! 37: * multiply to update the far-from-diagonal matrix entries.
! 38: * = 2: ZLAQR5 accumulates reflections, uses matrix-matrix
! 39: * multiply to update the far-from-diagonal matrix entries,
! 40: * and takes advantage of 2-by-2 block structure during
! 41: * matrix multiplies.
! 42: *
! 43: * N (input) integer scalar
! 44: * N is the order of the Hessenberg matrix H upon which this
! 45: * subroutine operates.
! 46: *
! 47: * KTOP (input) integer scalar
! 48: * KBOT (input) integer scalar
! 49: * These are the first and last rows and columns of an
! 50: * isolated diagonal block upon which the QR sweep is to be
! 51: * applied. It is assumed without a check that
! 52: * either KTOP = 1 or H(KTOP,KTOP-1) = 0
! 53: * and
! 54: * either KBOT = N or H(KBOT+1,KBOT) = 0.
! 55: *
! 56: * NSHFTS (input) integer scalar
! 57: * NSHFTS gives the number of simultaneous shifts. NSHFTS
! 58: * must be positive and even.
! 59: *
! 60: * S (input/output) COMPLEX*16 array of size (NSHFTS)
! 61: * S contains the shifts of origin that define the multi-
! 62: * shift QR sweep. On output S may be reordered.
! 63: *
! 64: * H (input/output) COMPLEX*16 array of size (LDH,N)
! 65: * On input H contains a Hessenberg matrix. On output a
! 66: * multi-shift QR sweep with shifts SR(J)+i*SI(J) is applied
! 67: * to the isolated diagonal block in rows and columns KTOP
! 68: * through KBOT.
! 69: *
! 70: * LDH (input) integer scalar
! 71: * LDH is the leading dimension of H just as declared in the
! 72: * calling procedure. LDH.GE.MAX(1,N).
! 73: *
! 74: * ILOZ (input) INTEGER
! 75: * IHIZ (input) INTEGER
! 76: * Specify the rows of Z to which transformations must be
! 77: * applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N
! 78: *
! 79: * Z (input/output) COMPLEX*16 array of size (LDZ,IHI)
! 80: * If WANTZ = .TRUE., then the QR Sweep unitary
! 81: * similarity transformation is accumulated into
! 82: * Z(ILOZ:IHIZ,ILO:IHI) from the right.
! 83: * If WANTZ = .FALSE., then Z is unreferenced.
! 84: *
! 85: * LDZ (input) integer scalar
! 86: * LDA is the leading dimension of Z just as declared in
! 87: * the calling procedure. LDZ.GE.N.
! 88: *
! 89: * V (workspace) COMPLEX*16 array of size (LDV,NSHFTS/2)
! 90: *
! 91: * LDV (input) integer scalar
! 92: * LDV is the leading dimension of V as declared in the
! 93: * calling procedure. LDV.GE.3.
! 94: *
! 95: * U (workspace) COMPLEX*16 array of size
! 96: * (LDU,3*NSHFTS-3)
! 97: *
! 98: * LDU (input) integer scalar
! 99: * LDU is the leading dimension of U just as declared in the
! 100: * in the calling subroutine. LDU.GE.3*NSHFTS-3.
! 101: *
! 102: * NH (input) integer scalar
! 103: * NH is the number of columns in array WH available for
! 104: * workspace. NH.GE.1.
! 105: *
! 106: * WH (workspace) COMPLEX*16 array of size (LDWH,NH)
! 107: *
! 108: * LDWH (input) integer scalar
! 109: * Leading dimension of WH just as declared in the
! 110: * calling procedure. LDWH.GE.3*NSHFTS-3.
! 111: *
! 112: * NV (input) integer scalar
! 113: * NV is the number of rows in WV agailable for workspace.
! 114: * NV.GE.1.
! 115: *
! 116: * WV (workspace) COMPLEX*16 array of size
! 117: * (LDWV,3*NSHFTS-3)
! 118: *
! 119: * LDWV (input) integer scalar
! 120: * LDWV is the leading dimension of WV as declared in the
! 121: * in the calling subroutine. LDWV.GE.NV.
! 122: *
! 123: * ================================================================
! 124: * Based on contributions by
! 125: * Karen Braman and Ralph Byers, Department of Mathematics,
! 126: * University of Kansas, USA
! 127: *
! 128: * ================================================================
! 129: * Reference:
! 130: *
! 131: * K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
! 132: * Algorithm Part I: Maintaining Well Focused Shifts, and
! 133: * Level 3 Performance, SIAM Journal of Matrix Analysis,
! 134: * volume 23, pages 929--947, 2002.
! 135: *
! 136: * ================================================================
! 137: * .. Parameters ..
! 138: COMPLEX*16 ZERO, ONE
! 139: PARAMETER ( ZERO = ( 0.0d0, 0.0d0 ),
! 140: $ ONE = ( 1.0d0, 0.0d0 ) )
! 141: DOUBLE PRECISION RZERO, RONE
! 142: PARAMETER ( RZERO = 0.0d0, RONE = 1.0d0 )
! 143: * ..
! 144: * .. Local Scalars ..
! 145: COMPLEX*16 ALPHA, BETA, CDUM, REFSUM
! 146: DOUBLE PRECISION H11, H12, H21, H22, SAFMAX, SAFMIN, SCL,
! 147: $ SMLNUM, TST1, TST2, ULP
! 148: INTEGER I2, I4, INCOL, J, J2, J4, JBOT, JCOL, JLEN,
! 149: $ JROW, JTOP, K, K1, KDU, KMS, KNZ, KRCOL, KZS,
! 150: $ M, M22, MBOT, MEND, MSTART, MTOP, NBMPS, NDCOL,
! 151: $ NS, NU
! 152: LOGICAL ACCUM, BLK22, BMP22
! 153: * ..
! 154: * .. External Functions ..
! 155: DOUBLE PRECISION DLAMCH
! 156: EXTERNAL DLAMCH
! 157: * ..
! 158: * .. Intrinsic Functions ..
! 159: *
! 160: INTRINSIC ABS, DBLE, DCONJG, DIMAG, MAX, MIN, MOD
! 161: * ..
! 162: * .. Local Arrays ..
! 163: COMPLEX*16 VT( 3 )
! 164: * ..
! 165: * .. External Subroutines ..
! 166: EXTERNAL DLABAD, ZGEMM, ZLACPY, ZLAQR1, ZLARFG, ZLASET,
! 167: $ ZTRMM
! 168: * ..
! 169: * .. Statement Functions ..
! 170: DOUBLE PRECISION CABS1
! 171: * ..
! 172: * .. Statement Function definitions ..
! 173: CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
! 174: * ..
! 175: * .. Executable Statements ..
! 176: *
! 177: * ==== If there are no shifts, then there is nothing to do. ====
! 178: *
! 179: IF( NSHFTS.LT.2 )
! 180: $ RETURN
! 181: *
! 182: * ==== If the active block is empty or 1-by-1, then there
! 183: * . is nothing to do. ====
! 184: *
! 185: IF( KTOP.GE.KBOT )
! 186: $ RETURN
! 187: *
! 188: * ==== NSHFTS is supposed to be even, but if it is odd,
! 189: * . then simply reduce it by one. ====
! 190: *
! 191: NS = NSHFTS - MOD( NSHFTS, 2 )
! 192: *
! 193: * ==== Machine constants for deflation ====
! 194: *
! 195: SAFMIN = DLAMCH( 'SAFE MINIMUM' )
! 196: SAFMAX = RONE / SAFMIN
! 197: CALL DLABAD( SAFMIN, SAFMAX )
! 198: ULP = DLAMCH( 'PRECISION' )
! 199: SMLNUM = SAFMIN*( DBLE( N ) / ULP )
! 200: *
! 201: * ==== Use accumulated reflections to update far-from-diagonal
! 202: * . entries ? ====
! 203: *
! 204: ACCUM = ( KACC22.EQ.1 ) .OR. ( KACC22.EQ.2 )
! 205: *
! 206: * ==== If so, exploit the 2-by-2 block structure? ====
! 207: *
! 208: BLK22 = ( NS.GT.2 ) .AND. ( KACC22.EQ.2 )
! 209: *
! 210: * ==== clear trash ====
! 211: *
! 212: IF( KTOP+2.LE.KBOT )
! 213: $ H( KTOP+2, KTOP ) = ZERO
! 214: *
! 215: * ==== NBMPS = number of 2-shift bulges in the chain ====
! 216: *
! 217: NBMPS = NS / 2
! 218: *
! 219: * ==== KDU = width of slab ====
! 220: *
! 221: KDU = 6*NBMPS - 3
! 222: *
! 223: * ==== Create and chase chains of NBMPS bulges ====
! 224: *
! 225: DO 210 INCOL = 3*( 1-NBMPS ) + KTOP - 1, KBOT - 2, 3*NBMPS - 2
! 226: NDCOL = INCOL + KDU
! 227: IF( ACCUM )
! 228: $ CALL ZLASET( 'ALL', KDU, KDU, ZERO, ONE, U, LDU )
! 229: *
! 230: * ==== Near-the-diagonal bulge chase. The following loop
! 231: * . performs the near-the-diagonal part of a small bulge
! 232: * . multi-shift QR sweep. Each 6*NBMPS-2 column diagonal
! 233: * . chunk extends from column INCOL to column NDCOL
! 234: * . (including both column INCOL and column NDCOL). The
! 235: * . following loop chases a 3*NBMPS column long chain of
! 236: * . NBMPS bulges 3*NBMPS-2 columns to the right. (INCOL
! 237: * . may be less than KTOP and and NDCOL may be greater than
! 238: * . KBOT indicating phantom columns from which to chase
! 239: * . bulges before they are actually introduced or to which
! 240: * . to chase bulges beyond column KBOT.) ====
! 241: *
! 242: DO 140 KRCOL = INCOL, MIN( INCOL+3*NBMPS-3, KBOT-2 )
! 243: *
! 244: * ==== Bulges number MTOP to MBOT are active double implicit
! 245: * . shift bulges. There may or may not also be small
! 246: * . 2-by-2 bulge, if there is room. The inactive bulges
! 247: * . (if any) must wait until the active bulges have moved
! 248: * . down the diagonal to make room. The phantom matrix
! 249: * . paradigm described above helps keep track. ====
! 250: *
! 251: MTOP = MAX( 1, ( ( KTOP-1 )-KRCOL+2 ) / 3+1 )
! 252: MBOT = MIN( NBMPS, ( KBOT-KRCOL ) / 3 )
! 253: M22 = MBOT + 1
! 254: BMP22 = ( MBOT.LT.NBMPS ) .AND. ( KRCOL+3*( M22-1 ) ).EQ.
! 255: $ ( KBOT-2 )
! 256: *
! 257: * ==== Generate reflections to chase the chain right
! 258: * . one column. (The minimum value of K is KTOP-1.) ====
! 259: *
! 260: DO 10 M = MTOP, MBOT
! 261: K = KRCOL + 3*( M-1 )
! 262: IF( K.EQ.KTOP-1 ) THEN
! 263: CALL ZLAQR1( 3, H( KTOP, KTOP ), LDH, S( 2*M-1 ),
! 264: $ S( 2*M ), V( 1, M ) )
! 265: ALPHA = V( 1, M )
! 266: CALL ZLARFG( 3, ALPHA, V( 2, M ), 1, V( 1, M ) )
! 267: ELSE
! 268: BETA = H( K+1, K )
! 269: V( 2, M ) = H( K+2, K )
! 270: V( 3, M ) = H( K+3, K )
! 271: CALL ZLARFG( 3, BETA, V( 2, M ), 1, V( 1, M ) )
! 272: *
! 273: * ==== A Bulge may collapse because of vigilant
! 274: * . deflation or destructive underflow. In the
! 275: * . underflow case, try the two-small-subdiagonals
! 276: * . trick to try to reinflate the bulge. ====
! 277: *
! 278: IF( H( K+3, K ).NE.ZERO .OR. H( K+3, K+1 ).NE.
! 279: $ ZERO .OR. H( K+3, K+2 ).EQ.ZERO ) THEN
! 280: *
! 281: * ==== Typical case: not collapsed (yet). ====
! 282: *
! 283: H( K+1, K ) = BETA
! 284: H( K+2, K ) = ZERO
! 285: H( K+3, K ) = ZERO
! 286: ELSE
! 287: *
! 288: * ==== Atypical case: collapsed. Attempt to
! 289: * . reintroduce ignoring H(K+1,K) and H(K+2,K).
! 290: * . If the fill resulting from the new
! 291: * . reflector is too large, then abandon it.
! 292: * . Otherwise, use the new one. ====
! 293: *
! 294: CALL ZLAQR1( 3, H( K+1, K+1 ), LDH, S( 2*M-1 ),
! 295: $ S( 2*M ), VT )
! 296: ALPHA = VT( 1 )
! 297: CALL ZLARFG( 3, ALPHA, VT( 2 ), 1, VT( 1 ) )
! 298: REFSUM = DCONJG( VT( 1 ) )*
! 299: $ ( H( K+1, K )+DCONJG( VT( 2 ) )*
! 300: $ H( K+2, K ) )
! 301: *
! 302: IF( CABS1( H( K+2, K )-REFSUM*VT( 2 ) )+
! 303: $ CABS1( REFSUM*VT( 3 ) ).GT.ULP*
! 304: $ ( CABS1( H( K, K ) )+CABS1( H( K+1,
! 305: $ K+1 ) )+CABS1( H( K+2, K+2 ) ) ) ) THEN
! 306: *
! 307: * ==== Starting a new bulge here would
! 308: * . create non-negligible fill. Use
! 309: * . the old one with trepidation. ====
! 310: *
! 311: H( K+1, K ) = BETA
! 312: H( K+2, K ) = ZERO
! 313: H( K+3, K ) = ZERO
! 314: ELSE
! 315: *
! 316: * ==== Stating a new bulge here would
! 317: * . create only negligible fill.
! 318: * . Replace the old reflector with
! 319: * . the new one. ====
! 320: *
! 321: H( K+1, K ) = H( K+1, K ) - REFSUM
! 322: H( K+2, K ) = ZERO
! 323: H( K+3, K ) = ZERO
! 324: V( 1, M ) = VT( 1 )
! 325: V( 2, M ) = VT( 2 )
! 326: V( 3, M ) = VT( 3 )
! 327: END IF
! 328: END IF
! 329: END IF
! 330: 10 CONTINUE
! 331: *
! 332: * ==== Generate a 2-by-2 reflection, if needed. ====
! 333: *
! 334: K = KRCOL + 3*( M22-1 )
! 335: IF( BMP22 ) THEN
! 336: IF( K.EQ.KTOP-1 ) THEN
! 337: CALL ZLAQR1( 2, H( K+1, K+1 ), LDH, S( 2*M22-1 ),
! 338: $ S( 2*M22 ), V( 1, M22 ) )
! 339: BETA = V( 1, M22 )
! 340: CALL ZLARFG( 2, BETA, V( 2, M22 ), 1, V( 1, M22 ) )
! 341: ELSE
! 342: BETA = H( K+1, K )
! 343: V( 2, M22 ) = H( K+2, K )
! 344: CALL ZLARFG( 2, BETA, V( 2, M22 ), 1, V( 1, M22 ) )
! 345: H( K+1, K ) = BETA
! 346: H( K+2, K ) = ZERO
! 347: END IF
! 348: END IF
! 349: *
! 350: * ==== Multiply H by reflections from the left ====
! 351: *
! 352: IF( ACCUM ) THEN
! 353: JBOT = MIN( NDCOL, KBOT )
! 354: ELSE IF( WANTT ) THEN
! 355: JBOT = N
! 356: ELSE
! 357: JBOT = KBOT
! 358: END IF
! 359: DO 30 J = MAX( KTOP, KRCOL ), JBOT
! 360: MEND = MIN( MBOT, ( J-KRCOL+2 ) / 3 )
! 361: DO 20 M = MTOP, MEND
! 362: K = KRCOL + 3*( M-1 )
! 363: REFSUM = DCONJG( V( 1, M ) )*
! 364: $ ( H( K+1, J )+DCONJG( V( 2, M ) )*
! 365: $ H( K+2, J )+DCONJG( V( 3, M ) )*H( K+3, J ) )
! 366: H( K+1, J ) = H( K+1, J ) - REFSUM
! 367: H( K+2, J ) = H( K+2, J ) - REFSUM*V( 2, M )
! 368: H( K+3, J ) = H( K+3, J ) - REFSUM*V( 3, M )
! 369: 20 CONTINUE
! 370: 30 CONTINUE
! 371: IF( BMP22 ) THEN
! 372: K = KRCOL + 3*( M22-1 )
! 373: DO 40 J = MAX( K+1, KTOP ), JBOT
! 374: REFSUM = DCONJG( V( 1, M22 ) )*
! 375: $ ( H( K+1, J )+DCONJG( V( 2, M22 ) )*
! 376: $ H( K+2, J ) )
! 377: H( K+1, J ) = H( K+1, J ) - REFSUM
! 378: H( K+2, J ) = H( K+2, J ) - REFSUM*V( 2, M22 )
! 379: 40 CONTINUE
! 380: END IF
! 381: *
! 382: * ==== Multiply H by reflections from the right.
! 383: * . Delay filling in the last row until the
! 384: * . vigilant deflation check is complete. ====
! 385: *
! 386: IF( ACCUM ) THEN
! 387: JTOP = MAX( KTOP, INCOL )
! 388: ELSE IF( WANTT ) THEN
! 389: JTOP = 1
! 390: ELSE
! 391: JTOP = KTOP
! 392: END IF
! 393: DO 80 M = MTOP, MBOT
! 394: IF( V( 1, M ).NE.ZERO ) THEN
! 395: K = KRCOL + 3*( M-1 )
! 396: DO 50 J = JTOP, MIN( KBOT, K+3 )
! 397: REFSUM = V( 1, M )*( H( J, K+1 )+V( 2, M )*
! 398: $ H( J, K+2 )+V( 3, M )*H( J, K+3 ) )
! 399: H( J, K+1 ) = H( J, K+1 ) - REFSUM
! 400: H( J, K+2 ) = H( J, K+2 ) -
! 401: $ REFSUM*DCONJG( V( 2, M ) )
! 402: H( J, K+3 ) = H( J, K+3 ) -
! 403: $ REFSUM*DCONJG( V( 3, M ) )
! 404: 50 CONTINUE
! 405: *
! 406: IF( ACCUM ) THEN
! 407: *
! 408: * ==== Accumulate U. (If necessary, update Z later
! 409: * . with with an efficient matrix-matrix
! 410: * . multiply.) ====
! 411: *
! 412: KMS = K - INCOL
! 413: DO 60 J = MAX( 1, KTOP-INCOL ), KDU
! 414: REFSUM = V( 1, M )*( U( J, KMS+1 )+V( 2, M )*
! 415: $ U( J, KMS+2 )+V( 3, M )*U( J, KMS+3 ) )
! 416: U( J, KMS+1 ) = U( J, KMS+1 ) - REFSUM
! 417: U( J, KMS+2 ) = U( J, KMS+2 ) -
! 418: $ REFSUM*DCONJG( V( 2, M ) )
! 419: U( J, KMS+3 ) = U( J, KMS+3 ) -
! 420: $ REFSUM*DCONJG( V( 3, M ) )
! 421: 60 CONTINUE
! 422: ELSE IF( WANTZ ) THEN
! 423: *
! 424: * ==== U is not accumulated, so update Z
! 425: * . now by multiplying by reflections
! 426: * . from the right. ====
! 427: *
! 428: DO 70 J = ILOZ, IHIZ
! 429: REFSUM = V( 1, M )*( Z( J, K+1 )+V( 2, M )*
! 430: $ Z( J, K+2 )+V( 3, M )*Z( J, K+3 ) )
! 431: Z( J, K+1 ) = Z( J, K+1 ) - REFSUM
! 432: Z( J, K+2 ) = Z( J, K+2 ) -
! 433: $ REFSUM*DCONJG( V( 2, M ) )
! 434: Z( J, K+3 ) = Z( J, K+3 ) -
! 435: $ REFSUM*DCONJG( V( 3, M ) )
! 436: 70 CONTINUE
! 437: END IF
! 438: END IF
! 439: 80 CONTINUE
! 440: *
! 441: * ==== Special case: 2-by-2 reflection (if needed) ====
! 442: *
! 443: K = KRCOL + 3*( M22-1 )
! 444: IF( BMP22 .AND. ( V( 1, M22 ).NE.ZERO ) ) THEN
! 445: DO 90 J = JTOP, MIN( KBOT, K+3 )
! 446: REFSUM = V( 1, M22 )*( H( J, K+1 )+V( 2, M22 )*
! 447: $ H( J, K+2 ) )
! 448: H( J, K+1 ) = H( J, K+1 ) - REFSUM
! 449: H( J, K+2 ) = H( J, K+2 ) -
! 450: $ REFSUM*DCONJG( V( 2, M22 ) )
! 451: 90 CONTINUE
! 452: *
! 453: IF( ACCUM ) THEN
! 454: KMS = K - INCOL
! 455: DO 100 J = MAX( 1, KTOP-INCOL ), KDU
! 456: REFSUM = V( 1, M22 )*( U( J, KMS+1 )+V( 2, M22 )*
! 457: $ U( J, KMS+2 ) )
! 458: U( J, KMS+1 ) = U( J, KMS+1 ) - REFSUM
! 459: U( J, KMS+2 ) = U( J, KMS+2 ) -
! 460: $ REFSUM*DCONJG( V( 2, M22 ) )
! 461: 100 CONTINUE
! 462: ELSE IF( WANTZ ) THEN
! 463: DO 110 J = ILOZ, IHIZ
! 464: REFSUM = V( 1, M22 )*( Z( J, K+1 )+V( 2, M22 )*
! 465: $ Z( J, K+2 ) )
! 466: Z( J, K+1 ) = Z( J, K+1 ) - REFSUM
! 467: Z( J, K+2 ) = Z( J, K+2 ) -
! 468: $ REFSUM*DCONJG( V( 2, M22 ) )
! 469: 110 CONTINUE
! 470: END IF
! 471: END IF
! 472: *
! 473: * ==== Vigilant deflation check ====
! 474: *
! 475: MSTART = MTOP
! 476: IF( KRCOL+3*( MSTART-1 ).LT.KTOP )
! 477: $ MSTART = MSTART + 1
! 478: MEND = MBOT
! 479: IF( BMP22 )
! 480: $ MEND = MEND + 1
! 481: IF( KRCOL.EQ.KBOT-2 )
! 482: $ MEND = MEND + 1
! 483: DO 120 M = MSTART, MEND
! 484: K = MIN( KBOT-1, KRCOL+3*( M-1 ) )
! 485: *
! 486: * ==== The following convergence test requires that
! 487: * . the tradition small-compared-to-nearby-diagonals
! 488: * . criterion and the Ahues & Tisseur (LAWN 122, 1997)
! 489: * . criteria both be satisfied. The latter improves
! 490: * . accuracy in some examples. Falling back on an
! 491: * . alternate convergence criterion when TST1 or TST2
! 492: * . is zero (as done here) is traditional but probably
! 493: * . unnecessary. ====
! 494: *
! 495: IF( H( K+1, K ).NE.ZERO ) THEN
! 496: TST1 = CABS1( H( K, K ) ) + CABS1( H( K+1, K+1 ) )
! 497: IF( TST1.EQ.RZERO ) THEN
! 498: IF( K.GE.KTOP+1 )
! 499: $ TST1 = TST1 + CABS1( H( K, K-1 ) )
! 500: IF( K.GE.KTOP+2 )
! 501: $ TST1 = TST1 + CABS1( H( K, K-2 ) )
! 502: IF( K.GE.KTOP+3 )
! 503: $ TST1 = TST1 + CABS1( H( K, K-3 ) )
! 504: IF( K.LE.KBOT-2 )
! 505: $ TST1 = TST1 + CABS1( H( K+2, K+1 ) )
! 506: IF( K.LE.KBOT-3 )
! 507: $ TST1 = TST1 + CABS1( H( K+3, K+1 ) )
! 508: IF( K.LE.KBOT-4 )
! 509: $ TST1 = TST1 + CABS1( H( K+4, K+1 ) )
! 510: END IF
! 511: IF( CABS1( H( K+1, K ) ).LE.MAX( SMLNUM, ULP*TST1 ) )
! 512: $ THEN
! 513: H12 = MAX( CABS1( H( K+1, K ) ),
! 514: $ CABS1( H( K, K+1 ) ) )
! 515: H21 = MIN( CABS1( H( K+1, K ) ),
! 516: $ CABS1( H( K, K+1 ) ) )
! 517: H11 = MAX( CABS1( H( K+1, K+1 ) ),
! 518: $ CABS1( H( K, K )-H( K+1, K+1 ) ) )
! 519: H22 = MIN( CABS1( H( K+1, K+1 ) ),
! 520: $ CABS1( H( K, K )-H( K+1, K+1 ) ) )
! 521: SCL = H11 + H12
! 522: TST2 = H22*( H11 / SCL )
! 523: *
! 524: IF( TST2.EQ.RZERO .OR. H21*( H12 / SCL ).LE.
! 525: $ MAX( SMLNUM, ULP*TST2 ) )H( K+1, K ) = ZERO
! 526: END IF
! 527: END IF
! 528: 120 CONTINUE
! 529: *
! 530: * ==== Fill in the last row of each bulge. ====
! 531: *
! 532: MEND = MIN( NBMPS, ( KBOT-KRCOL-1 ) / 3 )
! 533: DO 130 M = MTOP, MEND
! 534: K = KRCOL + 3*( M-1 )
! 535: REFSUM = V( 1, M )*V( 3, M )*H( K+4, K+3 )
! 536: H( K+4, K+1 ) = -REFSUM
! 537: H( K+4, K+2 ) = -REFSUM*DCONJG( V( 2, M ) )
! 538: H( K+4, K+3 ) = H( K+4, K+3 ) -
! 539: $ REFSUM*DCONJG( V( 3, M ) )
! 540: 130 CONTINUE
! 541: *
! 542: * ==== End of near-the-diagonal bulge chase. ====
! 543: *
! 544: 140 CONTINUE
! 545: *
! 546: * ==== Use U (if accumulated) to update far-from-diagonal
! 547: * . entries in H. If required, use U to update Z as
! 548: * . well. ====
! 549: *
! 550: IF( ACCUM ) THEN
! 551: IF( WANTT ) THEN
! 552: JTOP = 1
! 553: JBOT = N
! 554: ELSE
! 555: JTOP = KTOP
! 556: JBOT = KBOT
! 557: END IF
! 558: IF( ( .NOT.BLK22 ) .OR. ( INCOL.LT.KTOP ) .OR.
! 559: $ ( NDCOL.GT.KBOT ) .OR. ( NS.LE.2 ) ) THEN
! 560: *
! 561: * ==== Updates not exploiting the 2-by-2 block
! 562: * . structure of U. K1 and NU keep track of
! 563: * . the location and size of U in the special
! 564: * . cases of introducing bulges and chasing
! 565: * . bulges off the bottom. In these special
! 566: * . cases and in case the number of shifts
! 567: * . is NS = 2, there is no 2-by-2 block
! 568: * . structure to exploit. ====
! 569: *
! 570: K1 = MAX( 1, KTOP-INCOL )
! 571: NU = ( KDU-MAX( 0, NDCOL-KBOT ) ) - K1 + 1
! 572: *
! 573: * ==== Horizontal Multiply ====
! 574: *
! 575: DO 150 JCOL = MIN( NDCOL, KBOT ) + 1, JBOT, NH
! 576: JLEN = MIN( NH, JBOT-JCOL+1 )
! 577: CALL ZGEMM( 'C', 'N', NU, JLEN, NU, ONE, U( K1, K1 ),
! 578: $ LDU, H( INCOL+K1, JCOL ), LDH, ZERO, WH,
! 579: $ LDWH )
! 580: CALL ZLACPY( 'ALL', NU, JLEN, WH, LDWH,
! 581: $ H( INCOL+K1, JCOL ), LDH )
! 582: 150 CONTINUE
! 583: *
! 584: * ==== Vertical multiply ====
! 585: *
! 586: DO 160 JROW = JTOP, MAX( KTOP, INCOL ) - 1, NV
! 587: JLEN = MIN( NV, MAX( KTOP, INCOL )-JROW )
! 588: CALL ZGEMM( 'N', 'N', JLEN, NU, NU, ONE,
! 589: $ H( JROW, INCOL+K1 ), LDH, U( K1, K1 ),
! 590: $ LDU, ZERO, WV, LDWV )
! 591: CALL ZLACPY( 'ALL', JLEN, NU, WV, LDWV,
! 592: $ H( JROW, INCOL+K1 ), LDH )
! 593: 160 CONTINUE
! 594: *
! 595: * ==== Z multiply (also vertical) ====
! 596: *
! 597: IF( WANTZ ) THEN
! 598: DO 170 JROW = ILOZ, IHIZ, NV
! 599: JLEN = MIN( NV, IHIZ-JROW+1 )
! 600: CALL ZGEMM( 'N', 'N', JLEN, NU, NU, ONE,
! 601: $ Z( JROW, INCOL+K1 ), LDZ, U( K1, K1 ),
! 602: $ LDU, ZERO, WV, LDWV )
! 603: CALL ZLACPY( 'ALL', JLEN, NU, WV, LDWV,
! 604: $ Z( JROW, INCOL+K1 ), LDZ )
! 605: 170 CONTINUE
! 606: END IF
! 607: ELSE
! 608: *
! 609: * ==== Updates exploiting U's 2-by-2 block structure.
! 610: * . (I2, I4, J2, J4 are the last rows and columns
! 611: * . of the blocks.) ====
! 612: *
! 613: I2 = ( KDU+1 ) / 2
! 614: I4 = KDU
! 615: J2 = I4 - I2
! 616: J4 = KDU
! 617: *
! 618: * ==== KZS and KNZ deal with the band of zeros
! 619: * . along the diagonal of one of the triangular
! 620: * . blocks. ====
! 621: *
! 622: KZS = ( J4-J2 ) - ( NS+1 )
! 623: KNZ = NS + 1
! 624: *
! 625: * ==== Horizontal multiply ====
! 626: *
! 627: DO 180 JCOL = MIN( NDCOL, KBOT ) + 1, JBOT, NH
! 628: JLEN = MIN( NH, JBOT-JCOL+1 )
! 629: *
! 630: * ==== Copy bottom of H to top+KZS of scratch ====
! 631: * (The first KZS rows get multiplied by zero.) ====
! 632: *
! 633: CALL ZLACPY( 'ALL', KNZ, JLEN, H( INCOL+1+J2, JCOL ),
! 634: $ LDH, WH( KZS+1, 1 ), LDWH )
! 635: *
! 636: * ==== Multiply by U21' ====
! 637: *
! 638: CALL ZLASET( 'ALL', KZS, JLEN, ZERO, ZERO, WH, LDWH )
! 639: CALL ZTRMM( 'L', 'U', 'C', 'N', KNZ, JLEN, ONE,
! 640: $ U( J2+1, 1+KZS ), LDU, WH( KZS+1, 1 ),
! 641: $ LDWH )
! 642: *
! 643: * ==== Multiply top of H by U11' ====
! 644: *
! 645: CALL ZGEMM( 'C', 'N', I2, JLEN, J2, ONE, U, LDU,
! 646: $ H( INCOL+1, JCOL ), LDH, ONE, WH, LDWH )
! 647: *
! 648: * ==== Copy top of H to bottom of WH ====
! 649: *
! 650: CALL ZLACPY( 'ALL', J2, JLEN, H( INCOL+1, JCOL ), LDH,
! 651: $ WH( I2+1, 1 ), LDWH )
! 652: *
! 653: * ==== Multiply by U21' ====
! 654: *
! 655: CALL ZTRMM( 'L', 'L', 'C', 'N', J2, JLEN, ONE,
! 656: $ U( 1, I2+1 ), LDU, WH( I2+1, 1 ), LDWH )
! 657: *
! 658: * ==== Multiply by U22 ====
! 659: *
! 660: CALL ZGEMM( 'C', 'N', I4-I2, JLEN, J4-J2, ONE,
! 661: $ U( J2+1, I2+1 ), LDU,
! 662: $ H( INCOL+1+J2, JCOL ), LDH, ONE,
! 663: $ WH( I2+1, 1 ), LDWH )
! 664: *
! 665: * ==== Copy it back ====
! 666: *
! 667: CALL ZLACPY( 'ALL', KDU, JLEN, WH, LDWH,
! 668: $ H( INCOL+1, JCOL ), LDH )
! 669: 180 CONTINUE
! 670: *
! 671: * ==== Vertical multiply ====
! 672: *
! 673: DO 190 JROW = JTOP, MAX( INCOL, KTOP ) - 1, NV
! 674: JLEN = MIN( NV, MAX( INCOL, KTOP )-JROW )
! 675: *
! 676: * ==== Copy right of H to scratch (the first KZS
! 677: * . columns get multiplied by zero) ====
! 678: *
! 679: CALL ZLACPY( 'ALL', JLEN, KNZ, H( JROW, INCOL+1+J2 ),
! 680: $ LDH, WV( 1, 1+KZS ), LDWV )
! 681: *
! 682: * ==== Multiply by U21 ====
! 683: *
! 684: CALL ZLASET( 'ALL', JLEN, KZS, ZERO, ZERO, WV, LDWV )
! 685: CALL ZTRMM( 'R', 'U', 'N', 'N', JLEN, KNZ, ONE,
! 686: $ U( J2+1, 1+KZS ), LDU, WV( 1, 1+KZS ),
! 687: $ LDWV )
! 688: *
! 689: * ==== Multiply by U11 ====
! 690: *
! 691: CALL ZGEMM( 'N', 'N', JLEN, I2, J2, ONE,
! 692: $ H( JROW, INCOL+1 ), LDH, U, LDU, ONE, WV,
! 693: $ LDWV )
! 694: *
! 695: * ==== Copy left of H to right of scratch ====
! 696: *
! 697: CALL ZLACPY( 'ALL', JLEN, J2, H( JROW, INCOL+1 ), LDH,
! 698: $ WV( 1, 1+I2 ), LDWV )
! 699: *
! 700: * ==== Multiply by U21 ====
! 701: *
! 702: CALL ZTRMM( 'R', 'L', 'N', 'N', JLEN, I4-I2, ONE,
! 703: $ U( 1, I2+1 ), LDU, WV( 1, 1+I2 ), LDWV )
! 704: *
! 705: * ==== Multiply by U22 ====
! 706: *
! 707: CALL ZGEMM( 'N', 'N', JLEN, I4-I2, J4-J2, ONE,
! 708: $ H( JROW, INCOL+1+J2 ), LDH,
! 709: $ U( J2+1, I2+1 ), LDU, ONE, WV( 1, 1+I2 ),
! 710: $ LDWV )
! 711: *
! 712: * ==== Copy it back ====
! 713: *
! 714: CALL ZLACPY( 'ALL', JLEN, KDU, WV, LDWV,
! 715: $ H( JROW, INCOL+1 ), LDH )
! 716: 190 CONTINUE
! 717: *
! 718: * ==== Multiply Z (also vertical) ====
! 719: *
! 720: IF( WANTZ ) THEN
! 721: DO 200 JROW = ILOZ, IHIZ, NV
! 722: JLEN = MIN( NV, IHIZ-JROW+1 )
! 723: *
! 724: * ==== Copy right of Z to left of scratch (first
! 725: * . KZS columns get multiplied by zero) ====
! 726: *
! 727: CALL ZLACPY( 'ALL', JLEN, KNZ,
! 728: $ Z( JROW, INCOL+1+J2 ), LDZ,
! 729: $ WV( 1, 1+KZS ), LDWV )
! 730: *
! 731: * ==== Multiply by U12 ====
! 732: *
! 733: CALL ZLASET( 'ALL', JLEN, KZS, ZERO, ZERO, WV,
! 734: $ LDWV )
! 735: CALL ZTRMM( 'R', 'U', 'N', 'N', JLEN, KNZ, ONE,
! 736: $ U( J2+1, 1+KZS ), LDU, WV( 1, 1+KZS ),
! 737: $ LDWV )
! 738: *
! 739: * ==== Multiply by U11 ====
! 740: *
! 741: CALL ZGEMM( 'N', 'N', JLEN, I2, J2, ONE,
! 742: $ Z( JROW, INCOL+1 ), LDZ, U, LDU, ONE,
! 743: $ WV, LDWV )
! 744: *
! 745: * ==== Copy left of Z to right of scratch ====
! 746: *
! 747: CALL ZLACPY( 'ALL', JLEN, J2, Z( JROW, INCOL+1 ),
! 748: $ LDZ, WV( 1, 1+I2 ), LDWV )
! 749: *
! 750: * ==== Multiply by U21 ====
! 751: *
! 752: CALL ZTRMM( 'R', 'L', 'N', 'N', JLEN, I4-I2, ONE,
! 753: $ U( 1, I2+1 ), LDU, WV( 1, 1+I2 ),
! 754: $ LDWV )
! 755: *
! 756: * ==== Multiply by U22 ====
! 757: *
! 758: CALL ZGEMM( 'N', 'N', JLEN, I4-I2, J4-J2, ONE,
! 759: $ Z( JROW, INCOL+1+J2 ), LDZ,
! 760: $ U( J2+1, I2+1 ), LDU, ONE,
! 761: $ WV( 1, 1+I2 ), LDWV )
! 762: *
! 763: * ==== Copy the result back to Z ====
! 764: *
! 765: CALL ZLACPY( 'ALL', JLEN, KDU, WV, LDWV,
! 766: $ Z( JROW, INCOL+1 ), LDZ )
! 767: 200 CONTINUE
! 768: END IF
! 769: END IF
! 770: END IF
! 771: 210 CONTINUE
! 772: *
! 773: * ==== End of ZLAQR5 ====
! 774: *
! 775: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlaqr5.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack