0001:       SUBROUTINE ZHBGST( VECT, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, X,
0002:      $                   LDX, WORK, RWORK, INFO )
0003: *
0004: *  -- LAPACK routine (version 3.2) --
0005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
0006: *     November 2006
0007: *
0008: *     .. Scalar Arguments ..
0009:       CHARACTER          UPLO, VECT
0010:       INTEGER            INFO, KA, KB, LDAB, LDBB, LDX, N
0011: *     ..
0012: *     .. Array Arguments ..
0013:       DOUBLE PRECISION   RWORK( * )
0014:       COMPLEX*16         AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
0015:      $                   X( LDX, * )
0016: *     ..
0017: *
0018: *  Purpose
0019: *  =======
0020: *
0021: *  ZHBGST reduces a complex Hermitian-definite banded generalized
0022: *  eigenproblem  A*x = lambda*B*x  to standard form  C*y = lambda*y,
0023: *  such that C has the same bandwidth as A.
0024: *
0025: *  B must have been previously factorized as S**H*S by ZPBSTF, using a
0026: *  split Cholesky factorization. A is overwritten by C = X**H*A*X, where
0027: *  X = S**(-1)*Q and Q is a unitary matrix chosen to preserve the
0028: *  bandwidth of A.
0029: *
0030: *  Arguments
0031: *  =========
0032: *
0033: *  VECT    (input) CHARACTER*1
0034: *          = 'N':  do not form the transformation matrix X;
0035: *          = 'V':  form X.
0036: *
0037: *  UPLO    (input) CHARACTER*1
0038: *          = 'U':  Upper triangle of A is stored;
0039: *          = 'L':  Lower triangle of A is stored.
0040: *
0041: *  N       (input) INTEGER
0042: *          The order of the matrices A and B.  N >= 0.
0043: *
0044: *  KA      (input) INTEGER
0045: *          The number of superdiagonals of the matrix A if UPLO = 'U',
0046: *          or the number of subdiagonals if UPLO = 'L'.  KA >= 0.
0047: *
0048: *  KB      (input) INTEGER
0049: *          The number of superdiagonals of the matrix B if UPLO = 'U',
0050: *          or the number of subdiagonals if UPLO = 'L'.  KA >= KB >= 0.
0051: *
0052: *  AB      (input/output) COMPLEX*16 array, dimension (LDAB,N)
0053: *          On entry, the upper or lower triangle of the Hermitian band
0054: *          matrix A, stored in the first ka+1 rows of the array.  The
0055: *          j-th column of A is stored in the j-th column of the array AB
0056: *          as follows:
0057: *          if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;
0058: *          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+ka).
0059: *
0060: *          On exit, the transformed matrix X**H*A*X, stored in the same
0061: *          format as A.
0062: *
0063: *  LDAB    (input) INTEGER
0064: *          The leading dimension of the array AB.  LDAB >= KA+1.
0065: *
0066: *  BB      (input) COMPLEX*16 array, dimension (LDBB,N)
0067: *          The banded factor S from the split Cholesky factorization of
0068: *          B, as returned by ZPBSTF, stored in the first kb+1 rows of
0069: *          the array.
0070: *
0071: *  LDBB    (input) INTEGER
0072: *          The leading dimension of the array BB.  LDBB >= KB+1.
0073: *
0074: *  X       (output) COMPLEX*16 array, dimension (LDX,N)
0075: *          If VECT = 'V', the n-by-n matrix X.
0076: *          If VECT = 'N', the array X is not referenced.
0077: *
0078: *  LDX     (input) INTEGER
0079: *          The leading dimension of the array X.
0080: *          LDX >= max(1,N) if VECT = 'V'; LDX >= 1 otherwise.
0081: *
0082: *  WORK    (workspace) COMPLEX*16 array, dimension (N)
0083: *
0084: *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
0085: *
0086: *  INFO    (output) INTEGER
0087: *          = 0:  successful exit
0088: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
0089: *
0090: *  =====================================================================
0091: *
0092: *     .. Parameters ..
0093:       COMPLEX*16         CZERO, CONE
0094:       DOUBLE PRECISION   ONE
0095:       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
0096:      $                   CONE = ( 1.0D+0, 0.0D+0 ), ONE = 1.0D+0 )
0097: *     ..
0098: *     .. Local Scalars ..
0099:       LOGICAL            UPDATE, UPPER, WANTX
0100:       INTEGER            I, I0, I1, I2, INCA, J, J1, J1T, J2, J2T, K,
0101:      $                   KA1, KB1, KBT, L, M, NR, NRT, NX
0102:       DOUBLE PRECISION   BII
0103:       COMPLEX*16         RA, RA1, T
0104: *     ..
0105: *     .. External Functions ..
0106:       LOGICAL            LSAME
0107:       EXTERNAL           LSAME
0108: *     ..
0109: *     .. External Subroutines ..
0110:       EXTERNAL           XERBLA, ZDSCAL, ZGERC, ZGERU, ZLACGV, ZLAR2V,
0111:      $                   ZLARGV, ZLARTG, ZLARTV, ZLASET, ZROT
0112: *     ..
0113: *     .. Intrinsic Functions ..
0114:       INTRINSIC          DBLE, DCONJG, MAX, MIN
0115: *     ..
0116: *     .. Executable Statements ..
0117: *
0118: *     Test the input parameters
0119: *
0120:       WANTX = LSAME( VECT, 'V' )
0121:       UPPER = LSAME( UPLO, 'U' )
0122:       KA1 = KA + 1
0123:       KB1 = KB + 1
0124:       INFO = 0
0125:       IF( .NOT.WANTX .AND. .NOT.LSAME( VECT, 'N' ) ) THEN
0126:          INFO = -1
0127:       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
0128:          INFO = -2
0129:       ELSE IF( N.LT.0 ) THEN
0130:          INFO = -3
0131:       ELSE IF( KA.LT.0 ) THEN
0132:          INFO = -4
0133:       ELSE IF( KB.LT.0 .OR. KB.GT.KA ) THEN
0134:          INFO = -5
0135:       ELSE IF( LDAB.LT.KA+1 ) THEN
0136:          INFO = -7
0137:       ELSE IF( LDBB.LT.KB+1 ) THEN
0138:          INFO = -9
0139:       ELSE IF( LDX.LT.1 .OR. WANTX .AND. LDX.LT.MAX( 1, N ) ) THEN
0140:          INFO = -11
0141:       END IF
0142:       IF( INFO.NE.0 ) THEN
0143:          CALL XERBLA( 'ZHBGST', -INFO )
0144:          RETURN
0145:       END IF
0146: *
0147: *     Quick return if possible
0148: *
0149:       IF( N.EQ.0 )
0150:      $   RETURN
0151: *
0152:       INCA = LDAB*KA1
0153: *
0154: *     Initialize X to the unit matrix, if needed
0155: *
0156:       IF( WANTX )
0157:      $   CALL ZLASET( 'Full', N, N, CZERO, CONE, X, LDX )
0158: *
0159: *     Set M to the splitting point m. It must be the same value as is
0160: *     used in ZPBSTF. The chosen value allows the arrays WORK and RWORK
0161: *     to be of dimension (N).
0162: *
0163:       M = ( N+KB ) / 2
0164: *
0165: *     The routine works in two phases, corresponding to the two halves
0166: *     of the split Cholesky factorization of B as S**H*S where
0167: *
0168: *     S = ( U    )
0169: *         ( M  L )
0170: *
0171: *     with U upper triangular of order m, and L lower triangular of
0172: *     order n-m. S has the same bandwidth as B.
0173: *
0174: *     S is treated as a product of elementary matrices:
0175: *
0176: *     S = S(m)*S(m-1)*...*S(2)*S(1)*S(m+1)*S(m+2)*...*S(n-1)*S(n)
0177: *
0178: *     where S(i) is determined by the i-th row of S.
0179: *
0180: *     In phase 1, the index i takes the values n, n-1, ... , m+1;
0181: *     in phase 2, it takes the values 1, 2, ... , m.
0182: *
0183: *     For each value of i, the current matrix A is updated by forming
0184: *     inv(S(i))**H*A*inv(S(i)). This creates a triangular bulge outside
0185: *     the band of A. The bulge is then pushed down toward the bottom of
0186: *     A in phase 1, and up toward the top of A in phase 2, by applying
0187: *     plane rotations.
0188: *
0189: *     There are kb*(kb+1)/2 elements in the bulge, but at most 2*kb-1
0190: *     of them are linearly independent, so annihilating a bulge requires
0191: *     only 2*kb-1 plane rotations. The rotations are divided into a 1st
0192: *     set of kb-1 rotations, and a 2nd set of kb rotations.
0193: *
0194: *     Wherever possible, rotations are generated and applied in vector
0195: *     operations of length NR between the indices J1 and J2 (sometimes
0196: *     replaced by modified values NRT, J1T or J2T).
0197: *
0198: *     The real cosines and complex sines of the rotations are stored in
0199: *     the arrays RWORK and WORK, those of the 1st set in elements
0200: *     2:m-kb-1, and those of the 2nd set in elements m-kb+1:n.
0201: *
0202: *     The bulges are not formed explicitly; nonzero elements outside the
0203: *     band are created only when they are required for generating new
0204: *     rotations; they are stored in the array WORK, in positions where
0205: *     they are later overwritten by the sines of the rotations which
0206: *     annihilate them.
0207: *
0208: *     **************************** Phase 1 *****************************
0209: *
0210: *     The logical structure of this phase is:
0211: *
0212: *     UPDATE = .TRUE.
0213: *     DO I = N, M + 1, -1
0214: *        use S(i) to update A and create a new bulge
0215: *        apply rotations to push all bulges KA positions downward
0216: *     END DO
0217: *     UPDATE = .FALSE.
0218: *     DO I = M + KA + 1, N - 1
0219: *        apply rotations to push all bulges KA positions downward
0220: *     END DO
0221: *
0222: *     To avoid duplicating code, the two loops are merged.
0223: *
0224:       UPDATE = .TRUE.
0225:       I = N + 1
0226:    10 CONTINUE
0227:       IF( UPDATE ) THEN
0228:          I = I - 1
0229:          KBT = MIN( KB, I-1 )
0230:          I0 = I - 1
0231:          I1 = MIN( N, I+KA )
0232:          I2 = I - KBT + KA1
0233:          IF( I.LT.M+1 ) THEN
0234:             UPDATE = .FALSE.
0235:             I = I + 1
0236:             I0 = M
0237:             IF( KA.EQ.0 )
0238:      $         GO TO 480
0239:             GO TO 10
0240:          END IF
0241:       ELSE
0242:          I = I + KA
0243:          IF( I.GT.N-1 )
0244:      $      GO TO 480
0245:       END IF
0246: *
0247:       IF( UPPER ) THEN
0248: *
0249: *        Transform A, working with the upper triangle
0250: *
0251:          IF( UPDATE ) THEN
0252: *
0253: *           Form  inv(S(i))**H * A * inv(S(i))
0254: *
0255:             BII = DBLE( BB( KB1, I ) )
0256:             AB( KA1, I ) = ( DBLE( AB( KA1, I ) ) / BII ) / BII
0257:             DO 20 J = I + 1, I1
0258:                AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII
0259:    20       CONTINUE
0260:             DO 30 J = MAX( 1, I-KA ), I - 1
0261:                AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII
0262:    30       CONTINUE
0263:             DO 60 K = I - KBT, I - 1
0264:                DO 40 J = I - KBT, K
0265:                   AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
0266:      $                               BB( J-I+KB1, I )*
0267:      $                               DCONJG( AB( K-I+KA1, I ) ) -
0268:      $                               DCONJG( BB( K-I+KB1, I ) )*
0269:      $                               AB( J-I+KA1, I ) +
0270:      $                               DBLE( AB( KA1, I ) )*
0271:      $                               BB( J-I+KB1, I )*
0272:      $                               DCONJG( BB( K-I+KB1, I ) )
0273:    40          CONTINUE
0274:                DO 50 J = MAX( 1, I-KA ), I - KBT - 1
0275:                   AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
0276:      $                               DCONJG( BB( K-I+KB1, I ) )*
0277:      $                               AB( J-I+KA1, I )
0278:    50          CONTINUE
0279:    60       CONTINUE
0280:             DO 80 J = I, I1
0281:                DO 70 K = MAX( J-KA, I-KBT ), I - 1
0282:                   AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
0283:      $                               BB( K-I+KB1, I )*AB( I-J+KA1, J )
0284:    70          CONTINUE
0285:    80       CONTINUE
0286: *
0287:             IF( WANTX ) THEN
0288: *
0289: *              post-multiply X by inv(S(i))
0290: *
0291:                CALL ZDSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
0292:                IF( KBT.GT.0 )
0293:      $            CALL ZGERC( N-M, KBT, -CONE, X( M+1, I ), 1,
0294:      $                        BB( KB1-KBT, I ), 1, X( M+1, I-KBT ),
0295:      $                        LDX )
0296:             END IF
0297: *
0298: *           store a(i,i1) in RA1 for use in next loop over K
0299: *
0300:             RA1 = AB( I-I1+KA1, I1 )
0301:          END IF
0302: *
0303: *        Generate and apply vectors of rotations to chase all the
0304: *        existing bulges KA positions down toward the bottom of the
0305: *        band
0306: *
0307:          DO 130 K = 1, KB - 1
0308:             IF( UPDATE ) THEN
0309: *
0310: *              Determine the rotations which would annihilate the bulge
0311: *              which has in theory just been created
0312: *
0313:                IF( I-K+KA.LT.N .AND. I-K.GT.1 ) THEN
0314: *
0315: *                 generate rotation to annihilate a(i,i-k+ka+1)
0316: *
0317:                   CALL ZLARTG( AB( K+1, I-K+KA ), RA1,
0318:      $                         RWORK( I-K+KA-M ), WORK( I-K+KA-M ), RA )
0319: *
0320: *                 create nonzero element a(i-k,i-k+ka+1) outside the
0321: *                 band and store it in WORK(i-k)
0322: *
0323:                   T = -BB( KB1-K, I )*RA1
0324:                   WORK( I-K ) = RWORK( I-K+KA-M )*T -
0325:      $                          DCONJG( WORK( I-K+KA-M ) )*
0326:      $                          AB( 1, I-K+KA )
0327:                   AB( 1, I-K+KA ) = WORK( I-K+KA-M )*T +
0328:      $                              RWORK( I-K+KA-M )*AB( 1, I-K+KA )
0329:                   RA1 = RA
0330:                END IF
0331:             END IF
0332:             J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
0333:             NR = ( N-J2+KA ) / KA1
0334:             J1 = J2 + ( NR-1 )*KA1
0335:             IF( UPDATE ) THEN
0336:                J2T = MAX( J2, I+2*KA-K+1 )
0337:             ELSE
0338:                J2T = J2
0339:             END IF
0340:             NRT = ( N-J2T+KA ) / KA1
0341:             DO 90 J = J2T, J1, KA1
0342: *
0343: *              create nonzero element a(j-ka,j+1) outside the band
0344: *              and store it in WORK(j-m)
0345: *
0346:                WORK( J-M ) = WORK( J-M )*AB( 1, J+1 )
0347:                AB( 1, J+1 ) = RWORK( J-M )*AB( 1, J+1 )
0348:    90       CONTINUE
0349: *
0350: *           generate rotations in 1st set to annihilate elements which
0351: *           have been created outside the band
0352: *
0353:             IF( NRT.GT.0 )
0354:      $         CALL ZLARGV( NRT, AB( 1, J2T ), INCA, WORK( J2T-M ), KA1,
0355:      $                      RWORK( J2T-M ), KA1 )
0356:             IF( NR.GT.0 ) THEN
0357: *
0358: *              apply rotations in 1st set from the right
0359: *
0360:                DO 100 L = 1, KA - 1
0361:                   CALL ZLARTV( NR, AB( KA1-L, J2 ), INCA,
0362:      $                         AB( KA-L, J2+1 ), INCA, RWORK( J2-M ),
0363:      $                         WORK( J2-M ), KA1 )
0364:   100          CONTINUE
0365: *
0366: *              apply rotations in 1st set from both sides to diagonal
0367: *              blocks
0368: *
0369:                CALL ZLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ),
0370:      $                      AB( KA, J2+1 ), INCA, RWORK( J2-M ),
0371:      $                      WORK( J2-M ), KA1 )
0372: *
0373:                CALL ZLACGV( NR, WORK( J2-M ), KA1 )
0374:             END IF
0375: *
0376: *           start applying rotations in 1st set from the left
0377: *
0378:             DO 110 L = KA - 1, KB - K + 1, -1
0379:                NRT = ( N-J2+L ) / KA1
0380:                IF( NRT.GT.0 )
0381:      $            CALL ZLARTV( NRT, AB( L, J2+KA1-L ), INCA,
0382:      $                         AB( L+1, J2+KA1-L ), INCA, RWORK( J2-M ),
0383:      $                         WORK( J2-M ), KA1 )
0384:   110       CONTINUE
0385: *
0386:             IF( WANTX ) THEN
0387: *
0388: *              post-multiply X by product of rotations in 1st set
0389: *
0390:                DO 120 J = J2, J1, KA1
0391:                   CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
0392:      $                       RWORK( J-M ), DCONJG( WORK( J-M ) ) )
0393:   120          CONTINUE
0394:             END IF
0395:   130    CONTINUE
0396: *
0397:          IF( UPDATE ) THEN
0398:             IF( I2.LE.N .AND. KBT.GT.0 ) THEN
0399: *
0400: *              create nonzero element a(i-kbt,i-kbt+ka+1) outside the
0401: *              band and store it in WORK(i-kbt)
0402: *
0403:                WORK( I-KBT ) = -BB( KB1-KBT, I )*RA1
0404:             END IF
0405:          END IF
0406: *
0407:          DO 170 K = KB, 1, -1
0408:             IF( UPDATE ) THEN
0409:                J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1
0410:             ELSE
0411:                J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
0412:             END IF
0413: *
0414: *           finish applying rotations in 2nd set from the left
0415: *
0416:             DO 140 L = KB - K, 1, -1
0417:                NRT = ( N-J2+KA+L ) / KA1
0418:                IF( NRT.GT.0 )
0419:      $            CALL ZLARTV( NRT, AB( L, J2-L+1 ), INCA,
0420:      $                         AB( L+1, J2-L+1 ), INCA, RWORK( J2-KA ),
0421:      $                         WORK( J2-KA ), KA1 )
0422:   140       CONTINUE
0423:             NR = ( N-J2+KA ) / KA1
0424:             J1 = J2 + ( NR-1 )*KA1
0425:             DO 150 J = J1, J2, -KA1
0426:                WORK( J ) = WORK( J-KA )
0427:                RWORK( J ) = RWORK( J-KA )
0428:   150       CONTINUE
0429:             DO 160 J = J2, J1, KA1
0430: *
0431: *              create nonzero element a(j-ka,j+1) outside the band
0432: *              and store it in WORK(j)
0433: *
0434:                WORK( J ) = WORK( J )*AB( 1, J+1 )
0435:                AB( 1, J+1 ) = RWORK( J )*AB( 1, J+1 )
0436:   160       CONTINUE
0437:             IF( UPDATE ) THEN
0438:                IF( I-K.LT.N-KA .AND. K.LE.KBT )
0439:      $            WORK( I-K+KA ) = WORK( I-K )
0440:             END IF
0441:   170    CONTINUE
0442: *
0443:          DO 210 K = KB, 1, -1
0444:             J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
0445:             NR = ( N-J2+KA ) / KA1
0446:             J1 = J2 + ( NR-1 )*KA1
0447:             IF( NR.GT.0 ) THEN
0448: *
0449: *              generate rotations in 2nd set to annihilate elements
0450: *              which have been created outside the band
0451: *
0452:                CALL ZLARGV( NR, AB( 1, J2 ), INCA, WORK( J2 ), KA1,
0453:      $                      RWORK( J2 ), KA1 )
0454: *
0455: *              apply rotations in 2nd set from the right
0456: *
0457:                DO 180 L = 1, KA - 1
0458:                   CALL ZLARTV( NR, AB( KA1-L, J2 ), INCA,
0459:      $                         AB( KA-L, J2+1 ), INCA, RWORK( J2 ),
0460:      $                         WORK( J2 ), KA1 )
0461:   180          CONTINUE
0462: *
0463: *              apply rotations in 2nd set from both sides to diagonal
0464: *              blocks
0465: *
0466:                CALL ZLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ),
0467:      $                      AB( KA, J2+1 ), INCA, RWORK( J2 ),
0468:      $                      WORK( J2 ), KA1 )
0469: *
0470:                CALL ZLACGV( NR, WORK( J2 ), KA1 )
0471:             END IF
0472: *
0473: *           start applying rotations in 2nd set from the left
0474: *
0475:             DO 190 L = KA - 1, KB - K + 1, -1
0476:                NRT = ( N-J2+L ) / KA1
0477:                IF( NRT.GT.0 )
0478:      $            CALL ZLARTV( NRT, AB( L, J2+KA1-L ), INCA,
0479:      $                         AB( L+1, J2+KA1-L ), INCA, RWORK( J2 ),
0480:      $                         WORK( J2 ), KA1 )
0481:   190       CONTINUE
0482: *
0483:             IF( WANTX ) THEN
0484: *
0485: *              post-multiply X by product of rotations in 2nd set
0486: *
0487:                DO 200 J = J2, J1, KA1
0488:                   CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
0489:      $                       RWORK( J ), DCONJG( WORK( J ) ) )
0490:   200          CONTINUE
0491:             END IF
0492:   210    CONTINUE
0493: *
0494:          DO 230 K = 1, KB - 1
0495:             J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
0496: *
0497: *           finish applying rotations in 1st set from the left
0498: *
0499:             DO 220 L = KB - K, 1, -1
0500:                NRT = ( N-J2+L ) / KA1
0501:                IF( NRT.GT.0 )
0502:      $            CALL ZLARTV( NRT, AB( L, J2+KA1-L ), INCA,
0503:      $                         AB( L+1, J2+KA1-L ), INCA, RWORK( J2-M ),
0504:      $                         WORK( J2-M ), KA1 )
0505:   220       CONTINUE
0506:   230    CONTINUE
0507: *
0508:          IF( KB.GT.1 ) THEN
0509:             DO 240 J = N - 1, I2 + KA, -1
0510:                RWORK( J-M ) = RWORK( J-KA-M )
0511:                WORK( J-M ) = WORK( J-KA-M )
0512:   240       CONTINUE
0513:          END IF
0514: *
0515:       ELSE
0516: *
0517: *        Transform A, working with the lower triangle
0518: *
0519:          IF( UPDATE ) THEN
0520: *
0521: *           Form  inv(S(i))**H * A * inv(S(i))
0522: *
0523:             BII = DBLE( BB( 1, I ) )
0524:             AB( 1, I ) = ( DBLE( AB( 1, I ) ) / BII ) / BII
0525:             DO 250 J = I + 1, I1
0526:                AB( J-I+1, I ) = AB( J-I+1, I ) / BII
0527:   250       CONTINUE
0528:             DO 260 J = MAX( 1, I-KA ), I - 1
0529:                AB( I-J+1, J ) = AB( I-J+1, J ) / BII
0530:   260       CONTINUE
0531:             DO 290 K = I - KBT, I - 1
0532:                DO 270 J = I - KBT, K
0533:                   AB( K-J+1, J ) = AB( K-J+1, J ) -
0534:      $                             BB( I-J+1, J )*DCONJG( AB( I-K+1,
0535:      $                             K ) ) - DCONJG( BB( I-K+1, K ) )*
0536:      $                             AB( I-J+1, J ) + DBLE( AB( 1, I ) )*
0537:      $                             BB( I-J+1, J )*DCONJG( BB( I-K+1,
0538:      $                             K ) )
0539:   270          CONTINUE
0540:                DO 280 J = MAX( 1, I-KA ), I - KBT - 1
0541:                   AB( K-J+1, J ) = AB( K-J+1, J ) -
0542:      $                             DCONJG( BB( I-K+1, K ) )*
0543:      $                             AB( I-J+1, J )
0544:   280          CONTINUE
0545:   290       CONTINUE
0546:             DO 310 J = I, I1
0547:                DO 300 K = MAX( J-KA, I-KBT ), I - 1
0548:                   AB( J-K+1, K ) = AB( J-K+1, K ) -
0549:      $                             BB( I-K+1, K )*AB( J-I+1, I )
0550:   300          CONTINUE
0551:   310       CONTINUE
0552: *
0553:             IF( WANTX ) THEN
0554: *
0555: *              post-multiply X by inv(S(i))
0556: *
0557:                CALL ZDSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
0558:                IF( KBT.GT.0 )
0559:      $            CALL ZGERU( N-M, KBT, -CONE, X( M+1, I ), 1,
0560:      $                        BB( KBT+1, I-KBT ), LDBB-1,
0561:      $                        X( M+1, I-KBT ), LDX )
0562:             END IF
0563: *
0564: *           store a(i1,i) in RA1 for use in next loop over K
0565: *
0566:             RA1 = AB( I1-I+1, I )
0567:          END IF
0568: *
0569: *        Generate and apply vectors of rotations to chase all the
0570: *        existing bulges KA positions down toward the bottom of the
0571: *        band
0572: *
0573:          DO 360 K = 1, KB - 1
0574:             IF( UPDATE ) THEN
0575: *
0576: *              Determine the rotations which would annihilate the bulge
0577: *              which has in theory just been created
0578: *
0579:                IF( I-K+KA.LT.N .AND. I-K.GT.1 ) THEN
0580: *
0581: *                 generate rotation to annihilate a(i-k+ka+1,i)
0582: *
0583:                   CALL ZLARTG( AB( KA1-K, I ), RA1, RWORK( I-K+KA-M ),
0584:      $                         WORK( I-K+KA-M ), RA )
0585: *
0586: *                 create nonzero element a(i-k+ka+1,i-k) outside the
0587: *                 band and store it in WORK(i-k)
0588: *
0589:                   T = -BB( K+1, I-K )*RA1
0590:                   WORK( I-K ) = RWORK( I-K+KA-M )*T -
0591:      $                          DCONJG( WORK( I-K+KA-M ) )*
0592:      $                          AB( KA1, I-K )
0593:                   AB( KA1, I-K ) = WORK( I-K+KA-M )*T +
0594:      $                             RWORK( I-K+KA-M )*AB( KA1, I-K )
0595:                   RA1 = RA
0596:                END IF
0597:             END IF
0598:             J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
0599:             NR = ( N-J2+KA ) / KA1
0600:             J1 = J2 + ( NR-1 )*KA1
0601:             IF( UPDATE ) THEN
0602:                J2T = MAX( J2, I+2*KA-K+1 )
0603:             ELSE
0604:                J2T = J2
0605:             END IF
0606:             NRT = ( N-J2T+KA ) / KA1
0607:             DO 320 J = J2T, J1, KA1
0608: *
0609: *              create nonzero element a(j+1,j-ka) outside the band
0610: *              and store it in WORK(j-m)
0611: *
0612:                WORK( J-M ) = WORK( J-M )*AB( KA1, J-KA+1 )
0613:                AB( KA1, J-KA+1 ) = RWORK( J-M )*AB( KA1, J-KA+1 )
0614:   320       CONTINUE
0615: *
0616: *           generate rotations in 1st set to annihilate elements which
0617: *           have been created outside the band
0618: *
0619:             IF( NRT.GT.0 )
0620:      $         CALL ZLARGV( NRT, AB( KA1, J2T-KA ), INCA, WORK( J2T-M ),
0621:      $                      KA1, RWORK( J2T-M ), KA1 )
0622:             IF( NR.GT.0 ) THEN
0623: *
0624: *              apply rotations in 1st set from the left
0625: *
0626:                DO 330 L = 1, KA - 1
0627:                   CALL ZLARTV( NR, AB( L+1, J2-L ), INCA,
0628:      $                         AB( L+2, J2-L ), INCA, RWORK( J2-M ),
0629:      $                         WORK( J2-M ), KA1 )
0630:   330          CONTINUE
0631: *
0632: *              apply rotations in 1st set from both sides to diagonal
0633: *              blocks
0634: *
0635:                CALL ZLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
0636:      $                      INCA, RWORK( J2-M ), WORK( J2-M ), KA1 )
0637: *
0638:                CALL ZLACGV( NR, WORK( J2-M ), KA1 )
0639:             END IF
0640: *
0641: *           start applying rotations in 1st set from the right
0642: *
0643:             DO 340 L = KA - 1, KB - K + 1, -1
0644:                NRT = ( N-J2+L ) / KA1
0645:                IF( NRT.GT.0 )
0646:      $            CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
0647:      $                         AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ),
0648:      $                         WORK( J2-M ), KA1 )
0649:   340       CONTINUE
0650: *
0651:             IF( WANTX ) THEN
0652: *
0653: *              post-multiply X by product of rotations in 1st set
0654: *
0655:                DO 350 J = J2, J1, KA1
0656:                   CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
0657:      $                       RWORK( J-M ), WORK( J-M ) )
0658:   350          CONTINUE
0659:             END IF
0660:   360    CONTINUE
0661: *
0662:          IF( UPDATE ) THEN
0663:             IF( I2.LE.N .AND. KBT.GT.0 ) THEN
0664: *
0665: *              create nonzero element a(i-kbt+ka+1,i-kbt) outside the
0666: *              band and store it in WORK(i-kbt)
0667: *
0668:                WORK( I-KBT ) = -BB( KBT+1, I-KBT )*RA1
0669:             END IF
0670:          END IF
0671: *
0672:          DO 400 K = KB, 1, -1
0673:             IF( UPDATE ) THEN
0674:                J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1
0675:             ELSE
0676:                J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
0677:             END IF
0678: *
0679: *           finish applying rotations in 2nd set from the right
0680: *
0681:             DO 370 L = KB - K, 1, -1
0682:                NRT = ( N-J2+KA+L ) / KA1
0683:                IF( NRT.GT.0 )
0684:      $            CALL ZLARTV( NRT, AB( KA1-L+1, J2-KA ), INCA,
0685:      $                         AB( KA1-L, J2-KA+1 ), INCA,
0686:      $                         RWORK( J2-KA ), WORK( J2-KA ), KA1 )
0687:   370       CONTINUE
0688:             NR = ( N-J2+KA ) / KA1
0689:             J1 = J2 + ( NR-1 )*KA1
0690:             DO 380 J = J1, J2, -KA1
0691:                WORK( J ) = WORK( J-KA )
0692:                RWORK( J ) = RWORK( J-KA )
0693:   380       CONTINUE
0694:             DO 390 J = J2, J1, KA1
0695: *
0696: *              create nonzero element a(j+1,j-ka) outside the band
0697: *              and store it in WORK(j)
0698: *
0699:                WORK( J ) = WORK( J )*AB( KA1, J-KA+1 )
0700:                AB( KA1, J-KA+1 ) = RWORK( J )*AB( KA1, J-KA+1 )
0701:   390       CONTINUE
0702:             IF( UPDATE ) THEN
0703:                IF( I-K.LT.N-KA .AND. K.LE.KBT )
0704:      $            WORK( I-K+KA ) = WORK( I-K )
0705:             END IF
0706:   400    CONTINUE
0707: *
0708:          DO 440 K = KB, 1, -1
0709:             J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
0710:             NR = ( N-J2+KA ) / KA1
0711:             J1 = J2 + ( NR-1 )*KA1
0712:             IF( NR.GT.0 ) THEN
0713: *
0714: *              generate rotations in 2nd set to annihilate elements
0715: *              which have been created outside the band
0716: *
0717:                CALL ZLARGV( NR, AB( KA1, J2-KA ), INCA, WORK( J2 ), KA1,
0718:      $                      RWORK( J2 ), KA1 )
0719: *
0720: *              apply rotations in 2nd set from the left
0721: *
0722:                DO 410 L = 1, KA - 1
0723:                   CALL ZLARTV( NR, AB( L+1, J2-L ), INCA,
0724:      $                         AB( L+2, J2-L ), INCA, RWORK( J2 ),
0725:      $                         WORK( J2 ), KA1 )
0726:   410          CONTINUE
0727: *
0728: *              apply rotations in 2nd set from both sides to diagonal
0729: *              blocks
0730: *
0731:                CALL ZLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
0732:      $                      INCA, RWORK( J2 ), WORK( J2 ), KA1 )
0733: *
0734:                CALL ZLACGV( NR, WORK( J2 ), KA1 )
0735:             END IF
0736: *
0737: *           start applying rotations in 2nd set from the right
0738: *
0739:             DO 420 L = KA - 1, KB - K + 1, -1
0740:                NRT = ( N-J2+L ) / KA1
0741:                IF( NRT.GT.0 )
0742:      $            CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
0743:      $                         AB( KA1-L, J2+1 ), INCA, RWORK( J2 ),
0744:      $                         WORK( J2 ), KA1 )
0745:   420       CONTINUE
0746: *
0747:             IF( WANTX ) THEN
0748: *
0749: *              post-multiply X by product of rotations in 2nd set
0750: *
0751:                DO 430 J = J2, J1, KA1
0752:                   CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
0753:      $                       RWORK( J ), WORK( J ) )
0754:   430          CONTINUE
0755:             END IF
0756:   440    CONTINUE
0757: *
0758:          DO 460 K = 1, KB - 1
0759:             J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
0760: *
0761: *           finish applying rotations in 1st set from the right
0762: *
0763:             DO 450 L = KB - K, 1, -1
0764:                NRT = ( N-J2+L ) / KA1
0765:                IF( NRT.GT.0 )
0766:      $            CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
0767:      $                         AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ),
0768:      $                         WORK( J2-M ), KA1 )
0769:   450       CONTINUE
0770:   460    CONTINUE
0771: *
0772:          IF( KB.GT.1 ) THEN
0773:             DO 470 J = N - 1, I2 + KA, -1
0774:                RWORK( J-M ) = RWORK( J-KA-M )
0775:                WORK( J-M ) = WORK( J-KA-M )
0776:   470       CONTINUE
0777:          END IF
0778: *
0779:       END IF
0780: *
0781:       GO TO 10
0782: *
0783:   480 CONTINUE
0784: *
0785: *     **************************** Phase 2 *****************************
0786: *
0787: *     The logical structure of this phase is:
0788: *
0789: *     UPDATE = .TRUE.
0790: *     DO I = 1, M
0791: *        use S(i) to update A and create a new bulge
0792: *        apply rotations to push all bulges KA positions upward
0793: *     END DO
0794: *     UPDATE = .FALSE.
0795: *     DO I = M - KA - 1, 2, -1
0796: *        apply rotations to push all bulges KA positions upward
0797: *     END DO
0798: *
0799: *     To avoid duplicating code, the two loops are merged.
0800: *
0801:       UPDATE = .TRUE.
0802:       I = 0
0803:   490 CONTINUE
0804:       IF( UPDATE ) THEN
0805:          I = I + 1
0806:          KBT = MIN( KB, M-I )
0807:          I0 = I + 1
0808:          I1 = MAX( 1, I-KA )
0809:          I2 = I + KBT - KA1
0810:          IF( I.GT.M ) THEN
0811:             UPDATE = .FALSE.
0812:             I = I - 1
0813:             I0 = M + 1
0814:             IF( KA.EQ.0 )
0815:      $         RETURN
0816:             GO TO 490
0817:          END IF
0818:       ELSE
0819:          I = I - KA
0820:          IF( I.LT.2 )
0821:      $      RETURN
0822:       END IF
0823: *
0824:       IF( I.LT.M-KBT ) THEN
0825:          NX = M
0826:       ELSE
0827:          NX = N
0828:       END IF
0829: *
0830:       IF( UPPER ) THEN
0831: *
0832: *        Transform A, working with the upper triangle
0833: *
0834:          IF( UPDATE ) THEN
0835: *
0836: *           Form  inv(S(i))**H * A * inv(S(i))
0837: *
0838:             BII = DBLE( BB( KB1, I ) )
0839:             AB( KA1, I ) = ( DBLE( AB( KA1, I ) ) / BII ) / BII
0840:             DO 500 J = I1, I - 1
0841:                AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII
0842:   500       CONTINUE
0843:             DO 510 J = I + 1, MIN( N, I+KA )
0844:                AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII
0845:   510       CONTINUE
0846:             DO 540 K = I + 1, I + KBT
0847:                DO 520 J = K, I + KBT
0848:                   AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
0849:      $                               BB( I-J+KB1, J )*
0850:      $                               DCONJG( AB( I-K+KA1, K ) ) -
0851:      $                               DCONJG( BB( I-K+KB1, K ) )*
0852:      $                               AB( I-J+KA1, J ) +
0853:      $                               DBLE( AB( KA1, I ) )*
0854:      $                               BB( I-J+KB1, J )*
0855:      $                               DCONJG( BB( I-K+KB1, K ) )
0856:   520          CONTINUE
0857:                DO 530 J = I + KBT + 1, MIN( N, I+KA )
0858:                   AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
0859:      $                               DCONJG( BB( I-K+KB1, K ) )*
0860:      $                               AB( I-J+KA1, J )
0861:   530          CONTINUE
0862:   540       CONTINUE
0863:             DO 560 J = I1, I
0864:                DO 550 K = I + 1, MIN( J+KA, I+KBT )
0865:                   AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
0866:      $                               BB( I-K+KB1, K )*AB( J-I+KA1, I )
0867:   550          CONTINUE
0868:   560       CONTINUE
0869: *
0870:             IF( WANTX ) THEN
0871: *
0872: *              post-multiply X by inv(S(i))
0873: *
0874:                CALL ZDSCAL( NX, ONE / BII, X( 1, I ), 1 )
0875:                IF( KBT.GT.0 )
0876:      $            CALL ZGERU( NX, KBT, -CONE, X( 1, I ), 1,
0877:      $                        BB( KB, I+1 ), LDBB-1, X( 1, I+1 ), LDX )
0878:             END IF
0879: *
0880: *           store a(i1,i) in RA1 for use in next loop over K
0881: *
0882:             RA1 = AB( I1-I+KA1, I )
0883:          END IF
0884: *
0885: *        Generate and apply vectors of rotations to chase all the
0886: *        existing bulges KA positions up toward the top of the band
0887: *
0888:          DO 610 K = 1, KB - 1
0889:             IF( UPDATE ) THEN
0890: *
0891: *              Determine the rotations which would annihilate the bulge
0892: *              which has in theory just been created
0893: *
0894:                IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
0895: *
0896: *                 generate rotation to annihilate a(i+k-ka-1,i)
0897: *
0898:                   CALL ZLARTG( AB( K+1, I ), RA1, RWORK( I+K-KA ),
0899:      $                         WORK( I+K-KA ), RA )
0900: *
0901: *                 create nonzero element a(i+k-ka-1,i+k) outside the
0902: *                 band and store it in WORK(m-kb+i+k)
0903: *
0904:                   T = -BB( KB1-K, I+K )*RA1
0905:                   WORK( M-KB+I+K ) = RWORK( I+K-KA )*T -
0906:      $                               DCONJG( WORK( I+K-KA ) )*
0907:      $                               AB( 1, I+K )
0908:                   AB( 1, I+K ) = WORK( I+K-KA )*T +
0909:      $                           RWORK( I+K-KA )*AB( 1, I+K )
0910:                   RA1 = RA
0911:                END IF
0912:             END IF
0913:             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
0914:             NR = ( J2+KA-1 ) / KA1
0915:             J1 = J2 - ( NR-1 )*KA1
0916:             IF( UPDATE ) THEN
0917:                J2T = MIN( J2, I-2*KA+K-1 )
0918:             ELSE
0919:                J2T = J2
0920:             END IF
0921:             NRT = ( J2T+KA-1 ) / KA1
0922:             DO 570 J = J1, J2T, KA1
0923: *
0924: *              create nonzero element a(j-1,j+ka) outside the band
0925: *              and store it in WORK(j)
0926: *
0927:                WORK( J ) = WORK( J )*AB( 1, J+KA-1 )
0928:                AB( 1, J+KA-1 ) = RWORK( J )*AB( 1, J+KA-1 )
0929:   570       CONTINUE
0930: *
0931: *           generate rotations in 1st set to annihilate elements which
0932: *           have been created outside the band
0933: *
0934:             IF( NRT.GT.0 )
0935:      $         CALL ZLARGV( NRT, AB( 1, J1+KA ), INCA, WORK( J1 ), KA1,
0936:      $                      RWORK( J1 ), KA1 )
0937:             IF( NR.GT.0 ) THEN
0938: *
0939: *              apply rotations in 1st set from the left
0940: *
0941:                DO 580 L = 1, KA - 1
0942:                   CALL ZLARTV( NR, AB( KA1-L, J1+L ), INCA,
0943:      $                         AB( KA-L, J1+L ), INCA, RWORK( J1 ),
0944:      $                         WORK( J1 ), KA1 )
0945:   580          CONTINUE
0946: *
0947: *              apply rotations in 1st set from both sides to diagonal
0948: *              blocks
0949: *
0950:                CALL ZLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
0951:      $                      AB( KA, J1 ), INCA, RWORK( J1 ), WORK( J1 ),
0952:      $                      KA1 )
0953: *
0954:                CALL ZLACGV( NR, WORK( J1 ), KA1 )
0955:             END IF
0956: *
0957: *           start applying rotations in 1st set from the right
0958: *
0959:             DO 590 L = KA - 1, KB - K + 1, -1
0960:                NRT = ( J2+L-1 ) / KA1
0961:                J1T = J2 - ( NRT-1 )*KA1
0962:                IF( NRT.GT.0 )
0963:      $            CALL ZLARTV( NRT, AB( L, J1T ), INCA,
0964:      $                         AB( L+1, J1T-1 ), INCA, RWORK( J1T ),
0965:      $                         WORK( J1T ), KA1 )
0966:   590       CONTINUE
0967: *
0968:             IF( WANTX ) THEN
0969: *
0970: *              post-multiply X by product of rotations in 1st set
0971: *
0972:                DO 600 J = J1, J2, KA1
0973:                   CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
0974:      $                       RWORK( J ), WORK( J ) )
0975:   600          CONTINUE
0976:             END IF
0977:   610    CONTINUE
0978: *
0979:          IF( UPDATE ) THEN
0980:             IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
0981: *
0982: *              create nonzero element a(i+kbt-ka-1,i+kbt) outside the
0983: *              band and store it in WORK(m-kb+i+kbt)
0984: *
0985:                WORK( M-KB+I+KBT ) = -BB( KB1-KBT, I+KBT )*RA1
0986:             END IF
0987:          END IF
0988: *
0989:          DO 650 K = KB, 1, -1
0990:             IF( UPDATE ) THEN
0991:                J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1
0992:             ELSE
0993:                J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
0994:             END IF
0995: *
0996: *           finish applying rotations in 2nd set from the right
0997: *
0998:             DO 620 L = KB - K, 1, -1
0999:                NRT = ( J2+KA+L-1 ) / KA1
1000:                J1T = J2 - ( NRT-1 )*KA1
1001:                IF( NRT.GT.0 )
1002:      $            CALL ZLARTV( NRT, AB( L, J1T+KA ), INCA,
1003:      $                         AB( L+1, J1T+KA-1 ), INCA,
1004:      $                         RWORK( M-KB+J1T+KA ),
1005:      $                         WORK( M-KB+J1T+KA ), KA1 )
1006:   620       CONTINUE
1007:             NR = ( J2+KA-1 ) / KA1
1008:             J1 = J2 - ( NR-1 )*KA1
1009:             DO 630 J = J1, J2, KA1
1010:                WORK( M-KB+J ) = WORK( M-KB+J+KA )
1011:                RWORK( M-KB+J ) = RWORK( M-KB+J+KA )
1012:   630       CONTINUE
1013:             DO 640 J = J1, J2, KA1
1014: *
1015: *              create nonzero element a(j-1,j+ka) outside the band
1016: *              and store it in WORK(m-kb+j)
1017: *
1018:                WORK( M-KB+J ) = WORK( M-KB+J )*AB( 1, J+KA-1 )
1019:                AB( 1, J+KA-1 ) = RWORK( M-KB+J )*AB( 1, J+KA-1 )
1020:   640       CONTINUE
1021:             IF( UPDATE ) THEN
1022:                IF( I+K.GT.KA1 .AND. K.LE.KBT )
1023:      $            WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
1024:             END IF
1025:   650    CONTINUE
1026: *
1027:          DO 690 K = KB, 1, -1
1028:             J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1029:             NR = ( J2+KA-1 ) / KA1
1030:             J1 = J2 - ( NR-1 )*KA1
1031:             IF( NR.GT.0 ) THEN
1032: *
1033: *              generate rotations in 2nd set to annihilate elements
1034: *              which have been created outside the band
1035: *
1036:                CALL ZLARGV( NR, AB( 1, J1+KA ), INCA, WORK( M-KB+J1 ),
1037:      $                      KA1, RWORK( M-KB+J1 ), KA1 )
1038: *
1039: *              apply rotations in 2nd set from the left
1040: *
1041:                DO 660 L = 1, KA - 1
1042:                   CALL ZLARTV( NR, AB( KA1-L, J1+L ), INCA,
1043:      $                         AB( KA-L, J1+L ), INCA, RWORK( M-KB+J1 ),
1044:      $                         WORK( M-KB+J1 ), KA1 )
1045:   660          CONTINUE
1046: *
1047: *              apply rotations in 2nd set from both sides to diagonal
1048: *              blocks
1049: *
1050:                CALL ZLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
1051:      $                      AB( KA, J1 ), INCA, RWORK( M-KB+J1 ),
1052:      $                      WORK( M-KB+J1 ), KA1 )
1053: *
1054:                CALL ZLACGV( NR, WORK( M-KB+J1 ), KA1 )
1055:             END IF
1056: *
1057: *           start applying rotations in 2nd set from the right
1058: *
1059:             DO 670 L = KA - 1, KB - K + 1, -1
1060:                NRT = ( J2+L-1 ) / KA1
1061:                J1T = J2 - ( NRT-1 )*KA1
1062:                IF( NRT.GT.0 )
1063:      $            CALL ZLARTV( NRT, AB( L, J1T ), INCA,
1064:      $                         AB( L+1, J1T-1 ), INCA,
1065:      $                         RWORK( M-KB+J1T ), WORK( M-KB+J1T ),
1066:      $                         KA1 )
1067:   670       CONTINUE
1068: *
1069:             IF( WANTX ) THEN
1070: *
1071: *              post-multiply X by product of rotations in 2nd set
1072: *
1073:                DO 680 J = J1, J2, KA1
1074:                   CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1075:      $                       RWORK( M-KB+J ), WORK( M-KB+J ) )
1076:   680          CONTINUE
1077:             END IF
1078:   690    CONTINUE
1079: *
1080:          DO 710 K = 1, KB - 1
1081:             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
1082: *
1083: *           finish applying rotations in 1st set from the right
1084: *
1085:             DO 700 L = KB - K, 1, -1
1086:                NRT = ( J2+L-1 ) / KA1
1087:                J1T = J2 - ( NRT-1 )*KA1
1088:                IF( NRT.GT.0 )
1089:      $            CALL ZLARTV( NRT, AB( L, J1T ), INCA,
1090:      $                         AB( L+1, J1T-1 ), INCA, RWORK( J1T ),
1091:      $                         WORK( J1T ), KA1 )
1092:   700       CONTINUE
1093:   710    CONTINUE
1094: *
1095:          IF( KB.GT.1 ) THEN
1096:             DO 720 J = 2, I2 - KA
1097:                RWORK( J ) = RWORK( J+KA )
1098:                WORK( J ) = WORK( J+KA )
1099:   720       CONTINUE
1100:          END IF
1101: *
1102:       ELSE
1103: *
1104: *        Transform A, working with the lower triangle
1105: *
1106:          IF( UPDATE ) THEN
1107: *
1108: *           Form  inv(S(i))**H * A * inv(S(i))
1109: *
1110:             BII = DBLE( BB( 1, I ) )
1111:             AB( 1, I ) = ( DBLE( AB( 1, I ) ) / BII ) / BII
1112:             DO 730 J = I1, I - 1
1113:                AB( I-J+1, J ) = AB( I-J+1, J ) / BII
1114:   730       CONTINUE
1115:             DO 740 J = I + 1, MIN( N, I+KA )
1116:                AB( J-I+1, I ) = AB( J-I+1, I ) / BII
1117:   740       CONTINUE
1118:             DO 770 K = I + 1, I + KBT
1119:                DO 750 J = K, I + KBT
1120:                   AB( J-K+1, K ) = AB( J-K+1, K ) -
1121:      $                             BB( J-I+1, I )*DCONJG( AB( K-I+1,
1122:      $                             I ) ) - DCONJG( BB( K-I+1, I ) )*
1123:      $                             AB( J-I+1, I ) + DBLE( AB( 1, I ) )*
1124:      $                             BB( J-I+1, I )*DCONJG( BB( K-I+1,
1125:      $                             I ) )
1126:   750          CONTINUE
1127:                DO 760 J = I + KBT + 1, MIN( N, I+KA )
1128:                   AB( J-K+1, K ) = AB( J-K+1, K ) -
1129:      $                             DCONJG( BB( K-I+1, I ) )*
1130:      $                             AB( J-I+1, I )
1131:   760          CONTINUE
1132:   770       CONTINUE
1133:             DO 790 J = I1, I
1134:                DO 780 K = I + 1, MIN( J+KA, I+KBT )
1135:                   AB( K-J+1, J ) = AB( K-J+1, J ) -
1136:      $                             BB( K-I+1, I )*AB( I-J+1, J )
1137:   780          CONTINUE
1138:   790       CONTINUE
1139: *
1140:             IF( WANTX ) THEN
1141: *
1142: *              post-multiply X by inv(S(i))
1143: *
1144:                CALL ZDSCAL( NX, ONE / BII, X( 1, I ), 1 )
1145:                IF( KBT.GT.0 )
1146:      $            CALL ZGERC( NX, KBT, -CONE, X( 1, I ), 1, BB( 2, I ),
1147:      $                        1, X( 1, I+1 ), LDX )
1148:             END IF
1149: *
1150: *           store a(i,i1) in RA1 for use in next loop over K
1151: *
1152:             RA1 = AB( I-I1+1, I1 )
1153:          END IF
1154: *
1155: *        Generate and apply vectors of rotations to chase all the
1156: *        existing bulges KA positions up toward the top of the band
1157: *
1158:          DO 840 K = 1, KB - 1
1159:             IF( UPDATE ) THEN
1160: *
1161: *              Determine the rotations which would annihilate the bulge
1162: *              which has in theory just been created
1163: *
1164:                IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
1165: *
1166: *                 generate rotation to annihilate a(i,i+k-ka-1)
1167: *
1168:                   CALL ZLARTG( AB( KA1-K, I+K-KA ), RA1,
1169:      $                         RWORK( I+K-KA ), WORK( I+K-KA ), RA )
1170: *
1171: *                 create nonzero element a(i+k,i+k-ka-1) outside the
1172: *                 band and store it in WORK(m-kb+i+k)
1173: *
1174:                   T = -BB( K+1, I )*RA1
1175:                   WORK( M-KB+I+K ) = RWORK( I+K-KA )*T -
1176:      $                               DCONJG( WORK( I+K-KA ) )*
1177:      $                               AB( KA1, I+K-KA )
1178:                   AB( KA1, I+K-KA ) = WORK( I+K-KA )*T +
1179:      $                                RWORK( I+K-KA )*AB( KA1, I+K-KA )
1180:                   RA1 = RA
1181:                END IF
1182:             END IF
1183:             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
1184:             NR = ( J2+KA-1 ) / KA1
1185:             J1 = J2 - ( NR-1 )*KA1
1186:             IF( UPDATE ) THEN
1187:                J2T = MIN( J2, I-2*KA+K-1 )
1188:             ELSE
1189:                J2T = J2
1190:             END IF
1191:             NRT = ( J2T+KA-1 ) / KA1
1192:             DO 800 J = J1, J2T, KA1
1193: *
1194: *              create nonzero element a(j+ka,j-1) outside the band
1195: *              and store it in WORK(j)
1196: *
1197:                WORK( J ) = WORK( J )*AB( KA1, J-1 )
1198:                AB( KA1, J-1 ) = RWORK( J )*AB( KA1, J-1 )
1199:   800       CONTINUE
1200: *
1201: *           generate rotations in 1st set to annihilate elements which
1202: *           have been created outside the band
1203: *
1204:             IF( NRT.GT.0 )
1205:      $         CALL ZLARGV( NRT, AB( KA1, J1 ), INCA, WORK( J1 ), KA1,
1206:      $                      RWORK( J1 ), KA1 )
1207:             IF( NR.GT.0 ) THEN
1208: *
1209: *              apply rotations in 1st set from the right
1210: *
1211:                DO 810 L = 1, KA - 1
1212:                   CALL ZLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
1213:      $                         INCA, RWORK( J1 ), WORK( J1 ), KA1 )
1214:   810          CONTINUE
1215: *
1216: *              apply rotations in 1st set from both sides to diagonal
1217: *              blocks
1218: *
1219:                CALL ZLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
1220:      $                      AB( 2, J1-1 ), INCA, RWORK( J1 ),
1221:      $                      WORK( J1 ), KA1 )
1222: *
1223:                CALL ZLACGV( NR, WORK( J1 ), KA1 )
1224:             END IF
1225: *
1226: *           start applying rotations in 1st set from the left
1227: *
1228:             DO 820 L = KA - 1, KB - K + 1, -1
1229:                NRT = ( J2+L-1 ) / KA1
1230:                J1T = J2 - ( NRT-1 )*KA1
1231:                IF( NRT.GT.0 )
1232:      $            CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
1233:      $                         AB( KA1-L, J1T-KA1+L ), INCA,
1234:      $                         RWORK( J1T ), WORK( J1T ), KA1 )
1235:   820       CONTINUE
1236: *
1237:             IF( WANTX ) THEN
1238: *
1239: *              post-multiply X by product of rotations in 1st set
1240: *
1241:                DO 830 J = J1, J2, KA1
1242:                   CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1243:      $                       RWORK( J ), DCONJG( WORK( J ) ) )
1244:   830          CONTINUE
1245:             END IF
1246:   840    CONTINUE
1247: *
1248:          IF( UPDATE ) THEN
1249:             IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
1250: *
1251: *              create nonzero element a(i+kbt,i+kbt-ka-1) outside the
1252: *              band and store it in WORK(m-kb+i+kbt)
1253: *
1254:                WORK( M-KB+I+KBT ) = -BB( KBT+1, I )*RA1
1255:             END IF
1256:          END IF
1257: *
1258:          DO 880 K = KB, 1, -1
1259:             IF( UPDATE ) THEN
1260:                J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1
1261:             ELSE
1262:                J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1263:             END IF
1264: *
1265: *           finish applying rotations in 2nd set from the left
1266: *
1267:             DO 850 L = KB - K, 1, -1
1268:                NRT = ( J2+KA+L-1 ) / KA1
1269:                J1T = J2 - ( NRT-1 )*KA1
1270:                IF( NRT.GT.0 )
1271:      $            CALL ZLARTV( NRT, AB( KA1-L+1, J1T+L-1 ), INCA,
1272:      $                         AB( KA1-L, J1T+L-1 ), INCA,
1273:      $                         RWORK( M-KB+J1T+KA ),
1274:      $                         WORK( M-KB+J1T+KA ), KA1 )
1275:   850       CONTINUE
1276:             NR = ( J2+KA-1 ) / KA1
1277:             J1 = J2 - ( NR-1 )*KA1
1278:             DO 860 J = J1, J2, KA1
1279:                WORK( M-KB+J ) = WORK( M-KB+J+KA )
1280:                RWORK( M-KB+J ) = RWORK( M-KB+J+KA )
1281:   860       CONTINUE
1282:             DO 870 J = J1, J2, KA1
1283: *
1284: *              create nonzero element a(j+ka,j-1) outside the band
1285: *              and store it in WORK(m-kb+j)
1286: *
1287:                WORK( M-KB+J ) = WORK( M-KB+J )*AB( KA1, J-1 )
1288:                AB( KA1, J-1 ) = RWORK( M-KB+J )*AB( KA1, J-1 )
1289:   870       CONTINUE
1290:             IF( UPDATE ) THEN
1291:                IF( I+K.GT.KA1 .AND. K.LE.KBT )
1292:      $            WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
1293:             END IF
1294:   880    CONTINUE
1295: *
1296:          DO 920 K = KB, 1, -1
1297:             J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1298:             NR = ( J2+KA-1 ) / KA1
1299:             J1 = J2 - ( NR-1 )*KA1
1300:             IF( NR.GT.0 ) THEN
1301: *
1302: *              generate rotations in 2nd set to annihilate elements
1303: *              which have been created outside the band
1304: *
1305:                CALL ZLARGV( NR, AB( KA1, J1 ), INCA, WORK( M-KB+J1 ),
1306:      $                      KA1, RWORK( M-KB+J1 ), KA1 )
1307: *
1308: *              apply rotations in 2nd set from the right
1309: *
1310:                DO 890 L = 1, KA - 1
1311:                   CALL ZLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
1312:      $                         INCA, RWORK( M-KB+J1 ), WORK( M-KB+J1 ),
1313:      $                         KA1 )
1314:   890          CONTINUE
1315: *
1316: *              apply rotations in 2nd set from both sides to diagonal
1317: *              blocks
1318: *
1319:                CALL ZLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
1320:      $                      AB( 2, J1-1 ), INCA, RWORK( M-KB+J1 ),
1321:      $                      WORK( M-KB+J1 ), KA1 )
1322: *
1323:                CALL ZLACGV( NR, WORK( M-KB+J1 ), KA1 )
1324:             END IF
1325: *
1326: *           start applying rotations in 2nd set from the left
1327: *
1328:             DO 900 L = KA - 1, KB - K + 1, -1
1329:                NRT = ( J2+L-1 ) / KA1
1330:                J1T = J2 - ( NRT-1 )*KA1
1331:                IF( NRT.GT.0 )
1332:      $            CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
1333:      $                         AB( KA1-L, J1T-KA1+L ), INCA,
1334:      $                         RWORK( M-KB+J1T ), WORK( M-KB+J1T ),
1335:      $                         KA1 )
1336:   900       CONTINUE
1337: *
1338:             IF( WANTX ) THEN
1339: *
1340: *              post-multiply X by product of rotations in 2nd set
1341: *
1342:                DO 910 J = J1, J2, KA1
1343:                   CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1344:      $                       RWORK( M-KB+J ), DCONJG( WORK( M-KB+J ) ) )
1345:   910          CONTINUE
1346:             END IF
1347:   920    CONTINUE
1348: *
1349:          DO 940 K = 1, KB - 1
1350:             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
1351: *
1352: *           finish applying rotations in 1st set from the left
1353: *
1354:             DO 930 L = KB - K, 1, -1
1355:                NRT = ( J2+L-1 ) / KA1
1356:                J1T = J2 - ( NRT-1 )*KA1
1357:                IF( NRT.GT.0 )
1358:      $            CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
1359:      $                         AB( KA1-L, J1T-KA1+L ), INCA,
1360:      $                         RWORK( J1T ), WORK( J1T ), KA1 )
1361:   930       CONTINUE
1362:   940    CONTINUE
1363: *
1364:          IF( KB.GT.1 ) THEN
1365:             DO 950 J = 2, I2 - KA
1366:                RWORK( J ) = RWORK( J+KA )
1367:                WORK( J ) = WORK( J+KA )
1368:   950       CONTINUE
1369:          END IF
1370: *
1371:       END IF
1372: *
1373:       GO TO 490
1374: *
1375: *     End of ZHBGST
1376: *
1377:       END
1378: