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