![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
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.3.0) -- 6: * Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd.. 7: * November 2010 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 ) THEN 445: IF ( V( 1, M22 ).NE.ZERO ) THEN 446: DO 90 J = JTOP, MIN( KBOT, K+3 ) 447: REFSUM = V( 1, M22 )*( H( J, K+1 )+V( 2, M22 )* 448: $ H( J, K+2 ) ) 449: H( J, K+1 ) = H( J, K+1 ) - REFSUM 450: H( J, K+2 ) = H( J, K+2 ) - 451: $ REFSUM*DCONJG( V( 2, M22 ) ) 452: 90 CONTINUE 453: * 454: IF( ACCUM ) THEN 455: KMS = K - INCOL 456: DO 100 J = MAX( 1, KTOP-INCOL ), KDU 457: REFSUM = V( 1, M22 )*( U( J, KMS+1 )+ 458: $ V( 2, M22 )*U( J, KMS+2 ) ) 459: U( J, KMS+1 ) = U( J, KMS+1 ) - REFSUM 460: U( J, KMS+2 ) = U( J, KMS+2 ) - 461: $ REFSUM*DCONJG( V( 2, M22 ) ) 462: 100 CONTINUE 463: ELSE IF( WANTZ ) THEN 464: DO 110 J = ILOZ, IHIZ 465: REFSUM = V( 1, M22 )*( Z( J, K+1 )+V( 2, M22 )* 466: $ Z( J, K+2 ) ) 467: Z( J, K+1 ) = Z( J, K+1 ) - REFSUM 468: Z( J, K+2 ) = Z( J, K+2 ) - 469: $ REFSUM*DCONJG( V( 2, M22 ) ) 470: 110 CONTINUE 471: END IF 472: END IF 473: END IF 474: * 475: * ==== Vigilant deflation check ==== 476: * 477: MSTART = MTOP 478: IF( KRCOL+3*( MSTART-1 ).LT.KTOP ) 479: $ MSTART = MSTART + 1 480: MEND = MBOT 481: IF( BMP22 ) 482: $ MEND = MEND + 1 483: IF( KRCOL.EQ.KBOT-2 ) 484: $ MEND = MEND + 1 485: DO 120 M = MSTART, MEND 486: K = MIN( KBOT-1, KRCOL+3*( M-1 ) ) 487: * 488: * ==== The following convergence test requires that 489: * . the tradition small-compared-to-nearby-diagonals 490: * . criterion and the Ahues & Tisseur (LAWN 122, 1997) 491: * . criteria both be satisfied. The latter improves 492: * . accuracy in some examples. Falling back on an 493: * . alternate convergence criterion when TST1 or TST2 494: * . is zero (as done here) is traditional but probably 495: * . unnecessary. ==== 496: * 497: IF( H( K+1, K ).NE.ZERO ) THEN 498: TST1 = CABS1( H( K, K ) ) + CABS1( H( K+1, K+1 ) ) 499: IF( TST1.EQ.RZERO ) THEN 500: IF( K.GE.KTOP+1 ) 501: $ TST1 = TST1 + CABS1( H( K, K-1 ) ) 502: IF( K.GE.KTOP+2 ) 503: $ TST1 = TST1 + CABS1( H( K, K-2 ) ) 504: IF( K.GE.KTOP+3 ) 505: $ TST1 = TST1 + CABS1( H( K, K-3 ) ) 506: IF( K.LE.KBOT-2 ) 507: $ TST1 = TST1 + CABS1( H( K+2, K+1 ) ) 508: IF( K.LE.KBOT-3 ) 509: $ TST1 = TST1 + CABS1( H( K+3, K+1 ) ) 510: IF( K.LE.KBOT-4 ) 511: $ TST1 = TST1 + CABS1( H( K+4, K+1 ) ) 512: END IF 513: IF( CABS1( H( K+1, K ) ).LE.MAX( SMLNUM, ULP*TST1 ) ) 514: $ THEN 515: H12 = MAX( CABS1( H( K+1, K ) ), 516: $ CABS1( H( K, K+1 ) ) ) 517: H21 = MIN( CABS1( H( K+1, K ) ), 518: $ CABS1( H( K, K+1 ) ) ) 519: H11 = MAX( CABS1( H( K+1, K+1 ) ), 520: $ CABS1( H( K, K )-H( K+1, K+1 ) ) ) 521: H22 = MIN( CABS1( H( K+1, K+1 ) ), 522: $ CABS1( H( K, K )-H( K+1, K+1 ) ) ) 523: SCL = H11 + H12 524: TST2 = H22*( H11 / SCL ) 525: * 526: IF( TST2.EQ.RZERO .OR. H21*( H12 / SCL ).LE. 527: $ MAX( SMLNUM, ULP*TST2 ) )H( K+1, K ) = ZERO 528: END IF 529: END IF 530: 120 CONTINUE 531: * 532: * ==== Fill in the last row of each bulge. ==== 533: * 534: MEND = MIN( NBMPS, ( KBOT-KRCOL-1 ) / 3 ) 535: DO 130 M = MTOP, MEND 536: K = KRCOL + 3*( M-1 ) 537: REFSUM = V( 1, M )*V( 3, M )*H( K+4, K+3 ) 538: H( K+4, K+1 ) = -REFSUM 539: H( K+4, K+2 ) = -REFSUM*DCONJG( V( 2, M ) ) 540: H( K+4, K+3 ) = H( K+4, K+3 ) - 541: $ REFSUM*DCONJG( V( 3, M ) ) 542: 130 CONTINUE 543: * 544: * ==== End of near-the-diagonal bulge chase. ==== 545: * 546: 140 CONTINUE 547: * 548: * ==== Use U (if accumulated) to update far-from-diagonal 549: * . entries in H. If required, use U to update Z as 550: * . well. ==== 551: * 552: IF( ACCUM ) THEN 553: IF( WANTT ) THEN 554: JTOP = 1 555: JBOT = N 556: ELSE 557: JTOP = KTOP 558: JBOT = KBOT 559: END IF 560: IF( ( .NOT.BLK22 ) .OR. ( INCOL.LT.KTOP ) .OR. 561: $ ( NDCOL.GT.KBOT ) .OR. ( NS.LE.2 ) ) THEN 562: * 563: * ==== Updates not exploiting the 2-by-2 block 564: * . structure of U. K1 and NU keep track of 565: * . the location and size of U in the special 566: * . cases of introducing bulges and chasing 567: * . bulges off the bottom. In these special 568: * . cases and in case the number of shifts 569: * . is NS = 2, there is no 2-by-2 block 570: * . structure to exploit. ==== 571: * 572: K1 = MAX( 1, KTOP-INCOL ) 573: NU = ( KDU-MAX( 0, NDCOL-KBOT ) ) - K1 + 1 574: * 575: * ==== Horizontal Multiply ==== 576: * 577: DO 150 JCOL = MIN( NDCOL, KBOT ) + 1, JBOT, NH 578: JLEN = MIN( NH, JBOT-JCOL+1 ) 579: CALL ZGEMM( 'C', 'N', NU, JLEN, NU, ONE, U( K1, K1 ), 580: $ LDU, H( INCOL+K1, JCOL ), LDH, ZERO, WH, 581: $ LDWH ) 582: CALL ZLACPY( 'ALL', NU, JLEN, WH, LDWH, 583: $ H( INCOL+K1, JCOL ), LDH ) 584: 150 CONTINUE 585: * 586: * ==== Vertical multiply ==== 587: * 588: DO 160 JROW = JTOP, MAX( KTOP, INCOL ) - 1, NV 589: JLEN = MIN( NV, MAX( KTOP, INCOL )-JROW ) 590: CALL ZGEMM( 'N', 'N', JLEN, NU, NU, ONE, 591: $ H( JROW, INCOL+K1 ), LDH, U( K1, K1 ), 592: $ LDU, ZERO, WV, LDWV ) 593: CALL ZLACPY( 'ALL', JLEN, NU, WV, LDWV, 594: $ H( JROW, INCOL+K1 ), LDH ) 595: 160 CONTINUE 596: * 597: * ==== Z multiply (also vertical) ==== 598: * 599: IF( WANTZ ) THEN 600: DO 170 JROW = ILOZ, IHIZ, NV 601: JLEN = MIN( NV, IHIZ-JROW+1 ) 602: CALL ZGEMM( 'N', 'N', JLEN, NU, NU, ONE, 603: $ Z( JROW, INCOL+K1 ), LDZ, U( K1, K1 ), 604: $ LDU, ZERO, WV, LDWV ) 605: CALL ZLACPY( 'ALL', JLEN, NU, WV, LDWV, 606: $ Z( JROW, INCOL+K1 ), LDZ ) 607: 170 CONTINUE 608: END IF 609: ELSE 610: * 611: * ==== Updates exploiting U's 2-by-2 block structure. 612: * . (I2, I4, J2, J4 are the last rows and columns 613: * . of the blocks.) ==== 614: * 615: I2 = ( KDU+1 ) / 2 616: I4 = KDU 617: J2 = I4 - I2 618: J4 = KDU 619: * 620: * ==== KZS and KNZ deal with the band of zeros 621: * . along the diagonal of one of the triangular 622: * . blocks. ==== 623: * 624: KZS = ( J4-J2 ) - ( NS+1 ) 625: KNZ = NS + 1 626: * 627: * ==== Horizontal multiply ==== 628: * 629: DO 180 JCOL = MIN( NDCOL, KBOT ) + 1, JBOT, NH 630: JLEN = MIN( NH, JBOT-JCOL+1 ) 631: * 632: * ==== Copy bottom of H to top+KZS of scratch ==== 633: * (The first KZS rows get multiplied by zero.) ==== 634: * 635: CALL ZLACPY( 'ALL', KNZ, JLEN, H( INCOL+1+J2, JCOL ), 636: $ LDH, WH( KZS+1, 1 ), LDWH ) 637: * 638: * ==== Multiply by U21' ==== 639: * 640: CALL ZLASET( 'ALL', KZS, JLEN, ZERO, ZERO, WH, LDWH ) 641: CALL ZTRMM( 'L', 'U', 'C', 'N', KNZ, JLEN, ONE, 642: $ U( J2+1, 1+KZS ), LDU, WH( KZS+1, 1 ), 643: $ LDWH ) 644: * 645: * ==== Multiply top of H by U11' ==== 646: * 647: CALL ZGEMM( 'C', 'N', I2, JLEN, J2, ONE, U, LDU, 648: $ H( INCOL+1, JCOL ), LDH, ONE, WH, LDWH ) 649: * 650: * ==== Copy top of H to bottom of WH ==== 651: * 652: CALL ZLACPY( 'ALL', J2, JLEN, H( INCOL+1, JCOL ), LDH, 653: $ WH( I2+1, 1 ), LDWH ) 654: * 655: * ==== Multiply by U21' ==== 656: * 657: CALL ZTRMM( 'L', 'L', 'C', 'N', J2, JLEN, ONE, 658: $ U( 1, I2+1 ), LDU, WH( I2+1, 1 ), LDWH ) 659: * 660: * ==== Multiply by U22 ==== 661: * 662: CALL ZGEMM( 'C', 'N', I4-I2, JLEN, J4-J2, ONE, 663: $ U( J2+1, I2+1 ), LDU, 664: $ H( INCOL+1+J2, JCOL ), LDH, ONE, 665: $ WH( I2+1, 1 ), LDWH ) 666: * 667: * ==== Copy it back ==== 668: * 669: CALL ZLACPY( 'ALL', KDU, JLEN, WH, LDWH, 670: $ H( INCOL+1, JCOL ), LDH ) 671: 180 CONTINUE 672: * 673: * ==== Vertical multiply ==== 674: * 675: DO 190 JROW = JTOP, MAX( INCOL, KTOP ) - 1, NV 676: JLEN = MIN( NV, MAX( INCOL, KTOP )-JROW ) 677: * 678: * ==== Copy right of H to scratch (the first KZS 679: * . columns get multiplied by zero) ==== 680: * 681: CALL ZLACPY( 'ALL', JLEN, KNZ, H( JROW, INCOL+1+J2 ), 682: $ LDH, WV( 1, 1+KZS ), LDWV ) 683: * 684: * ==== Multiply by U21 ==== 685: * 686: CALL ZLASET( 'ALL', JLEN, KZS, ZERO, ZERO, WV, LDWV ) 687: CALL ZTRMM( 'R', 'U', 'N', 'N', JLEN, KNZ, ONE, 688: $ U( J2+1, 1+KZS ), LDU, WV( 1, 1+KZS ), 689: $ LDWV ) 690: * 691: * ==== Multiply by U11 ==== 692: * 693: CALL ZGEMM( 'N', 'N', JLEN, I2, J2, ONE, 694: $ H( JROW, INCOL+1 ), LDH, U, LDU, ONE, WV, 695: $ LDWV ) 696: * 697: * ==== Copy left of H to right of scratch ==== 698: * 699: CALL ZLACPY( 'ALL', JLEN, J2, H( JROW, INCOL+1 ), LDH, 700: $ WV( 1, 1+I2 ), LDWV ) 701: * 702: * ==== Multiply by U21 ==== 703: * 704: CALL ZTRMM( 'R', 'L', 'N', 'N', JLEN, I4-I2, ONE, 705: $ U( 1, I2+1 ), LDU, WV( 1, 1+I2 ), LDWV ) 706: * 707: * ==== Multiply by U22 ==== 708: * 709: CALL ZGEMM( 'N', 'N', JLEN, I4-I2, J4-J2, ONE, 710: $ H( JROW, INCOL+1+J2 ), LDH, 711: $ U( J2+1, I2+1 ), LDU, ONE, WV( 1, 1+I2 ), 712: $ LDWV ) 713: * 714: * ==== Copy it back ==== 715: * 716: CALL ZLACPY( 'ALL', JLEN, KDU, WV, LDWV, 717: $ H( JROW, INCOL+1 ), LDH ) 718: 190 CONTINUE 719: * 720: * ==== Multiply Z (also vertical) ==== 721: * 722: IF( WANTZ ) THEN 723: DO 200 JROW = ILOZ, IHIZ, NV 724: JLEN = MIN( NV, IHIZ-JROW+1 ) 725: * 726: * ==== Copy right of Z to left of scratch (first 727: * . KZS columns get multiplied by zero) ==== 728: * 729: CALL ZLACPY( 'ALL', JLEN, KNZ, 730: $ Z( JROW, INCOL+1+J2 ), LDZ, 731: $ WV( 1, 1+KZS ), LDWV ) 732: * 733: * ==== Multiply by U12 ==== 734: * 735: CALL ZLASET( 'ALL', JLEN, KZS, ZERO, ZERO, WV, 736: $ LDWV ) 737: CALL ZTRMM( 'R', 'U', 'N', 'N', JLEN, KNZ, ONE, 738: $ U( J2+1, 1+KZS ), LDU, WV( 1, 1+KZS ), 739: $ LDWV ) 740: * 741: * ==== Multiply by U11 ==== 742: * 743: CALL ZGEMM( 'N', 'N', JLEN, I2, J2, ONE, 744: $ Z( JROW, INCOL+1 ), LDZ, U, LDU, ONE, 745: $ WV, LDWV ) 746: * 747: * ==== Copy left of Z to right of scratch ==== 748: * 749: CALL ZLACPY( 'ALL', JLEN, J2, Z( JROW, INCOL+1 ), 750: $ LDZ, WV( 1, 1+I2 ), LDWV ) 751: * 752: * ==== Multiply by U21 ==== 753: * 754: CALL ZTRMM( 'R', 'L', 'N', 'N', JLEN, I4-I2, ONE, 755: $ U( 1, I2+1 ), LDU, WV( 1, 1+I2 ), 756: $ LDWV ) 757: * 758: * ==== Multiply by U22 ==== 759: * 760: CALL ZGEMM( 'N', 'N', JLEN, I4-I2, J4-J2, ONE, 761: $ Z( JROW, INCOL+1+J2 ), LDZ, 762: $ U( J2+1, I2+1 ), LDU, ONE, 763: $ WV( 1, 1+I2 ), LDWV ) 764: * 765: * ==== Copy the result back to Z ==== 766: * 767: CALL ZLACPY( 'ALL', JLEN, KDU, WV, LDWV, 768: $ Z( JROW, INCOL+1 ), LDZ ) 769: 200 CONTINUE 770: END IF 771: END IF 772: END IF 773: 210 CONTINUE 774: * 775: * ==== End of ZLAQR5 ==== 776: * 777: END