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