Annotation of rpl/lapack/lapack/zhbgst.f, revision 1.1
1.1 ! bertrand 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) --
! 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 7: * November 2006
! 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, I2 + 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, I2 + 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhbgst.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack