![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.2.2.
1: SUBROUTINE ZHBGST( VECT, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, X, 2: $ LDX, WORK, RWORK, INFO ) 3: * 4: * -- LAPACK routine (version 3.2.2) -- 5: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 7: * June 2010 8: * 9: * .. Scalar Arguments .. 10: CHARACTER UPLO, VECT 11: INTEGER INFO, KA, KB, LDAB, LDBB, LDX, N 12: * .. 13: * .. Array Arguments .. 14: DOUBLE PRECISION RWORK( * ) 15: COMPLEX*16 AB( LDAB, * ), BB( LDBB, * ), WORK( * ), 16: $ X( LDX, * ) 17: * .. 18: * 19: * Purpose 20: * ======= 21: * 22: * ZHBGST reduces a complex Hermitian-definite banded generalized 23: * eigenproblem A*x = lambda*B*x to standard form C*y = lambda*y, 24: * such that C has the same bandwidth as A. 25: * 26: * B must have been previously factorized as S**H*S by ZPBSTF, using a 27: * split Cholesky factorization. A is overwritten by C = X**H*A*X, where 28: * X = S**(-1)*Q and Q is a unitary matrix chosen to preserve the 29: * bandwidth of A. 30: * 31: * Arguments 32: * ========= 33: * 34: * VECT (input) CHARACTER*1 35: * = 'N': do not form the transformation matrix X; 36: * = 'V': form X. 37: * 38: * UPLO (input) CHARACTER*1 39: * = 'U': Upper triangle of A is stored; 40: * = 'L': Lower triangle of A is stored. 41: * 42: * N (input) INTEGER 43: * The order of the matrices A and B. N >= 0. 44: * 45: * KA (input) INTEGER 46: * The number of superdiagonals of the matrix A if UPLO = 'U', 47: * or the number of subdiagonals if UPLO = 'L'. KA >= 0. 48: * 49: * KB (input) INTEGER 50: * The number of superdiagonals of the matrix B if UPLO = 'U', 51: * or the number of subdiagonals if UPLO = 'L'. KA >= KB >= 0. 52: * 53: * AB (input/output) COMPLEX*16 array, dimension (LDAB,N) 54: * On entry, the upper or lower triangle of the Hermitian band 55: * matrix A, stored in the first ka+1 rows of the array. The 56: * j-th column of A is stored in the j-th column of the array AB 57: * as follows: 58: * if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j; 59: * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+ka). 60: * 61: * On exit, the transformed matrix X**H*A*X, stored in the same 62: * format as A. 63: * 64: * LDAB (input) INTEGER 65: * The leading dimension of the array AB. LDAB >= KA+1. 66: * 67: * BB (input) COMPLEX*16 array, dimension (LDBB,N) 68: * The banded factor S from the split Cholesky factorization of 69: * B, as returned by ZPBSTF, stored in the first kb+1 rows of 70: * the array. 71: * 72: * LDBB (input) INTEGER 73: * The leading dimension of the array BB. LDBB >= KB+1. 74: * 75: * X (output) COMPLEX*16 array, dimension (LDX,N) 76: * If VECT = 'V', the n-by-n matrix X. 77: * If VECT = 'N', the array X is not referenced. 78: * 79: * LDX (input) INTEGER 80: * The leading dimension of the array X. 81: * LDX >= max(1,N) if VECT = 'V'; LDX >= 1 otherwise. 82: * 83: * WORK (workspace) COMPLEX*16 array, dimension (N) 84: * 85: * RWORK (workspace) DOUBLE PRECISION array, dimension (N) 86: * 87: * INFO (output) INTEGER 88: * = 0: successful exit 89: * < 0: if INFO = -i, the i-th argument had an illegal value. 90: * 91: * ===================================================================== 92: * 93: * .. Parameters .. 94: COMPLEX*16 CZERO, CONE 95: DOUBLE PRECISION ONE 96: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 97: $ CONE = ( 1.0D+0, 0.0D+0 ), ONE = 1.0D+0 ) 98: * .. 99: * .. Local Scalars .. 100: LOGICAL UPDATE, UPPER, WANTX 101: INTEGER I, I0, I1, I2, INCA, J, J1, J1T, J2, J2T, K, 102: $ KA1, KB1, KBT, L, M, NR, NRT, NX 103: DOUBLE PRECISION BII 104: COMPLEX*16 RA, RA1, T 105: * .. 106: * .. External Functions .. 107: LOGICAL LSAME 108: EXTERNAL LSAME 109: * .. 110: * .. External Subroutines .. 111: EXTERNAL XERBLA, ZDSCAL, ZGERC, ZGERU, ZLACGV, ZLAR2V, 112: $ ZLARGV, ZLARTG, ZLARTV, ZLASET, ZROT 113: * .. 114: * .. Intrinsic Functions .. 115: INTRINSIC DBLE, DCONJG, MAX, MIN 116: * .. 117: * .. Executable Statements .. 118: * 119: * Test the input parameters 120: * 121: WANTX = LSAME( VECT, 'V' ) 122: UPPER = LSAME( UPLO, 'U' ) 123: KA1 = KA + 1 124: KB1 = KB + 1 125: INFO = 0 126: IF( .NOT.WANTX .AND. .NOT.LSAME( VECT, 'N' ) ) THEN 127: INFO = -1 128: ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 129: INFO = -2 130: ELSE IF( N.LT.0 ) THEN 131: INFO = -3 132: ELSE IF( KA.LT.0 ) THEN 133: INFO = -4 134: ELSE IF( KB.LT.0 .OR. KB.GT.KA ) THEN 135: INFO = -5 136: ELSE IF( LDAB.LT.KA+1 ) THEN 137: INFO = -7 138: ELSE IF( LDBB.LT.KB+1 ) THEN 139: INFO = -9 140: ELSE IF( LDX.LT.1 .OR. WANTX .AND. LDX.LT.MAX( 1, N ) ) THEN 141: INFO = -11 142: END IF 143: IF( INFO.NE.0 ) THEN 144: CALL XERBLA( 'ZHBGST', -INFO ) 145: RETURN 146: END IF 147: * 148: * Quick return if possible 149: * 150: IF( N.EQ.0 ) 151: $ RETURN 152: * 153: INCA = LDAB*KA1 154: * 155: * Initialize X to the unit matrix, if needed 156: * 157: IF( WANTX ) 158: $ CALL ZLASET( 'Full', N, N, CZERO, CONE, X, LDX ) 159: * 160: * Set M to the splitting point m. It must be the same value as is 161: * used in ZPBSTF. The chosen value allows the arrays WORK and RWORK 162: * to be of dimension (N). 163: * 164: M = ( N+KB ) / 2 165: * 166: * The routine works in two phases, corresponding to the two halves 167: * of the split Cholesky factorization of B as S**H*S where 168: * 169: * S = ( U ) 170: * ( M L ) 171: * 172: * with U upper triangular of order m, and L lower triangular of 173: * order n-m. S has the same bandwidth as B. 174: * 175: * S is treated as a product of elementary matrices: 176: * 177: * S = S(m)*S(m-1)*...*S(2)*S(1)*S(m+1)*S(m+2)*...*S(n-1)*S(n) 178: * 179: * where S(i) is determined by the i-th row of S. 180: * 181: * In phase 1, the index i takes the values n, n-1, ... , m+1; 182: * in phase 2, it takes the values 1, 2, ... , m. 183: * 184: * For each value of i, the current matrix A is updated by forming 185: * inv(S(i))**H*A*inv(S(i)). This creates a triangular bulge outside 186: * the band of A. The bulge is then pushed down toward the bottom of 187: * A in phase 1, and up toward the top of A in phase 2, by applying 188: * plane rotations. 189: * 190: * There are kb*(kb+1)/2 elements in the bulge, but at most 2*kb-1 191: * of them are linearly independent, so annihilating a bulge requires 192: * only 2*kb-1 plane rotations. The rotations are divided into a 1st 193: * set of kb-1 rotations, and a 2nd set of kb rotations. 194: * 195: * Wherever possible, rotations are generated and applied in vector 196: * operations of length NR between the indices J1 and J2 (sometimes 197: * replaced by modified values NRT, J1T or J2T). 198: * 199: * The real cosines and complex sines of the rotations are stored in 200: * the arrays RWORK and WORK, those of the 1st set in elements 201: * 2:m-kb-1, and those of the 2nd set in elements m-kb+1:n. 202: * 203: * The bulges are not formed explicitly; nonzero elements outside the 204: * band are created only when they are required for generating new 205: * rotations; they are stored in the array WORK, in positions where 206: * they are later overwritten by the sines of the rotations which 207: * annihilate them. 208: * 209: * **************************** Phase 1 ***************************** 210: * 211: * The logical structure of this phase is: 212: * 213: * UPDATE = .TRUE. 214: * DO I = N, M + 1, -1 215: * use S(i) to update A and create a new bulge 216: * apply rotations to push all bulges KA positions downward 217: * END DO 218: * UPDATE = .FALSE. 219: * DO I = M + KA + 1, N - 1 220: * apply rotations to push all bulges KA positions downward 221: * END DO 222: * 223: * To avoid duplicating code, the two loops are merged. 224: * 225: UPDATE = .TRUE. 226: I = N + 1 227: 10 CONTINUE 228: IF( UPDATE ) THEN 229: I = I - 1 230: KBT = MIN( KB, I-1 ) 231: I0 = I - 1 232: I1 = MIN( N, I+KA ) 233: I2 = I - KBT + KA1 234: IF( I.LT.M+1 ) THEN 235: UPDATE = .FALSE. 236: I = I + 1 237: I0 = M 238: IF( KA.EQ.0 ) 239: $ GO TO 480 240: GO TO 10 241: END IF 242: ELSE 243: I = I + KA 244: IF( I.GT.N-1 ) 245: $ GO TO 480 246: END IF 247: * 248: IF( UPPER ) THEN 249: * 250: * Transform A, working with the upper triangle 251: * 252: IF( UPDATE ) THEN 253: * 254: * Form inv(S(i))**H * A * inv(S(i)) 255: * 256: BII = DBLE( BB( KB1, I ) ) 257: AB( KA1, I ) = ( DBLE( AB( KA1, I ) ) / BII ) / BII 258: DO 20 J = I + 1, I1 259: AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII 260: 20 CONTINUE 261: DO 30 J = MAX( 1, I-KA ), I - 1 262: AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII 263: 30 CONTINUE 264: DO 60 K = I - KBT, I - 1 265: DO 40 J = I - KBT, K 266: AB( J-K+KA1, K ) = AB( J-K+KA1, K ) - 267: $ BB( J-I+KB1, I )* 268: $ DCONJG( AB( K-I+KA1, I ) ) - 269: $ DCONJG( BB( K-I+KB1, I ) )* 270: $ AB( J-I+KA1, I ) + 271: $ DBLE( AB( KA1, I ) )* 272: $ BB( J-I+KB1, I )* 273: $ DCONJG( BB( K-I+KB1, I ) ) 274: 40 CONTINUE 275: DO 50 J = MAX( 1, I-KA ), I - KBT - 1 276: AB( J-K+KA1, K ) = AB( J-K+KA1, K ) - 277: $ DCONJG( BB( K-I+KB1, I ) )* 278: $ AB( J-I+KA1, I ) 279: 50 CONTINUE 280: 60 CONTINUE 281: DO 80 J = I, I1 282: DO 70 K = MAX( J-KA, I-KBT ), I - 1 283: AB( K-J+KA1, J ) = AB( K-J+KA1, J ) - 284: $ BB( K-I+KB1, I )*AB( I-J+KA1, J ) 285: 70 CONTINUE 286: 80 CONTINUE 287: * 288: IF( WANTX ) THEN 289: * 290: * post-multiply X by inv(S(i)) 291: * 292: CALL ZDSCAL( N-M, ONE / BII, X( M+1, I ), 1 ) 293: IF( KBT.GT.0 ) 294: $ CALL ZGERC( N-M, KBT, -CONE, X( M+1, I ), 1, 295: $ BB( KB1-KBT, I ), 1, X( M+1, I-KBT ), 296: $ LDX ) 297: END IF 298: * 299: * store a(i,i1) in RA1 for use in next loop over K 300: * 301: RA1 = AB( I-I1+KA1, I1 ) 302: END IF 303: * 304: * Generate and apply vectors of rotations to chase all the 305: * existing bulges KA positions down toward the bottom of the 306: * band 307: * 308: DO 130 K = 1, KB - 1 309: IF( UPDATE ) THEN 310: * 311: * Determine the rotations which would annihilate the bulge 312: * which has in theory just been created 313: * 314: IF( I-K+KA.LT.N .AND. I-K.GT.1 ) THEN 315: * 316: * generate rotation to annihilate a(i,i-k+ka+1) 317: * 318: CALL ZLARTG( AB( K+1, I-K+KA ), RA1, 319: $ RWORK( I-K+KA-M ), WORK( I-K+KA-M ), RA ) 320: * 321: * create nonzero element a(i-k,i-k+ka+1) outside the 322: * band and store it in WORK(i-k) 323: * 324: T = -BB( KB1-K, I )*RA1 325: WORK( I-K ) = RWORK( I-K+KA-M )*T - 326: $ DCONJG( WORK( I-K+KA-M ) )* 327: $ AB( 1, I-K+KA ) 328: AB( 1, I-K+KA ) = WORK( I-K+KA-M )*T + 329: $ RWORK( I-K+KA-M )*AB( 1, I-K+KA ) 330: RA1 = RA 331: END IF 332: END IF 333: J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1 334: NR = ( N-J2+KA ) / KA1 335: J1 = J2 + ( NR-1 )*KA1 336: IF( UPDATE ) THEN 337: J2T = MAX( J2, I+2*KA-K+1 ) 338: ELSE 339: J2T = J2 340: END IF 341: NRT = ( N-J2T+KA ) / KA1 342: DO 90 J = J2T, J1, KA1 343: * 344: * create nonzero element a(j-ka,j+1) outside the band 345: * and store it in WORK(j-m) 346: * 347: WORK( J-M ) = WORK( J-M )*AB( 1, J+1 ) 348: AB( 1, J+1 ) = RWORK( J-M )*AB( 1, J+1 ) 349: 90 CONTINUE 350: * 351: * generate rotations in 1st set to annihilate elements which 352: * have been created outside the band 353: * 354: IF( NRT.GT.0 ) 355: $ CALL ZLARGV( NRT, AB( 1, J2T ), INCA, WORK( J2T-M ), KA1, 356: $ RWORK( J2T-M ), KA1 ) 357: IF( NR.GT.0 ) THEN 358: * 359: * apply rotations in 1st set from the right 360: * 361: DO 100 L = 1, KA - 1 362: CALL ZLARTV( NR, AB( KA1-L, J2 ), INCA, 363: $ AB( KA-L, J2+1 ), INCA, RWORK( J2-M ), 364: $ WORK( J2-M ), KA1 ) 365: 100 CONTINUE 366: * 367: * apply rotations in 1st set from both sides to diagonal 368: * blocks 369: * 370: CALL ZLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ), 371: $ AB( KA, J2+1 ), INCA, RWORK( J2-M ), 372: $ WORK( J2-M ), KA1 ) 373: * 374: CALL ZLACGV( NR, WORK( J2-M ), KA1 ) 375: END IF 376: * 377: * start applying rotations in 1st set from the left 378: * 379: DO 110 L = KA - 1, KB - K + 1, -1 380: NRT = ( N-J2+L ) / KA1 381: IF( NRT.GT.0 ) 382: $ CALL ZLARTV( NRT, AB( L, J2+KA1-L ), INCA, 383: $ AB( L+1, J2+KA1-L ), INCA, RWORK( J2-M ), 384: $ WORK( J2-M ), KA1 ) 385: 110 CONTINUE 386: * 387: IF( WANTX ) THEN 388: * 389: * post-multiply X by product of rotations in 1st set 390: * 391: DO 120 J = J2, J1, KA1 392: CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1, 393: $ RWORK( J-M ), DCONJG( WORK( J-M ) ) ) 394: 120 CONTINUE 395: END IF 396: 130 CONTINUE 397: * 398: IF( UPDATE ) THEN 399: IF( I2.LE.N .AND. KBT.GT.0 ) THEN 400: * 401: * create nonzero element a(i-kbt,i-kbt+ka+1) outside the 402: * band and store it in WORK(i-kbt) 403: * 404: WORK( I-KBT ) = -BB( KB1-KBT, I )*RA1 405: END IF 406: END IF 407: * 408: DO 170 K = KB, 1, -1 409: IF( UPDATE ) THEN 410: J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1 411: ELSE 412: J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1 413: END IF 414: * 415: * finish applying rotations in 2nd set from the left 416: * 417: DO 140 L = KB - K, 1, -1 418: NRT = ( N-J2+KA+L ) / KA1 419: IF( NRT.GT.0 ) 420: $ CALL ZLARTV( NRT, AB( L, J2-L+1 ), INCA, 421: $ AB( L+1, J2-L+1 ), INCA, RWORK( J2-KA ), 422: $ WORK( J2-KA ), KA1 ) 423: 140 CONTINUE 424: NR = ( N-J2+KA ) / KA1 425: J1 = J2 + ( NR-1 )*KA1 426: DO 150 J = J1, J2, -KA1 427: WORK( J ) = WORK( J-KA ) 428: RWORK( J ) = RWORK( J-KA ) 429: 150 CONTINUE 430: DO 160 J = J2, J1, KA1 431: * 432: * create nonzero element a(j-ka,j+1) outside the band 433: * and store it in WORK(j) 434: * 435: WORK( J ) = WORK( J )*AB( 1, J+1 ) 436: AB( 1, J+1 ) = RWORK( J )*AB( 1, J+1 ) 437: 160 CONTINUE 438: IF( UPDATE ) THEN 439: IF( I-K.LT.N-KA .AND. K.LE.KBT ) 440: $ WORK( I-K+KA ) = WORK( I-K ) 441: END IF 442: 170 CONTINUE 443: * 444: DO 210 K = KB, 1, -1 445: J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1 446: NR = ( N-J2+KA ) / KA1 447: J1 = J2 + ( NR-1 )*KA1 448: IF( NR.GT.0 ) THEN 449: * 450: * generate rotations in 2nd set to annihilate elements 451: * which have been created outside the band 452: * 453: CALL ZLARGV( NR, AB( 1, J2 ), INCA, WORK( J2 ), KA1, 454: $ RWORK( J2 ), KA1 ) 455: * 456: * apply rotations in 2nd set from the right 457: * 458: DO 180 L = 1, KA - 1 459: CALL ZLARTV( NR, AB( KA1-L, J2 ), INCA, 460: $ AB( KA-L, J2+1 ), INCA, RWORK( J2 ), 461: $ WORK( J2 ), KA1 ) 462: 180 CONTINUE 463: * 464: * apply rotations in 2nd set from both sides to diagonal 465: * blocks 466: * 467: CALL ZLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ), 468: $ AB( KA, J2+1 ), INCA, RWORK( J2 ), 469: $ WORK( J2 ), KA1 ) 470: * 471: CALL ZLACGV( NR, WORK( J2 ), KA1 ) 472: END IF 473: * 474: * start applying rotations in 2nd set from the left 475: * 476: DO 190 L = KA - 1, KB - K + 1, -1 477: NRT = ( N-J2+L ) / KA1 478: IF( NRT.GT.0 ) 479: $ CALL ZLARTV( NRT, AB( L, J2+KA1-L ), INCA, 480: $ AB( L+1, J2+KA1-L ), INCA, RWORK( J2 ), 481: $ WORK( J2 ), KA1 ) 482: 190 CONTINUE 483: * 484: IF( WANTX ) THEN 485: * 486: * post-multiply X by product of rotations in 2nd set 487: * 488: DO 200 J = J2, J1, KA1 489: CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1, 490: $ RWORK( J ), DCONJG( WORK( J ) ) ) 491: 200 CONTINUE 492: END IF 493: 210 CONTINUE 494: * 495: DO 230 K = 1, KB - 1 496: J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1 497: * 498: * finish applying rotations in 1st set from the left 499: * 500: DO 220 L = KB - K, 1, -1 501: NRT = ( N-J2+L ) / KA1 502: IF( NRT.GT.0 ) 503: $ CALL ZLARTV( NRT, AB( L, J2+KA1-L ), INCA, 504: $ AB( L+1, J2+KA1-L ), INCA, RWORK( J2-M ), 505: $ WORK( J2-M ), KA1 ) 506: 220 CONTINUE 507: 230 CONTINUE 508: * 509: IF( KB.GT.1 ) THEN 510: DO 240 J = N - 1, J2 + KA, -1 511: RWORK( J-M ) = RWORK( J-KA-M ) 512: WORK( J-M ) = WORK( J-KA-M ) 513: 240 CONTINUE 514: END IF 515: * 516: ELSE 517: * 518: * Transform A, working with the lower triangle 519: * 520: IF( UPDATE ) THEN 521: * 522: * Form inv(S(i))**H * A * inv(S(i)) 523: * 524: BII = DBLE( BB( 1, I ) ) 525: AB( 1, I ) = ( DBLE( AB( 1, I ) ) / BII ) / BII 526: DO 250 J = I + 1, I1 527: AB( J-I+1, I ) = AB( J-I+1, I ) / BII 528: 250 CONTINUE 529: DO 260 J = MAX( 1, I-KA ), I - 1 530: AB( I-J+1, J ) = AB( I-J+1, J ) / BII 531: 260 CONTINUE 532: DO 290 K = I - KBT, I - 1 533: DO 270 J = I - KBT, K 534: AB( K-J+1, J ) = AB( K-J+1, J ) - 535: $ BB( I-J+1, J )*DCONJG( AB( I-K+1, 536: $ K ) ) - DCONJG( BB( I-K+1, K ) )* 537: $ AB( I-J+1, J ) + DBLE( AB( 1, I ) )* 538: $ BB( I-J+1, J )*DCONJG( BB( I-K+1, 539: $ K ) ) 540: 270 CONTINUE 541: DO 280 J = MAX( 1, I-KA ), I - KBT - 1 542: AB( K-J+1, J ) = AB( K-J+1, J ) - 543: $ DCONJG( BB( I-K+1, K ) )* 544: $ AB( I-J+1, J ) 545: 280 CONTINUE 546: 290 CONTINUE 547: DO 310 J = I, I1 548: DO 300 K = MAX( J-KA, I-KBT ), I - 1 549: AB( J-K+1, K ) = AB( J-K+1, K ) - 550: $ BB( I-K+1, K )*AB( J-I+1, I ) 551: 300 CONTINUE 552: 310 CONTINUE 553: * 554: IF( WANTX ) THEN 555: * 556: * post-multiply X by inv(S(i)) 557: * 558: CALL ZDSCAL( N-M, ONE / BII, X( M+1, I ), 1 ) 559: IF( KBT.GT.0 ) 560: $ CALL ZGERU( N-M, KBT, -CONE, X( M+1, I ), 1, 561: $ BB( KBT+1, I-KBT ), LDBB-1, 562: $ X( M+1, I-KBT ), LDX ) 563: END IF 564: * 565: * store a(i1,i) in RA1 for use in next loop over K 566: * 567: RA1 = AB( I1-I+1, I ) 568: END IF 569: * 570: * Generate and apply vectors of rotations to chase all the 571: * existing bulges KA positions down toward the bottom of the 572: * band 573: * 574: DO 360 K = 1, KB - 1 575: IF( UPDATE ) THEN 576: * 577: * Determine the rotations which would annihilate the bulge 578: * which has in theory just been created 579: * 580: IF( I-K+KA.LT.N .AND. I-K.GT.1 ) THEN 581: * 582: * generate rotation to annihilate a(i-k+ka+1,i) 583: * 584: CALL ZLARTG( AB( KA1-K, I ), RA1, RWORK( I-K+KA-M ), 585: $ WORK( I-K+KA-M ), RA ) 586: * 587: * create nonzero element a(i-k+ka+1,i-k) outside the 588: * band and store it in WORK(i-k) 589: * 590: T = -BB( K+1, I-K )*RA1 591: WORK( I-K ) = RWORK( I-K+KA-M )*T - 592: $ DCONJG( WORK( I-K+KA-M ) )* 593: $ AB( KA1, I-K ) 594: AB( KA1, I-K ) = WORK( I-K+KA-M )*T + 595: $ RWORK( I-K+KA-M )*AB( KA1, I-K ) 596: RA1 = RA 597: END IF 598: END IF 599: J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1 600: NR = ( N-J2+KA ) / KA1 601: J1 = J2 + ( NR-1 )*KA1 602: IF( UPDATE ) THEN 603: J2T = MAX( J2, I+2*KA-K+1 ) 604: ELSE 605: J2T = J2 606: END IF 607: NRT = ( N-J2T+KA ) / KA1 608: DO 320 J = J2T, J1, KA1 609: * 610: * create nonzero element a(j+1,j-ka) outside the band 611: * and store it in WORK(j-m) 612: * 613: WORK( J-M ) = WORK( J-M )*AB( KA1, J-KA+1 ) 614: AB( KA1, J-KA+1 ) = RWORK( J-M )*AB( KA1, J-KA+1 ) 615: 320 CONTINUE 616: * 617: * generate rotations in 1st set to annihilate elements which 618: * have been created outside the band 619: * 620: IF( NRT.GT.0 ) 621: $ CALL ZLARGV( NRT, AB( KA1, J2T-KA ), INCA, WORK( J2T-M ), 622: $ KA1, RWORK( J2T-M ), KA1 ) 623: IF( NR.GT.0 ) THEN 624: * 625: * apply rotations in 1st set from the left 626: * 627: DO 330 L = 1, KA - 1 628: CALL ZLARTV( NR, AB( L+1, J2-L ), INCA, 629: $ AB( L+2, J2-L ), INCA, RWORK( J2-M ), 630: $ WORK( J2-M ), KA1 ) 631: 330 CONTINUE 632: * 633: * apply rotations in 1st set from both sides to diagonal 634: * blocks 635: * 636: CALL ZLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ), 637: $ INCA, RWORK( J2-M ), WORK( J2-M ), KA1 ) 638: * 639: CALL ZLACGV( NR, WORK( J2-M ), KA1 ) 640: END IF 641: * 642: * start applying rotations in 1st set from the right 643: * 644: DO 340 L = KA - 1, KB - K + 1, -1 645: NRT = ( N-J2+L ) / KA1 646: IF( NRT.GT.0 ) 647: $ CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA, 648: $ AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ), 649: $ WORK( J2-M ), KA1 ) 650: 340 CONTINUE 651: * 652: IF( WANTX ) THEN 653: * 654: * post-multiply X by product of rotations in 1st set 655: * 656: DO 350 J = J2, J1, KA1 657: CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1, 658: $ RWORK( J-M ), WORK( J-M ) ) 659: 350 CONTINUE 660: END IF 661: 360 CONTINUE 662: * 663: IF( UPDATE ) THEN 664: IF( I2.LE.N .AND. KBT.GT.0 ) THEN 665: * 666: * create nonzero element a(i-kbt+ka+1,i-kbt) outside the 667: * band and store it in WORK(i-kbt) 668: * 669: WORK( I-KBT ) = -BB( KBT+1, I-KBT )*RA1 670: END IF 671: END IF 672: * 673: DO 400 K = KB, 1, -1 674: IF( UPDATE ) THEN 675: J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1 676: ELSE 677: J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1 678: END IF 679: * 680: * finish applying rotations in 2nd set from the right 681: * 682: DO 370 L = KB - K, 1, -1 683: NRT = ( N-J2+KA+L ) / KA1 684: IF( NRT.GT.0 ) 685: $ CALL ZLARTV( NRT, AB( KA1-L+1, J2-KA ), INCA, 686: $ AB( KA1-L, J2-KA+1 ), INCA, 687: $ RWORK( J2-KA ), WORK( J2-KA ), KA1 ) 688: 370 CONTINUE 689: NR = ( N-J2+KA ) / KA1 690: J1 = J2 + ( NR-1 )*KA1 691: DO 380 J = J1, J2, -KA1 692: WORK( J ) = WORK( J-KA ) 693: RWORK( J ) = RWORK( J-KA ) 694: 380 CONTINUE 695: DO 390 J = J2, J1, KA1 696: * 697: * create nonzero element a(j+1,j-ka) outside the band 698: * and store it in WORK(j) 699: * 700: WORK( J ) = WORK( J )*AB( KA1, J-KA+1 ) 701: AB( KA1, J-KA+1 ) = RWORK( J )*AB( KA1, J-KA+1 ) 702: 390 CONTINUE 703: IF( UPDATE ) THEN 704: IF( I-K.LT.N-KA .AND. K.LE.KBT ) 705: $ WORK( I-K+KA ) = WORK( I-K ) 706: END IF 707: 400 CONTINUE 708: * 709: DO 440 K = KB, 1, -1 710: J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1 711: NR = ( N-J2+KA ) / KA1 712: J1 = J2 + ( NR-1 )*KA1 713: IF( NR.GT.0 ) THEN 714: * 715: * generate rotations in 2nd set to annihilate elements 716: * which have been created outside the band 717: * 718: CALL ZLARGV( NR, AB( KA1, J2-KA ), INCA, WORK( J2 ), KA1, 719: $ RWORK( J2 ), KA1 ) 720: * 721: * apply rotations in 2nd set from the left 722: * 723: DO 410 L = 1, KA - 1 724: CALL ZLARTV( NR, AB( L+1, J2-L ), INCA, 725: $ AB( L+2, J2-L ), INCA, RWORK( J2 ), 726: $ WORK( J2 ), KA1 ) 727: 410 CONTINUE 728: * 729: * apply rotations in 2nd set from both sides to diagonal 730: * blocks 731: * 732: CALL ZLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ), 733: $ INCA, RWORK( J2 ), WORK( J2 ), KA1 ) 734: * 735: CALL ZLACGV( NR, WORK( J2 ), KA1 ) 736: END IF 737: * 738: * start applying rotations in 2nd set from the right 739: * 740: DO 420 L = KA - 1, KB - K + 1, -1 741: NRT = ( N-J2+L ) / KA1 742: IF( NRT.GT.0 ) 743: $ CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA, 744: $ AB( KA1-L, J2+1 ), INCA, RWORK( J2 ), 745: $ WORK( J2 ), KA1 ) 746: 420 CONTINUE 747: * 748: IF( WANTX ) THEN 749: * 750: * post-multiply X by product of rotations in 2nd set 751: * 752: DO 430 J = J2, J1, KA1 753: CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1, 754: $ RWORK( J ), WORK( J ) ) 755: 430 CONTINUE 756: END IF 757: 440 CONTINUE 758: * 759: DO 460 K = 1, KB - 1 760: J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1 761: * 762: * finish applying rotations in 1st set from the right 763: * 764: DO 450 L = KB - K, 1, -1 765: NRT = ( N-J2+L ) / KA1 766: IF( NRT.GT.0 ) 767: $ CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA, 768: $ AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ), 769: $ WORK( J2-M ), KA1 ) 770: 450 CONTINUE 771: 460 CONTINUE 772: * 773: IF( KB.GT.1 ) THEN 774: DO 470 J = N - 1, J2 + KA, -1 775: RWORK( J-M ) = RWORK( J-KA-M ) 776: WORK( J-M ) = WORK( J-KA-M ) 777: 470 CONTINUE 778: END IF 779: * 780: END IF 781: * 782: GO TO 10 783: * 784: 480 CONTINUE 785: * 786: * **************************** Phase 2 ***************************** 787: * 788: * The logical structure of this phase is: 789: * 790: * UPDATE = .TRUE. 791: * DO I = 1, M 792: * use S(i) to update A and create a new bulge 793: * apply rotations to push all bulges KA positions upward 794: * END DO 795: * UPDATE = .FALSE. 796: * DO I = M - KA - 1, 2, -1 797: * apply rotations to push all bulges KA positions upward 798: * END DO 799: * 800: * To avoid duplicating code, the two loops are merged. 801: * 802: UPDATE = .TRUE. 803: I = 0 804: 490 CONTINUE 805: IF( UPDATE ) THEN 806: I = I + 1 807: KBT = MIN( KB, M-I ) 808: I0 = I + 1 809: I1 = MAX( 1, I-KA ) 810: I2 = I + KBT - KA1 811: IF( I.GT.M ) THEN 812: UPDATE = .FALSE. 813: I = I - 1 814: I0 = M + 1 815: IF( KA.EQ.0 ) 816: $ RETURN 817: GO TO 490 818: END IF 819: ELSE 820: I = I - KA 821: IF( I.LT.2 ) 822: $ RETURN 823: END IF 824: * 825: IF( I.LT.M-KBT ) THEN 826: NX = M 827: ELSE 828: NX = N 829: END IF 830: * 831: IF( UPPER ) THEN 832: * 833: * Transform A, working with the upper triangle 834: * 835: IF( UPDATE ) THEN 836: * 837: * Form inv(S(i))**H * A * inv(S(i)) 838: * 839: BII = DBLE( BB( KB1, I ) ) 840: AB( KA1, I ) = ( DBLE( AB( KA1, I ) ) / BII ) / BII 841: DO 500 J = I1, I - 1 842: AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII 843: 500 CONTINUE 844: DO 510 J = I + 1, MIN( N, I+KA ) 845: AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII 846: 510 CONTINUE 847: DO 540 K = I + 1, I + KBT 848: DO 520 J = K, I + KBT 849: AB( K-J+KA1, J ) = AB( K-J+KA1, J ) - 850: $ BB( I-J+KB1, J )* 851: $ DCONJG( AB( I-K+KA1, K ) ) - 852: $ DCONJG( BB( I-K+KB1, K ) )* 853: $ AB( I-J+KA1, J ) + 854: $ DBLE( AB( KA1, I ) )* 855: $ BB( I-J+KB1, J )* 856: $ DCONJG( BB( I-K+KB1, K ) ) 857: 520 CONTINUE 858: DO 530 J = I + KBT + 1, MIN( N, I+KA ) 859: AB( K-J+KA1, J ) = AB( K-J+KA1, J ) - 860: $ DCONJG( BB( I-K+KB1, K ) )* 861: $ AB( I-J+KA1, J ) 862: 530 CONTINUE 863: 540 CONTINUE 864: DO 560 J = I1, I 865: DO 550 K = I + 1, MIN( J+KA, I+KBT ) 866: AB( J-K+KA1, K ) = AB( J-K+KA1, K ) - 867: $ BB( I-K+KB1, K )*AB( J-I+KA1, I ) 868: 550 CONTINUE 869: 560 CONTINUE 870: * 871: IF( WANTX ) THEN 872: * 873: * post-multiply X by inv(S(i)) 874: * 875: CALL ZDSCAL( NX, ONE / BII, X( 1, I ), 1 ) 876: IF( KBT.GT.0 ) 877: $ CALL ZGERU( NX, KBT, -CONE, X( 1, I ), 1, 878: $ BB( KB, I+1 ), LDBB-1, X( 1, I+1 ), LDX ) 879: END IF 880: * 881: * store a(i1,i) in RA1 for use in next loop over K 882: * 883: RA1 = AB( I1-I+KA1, I ) 884: END IF 885: * 886: * Generate and apply vectors of rotations to chase all the 887: * existing bulges KA positions up toward the top of the band 888: * 889: DO 610 K = 1, KB - 1 890: IF( UPDATE ) THEN 891: * 892: * Determine the rotations which would annihilate the bulge 893: * which has in theory just been created 894: * 895: IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN 896: * 897: * generate rotation to annihilate a(i+k-ka-1,i) 898: * 899: CALL ZLARTG( AB( K+1, I ), RA1, RWORK( I+K-KA ), 900: $ WORK( I+K-KA ), RA ) 901: * 902: * create nonzero element a(i+k-ka-1,i+k) outside the 903: * band and store it in WORK(m-kb+i+k) 904: * 905: T = -BB( KB1-K, I+K )*RA1 906: WORK( M-KB+I+K ) = RWORK( I+K-KA )*T - 907: $ DCONJG( WORK( I+K-KA ) )* 908: $ AB( 1, I+K ) 909: AB( 1, I+K ) = WORK( I+K-KA )*T + 910: $ RWORK( I+K-KA )*AB( 1, I+K ) 911: RA1 = RA 912: END IF 913: END IF 914: J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1 915: NR = ( J2+KA-1 ) / KA1 916: J1 = J2 - ( NR-1 )*KA1 917: IF( UPDATE ) THEN 918: J2T = MIN( J2, I-2*KA+K-1 ) 919: ELSE 920: J2T = J2 921: END IF 922: NRT = ( J2T+KA-1 ) / KA1 923: DO 570 J = J1, J2T, KA1 924: * 925: * create nonzero element a(j-1,j+ka) outside the band 926: * and store it in WORK(j) 927: * 928: WORK( J ) = WORK( J )*AB( 1, J+KA-1 ) 929: AB( 1, J+KA-1 ) = RWORK( J )*AB( 1, J+KA-1 ) 930: 570 CONTINUE 931: * 932: * generate rotations in 1st set to annihilate elements which 933: * have been created outside the band 934: * 935: IF( NRT.GT.0 ) 936: $ CALL ZLARGV( NRT, AB( 1, J1+KA ), INCA, WORK( J1 ), KA1, 937: $ RWORK( J1 ), KA1 ) 938: IF( NR.GT.0 ) THEN 939: * 940: * apply rotations in 1st set from the left 941: * 942: DO 580 L = 1, KA - 1 943: CALL ZLARTV( NR, AB( KA1-L, J1+L ), INCA, 944: $ AB( KA-L, J1+L ), INCA, RWORK( J1 ), 945: $ WORK( J1 ), KA1 ) 946: 580 CONTINUE 947: * 948: * apply rotations in 1st set from both sides to diagonal 949: * blocks 950: * 951: CALL ZLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ), 952: $ AB( KA, J1 ), INCA, RWORK( J1 ), WORK( J1 ), 953: $ KA1 ) 954: * 955: CALL ZLACGV( NR, WORK( J1 ), KA1 ) 956: END IF 957: * 958: * start applying rotations in 1st set from the right 959: * 960: DO 590 L = KA - 1, KB - K + 1, -1 961: NRT = ( J2+L-1 ) / KA1 962: J1T = J2 - ( NRT-1 )*KA1 963: IF( NRT.GT.0 ) 964: $ CALL ZLARTV( NRT, AB( L, J1T ), INCA, 965: $ AB( L+1, J1T-1 ), INCA, RWORK( J1T ), 966: $ WORK( J1T ), KA1 ) 967: 590 CONTINUE 968: * 969: IF( WANTX ) THEN 970: * 971: * post-multiply X by product of rotations in 1st set 972: * 973: DO 600 J = J1, J2, KA1 974: CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1, 975: $ RWORK( J ), WORK( J ) ) 976: 600 CONTINUE 977: END IF 978: 610 CONTINUE 979: * 980: IF( UPDATE ) THEN 981: IF( I2.GT.0 .AND. KBT.GT.0 ) THEN 982: * 983: * create nonzero element a(i+kbt-ka-1,i+kbt) outside the 984: * band and store it in WORK(m-kb+i+kbt) 985: * 986: WORK( M-KB+I+KBT ) = -BB( KB1-KBT, I+KBT )*RA1 987: END IF 988: END IF 989: * 990: DO 650 K = KB, 1, -1 991: IF( UPDATE ) THEN 992: J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1 993: ELSE 994: J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1 995: END IF 996: * 997: * finish applying rotations in 2nd set from the right 998: * 999: DO 620 L = KB - K, 1, -1 1000: NRT = ( J2+KA+L-1 ) / KA1 1001: J1T = J2 - ( NRT-1 )*KA1 1002: IF( NRT.GT.0 ) 1003: $ CALL ZLARTV( NRT, AB( L, J1T+KA ), INCA, 1004: $ AB( L+1, J1T+KA-1 ), INCA, 1005: $ RWORK( M-KB+J1T+KA ), 1006: $ WORK( M-KB+J1T+KA ), KA1 ) 1007: 620 CONTINUE 1008: NR = ( J2+KA-1 ) / KA1 1009: J1 = J2 - ( NR-1 )*KA1 1010: DO 630 J = J1, J2, KA1 1011: WORK( M-KB+J ) = WORK( M-KB+J+KA ) 1012: RWORK( M-KB+J ) = RWORK( M-KB+J+KA ) 1013: 630 CONTINUE 1014: DO 640 J = J1, J2, KA1 1015: * 1016: * create nonzero element a(j-1,j+ka) outside the band 1017: * and store it in WORK(m-kb+j) 1018: * 1019: WORK( M-KB+J ) = WORK( M-KB+J )*AB( 1, J+KA-1 ) 1020: AB( 1, J+KA-1 ) = RWORK( M-KB+J )*AB( 1, J+KA-1 ) 1021: 640 CONTINUE 1022: IF( UPDATE ) THEN 1023: IF( I+K.GT.KA1 .AND. K.LE.KBT ) 1024: $ WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K ) 1025: END IF 1026: 650 CONTINUE 1027: * 1028: DO 690 K = KB, 1, -1 1029: J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1 1030: NR = ( J2+KA-1 ) / KA1 1031: J1 = J2 - ( NR-1 )*KA1 1032: IF( NR.GT.0 ) THEN 1033: * 1034: * generate rotations in 2nd set to annihilate elements 1035: * which have been created outside the band 1036: * 1037: CALL ZLARGV( NR, AB( 1, J1+KA ), INCA, WORK( M-KB+J1 ), 1038: $ KA1, RWORK( M-KB+J1 ), KA1 ) 1039: * 1040: * apply rotations in 2nd set from the left 1041: * 1042: DO 660 L = 1, KA - 1 1043: CALL ZLARTV( NR, AB( KA1-L, J1+L ), INCA, 1044: $ AB( KA-L, J1+L ), INCA, RWORK( M-KB+J1 ), 1045: $ WORK( M-KB+J1 ), KA1 ) 1046: 660 CONTINUE 1047: * 1048: * apply rotations in 2nd set from both sides to diagonal 1049: * blocks 1050: * 1051: CALL ZLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ), 1052: $ AB( KA, J1 ), INCA, RWORK( M-KB+J1 ), 1053: $ WORK( M-KB+J1 ), KA1 ) 1054: * 1055: CALL ZLACGV( NR, WORK( M-KB+J1 ), KA1 ) 1056: END IF 1057: * 1058: * start applying rotations in 2nd set from the right 1059: * 1060: DO 670 L = KA - 1, KB - K + 1, -1 1061: NRT = ( J2+L-1 ) / KA1 1062: J1T = J2 - ( NRT-1 )*KA1 1063: IF( NRT.GT.0 ) 1064: $ CALL ZLARTV( NRT, AB( L, J1T ), INCA, 1065: $ AB( L+1, J1T-1 ), INCA, 1066: $ RWORK( M-KB+J1T ), WORK( M-KB+J1T ), 1067: $ KA1 ) 1068: 670 CONTINUE 1069: * 1070: IF( WANTX ) THEN 1071: * 1072: * post-multiply X by product of rotations in 2nd set 1073: * 1074: DO 680 J = J1, J2, KA1 1075: CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1, 1076: $ RWORK( M-KB+J ), WORK( M-KB+J ) ) 1077: 680 CONTINUE 1078: END IF 1079: 690 CONTINUE 1080: * 1081: DO 710 K = 1, KB - 1 1082: J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1 1083: * 1084: * finish applying rotations in 1st set from the right 1085: * 1086: DO 700 L = KB - K, 1, -1 1087: NRT = ( J2+L-1 ) / KA1 1088: J1T = J2 - ( NRT-1 )*KA1 1089: IF( NRT.GT.0 ) 1090: $ CALL ZLARTV( NRT, AB( L, J1T ), INCA, 1091: $ AB( L+1, J1T-1 ), INCA, RWORK( J1T ), 1092: $ WORK( J1T ), KA1 ) 1093: 700 CONTINUE 1094: 710 CONTINUE 1095: * 1096: IF( KB.GT.1 ) THEN 1097: DO 720 J = 2, I2 - KA 1098: RWORK( J ) = RWORK( J+KA ) 1099: WORK( J ) = WORK( J+KA ) 1100: 720 CONTINUE 1101: END IF 1102: * 1103: ELSE 1104: * 1105: * Transform A, working with the lower triangle 1106: * 1107: IF( UPDATE ) THEN 1108: * 1109: * Form inv(S(i))**H * A * inv(S(i)) 1110: * 1111: BII = DBLE( BB( 1, I ) ) 1112: AB( 1, I ) = ( DBLE( AB( 1, I ) ) / BII ) / BII 1113: DO 730 J = I1, I - 1 1114: AB( I-J+1, J ) = AB( I-J+1, J ) / BII 1115: 730 CONTINUE 1116: DO 740 J = I + 1, MIN( N, I+KA ) 1117: AB( J-I+1, I ) = AB( J-I+1, I ) / BII 1118: 740 CONTINUE 1119: DO 770 K = I + 1, I + KBT 1120: DO 750 J = K, I + KBT 1121: AB( J-K+1, K ) = AB( J-K+1, K ) - 1122: $ BB( J-I+1, I )*DCONJG( AB( K-I+1, 1123: $ I ) ) - DCONJG( BB( K-I+1, I ) )* 1124: $ AB( J-I+1, I ) + DBLE( AB( 1, I ) )* 1125: $ BB( J-I+1, I )*DCONJG( BB( K-I+1, 1126: $ I ) ) 1127: 750 CONTINUE 1128: DO 760 J = I + KBT + 1, MIN( N, I+KA ) 1129: AB( J-K+1, K ) = AB( J-K+1, K ) - 1130: $ DCONJG( BB( K-I+1, I ) )* 1131: $ AB( J-I+1, I ) 1132: 760 CONTINUE 1133: 770 CONTINUE 1134: DO 790 J = I1, I 1135: DO 780 K = I + 1, MIN( J+KA, I+KBT ) 1136: AB( K-J+1, J ) = AB( K-J+1, J ) - 1137: $ BB( K-I+1, I )*AB( I-J+1, J ) 1138: 780 CONTINUE 1139: 790 CONTINUE 1140: * 1141: IF( WANTX ) THEN 1142: * 1143: * post-multiply X by inv(S(i)) 1144: * 1145: CALL ZDSCAL( NX, ONE / BII, X( 1, I ), 1 ) 1146: IF( KBT.GT.0 ) 1147: $ CALL ZGERC( NX, KBT, -CONE, X( 1, I ), 1, BB( 2, I ), 1148: $ 1, X( 1, I+1 ), LDX ) 1149: END IF 1150: * 1151: * store a(i,i1) in RA1 for use in next loop over K 1152: * 1153: RA1 = AB( I-I1+1, I1 ) 1154: END IF 1155: * 1156: * Generate and apply vectors of rotations to chase all the 1157: * existing bulges KA positions up toward the top of the band 1158: * 1159: DO 840 K = 1, KB - 1 1160: IF( UPDATE ) THEN 1161: * 1162: * Determine the rotations which would annihilate the bulge 1163: * which has in theory just been created 1164: * 1165: IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN 1166: * 1167: * generate rotation to annihilate a(i,i+k-ka-1) 1168: * 1169: CALL ZLARTG( AB( KA1-K, I+K-KA ), RA1, 1170: $ RWORK( I+K-KA ), WORK( I+K-KA ), RA ) 1171: * 1172: * create nonzero element a(i+k,i+k-ka-1) outside the 1173: * band and store it in WORK(m-kb+i+k) 1174: * 1175: T = -BB( K+1, I )*RA1 1176: WORK( M-KB+I+K ) = RWORK( I+K-KA )*T - 1177: $ DCONJG( WORK( I+K-KA ) )* 1178: $ AB( KA1, I+K-KA ) 1179: AB( KA1, I+K-KA ) = WORK( I+K-KA )*T + 1180: $ RWORK( I+K-KA )*AB( KA1, I+K-KA ) 1181: RA1 = RA 1182: END IF 1183: END IF 1184: J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1 1185: NR = ( J2+KA-1 ) / KA1 1186: J1 = J2 - ( NR-1 )*KA1 1187: IF( UPDATE ) THEN 1188: J2T = MIN( J2, I-2*KA+K-1 ) 1189: ELSE 1190: J2T = J2 1191: END IF 1192: NRT = ( J2T+KA-1 ) / KA1 1193: DO 800 J = J1, J2T, KA1 1194: * 1195: * create nonzero element a(j+ka,j-1) outside the band 1196: * and store it in WORK(j) 1197: * 1198: WORK( J ) = WORK( J )*AB( KA1, J-1 ) 1199: AB( KA1, J-1 ) = RWORK( J )*AB( KA1, J-1 ) 1200: 800 CONTINUE 1201: * 1202: * generate rotations in 1st set to annihilate elements which 1203: * have been created outside the band 1204: * 1205: IF( NRT.GT.0 ) 1206: $ CALL ZLARGV( NRT, AB( KA1, J1 ), INCA, WORK( J1 ), KA1, 1207: $ RWORK( J1 ), KA1 ) 1208: IF( NR.GT.0 ) THEN 1209: * 1210: * apply rotations in 1st set from the right 1211: * 1212: DO 810 L = 1, KA - 1 1213: CALL ZLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ), 1214: $ INCA, RWORK( J1 ), WORK( J1 ), KA1 ) 1215: 810 CONTINUE 1216: * 1217: * apply rotations in 1st set from both sides to diagonal 1218: * blocks 1219: * 1220: CALL ZLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ), 1221: $ AB( 2, J1-1 ), INCA, RWORK( J1 ), 1222: $ WORK( J1 ), KA1 ) 1223: * 1224: CALL ZLACGV( NR, WORK( J1 ), KA1 ) 1225: END IF 1226: * 1227: * start applying rotations in 1st set from the left 1228: * 1229: DO 820 L = KA - 1, KB - K + 1, -1 1230: NRT = ( J2+L-1 ) / KA1 1231: J1T = J2 - ( NRT-1 )*KA1 1232: IF( NRT.GT.0 ) 1233: $ CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA, 1234: $ AB( KA1-L, J1T-KA1+L ), INCA, 1235: $ RWORK( J1T ), WORK( J1T ), KA1 ) 1236: 820 CONTINUE 1237: * 1238: IF( WANTX ) THEN 1239: * 1240: * post-multiply X by product of rotations in 1st set 1241: * 1242: DO 830 J = J1, J2, KA1 1243: CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1, 1244: $ RWORK( J ), DCONJG( WORK( J ) ) ) 1245: 830 CONTINUE 1246: END IF 1247: 840 CONTINUE 1248: * 1249: IF( UPDATE ) THEN 1250: IF( I2.GT.0 .AND. KBT.GT.0 ) THEN 1251: * 1252: * create nonzero element a(i+kbt,i+kbt-ka-1) outside the 1253: * band and store it in WORK(m-kb+i+kbt) 1254: * 1255: WORK( M-KB+I+KBT ) = -BB( KBT+1, I )*RA1 1256: END IF 1257: END IF 1258: * 1259: DO 880 K = KB, 1, -1 1260: IF( UPDATE ) THEN 1261: J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1 1262: ELSE 1263: J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1 1264: END IF 1265: * 1266: * finish applying rotations in 2nd set from the left 1267: * 1268: DO 850 L = KB - K, 1, -1 1269: NRT = ( J2+KA+L-1 ) / KA1 1270: J1T = J2 - ( NRT-1 )*KA1 1271: IF( NRT.GT.0 ) 1272: $ CALL ZLARTV( NRT, AB( KA1-L+1, J1T+L-1 ), INCA, 1273: $ AB( KA1-L, J1T+L-1 ), INCA, 1274: $ RWORK( M-KB+J1T+KA ), 1275: $ WORK( M-KB+J1T+KA ), KA1 ) 1276: 850 CONTINUE 1277: NR = ( J2+KA-1 ) / KA1 1278: J1 = J2 - ( NR-1 )*KA1 1279: DO 860 J = J1, J2, KA1 1280: WORK( M-KB+J ) = WORK( M-KB+J+KA ) 1281: RWORK( M-KB+J ) = RWORK( M-KB+J+KA ) 1282: 860 CONTINUE 1283: DO 870 J = J1, J2, KA1 1284: * 1285: * create nonzero element a(j+ka,j-1) outside the band 1286: * and store it in WORK(m-kb+j) 1287: * 1288: WORK( M-KB+J ) = WORK( M-KB+J )*AB( KA1, J-1 ) 1289: AB( KA1, J-1 ) = RWORK( M-KB+J )*AB( KA1, J-1 ) 1290: 870 CONTINUE 1291: IF( UPDATE ) THEN 1292: IF( I+K.GT.KA1 .AND. K.LE.KBT ) 1293: $ WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K ) 1294: END IF 1295: 880 CONTINUE 1296: * 1297: DO 920 K = KB, 1, -1 1298: J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1 1299: NR = ( J2+KA-1 ) / KA1 1300: J1 = J2 - ( NR-1 )*KA1 1301: IF( NR.GT.0 ) THEN 1302: * 1303: * generate rotations in 2nd set to annihilate elements 1304: * which have been created outside the band 1305: * 1306: CALL ZLARGV( NR, AB( KA1, J1 ), INCA, WORK( M-KB+J1 ), 1307: $ KA1, RWORK( M-KB+J1 ), KA1 ) 1308: * 1309: * apply rotations in 2nd set from the right 1310: * 1311: DO 890 L = 1, KA - 1 1312: CALL ZLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ), 1313: $ INCA, RWORK( M-KB+J1 ), WORK( M-KB+J1 ), 1314: $ KA1 ) 1315: 890 CONTINUE 1316: * 1317: * apply rotations in 2nd set from both sides to diagonal 1318: * blocks 1319: * 1320: CALL ZLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ), 1321: $ AB( 2, J1-1 ), INCA, RWORK( M-KB+J1 ), 1322: $ WORK( M-KB+J1 ), KA1 ) 1323: * 1324: CALL ZLACGV( NR, WORK( M-KB+J1 ), KA1 ) 1325: END IF 1326: * 1327: * start applying rotations in 2nd set from the left 1328: * 1329: DO 900 L = KA - 1, KB - K + 1, -1 1330: NRT = ( J2+L-1 ) / KA1 1331: J1T = J2 - ( NRT-1 )*KA1 1332: IF( NRT.GT.0 ) 1333: $ CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA, 1334: $ AB( KA1-L, J1T-KA1+L ), INCA, 1335: $ RWORK( M-KB+J1T ), WORK( M-KB+J1T ), 1336: $ KA1 ) 1337: 900 CONTINUE 1338: * 1339: IF( WANTX ) THEN 1340: * 1341: * post-multiply X by product of rotations in 2nd set 1342: * 1343: DO 910 J = J1, J2, KA1 1344: CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1, 1345: $ RWORK( M-KB+J ), DCONJG( WORK( M-KB+J ) ) ) 1346: 910 CONTINUE 1347: END IF 1348: 920 CONTINUE 1349: * 1350: DO 940 K = 1, KB - 1 1351: J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1 1352: * 1353: * finish applying rotations in 1st set from the left 1354: * 1355: DO 930 L = KB - K, 1, -1 1356: NRT = ( J2+L-1 ) / KA1 1357: J1T = J2 - ( NRT-1 )*KA1 1358: IF( NRT.GT.0 ) 1359: $ CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA, 1360: $ AB( KA1-L, J1T-KA1+L ), INCA, 1361: $ RWORK( J1T ), WORK( J1T ), KA1 ) 1362: 930 CONTINUE 1363: 940 CONTINUE 1364: * 1365: IF( KB.GT.1 ) THEN 1366: DO 950 J = 2, I2 - KA 1367: RWORK( J ) = RWORK( J+KA ) 1368: WORK( J ) = WORK( J+KA ) 1369: 950 CONTINUE 1370: END IF 1371: * 1372: END IF 1373: * 1374: GO TO 490 1375: * 1376: * End of ZHBGST 1377: * 1378: END