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