0001:       SUBROUTINE CHBGST( 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:       REAL               RWORK( * )
0014:       COMPLEX            AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
0015:      $                   X( LDX, * )
0016: *     ..
0017: *
0018: *  Purpose
0019: *  =======
0020: *
0021: *  CHBGST 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 CPBSTF, 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 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 array, dimension (LDBB,N)
0067: *          The banded factor S from the split Cholesky factorization of
0068: *          B, as returned by CPBSTF, 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 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 array, dimension (N)
0083: *
0084: *  RWORK   (workspace) REAL 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            CZERO, CONE
0094:       REAL               ONE
0095:       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
0096:      $                   CONE = ( 1.0E+0, 0.0E+0 ), ONE = 1.0E+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:       REAL               BII
0103:       COMPLEX            RA, RA1, T
0104: *     ..
0105: *     .. External Functions ..
0106:       LOGICAL            LSAME
0107:       EXTERNAL           LSAME
0108: *     ..
0109: *     .. External Subroutines ..
0110:       EXTERNAL           CGERC, CGERU, CLACGV, CLAR2V, CLARGV, CLARTG,
0111:      $                   CLARTV, CLASET, CROT, CSSCAL, XERBLA
0112: *     ..
0113: *     .. Intrinsic Functions ..
0114:       INTRINSIC          CONJG, MAX, MIN, REAL
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( 'CHBGST', -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 CLASET( '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 CPBSTF. 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 = REAL( BB( KB1, I ) )
0256:             AB( KA1, I ) = ( REAL( 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:      $                               CONJG( AB( K-I+KA1, I ) ) -
0268:      $                               CONJG( BB( K-I+KB1, I ) )*
0269:      $                               AB( J-I+KA1, I ) +
0270:      $                               REAL( AB( KA1, I ) )*
0271:      $                               BB( J-I+KB1, I )*
0272:      $                               CONJG( 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:      $                               CONJG( 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 CSSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
0292:                IF( KBT.GT.0 )
0293:      $            CALL CGERC( 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 CLARTG( 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:      $                          CONJG( 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 CLARGV( 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 CLARTV( 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 CLAR2V( 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 CLACGV( 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 CLARTV( 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 CROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
0392:      $                       RWORK( J-M ), CONJG( 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 CLARTV( 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 CLARGV( 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 CLARTV( 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 CLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ),
0467:      $                      AB( KA, J2+1 ), INCA, RWORK( J2 ),
0468:      $                      WORK( J2 ), KA1 )
0469: *
0470:                CALL CLACGV( 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 CLARTV( 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 CROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
0489:      $                       RWORK( J ), CONJG( 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 CLARTV( 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 = REAL( BB( 1, I ) )
0524:             AB( 1, I ) = ( REAL( 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 )*CONJG( AB( I-K+1,
0535:      $                             K ) ) - CONJG( BB( I-K+1, K ) )*
0536:      $                             AB( I-J+1, J ) + REAL( AB( 1, I ) )*
0537:      $                             BB( I-J+1, J )*CONJG( 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:      $                             CONJG( 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 CSSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
0558:                IF( KBT.GT.0 )
0559:      $            CALL CGERU( 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 CLARTG( 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:      $                          CONJG( WORK( I-K+KA-M ) )*AB( KA1, I-K )
0592:                   AB( KA1, I-K ) = WORK( I-K+KA-M )*T +
0593:      $                             RWORK( I-K+KA-M )*AB( KA1, I-K )
0594:                   RA1 = RA
0595:                END IF
0596:             END IF
0597:             J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
0598:             NR = ( N-J2+KA ) / KA1
0599:             J1 = J2 + ( NR-1 )*KA1
0600:             IF( UPDATE ) THEN
0601:                J2T = MAX( J2, I+2*KA-K+1 )
0602:             ELSE
0603:                J2T = J2
0604:             END IF
0605:             NRT = ( N-J2T+KA ) / KA1
0606:             DO 320 J = J2T, J1, KA1
0607: *
0608: *              create nonzero element a(j+1,j-ka) outside the band
0609: *              and store it in WORK(j-m)
0610: *
0611:                WORK( J-M ) = WORK( J-M )*AB( KA1, J-KA+1 )
0612:                AB( KA1, J-KA+1 ) = RWORK( J-M )*AB( KA1, J-KA+1 )
0613:   320       CONTINUE
0614: *
0615: *           generate rotations in 1st set to annihilate elements which
0616: *           have been created outside the band
0617: *
0618:             IF( NRT.GT.0 )
0619:      $         CALL CLARGV( NRT, AB( KA1, J2T-KA ), INCA, WORK( J2T-M ),
0620:      $                      KA1, RWORK( J2T-M ), KA1 )
0621:             IF( NR.GT.0 ) THEN
0622: *
0623: *              apply rotations in 1st set from the left
0624: *
0625:                DO 330 L = 1, KA - 1
0626:                   CALL CLARTV( NR, AB( L+1, J2-L ), INCA,
0627:      $                         AB( L+2, J2-L ), INCA, RWORK( J2-M ),
0628:      $                         WORK( J2-M ), KA1 )
0629:   330          CONTINUE
0630: *
0631: *              apply rotations in 1st set from both sides to diagonal
0632: *              blocks
0633: *
0634:                CALL CLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
0635:      $                      INCA, RWORK( J2-M ), WORK( J2-M ), KA1 )
0636: *
0637:                CALL CLACGV( NR, WORK( J2-M ), KA1 )
0638:             END IF
0639: *
0640: *           start applying rotations in 1st set from the right
0641: *
0642:             DO 340 L = KA - 1, KB - K + 1, -1
0643:                NRT = ( N-J2+L ) / KA1
0644:                IF( NRT.GT.0 )
0645:      $            CALL CLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
0646:      $                         AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ),
0647:      $                         WORK( J2-M ), KA1 )
0648:   340       CONTINUE
0649: *
0650:             IF( WANTX ) THEN
0651: *
0652: *              post-multiply X by product of rotations in 1st set
0653: *
0654:                DO 350 J = J2, J1, KA1
0655:                   CALL CROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
0656:      $                       RWORK( J-M ), WORK( J-M ) )
0657:   350          CONTINUE
0658:             END IF
0659:   360    CONTINUE
0660: *
0661:          IF( UPDATE ) THEN
0662:             IF( I2.LE.N .AND. KBT.GT.0 ) THEN
0663: *
0664: *              create nonzero element a(i-kbt+ka+1,i-kbt) outside the
0665: *              band and store it in WORK(i-kbt)
0666: *
0667:                WORK( I-KBT ) = -BB( KBT+1, I-KBT )*RA1
0668:             END IF
0669:          END IF
0670: *
0671:          DO 400 K = KB, 1, -1
0672:             IF( UPDATE ) THEN
0673:                J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1
0674:             ELSE
0675:                J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
0676:             END IF
0677: *
0678: *           finish applying rotations in 2nd set from the right
0679: *
0680:             DO 370 L = KB - K, 1, -1
0681:                NRT = ( N-J2+KA+L ) / KA1
0682:                IF( NRT.GT.0 )
0683:      $            CALL CLARTV( NRT, AB( KA1-L+1, J2-KA ), INCA,
0684:      $                         AB( KA1-L, J2-KA+1 ), INCA,
0685:      $                         RWORK( J2-KA ), WORK( J2-KA ), KA1 )
0686:   370       CONTINUE
0687:             NR = ( N-J2+KA ) / KA1
0688:             J1 = J2 + ( NR-1 )*KA1
0689:             DO 380 J = J1, J2, -KA1
0690:                WORK( J ) = WORK( J-KA )
0691:                RWORK( J ) = RWORK( J-KA )
0692:   380       CONTINUE
0693:             DO 390 J = J2, J1, KA1
0694: *
0695: *              create nonzero element a(j+1,j-ka) outside the band
0696: *              and store it in WORK(j)
0697: *
0698:                WORK( J ) = WORK( J )*AB( KA1, J-KA+1 )
0699:                AB( KA1, J-KA+1 ) = RWORK( J )*AB( KA1, J-KA+1 )
0700:   390       CONTINUE
0701:             IF( UPDATE ) THEN
0702:                IF( I-K.LT.N-KA .AND. K.LE.KBT )
0703:      $            WORK( I-K+KA ) = WORK( I-K )
0704:             END IF
0705:   400    CONTINUE
0706: *
0707:          DO 440 K = KB, 1, -1
0708:             J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
0709:             NR = ( N-J2+KA ) / KA1
0710:             J1 = J2 + ( NR-1 )*KA1
0711:             IF( NR.GT.0 ) THEN
0712: *
0713: *              generate rotations in 2nd set to annihilate elements
0714: *              which have been created outside the band
0715: *
0716:                CALL CLARGV( NR, AB( KA1, J2-KA ), INCA, WORK( J2 ), KA1,
0717:      $                      RWORK( J2 ), KA1 )
0718: *
0719: *              apply rotations in 2nd set from the left
0720: *
0721:                DO 410 L = 1, KA - 1
0722:                   CALL CLARTV( NR, AB( L+1, J2-L ), INCA,
0723:      $                         AB( L+2, J2-L ), INCA, RWORK( J2 ),
0724:      $                         WORK( J2 ), KA1 )
0725:   410          CONTINUE
0726: *
0727: *              apply rotations in 2nd set from both sides to diagonal
0728: *              blocks
0729: *
0730:                CALL CLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
0731:      $                      INCA, RWORK( J2 ), WORK( J2 ), KA1 )
0732: *
0733:                CALL CLACGV( NR, WORK( J2 ), KA1 )
0734:             END IF
0735: *
0736: *           start applying rotations in 2nd set from the right
0737: *
0738:             DO 420 L = KA - 1, KB - K + 1, -1
0739:                NRT = ( N-J2+L ) / KA1
0740:                IF( NRT.GT.0 )
0741:      $            CALL CLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
0742:      $                         AB( KA1-L, J2+1 ), INCA, RWORK( J2 ),
0743:      $                         WORK( J2 ), KA1 )
0744:   420       CONTINUE
0745: *
0746:             IF( WANTX ) THEN
0747: *
0748: *              post-multiply X by product of rotations in 2nd set
0749: *
0750:                DO 430 J = J2, J1, KA1
0751:                   CALL CROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
0752:      $                       RWORK( J ), WORK( J ) )
0753:   430          CONTINUE
0754:             END IF
0755:   440    CONTINUE
0756: *
0757:          DO 460 K = 1, KB - 1
0758:             J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
0759: *
0760: *           finish applying rotations in 1st set from the right
0761: *
0762:             DO 450 L = KB - K, 1, -1
0763:                NRT = ( N-J2+L ) / KA1
0764:                IF( NRT.GT.0 )
0765:      $            CALL CLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
0766:      $                         AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ),
0767:      $                         WORK( J2-M ), KA1 )
0768:   450       CONTINUE
0769:   460    CONTINUE
0770: *
0771:          IF( KB.GT.1 ) THEN
0772:             DO 470 J = N - 1, I2 + KA, -1
0773:                RWORK( J-M ) = RWORK( J-KA-M )
0774:                WORK( J-M ) = WORK( J-KA-M )
0775:   470       CONTINUE
0776:          END IF
0777: *
0778:       END IF
0779: *
0780:       GO TO 10
0781: *
0782:   480 CONTINUE
0783: *
0784: *     **************************** Phase 2 *****************************
0785: *
0786: *     The logical structure of this phase is:
0787: *
0788: *     UPDATE = .TRUE.
0789: *     DO I = 1, M
0790: *        use S(i) to update A and create a new bulge
0791: *        apply rotations to push all bulges KA positions upward
0792: *     END DO
0793: *     UPDATE = .FALSE.
0794: *     DO I = M - KA - 1, 2, -1
0795: *        apply rotations to push all bulges KA positions upward
0796: *     END DO
0797: *
0798: *     To avoid duplicating code, the two loops are merged.
0799: *
0800:       UPDATE = .TRUE.
0801:       I = 0
0802:   490 CONTINUE
0803:       IF( UPDATE ) THEN
0804:          I = I + 1
0805:          KBT = MIN( KB, M-I )
0806:          I0 = I + 1
0807:          I1 = MAX( 1, I-KA )
0808:          I2 = I + KBT - KA1
0809:          IF( I.GT.M ) THEN
0810:             UPDATE = .FALSE.
0811:             I = I - 1
0812:             I0 = M + 1
0813:             IF( KA.EQ.0 )
0814:      $         RETURN
0815:             GO TO 490
0816:          END IF
0817:       ELSE
0818:          I = I - KA
0819:          IF( I.LT.2 )
0820:      $      RETURN
0821:       END IF
0822: *
0823:       IF( I.LT.M-KBT ) THEN
0824:          NX = M
0825:       ELSE
0826:          NX = N
0827:       END IF
0828: *
0829:       IF( UPPER ) THEN
0830: *
0831: *        Transform A, working with the upper triangle
0832: *
0833:          IF( UPDATE ) THEN
0834: *
0835: *           Form  inv(S(i))**H * A * inv(S(i))
0836: *
0837:             BII = REAL( BB( KB1, I ) )
0838:             AB( KA1, I ) = ( REAL( AB( KA1, I ) ) / BII ) / BII
0839:             DO 500 J = I1, I - 1
0840:                AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII
0841:   500       CONTINUE
0842:             DO 510 J = I + 1, MIN( N, I+KA )
0843:                AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII
0844:   510       CONTINUE
0845:             DO 540 K = I + 1, I + KBT
0846:                DO 520 J = K, I + KBT
0847:                   AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
0848:      $                               BB( I-J+KB1, J )*
0849:      $                               CONJG( AB( I-K+KA1, K ) ) -
0850:      $                               CONJG( BB( I-K+KB1, K ) )*
0851:      $                               AB( I-J+KA1, J ) +
0852:      $                               REAL( AB( KA1, I ) )*
0853:      $                               BB( I-J+KB1, J )*
0854:      $                               CONJG( BB( I-K+KB1, K ) )
0855:   520          CONTINUE
0856:                DO 530 J = I + KBT + 1, MIN( N, I+KA )
0857:                   AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
0858:      $                               CONJG( BB( I-K+KB1, K ) )*
0859:      $                               AB( I-J+KA1, J )
0860:   530          CONTINUE
0861:   540       CONTINUE
0862:             DO 560 J = I1, I
0863:                DO 550 K = I + 1, MIN( J+KA, I+KBT )
0864:                   AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
0865:      $                               BB( I-K+KB1, K )*AB( J-I+KA1, I )
0866:   550          CONTINUE
0867:   560       CONTINUE
0868: *
0869:             IF( WANTX ) THEN
0870: *
0871: *              post-multiply X by inv(S(i))
0872: *
0873:                CALL CSSCAL( NX, ONE / BII, X( 1, I ), 1 )
0874:                IF( KBT.GT.0 )
0875:      $            CALL CGERU( NX, KBT, -CONE, X( 1, I ), 1,
0876:      $                        BB( KB, I+1 ), LDBB-1, X( 1, I+1 ), LDX )
0877:             END IF
0878: *
0879: *           store a(i1,i) in RA1 for use in next loop over K
0880: *
0881:             RA1 = AB( I1-I+KA1, I )
0882:          END IF
0883: *
0884: *        Generate and apply vectors of rotations to chase all the
0885: *        existing bulges KA positions up toward the top of the band
0886: *
0887:          DO 610 K = 1, KB - 1
0888:             IF( UPDATE ) THEN
0889: *
0890: *              Determine the rotations which would annihilate the bulge
0891: *              which has in theory just been created
0892: *
0893:                IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
0894: *
0895: *                 generate rotation to annihilate a(i+k-ka-1,i)
0896: *
0897:                   CALL CLARTG( AB( K+1, I ), RA1, RWORK( I+K-KA ),
0898:      $                         WORK( I+K-KA ), RA )
0899: *
0900: *                 create nonzero element a(i+k-ka-1,i+k) outside the
0901: *                 band and store it in WORK(m-kb+i+k)
0902: *
0903:                   T = -BB( KB1-K, I+K )*RA1
0904:                   WORK( M-KB+I+K ) = RWORK( I+K-KA )*T -
0905:      $                               CONJG( WORK( I+K-KA ) )*
0906:      $                               AB( 1, I+K )
0907:                   AB( 1, I+K ) = WORK( I+K-KA )*T +
0908:      $                           RWORK( I+K-KA )*AB( 1, I+K )
0909:                   RA1 = RA
0910:                END IF
0911:             END IF
0912:             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
0913:             NR = ( J2+KA-1 ) / KA1
0914:             J1 = J2 - ( NR-1 )*KA1
0915:             IF( UPDATE ) THEN
0916:                J2T = MIN( J2, I-2*KA+K-1 )
0917:             ELSE
0918:                J2T = J2
0919:             END IF
0920:             NRT = ( J2T+KA-1 ) / KA1
0921:             DO 570 J = J1, J2T, KA1
0922: *
0923: *              create nonzero element a(j-1,j+ka) outside the band
0924: *              and store it in WORK(j)
0925: *
0926:                WORK( J ) = WORK( J )*AB( 1, J+KA-1 )
0927:                AB( 1, J+KA-1 ) = RWORK( J )*AB( 1, J+KA-1 )
0928:   570       CONTINUE
0929: *
0930: *           generate rotations in 1st set to annihilate elements which
0931: *           have been created outside the band
0932: *
0933:             IF( NRT.GT.0 )
0934:      $         CALL CLARGV( NRT, AB( 1, J1+KA ), INCA, WORK( J1 ), KA1,
0935:      $                      RWORK( J1 ), KA1 )
0936:             IF( NR.GT.0 ) THEN
0937: *
0938: *              apply rotations in 1st set from the left
0939: *
0940:                DO 580 L = 1, KA - 1
0941:                   CALL CLARTV( NR, AB( KA1-L, J1+L ), INCA,
0942:      $                         AB( KA-L, J1+L ), INCA, RWORK( J1 ),
0943:      $                         WORK( J1 ), KA1 )
0944:   580          CONTINUE
0945: *
0946: *              apply rotations in 1st set from both sides to diagonal
0947: *              blocks
0948: *
0949:                CALL CLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
0950:      $                      AB( KA, J1 ), INCA, RWORK( J1 ), WORK( J1 ),
0951:      $                      KA1 )
0952: *
0953:                CALL CLACGV( NR, WORK( J1 ), KA1 )
0954:             END IF
0955: *
0956: *           start applying rotations in 1st set from the right
0957: *
0958:             DO 590 L = KA - 1, KB - K + 1, -1
0959:                NRT = ( J2+L-1 ) / KA1
0960:                J1T = J2 - ( NRT-1 )*KA1
0961:                IF( NRT.GT.0 )
0962:      $            CALL CLARTV( NRT, AB( L, J1T ), INCA,
0963:      $                         AB( L+1, J1T-1 ), INCA, RWORK( J1T ),
0964:      $                         WORK( J1T ), KA1 )
0965:   590       CONTINUE
0966: *
0967:             IF( WANTX ) THEN
0968: *
0969: *              post-multiply X by product of rotations in 1st set
0970: *
0971:                DO 600 J = J1, J2, KA1
0972:                   CALL CROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
0973:      $                       RWORK( J ), WORK( J ) )
0974:   600          CONTINUE
0975:             END IF
0976:   610    CONTINUE
0977: *
0978:          IF( UPDATE ) THEN
0979:             IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
0980: *
0981: *              create nonzero element a(i+kbt-ka-1,i+kbt) outside the
0982: *              band and store it in WORK(m-kb+i+kbt)
0983: *
0984:                WORK( M-KB+I+KBT ) = -BB( KB1-KBT, I+KBT )*RA1
0985:             END IF
0986:          END IF
0987: *
0988:          DO 650 K = KB, 1, -1
0989:             IF( UPDATE ) THEN
0990:                J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1
0991:             ELSE
0992:                J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
0993:             END IF
0994: *
0995: *           finish applying rotations in 2nd set from the right
0996: *
0997:             DO 620 L = KB - K, 1, -1
0998:                NRT = ( J2+KA+L-1 ) / KA1
0999:                J1T = J2 - ( NRT-1 )*KA1
1000:                IF( NRT.GT.0 )
1001:      $            CALL CLARTV( NRT, AB( L, J1T+KA ), INCA,
1002:      $                         AB( L+1, J1T+KA-1 ), INCA,
1003:      $                         RWORK( M-KB+J1T+KA ),
1004:      $                         WORK( M-KB+J1T+KA ), KA1 )
1005:   620       CONTINUE
1006:             NR = ( J2+KA-1 ) / KA1
1007:             J1 = J2 - ( NR-1 )*KA1
1008:             DO 630 J = J1, J2, KA1
1009:                WORK( M-KB+J ) = WORK( M-KB+J+KA )
1010:                RWORK( M-KB+J ) = RWORK( M-KB+J+KA )
1011:   630       CONTINUE
1012:             DO 640 J = J1, J2, KA1
1013: *
1014: *              create nonzero element a(j-1,j+ka) outside the band
1015: *              and store it in WORK(m-kb+j)
1016: *
1017:                WORK( M-KB+J ) = WORK( M-KB+J )*AB( 1, J+KA-1 )
1018:                AB( 1, J+KA-1 ) = RWORK( M-KB+J )*AB( 1, J+KA-1 )
1019:   640       CONTINUE
1020:             IF( UPDATE ) THEN
1021:                IF( I+K.GT.KA1 .AND. K.LE.KBT )
1022:      $            WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
1023:             END IF
1024:   650    CONTINUE
1025: *
1026:          DO 690 K = KB, 1, -1
1027:             J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1028:             NR = ( J2+KA-1 ) / KA1
1029:             J1 = J2 - ( NR-1 )*KA1
1030:             IF( NR.GT.0 ) THEN
1031: *
1032: *              generate rotations in 2nd set to annihilate elements
1033: *              which have been created outside the band
1034: *
1035:                CALL CLARGV( NR, AB( 1, J1+KA ), INCA, WORK( M-KB+J1 ),
1036:      $                      KA1, RWORK( M-KB+J1 ), KA1 )
1037: *
1038: *              apply rotations in 2nd set from the left
1039: *
1040:                DO 660 L = 1, KA - 1
1041:                   CALL CLARTV( NR, AB( KA1-L, J1+L ), INCA,
1042:      $                         AB( KA-L, J1+L ), INCA, RWORK( M-KB+J1 ),
1043:      $                         WORK( M-KB+J1 ), KA1 )
1044:   660          CONTINUE
1045: *
1046: *              apply rotations in 2nd set from both sides to diagonal
1047: *              blocks
1048: *
1049:                CALL CLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
1050:      $                      AB( KA, J1 ), INCA, RWORK( M-KB+J1 ),
1051:      $                      WORK( M-KB+J1 ), KA1 )
1052: *
1053:                CALL CLACGV( NR, WORK( M-KB+J1 ), KA1 )
1054:             END IF
1055: *
1056: *           start applying rotations in 2nd set from the right
1057: *
1058:             DO 670 L = KA - 1, KB - K + 1, -1
1059:                NRT = ( J2+L-1 ) / KA1
1060:                J1T = J2 - ( NRT-1 )*KA1
1061:                IF( NRT.GT.0 )
1062:      $            CALL CLARTV( NRT, AB( L, J1T ), INCA,
1063:      $                         AB( L+1, J1T-1 ), INCA,
1064:      $                         RWORK( M-KB+J1T ), WORK( M-KB+J1T ),
1065:      $                         KA1 )
1066:   670       CONTINUE
1067: *
1068:             IF( WANTX ) THEN
1069: *
1070: *              post-multiply X by product of rotations in 2nd set
1071: *
1072:                DO 680 J = J1, J2, KA1
1073:                   CALL CROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1074:      $                       RWORK( M-KB+J ), WORK( M-KB+J ) )
1075:   680          CONTINUE
1076:             END IF
1077:   690    CONTINUE
1078: *
1079:          DO 710 K = 1, KB - 1
1080:             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
1081: *
1082: *           finish applying rotations in 1st set from the right
1083: *
1084:             DO 700 L = KB - K, 1, -1
1085:                NRT = ( J2+L-1 ) / KA1
1086:                J1T = J2 - ( NRT-1 )*KA1
1087:                IF( NRT.GT.0 )
1088:      $            CALL CLARTV( NRT, AB( L, J1T ), INCA,
1089:      $                         AB( L+1, J1T-1 ), INCA, RWORK( J1T ),
1090:      $                         WORK( J1T ), KA1 )
1091:   700       CONTINUE
1092:   710    CONTINUE
1093: *
1094:          IF( KB.GT.1 ) THEN
1095:             DO 720 J = 2, I2 - KA
1096:                RWORK( J ) = RWORK( J+KA )
1097:                WORK( J ) = WORK( J+KA )
1098:   720       CONTINUE
1099:          END IF
1100: *
1101:       ELSE
1102: *
1103: *        Transform A, working with the lower triangle
1104: *
1105:          IF( UPDATE ) THEN
1106: *
1107: *           Form  inv(S(i))**H * A * inv(S(i))
1108: *
1109:             BII = REAL( BB( 1, I ) )
1110:             AB( 1, I ) = ( REAL( AB( 1, I ) ) / BII ) / BII
1111:             DO 730 J = I1, I - 1
1112:                AB( I-J+1, J ) = AB( I-J+1, J ) / BII
1113:   730       CONTINUE
1114:             DO 740 J = I + 1, MIN( N, I+KA )
1115:                AB( J-I+1, I ) = AB( J-I+1, I ) / BII
1116:   740       CONTINUE
1117:             DO 770 K = I + 1, I + KBT
1118:                DO 750 J = K, I + KBT
1119:                   AB( J-K+1, K ) = AB( J-K+1, K ) -
1120:      $                             BB( J-I+1, I )*CONJG( AB( K-I+1,
1121:      $                             I ) ) - CONJG( BB( K-I+1, I ) )*
1122:      $                             AB( J-I+1, I ) + REAL( AB( 1, I ) )*
1123:      $                             BB( J-I+1, I )*CONJG( BB( K-I+1,
1124:      $                             I ) )
1125:   750          CONTINUE
1126:                DO 760 J = I + KBT + 1, MIN( N, I+KA )
1127:                   AB( J-K+1, K ) = AB( J-K+1, K ) -
1128:      $                             CONJG( BB( K-I+1, I ) )*
1129:      $                             AB( J-I+1, I )
1130:   760          CONTINUE
1131:   770       CONTINUE
1132:             DO 790 J = I1, I
1133:                DO 780 K = I + 1, MIN( J+KA, I+KBT )
1134:                   AB( K-J+1, J ) = AB( K-J+1, J ) -
1135:      $                             BB( K-I+1, I )*AB( I-J+1, J )
1136:   780          CONTINUE
1137:   790       CONTINUE
1138: *
1139:             IF( WANTX ) THEN
1140: *
1141: *              post-multiply X by inv(S(i))
1142: *
1143:                CALL CSSCAL( NX, ONE / BII, X( 1, I ), 1 )
1144:                IF( KBT.GT.0 )
1145:      $            CALL CGERC( NX, KBT, -CONE, X( 1, I ), 1, BB( 2, I ),
1146:      $                        1, X( 1, I+1 ), LDX )
1147:             END IF
1148: *
1149: *           store a(i,i1) in RA1 for use in next loop over K
1150: *
1151:             RA1 = AB( I-I1+1, I1 )
1152:          END IF
1153: *
1154: *        Generate and apply vectors of rotations to chase all the
1155: *        existing bulges KA positions up toward the top of the band
1156: *
1157:          DO 840 K = 1, KB - 1
1158:             IF( UPDATE ) THEN
1159: *
1160: *              Determine the rotations which would annihilate the bulge
1161: *              which has in theory just been created
1162: *
1163:                IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
1164: *
1165: *                 generate rotation to annihilate a(i,i+k-ka-1)
1166: *
1167:                   CALL CLARTG( AB( KA1-K, I+K-KA ), RA1,
1168:      $                         RWORK( I+K-KA ), WORK( I+K-KA ), RA )
1169: *
1170: *                 create nonzero element a(i+k,i+k-ka-1) outside the
1171: *                 band and store it in WORK(m-kb+i+k)
1172: *
1173:                   T = -BB( K+1, I )*RA1
1174:                   WORK( M-KB+I+K ) = RWORK( I+K-KA )*T -
1175:      $                               CONJG( WORK( I+K-KA ) )*
1176:      $                               AB( KA1, I+K-KA )
1177:                   AB( KA1, I+K-KA ) = WORK( I+K-KA )*T +
1178:      $                                RWORK( I+K-KA )*AB( KA1, I+K-KA )
1179:                   RA1 = RA
1180:                END IF
1181:             END IF
1182:             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
1183:             NR = ( J2+KA-1 ) / KA1
1184:             J1 = J2 - ( NR-1 )*KA1
1185:             IF( UPDATE ) THEN
1186:                J2T = MIN( J2, I-2*KA+K-1 )
1187:             ELSE
1188:                J2T = J2
1189:             END IF
1190:             NRT = ( J2T+KA-1 ) / KA1
1191:             DO 800 J = J1, J2T, KA1
1192: *
1193: *              create nonzero element a(j+ka,j-1) outside the band
1194: *              and store it in WORK(j)
1195: *
1196:                WORK( J ) = WORK( J )*AB( KA1, J-1 )
1197:                AB( KA1, J-1 ) = RWORK( J )*AB( KA1, J-1 )
1198:   800       CONTINUE
1199: *
1200: *           generate rotations in 1st set to annihilate elements which
1201: *           have been created outside the band
1202: *
1203:             IF( NRT.GT.0 )
1204:      $         CALL CLARGV( NRT, AB( KA1, J1 ), INCA, WORK( J1 ), KA1,
1205:      $                      RWORK( J1 ), KA1 )
1206:             IF( NR.GT.0 ) THEN
1207: *
1208: *              apply rotations in 1st set from the right
1209: *
1210:                DO 810 L = 1, KA - 1
1211:                   CALL CLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
1212:      $                         INCA, RWORK( J1 ), WORK( J1 ), KA1 )
1213:   810          CONTINUE
1214: *
1215: *              apply rotations in 1st set from both sides to diagonal
1216: *              blocks
1217: *
1218:                CALL CLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
1219:      $                      AB( 2, J1-1 ), INCA, RWORK( J1 ),
1220:      $                      WORK( J1 ), KA1 )
1221: *
1222:                CALL CLACGV( NR, WORK( J1 ), KA1 )
1223:             END IF
1224: *
1225: *           start applying rotations in 1st set from the left
1226: *
1227:             DO 820 L = KA - 1, KB - K + 1, -1
1228:                NRT = ( J2+L-1 ) / KA1
1229:                J1T = J2 - ( NRT-1 )*KA1
1230:                IF( NRT.GT.0 )
1231:      $            CALL CLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
1232:      $                         AB( KA1-L, J1T-KA1+L ), INCA,
1233:      $                         RWORK( J1T ), WORK( J1T ), KA1 )
1234:   820       CONTINUE
1235: *
1236:             IF( WANTX ) THEN
1237: *
1238: *              post-multiply X by product of rotations in 1st set
1239: *
1240:                DO 830 J = J1, J2, KA1
1241:                   CALL CROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1242:      $                       RWORK( J ), CONJG( WORK( J ) ) )
1243:   830          CONTINUE
1244:             END IF
1245:   840    CONTINUE
1246: *
1247:          IF( UPDATE ) THEN
1248:             IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
1249: *
1250: *              create nonzero element a(i+kbt,i+kbt-ka-1) outside the
1251: *              band and store it in WORK(m-kb+i+kbt)
1252: *
1253:                WORK( M-KB+I+KBT ) = -BB( KBT+1, I )*RA1
1254:             END IF
1255:          END IF
1256: *
1257:          DO 880 K = KB, 1, -1
1258:             IF( UPDATE ) THEN
1259:                J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1
1260:             ELSE
1261:                J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1262:             END IF
1263: *
1264: *           finish applying rotations in 2nd set from the left
1265: *
1266:             DO 850 L = KB - K, 1, -1
1267:                NRT = ( J2+KA+L-1 ) / KA1
1268:                J1T = J2 - ( NRT-1 )*KA1
1269:                IF( NRT.GT.0 )
1270:      $            CALL CLARTV( NRT, AB( KA1-L+1, J1T+L-1 ), INCA,
1271:      $                         AB( KA1-L, J1T+L-1 ), INCA,
1272:      $                         RWORK( M-KB+J1T+KA ),
1273:      $                         WORK( M-KB+J1T+KA ), KA1 )
1274:   850       CONTINUE
1275:             NR = ( J2+KA-1 ) / KA1
1276:             J1 = J2 - ( NR-1 )*KA1
1277:             DO 860 J = J1, J2, KA1
1278:                WORK( M-KB+J ) = WORK( M-KB+J+KA )
1279:                RWORK( M-KB+J ) = RWORK( M-KB+J+KA )
1280:   860       CONTINUE
1281:             DO 870 J = J1, J2, KA1
1282: *
1283: *              create nonzero element a(j+ka,j-1) outside the band
1284: *              and store it in WORK(m-kb+j)
1285: *
1286:                WORK( M-KB+J ) = WORK( M-KB+J )*AB( KA1, J-1 )
1287:                AB( KA1, J-1 ) = RWORK( M-KB+J )*AB( KA1, J-1 )
1288:   870       CONTINUE
1289:             IF( UPDATE ) THEN
1290:                IF( I+K.GT.KA1 .AND. K.LE.KBT )
1291:      $            WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
1292:             END IF
1293:   880    CONTINUE
1294: *
1295:          DO 920 K = KB, 1, -1
1296:             J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1297:             NR = ( J2+KA-1 ) / KA1
1298:             J1 = J2 - ( NR-1 )*KA1
1299:             IF( NR.GT.0 ) THEN
1300: *
1301: *              generate rotations in 2nd set to annihilate elements
1302: *              which have been created outside the band
1303: *
1304:                CALL CLARGV( NR, AB( KA1, J1 ), INCA, WORK( M-KB+J1 ),
1305:      $                      KA1, RWORK( M-KB+J1 ), KA1 )
1306: *
1307: *              apply rotations in 2nd set from the right
1308: *
1309:                DO 890 L = 1, KA - 1
1310:                   CALL CLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
1311:      $                         INCA, RWORK( M-KB+J1 ), WORK( M-KB+J1 ),
1312:      $                         KA1 )
1313:   890          CONTINUE
1314: *
1315: *              apply rotations in 2nd set from both sides to diagonal
1316: *              blocks
1317: *
1318:                CALL CLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
1319:      $                      AB( 2, J1-1 ), INCA, RWORK( M-KB+J1 ),
1320:      $                      WORK( M-KB+J1 ), KA1 )
1321: *
1322:                CALL CLACGV( NR, WORK( M-KB+J1 ), KA1 )
1323:             END IF
1324: *
1325: *           start applying rotations in 2nd set from the left
1326: *
1327:             DO 900 L = KA - 1, KB - K + 1, -1
1328:                NRT = ( J2+L-1 ) / KA1
1329:                J1T = J2 - ( NRT-1 )*KA1
1330:                IF( NRT.GT.0 )
1331:      $            CALL CLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
1332:      $                         AB( KA1-L, J1T-KA1+L ), INCA,
1333:      $                         RWORK( M-KB+J1T ), WORK( M-KB+J1T ),
1334:      $                         KA1 )
1335:   900       CONTINUE
1336: *
1337:             IF( WANTX ) THEN
1338: *
1339: *              post-multiply X by product of rotations in 2nd set
1340: *
1341:                DO 910 J = J1, J2, KA1
1342:                   CALL CROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1343:      $                       RWORK( M-KB+J ), CONJG( WORK( M-KB+J ) ) )
1344:   910          CONTINUE
1345:             END IF
1346:   920    CONTINUE
1347: *
1348:          DO 940 K = 1, KB - 1
1349:             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
1350: *
1351: *           finish applying rotations in 1st set from the left
1352: *
1353:             DO 930 L = KB - K, 1, -1
1354:                NRT = ( J2+L-1 ) / KA1
1355:                J1T = J2 - ( NRT-1 )*KA1
1356:                IF( NRT.GT.0 )
1357:      $            CALL CLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
1358:      $                         AB( KA1-L, J1T-KA1+L ), INCA,
1359:      $                         RWORK( J1T ), WORK( J1T ), KA1 )
1360:   930       CONTINUE
1361:   940    CONTINUE
1362: *
1363:          IF( KB.GT.1 ) THEN
1364:             DO 950 J = 2, I2 - KA
1365:                RWORK( J ) = RWORK( J+KA )
1366:                WORK( J ) = WORK( J+KA )
1367:   950       CONTINUE
1368:          END IF
1369: *
1370:       END IF
1371: *
1372:       GO TO 490
1373: *
1374: *     End of CHBGST
1375: *
1376:       END
1377: