SRC/ssbgst.f

   1 *> \brief \b SSBGST
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SSBGST + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssbgst.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssbgst.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssbgst.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SSBGST( VECT, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, X,
  22 *                          LDX, WORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          UPLO, VECT
  26 *       INTEGER            INFO, KA, KB, LDAB, LDBB, LDX, N
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       REAL               AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
  30 *      $                   X( LDX, * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> SSBGST reduces a real symmetric-definite banded generalized
  40 *> eigenproblem  A*x = lambda*B*x  to standard form  C*y = lambda*y,
  41 *> such that C has the same bandwidth as A.
  42 *>
  43 *> B must have been previously factorized as S**T*S by SPBSTF, using a
  44 *> split Cholesky factorization. A is overwritten by C = X**T*A*X, where
  45 *> X = S**(-1)*Q and Q is an orthogonal matrix chosen to preserve the
  46 *> bandwidth of A.
  47 *> \endverbatim
  48 *
  49 *  Arguments:
  50 *  ==========
  51 *
  52 *> \param[in] VECT
  53 *> \verbatim
  54 *>          VECT is CHARACTER*1
  55 *>          = 'N':  do not form the transformation matrix X;
  56 *>          = 'V':  form X.
  57 *> \endverbatim
  58 *>
  59 *> \param[in] UPLO
  60 *> \verbatim
  61 *>          UPLO is CHARACTER*1
  62 *>          = 'U':  Upper triangle of A is stored;
  63 *>          = 'L':  Lower triangle of A is stored.
  64 *> \endverbatim
  65 *>
  66 *> \param[in] N
  67 *> \verbatim
  68 *>          N is INTEGER
  69 *>          The order of the matrices A and B.  N >= 0.
  70 *> \endverbatim
  71 *>
  72 *> \param[in] KA
  73 *> \verbatim
  74 *>          KA is INTEGER
  75 *>          The number of superdiagonals of the matrix A if UPLO = 'U',
  76 *>          or the number of subdiagonals if UPLO = 'L'.  KA >= 0.
  77 *> \endverbatim
  78 *>
  79 *> \param[in] KB
  80 *> \verbatim
  81 *>          KB is INTEGER
  82 *>          The number of superdiagonals of the matrix B if UPLO = 'U',
  83 *>          or the number of subdiagonals if UPLO = 'L'.  KA >= KB >= 0.
  84 *> \endverbatim
  85 *>
  86 *> \param[in,out] AB
  87 *> \verbatim
  88 *>          AB is REAL array, dimension (LDAB,N)
  89 *>          On entry, the upper or lower triangle of the symmetric band
  90 *>          matrix A, stored in the first ka+1 rows of the array.  The
  91 *>          j-th column of A is stored in the j-th column of the array AB
  92 *>          as follows:
  93 *>          if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;
  94 *>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+ka).
  95 *>
  96 *>          On exit, the transformed matrix X**T*A*X, stored in the same
  97 *>          format as A.
  98 *> \endverbatim
  99 *>
 100 *> \param[in] LDAB
 101 *> \verbatim
 102 *>          LDAB is INTEGER
 103 *>          The leading dimension of the array AB.  LDAB >= KA+1.
 104 *> \endverbatim
 105 *>
 106 *> \param[in] BB
 107 *> \verbatim
 108 *>          BB is REAL array, dimension (LDBB,N)
 109 *>          The banded factor S from the split Cholesky factorization of
 110 *>          B, as returned by SPBSTF, stored in the first KB+1 rows of
 111 *>          the array.
 112 *> \endverbatim
 113 *>
 114 *> \param[in] LDBB
 115 *> \verbatim
 116 *>          LDBB is INTEGER
 117 *>          The leading dimension of the array BB.  LDBB >= KB+1.
 118 *> \endverbatim
 119 *>
 120 *> \param[out] X
 121 *> \verbatim
 122 *>          X is REAL array, dimension (LDX,N)
 123 *>          If VECT = 'V', the n-by-n matrix X.
 124 *>          If VECT = 'N', the array X is not referenced.
 125 *> \endverbatim
 126 *>
 127 *> \param[in] LDX
 128 *> \verbatim
 129 *>          LDX is INTEGER
 130 *>          The leading dimension of the array X.
 131 *>          LDX >= max(1,N) if VECT = 'V'; LDX >= 1 otherwise.
 132 *> \endverbatim
 133 *>
 134 *> \param[out] WORK
 135 *> \verbatim
 136 *>          WORK is REAL array, dimension (2*N)
 137 *> \endverbatim
 138 *>
 139 *> \param[out] INFO
 140 *> \verbatim
 141 *>          INFO is INTEGER
 142 *>          = 0:  successful exit
 143 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 144 *> \endverbatim
 145 *
 146 *  Authors:
 147 *  ========
 148 *
 149 *> \author Univ. of Tennessee
 150 *> \author Univ. of California Berkeley
 151 *> \author Univ. of Colorado Denver
 152 *> \author NAG Ltd.
 153 *
 154 *> \date November 2011
 155 *
 156 *> \ingroup realOTHERcomputational
 157 *
 158 *  =====================================================================
 159       SUBROUTINE SSBGST( VECT, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, X,
 160      $                   LDX, WORK, INFO )
 161 *
 162 *  -- LAPACK computational routine (version 3.4.0) --
 163 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 164 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 165 *     November 2011
 166 *
 167 *     .. Scalar Arguments ..
 168       CHARACTER          UPLO, VECT
 169       INTEGER            INFO, KA, KB, LDAB, LDBB, LDX, N
 170 *     ..
 171 *     .. Array Arguments ..
 172       REAL               AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
 173      $                   X( LDX, * )
 174 *     ..
 175 *
 176 *  =====================================================================
 177 *
 178 *     .. Parameters ..
 179       REAL               ZERO, ONE
 180       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 181 *     ..
 182 *     .. Local Scalars ..
 183       LOGICAL            UPDATE, UPPER, WANTX
 184       INTEGER            I, I0, I1, I2, INCA, J, J1, J1T, J2, J2T, K,
 185      $                   KA1, KB1, KBT, L, M, NR, NRT, NX
 186       REAL               BII, RA, RA1, T
 187 *     ..
 188 *     .. External Functions ..
 189       LOGICAL            LSAME
 190       EXTERNAL           LSAME
 191 *     ..
 192 *     .. External Subroutines ..
 193       EXTERNAL           SGER, SLAR2V, SLARGV, SLARTG, SLARTV, SLASET,
 194      $                   SROT, SSCAL, XERBLA
 195 *     ..
 196 *     .. Intrinsic Functions ..
 197       INTRINSIC          MAX, MIN
 198 *     ..
 199 *     .. Executable Statements ..
 200 *
 201 *     Test the input parameters
 202 *
 203       WANTX = LSAME( VECT, 'V' )
 204       UPPER = LSAME( UPLO, 'U' )
 205       KA1 = KA + 1
 206       KB1 = KB + 1
 207       INFO = 0
 208       IF( .NOT.WANTX .AND. .NOT.LSAME( VECT, 'N' ) ) THEN
 209          INFO = -1
 210       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 211          INFO = -2
 212       ELSE IF( N.LT.0 ) THEN
 213          INFO = -3
 214       ELSE IF( KA.LT.0 ) THEN
 215          INFO = -4
 216       ELSE IF( KB.LT.0 .OR. KB.GT.KA ) THEN
 217          INFO = -5
 218       ELSE IF( LDAB.LT.KA+1 ) THEN
 219          INFO = -7
 220       ELSE IF( LDBB.LT.KB+1 ) THEN
 221          INFO = -9
 222       ELSE IF( LDX.LT.1 .OR. WANTX .AND. LDX.LT.MAX( 1, N ) ) THEN
 223          INFO = -11
 224       END IF
 225       IF( INFO.NE.0 ) THEN
 226          CALL XERBLA( 'SSBGST', -INFO )
 227          RETURN
 228       END IF
 229 *
 230 *     Quick return if possible
 231 *
 232       IF( N.EQ.0 )
 233      $   RETURN
 234 *
 235       INCA = LDAB*KA1
 236 *
 237 *     Initialize X to the unit matrix, if needed
 238 *
 239       IF( WANTX )
 240      $   CALL SLASET( 'Full', N, N, ZERO, ONE, X, LDX )
 241 *
 242 *     Set M to the splitting point m. It must be the same value as is
 243 *     used in SPBSTF. The chosen value allows the arrays WORK and RWORK
 244 *     to be of dimension (N).
 245 *
 246       M = ( N+KB ) / 2
 247 *
 248 *     The routine works in two phases, corresponding to the two halves
 249 *     of the split Cholesky factorization of B as S**T*S where
 250 *
 251 *     S = ( U    )
 252 *         ( M  L )
 253 *
 254 *     with U upper triangular of order m, and L lower triangular of
 255 *     order n-m. S has the same bandwidth as B.
 256 *
 257 *     S is treated as a product of elementary matrices:
 258 *
 259 *     S = S(m)*S(m-1)*...*S(2)*S(1)*S(m+1)*S(m+2)*...*S(n-1)*S(n)
 260 *
 261 *     where S(i) is determined by the i-th row of S.
 262 *
 263 *     In phase 1, the index i takes the values n, n-1, ... , m+1;
 264 *     in phase 2, it takes the values 1, 2, ... , m.
 265 *
 266 *     For each value of i, the current matrix A is updated by forming
 267 *     inv(S(i))**T*A*inv(S(i)). This creates a triangular bulge outside
 268 *     the band of A. The bulge is then pushed down toward the bottom of
 269 *     A in phase 1, and up toward the top of A in phase 2, by applying
 270 *     plane rotations.
 271 *
 272 *     There are kb*(kb+1)/2 elements in the bulge, but at most 2*kb-1
 273 *     of them are linearly independent, so annihilating a bulge requires
 274 *     only 2*kb-1 plane rotations. The rotations are divided into a 1st
 275 *     set of kb-1 rotations, and a 2nd set of kb rotations.
 276 *
 277 *     Wherever possible, rotations are generated and applied in vector
 278 *     operations of length NR between the indices J1 and J2 (sometimes
 279 *     replaced by modified values NRT, J1T or J2T).
 280 *
 281 *     The cosines and sines of the rotations are stored in the array
 282 *     WORK. The cosines of the 1st set of rotations are stored in
 283 *     elements n+2:n+m-kb-1 and the sines of the 1st set in elements
 284 *     2:m-kb-1; the cosines of the 2nd set are stored in elements
 285 *     n+m-kb+1:2*n and the sines of the second set in elements m-kb+1:n.
 286 *
 287 *     The bulges are not formed explicitly; nonzero elements outside the
 288 *     band are created only when they are required for generating new
 289 *     rotations; they are stored in the array WORK, in positions where
 290 *     they are later overwritten by the sines of the rotations which
 291 *     annihilate them.
 292 *
 293 *     **************************** Phase 1 *****************************
 294 *
 295 *     The logical structure of this phase is:
 296 *
 297 *     UPDATE = .TRUE.
 298 *     DO I = N, M + 1, -1
 299 *        use S(i) to update A and create a new bulge
 300 *        apply rotations to push all bulges KA positions downward
 301 *     END DO
 302 *     UPDATE = .FALSE.
 303 *     DO I = M + KA + 1, N - 1
 304 *        apply rotations to push all bulges KA positions downward
 305 *     END DO
 306 *
 307 *     To avoid duplicating code, the two loops are merged.
 308 *
 309       UPDATE = .TRUE.
 310       I = N + 1
 311    10 CONTINUE
 312       IF( UPDATE ) THEN
 313          I = I - 1
 314          KBT = MIN( KB, I-1 )
 315          I0 = I - 1
 316          I1 = MIN( N, I+KA )
 317          I2 = I - KBT + KA1
 318          IF( I.LT.M+1 ) THEN
 319             UPDATE = .FALSE.
 320             I = I + 1
 321             I0 = M
 322             IF( KA.EQ.0 )
 323      $         GO TO 480
 324             GO TO 10
 325          END IF
 326       ELSE
 327          I = I + KA
 328          IF( I.GT.N-1 )
 329      $      GO TO 480
 330       END IF
 331 *
 332       IF( UPPER ) THEN
 333 *
 334 *        Transform A, working with the upper triangle
 335 *
 336          IF( UPDATE ) THEN
 337 *
 338 *           Form  inv(S(i))**T * A * inv(S(i))
 339 *
 340             BII = BB( KB1, I )
 341             DO 20 J = I, I1
 342                AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII
 343    20       CONTINUE
 344             DO 30 J = MAX( 1, I-KA ), I
 345                AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII
 346    30       CONTINUE
 347             DO 60 K = I - KBT, I - 1
 348                DO 40 J = I - KBT, K
 349                   AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
 350      $                               BB( J-I+KB1, I )*AB( K-I+KA1, I ) -
 351      $                               BB( K-I+KB1, I )*AB( J-I+KA1, I ) +
 352      $                               AB( KA1, I )*BB( J-I+KB1, I )*
 353      $                               BB( K-I+KB1, I )
 354    40          CONTINUE
 355                DO 50 J = MAX( 1, I-KA ), I - KBT - 1
 356                   AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
 357      $                               BB( K-I+KB1, I )*AB( J-I+KA1, I )
 358    50          CONTINUE
 359    60       CONTINUE
 360             DO 80 J = I, I1
 361                DO 70 K = MAX( J-KA, I-KBT ), I - 1
 362                   AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
 363      $                               BB( K-I+KB1, I )*AB( I-J+KA1, J )
 364    70          CONTINUE
 365    80       CONTINUE
 366 *
 367             IF( WANTX ) THEN
 368 *
 369 *              post-multiply X by inv(S(i))
 370 *
 371                CALL SSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
 372                IF( KBT.GT.0 )
 373      $            CALL SGER( N-M, KBT, -ONE, X( M+1, I ), 1,
 374      $                       BB( KB1-KBT, I ), 1, X( M+1, I-KBT ), LDX )
 375             END IF
 376 *
 377 *           store a(i,i1) in RA1 for use in next loop over K
 378 *
 379             RA1 = AB( I-I1+KA1, I1 )
 380          END IF
 381 *
 382 *        Generate and apply vectors of rotations to chase all the
 383 *        existing bulges KA positions down toward the bottom of the
 384 *        band
 385 *
 386          DO 130 K = 1, KB - 1
 387             IF( UPDATE ) THEN
 388 *
 389 *              Determine the rotations which would annihilate the bulge
 390 *              which has in theory just been created
 391 *
 392                IF( I-K+KA.LT.N .AND. I-K.GT.1 ) THEN
 393 *
 394 *                 generate rotation to annihilate a(i,i-k+ka+1)
 395 *
 396                   CALL SLARTG( AB( K+1, I-K+KA ), RA1,
 397      $                         WORK( N+I-K+KA-M ), WORK( I-K+KA-M ),
 398      $                         RA )
 399 *
 400 *                 create nonzero element a(i-k,i-k+ka+1) outside the
 401 *                 band and store it in WORK(i-k)
 402 *
 403                   T = -BB( KB1-K, I )*RA1
 404                   WORK( I-K ) = WORK( N+I-K+KA-M )*T -
 405      $                          WORK( I-K+KA-M )*AB( 1, I-K+KA )
 406                   AB( 1, I-K+KA ) = WORK( I-K+KA-M )*T +
 407      $                              WORK( N+I-K+KA-M )*AB( 1, I-K+KA )
 408                   RA1 = RA
 409                END IF
 410             END IF
 411             J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
 412             NR = ( N-J2+KA ) / KA1
 413             J1 = J2 + ( NR-1 )*KA1
 414             IF( UPDATE ) THEN
 415                J2T = MAX( J2, I+2*KA-K+1 )
 416             ELSE
 417                J2T = J2
 418             END IF
 419             NRT = ( N-J2T+KA ) / KA1
 420             DO 90 J = J2T, J1, KA1
 421 *
 422 *              create nonzero element a(j-ka,j+1) outside the band
 423 *              and store it in WORK(j-m)
 424 *
 425                WORK( J-M ) = WORK( J-M )*AB( 1, J+1 )
 426                AB( 1, J+1 ) = WORK( N+J-M )*AB( 1, J+1 )
 427    90       CONTINUE
 428 *
 429 *           generate rotations in 1st set to annihilate elements which
 430 *           have been created outside the band
 431 *
 432             IF( NRT.GT.0 )
 433      $         CALL SLARGV( NRT, AB( 1, J2T ), INCA, WORK( J2T-M ), KA1,
 434      $                      WORK( N+J2T-M ), KA1 )
 435             IF( NR.GT.0 ) THEN
 436 *
 437 *              apply rotations in 1st set from the right
 438 *
 439                DO 100 L = 1, KA - 1
 440                   CALL SLARTV( NR, AB( KA1-L, J2 ), INCA,
 441      $                         AB( KA-L, J2+1 ), INCA, WORK( N+J2-M ),
 442      $                         WORK( J2-M ), KA1 )
 443   100          CONTINUE
 444 *
 445 *              apply rotations in 1st set from both sides to diagonal
 446 *              blocks
 447 *
 448                CALL SLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ),
 449      $                      AB( KA, J2+1 ), INCA, WORK( N+J2-M ),
 450      $                      WORK( J2-M ), KA1 )
 451 *
 452             END IF
 453 *
 454 *           start applying rotations in 1st set from the left
 455 *
 456             DO 110 L = KA - 1, KB - K + 1, -1
 457                NRT = ( N-J2+L ) / KA1
 458                IF( NRT.GT.0 )
 459      $            CALL SLARTV( NRT, AB( L, J2+KA1-L ), INCA,
 460      $                         AB( L+1, J2+KA1-L ), INCA,
 461      $                         WORK( N+J2-M ), WORK( J2-M ), KA1 )
 462   110       CONTINUE
 463 *
 464             IF( WANTX ) THEN
 465 *
 466 *              post-multiply X by product of rotations in 1st set
 467 *
 468                DO 120 J = J2, J1, KA1
 469                   CALL SROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
 470      $                       WORK( N+J-M ), WORK( J-M ) )
 471   120          CONTINUE
 472             END IF
 473   130    CONTINUE
 474 *
 475          IF( UPDATE ) THEN
 476             IF( I2.LE.N .AND. KBT.GT.0 ) THEN
 477 *
 478 *              create nonzero element a(i-kbt,i-kbt+ka+1) outside the
 479 *              band and store it in WORK(i-kbt)
 480 *
 481                WORK( I-KBT ) = -BB( KB1-KBT, I )*RA1
 482             END IF
 483          END IF
 484 *
 485          DO 170 K = KB, 1, -1
 486             IF( UPDATE ) THEN
 487                J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1
 488             ELSE
 489                J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
 490             END IF
 491 *
 492 *           finish applying rotations in 2nd set from the left
 493 *
 494             DO 140 L = KB - K, 1, -1
 495                NRT = ( N-J2+KA+L ) / KA1
 496                IF( NRT.GT.0 )
 497      $            CALL SLARTV( NRT, AB( L, J2-L+1 ), INCA,
 498      $                         AB( L+1, J2-L+1 ), INCA, WORK( N+J2-KA ),
 499      $                         WORK( J2-KA ), KA1 )
 500   140       CONTINUE
 501             NR = ( N-J2+KA ) / KA1
 502             J1 = J2 + ( NR-1 )*KA1
 503             DO 150 J = J1, J2, -KA1
 504                WORK( J ) = WORK( J-KA )
 505                WORK( N+J ) = WORK( N+J-KA )
 506   150       CONTINUE
 507             DO 160 J = J2, J1, KA1
 508 *
 509 *              create nonzero element a(j-ka,j+1) outside the band
 510 *              and store it in WORK(j)
 511 *
 512                WORK( J ) = WORK( J )*AB( 1, J+1 )
 513                AB( 1, J+1 ) = WORK( N+J )*AB( 1, J+1 )
 514   160       CONTINUE
 515             IF( UPDATE ) THEN
 516                IF( I-K.LT.N-KA .AND. K.LE.KBT )
 517      $            WORK( I-K+KA ) = WORK( I-K )
 518             END IF
 519   170    CONTINUE
 520 *
 521          DO 210 K = KB, 1, -1
 522             J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
 523             NR = ( N-J2+KA ) / KA1
 524             J1 = J2 + ( NR-1 )*KA1
 525             IF( NR.GT.0 ) THEN
 526 *
 527 *              generate rotations in 2nd set to annihilate elements
 528 *              which have been created outside the band
 529 *
 530                CALL SLARGV( NR, AB( 1, J2 ), INCA, WORK( J2 ), KA1,
 531      $                      WORK( N+J2 ), KA1 )
 532 *
 533 *              apply rotations in 2nd set from the right
 534 *
 535                DO 180 L = 1, KA - 1
 536                   CALL SLARTV( NR, AB( KA1-L, J2 ), INCA,
 537      $                         AB( KA-L, J2+1 ), INCA, WORK( N+J2 ),
 538      $                         WORK( J2 ), KA1 )
 539   180          CONTINUE
 540 *
 541 *              apply rotations in 2nd set from both sides to diagonal
 542 *              blocks
 543 *
 544                CALL SLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ),
 545      $                      AB( KA, J2+1 ), INCA, WORK( N+J2 ),
 546      $                      WORK( J2 ), KA1 )
 547 *
 548             END IF
 549 *
 550 *           start applying rotations in 2nd set from the left
 551 *
 552             DO 190 L = KA - 1, KB - K + 1, -1
 553                NRT = ( N-J2+L ) / KA1
 554                IF( NRT.GT.0 )
 555      $            CALL SLARTV( NRT, AB( L, J2+KA1-L ), INCA,
 556      $                         AB( L+1, J2+KA1-L ), INCA, WORK( N+J2 ),
 557      $                         WORK( J2 ), KA1 )
 558   190       CONTINUE
 559 *
 560             IF( WANTX ) THEN
 561 *
 562 *              post-multiply X by product of rotations in 2nd set
 563 *
 564                DO 200 J = J2, J1, KA1
 565                   CALL SROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
 566      $                       WORK( N+J ), WORK( J ) )
 567   200          CONTINUE
 568             END IF
 569   210    CONTINUE
 570 *
 571          DO 230 K = 1, KB - 1
 572             J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
 573 *
 574 *           finish applying rotations in 1st set from the left
 575 *
 576             DO 220 L = KB - K, 1, -1
 577                NRT = ( N-J2+L ) / KA1
 578                IF( NRT.GT.0 )
 579      $            CALL SLARTV( NRT, AB( L, J2+KA1-L ), INCA,
 580      $                         AB( L+1, J2+KA1-L ), INCA,
 581      $                         WORK( N+J2-M ), WORK( J2-M ), KA1 )
 582   220       CONTINUE
 583   230    CONTINUE
 584 *
 585          IF( KB.GT.1 ) THEN
 586             DO 240 J = N - 1, I - KB + 2*KA + 1, -1
 587                WORK( N+J-M ) = WORK( N+J-KA-M )
 588                WORK( J-M ) = WORK( J-KA-M )
 589   240       CONTINUE
 590          END IF
 591 *
 592       ELSE
 593 *
 594 *        Transform A, working with the lower triangle
 595 *
 596          IF( UPDATE ) THEN
 597 *
 598 *           Form  inv(S(i))**T * A * inv(S(i))
 599 *
 600             BII = BB( 1, I )
 601             DO 250 J = I, I1
 602                AB( J-I+1, I ) = AB( J-I+1, I ) / BII
 603   250       CONTINUE
 604             DO 260 J = MAX( 1, I-KA ), I
 605                AB( I-J+1, J ) = AB( I-J+1, J ) / BII
 606   260       CONTINUE
 607             DO 290 K = I - KBT, I - 1
 608                DO 270 J = I - KBT, K
 609                   AB( K-J+1, J ) = AB( K-J+1, J ) -
 610      $                             BB( I-J+1, J )*AB( I-K+1, K ) -
 611      $                             BB( I-K+1, K )*AB( I-J+1, J ) +
 612      $                             AB( 1, I )*BB( I-J+1, J )*
 613      $                             BB( I-K+1, K )
 614   270          CONTINUE
 615                DO 280 J = MAX( 1, I-KA ), I - KBT - 1
 616                   AB( K-J+1, J ) = AB( K-J+1, J ) -
 617      $                             BB( I-K+1, K )*AB( I-J+1, J )
 618   280          CONTINUE
 619   290       CONTINUE
 620             DO 310 J = I, I1
 621                DO 300 K = MAX( J-KA, I-KBT ), I - 1
 622                   AB( J-K+1, K ) = AB( J-K+1, K ) -
 623      $                             BB( I-K+1, K )*AB( J-I+1, I )
 624   300          CONTINUE
 625   310       CONTINUE
 626 *
 627             IF( WANTX ) THEN
 628 *
 629 *              post-multiply X by inv(S(i))
 630 *
 631                CALL SSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
 632                IF( KBT.GT.0 )
 633      $            CALL SGER( N-M, KBT, -ONE, X( M+1, I ), 1,
 634      $                       BB( KBT+1, I-KBT ), LDBB-1,
 635      $                       X( M+1, I-KBT ), LDX )
 636             END IF
 637 *
 638 *           store a(i1,i) in RA1 for use in next loop over K
 639 *
 640             RA1 = AB( I1-I+1, I )
 641          END IF
 642 *
 643 *        Generate and apply vectors of rotations to chase all the
 644 *        existing bulges KA positions down toward the bottom of the
 645 *        band
 646 *
 647          DO 360 K = 1, KB - 1
 648             IF( UPDATE ) THEN
 649 *
 650 *              Determine the rotations which would annihilate the bulge
 651 *              which has in theory just been created
 652 *
 653                IF( I-K+KA.LT.N .AND. I-K.GT.1 ) THEN
 654 *
 655 *                 generate rotation to annihilate a(i-k+ka+1,i)
 656 *
 657                   CALL SLARTG( AB( KA1-K, I ), RA1, WORK( N+I-K+KA-M ),
 658      $                         WORK( I-K+KA-M ), RA )
 659 *
 660 *                 create nonzero element a(i-k+ka+1,i-k) outside the
 661 *                 band and store it in WORK(i-k)
 662 *
 663                   T = -BB( K+1, I-K )*RA1
 664                   WORK( I-K ) = WORK( N+I-K+KA-M )*T -
 665      $                          WORK( I-K+KA-M )*AB( KA1, I-K )
 666                   AB( KA1, I-K ) = WORK( I-K+KA-M )*T +
 667      $                             WORK( N+I-K+KA-M )*AB( KA1, I-K )
 668                   RA1 = RA
 669                END IF
 670             END IF
 671             J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
 672             NR = ( N-J2+KA ) / KA1
 673             J1 = J2 + ( NR-1 )*KA1
 674             IF( UPDATE ) THEN
 675                J2T = MAX( J2, I+2*KA-K+1 )
 676             ELSE
 677                J2T = J2
 678             END IF
 679             NRT = ( N-J2T+KA ) / KA1
 680             DO 320 J = J2T, J1, KA1
 681 *
 682 *              create nonzero element a(j+1,j-ka) outside the band
 683 *              and store it in WORK(j-m)
 684 *
 685                WORK( J-M ) = WORK( J-M )*AB( KA1, J-KA+1 )
 686                AB( KA1, J-KA+1 ) = WORK( N+J-M )*AB( KA1, J-KA+1 )
 687   320       CONTINUE
 688 *
 689 *           generate rotations in 1st set to annihilate elements which
 690 *           have been created outside the band
 691 *
 692             IF( NRT.GT.0 )
 693      $         CALL SLARGV( NRT, AB( KA1, J2T-KA ), INCA, WORK( J2T-M ),
 694      $                      KA1, WORK( N+J2T-M ), KA1 )
 695             IF( NR.GT.0 ) THEN
 696 *
 697 *              apply rotations in 1st set from the left
 698 *
 699                DO 330 L = 1, KA - 1
 700                   CALL SLARTV( NR, AB( L+1, J2-L ), INCA,
 701      $                         AB( L+2, J2-L ), INCA, WORK( N+J2-M ),
 702      $                         WORK( J2-M ), KA1 )
 703   330          CONTINUE
 704 *
 705 *              apply rotations in 1st set from both sides to diagonal
 706 *              blocks
 707 *
 708                CALL SLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
 709      $                      INCA, WORK( N+J2-M ), WORK( J2-M ), KA1 )
 710 *
 711             END IF
 712 *
 713 *           start applying rotations in 1st set from the right
 714 *
 715             DO 340 L = KA - 1, KB - K + 1, -1
 716                NRT = ( N-J2+L ) / KA1
 717                IF( NRT.GT.0 )
 718      $            CALL SLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
 719      $                         AB( KA1-L, J2+1 ), INCA, WORK( N+J2-M ),
 720      $                         WORK( J2-M ), KA1 )
 721   340       CONTINUE
 722 *
 723             IF( WANTX ) THEN
 724 *
 725 *              post-multiply X by product of rotations in 1st set
 726 *
 727                DO 350 J = J2, J1, KA1
 728                   CALL SROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
 729      $                       WORK( N+J-M ), WORK( J-M ) )
 730   350          CONTINUE
 731             END IF
 732   360    CONTINUE
 733 *
 734          IF( UPDATE ) THEN
 735             IF( I2.LE.N .AND. KBT.GT.0 ) THEN
 736 *
 737 *              create nonzero element a(i-kbt+ka+1,i-kbt) outside the
 738 *              band and store it in WORK(i-kbt)
 739 *
 740                WORK( I-KBT ) = -BB( KBT+1, I-KBT )*RA1
 741             END IF
 742          END IF
 743 *
 744          DO 400 K = KB, 1, -1
 745             IF( UPDATE ) THEN
 746                J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1
 747             ELSE
 748                J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
 749             END IF
 750 *
 751 *           finish applying rotations in 2nd set from the right
 752 *
 753             DO 370 L = KB - K, 1, -1
 754                NRT = ( N-J2+KA+L ) / KA1
 755                IF( NRT.GT.0 )
 756      $            CALL SLARTV( NRT, AB( KA1-L+1, J2-KA ), INCA,
 757      $                         AB( KA1-L, J2-KA+1 ), INCA,
 758      $                         WORK( N+J2-KA ), WORK( J2-KA ), KA1 )
 759   370       CONTINUE
 760             NR = ( N-J2+KA ) / KA1
 761             J1 = J2 + ( NR-1 )*KA1
 762             DO 380 J = J1, J2, -KA1
 763                WORK( J ) = WORK( J-KA )
 764                WORK( N+J ) = WORK( N+J-KA )
 765   380       CONTINUE
 766             DO 390 J = J2, J1, KA1
 767 *
 768 *              create nonzero element a(j+1,j-ka) outside the band
 769 *              and store it in WORK(j)
 770 *
 771                WORK( J ) = WORK( J )*AB( KA1, J-KA+1 )
 772                AB( KA1, J-KA+1 ) = WORK( N+J )*AB( KA1, J-KA+1 )
 773   390       CONTINUE
 774             IF( UPDATE ) THEN
 775                IF( I-K.LT.N-KA .AND. K.LE.KBT )
 776      $            WORK( I-K+KA ) = WORK( I-K )
 777             END IF
 778   400    CONTINUE
 779 *
 780          DO 440 K = KB, 1, -1
 781             J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
 782             NR = ( N-J2+KA ) / KA1
 783             J1 = J2 + ( NR-1 )*KA1
 784             IF( NR.GT.0 ) THEN
 785 *
 786 *              generate rotations in 2nd set to annihilate elements
 787 *              which have been created outside the band
 788 *
 789                CALL SLARGV( NR, AB( KA1, J2-KA ), INCA, WORK( J2 ), KA1,
 790      $                      WORK( N+J2 ), KA1 )
 791 *
 792 *              apply rotations in 2nd set from the left
 793 *
 794                DO 410 L = 1, KA - 1
 795                   CALL SLARTV( NR, AB( L+1, J2-L ), INCA,
 796      $                         AB( L+2, J2-L ), INCA, WORK( N+J2 ),
 797      $                         WORK( J2 ), KA1 )
 798   410          CONTINUE
 799 *
 800 *              apply rotations in 2nd set from both sides to diagonal
 801 *              blocks
 802 *
 803                CALL SLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
 804      $                      INCA, WORK( N+J2 ), WORK( J2 ), KA1 )
 805 *
 806             END IF
 807 *
 808 *           start applying rotations in 2nd set from the right
 809 *
 810             DO 420 L = KA - 1, KB - K + 1, -1
 811                NRT = ( N-J2+L ) / KA1
 812                IF( NRT.GT.0 )
 813      $            CALL SLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
 814      $                         AB( KA1-L, J2+1 ), INCA, WORK( N+J2 ),
 815      $                         WORK( J2 ), KA1 )
 816   420       CONTINUE
 817 *
 818             IF( WANTX ) THEN
 819 *
 820 *              post-multiply X by product of rotations in 2nd set
 821 *
 822                DO 430 J = J2, J1, KA1
 823                   CALL SROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
 824      $                       WORK( N+J ), WORK( J ) )
 825   430          CONTINUE
 826             END IF
 827   440    CONTINUE
 828 *
 829          DO 460 K = 1, KB - 1
 830             J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
 831 *
 832 *           finish applying rotations in 1st set from the right
 833 *
 834             DO 450 L = KB - K, 1, -1
 835                NRT = ( N-J2+L ) / KA1
 836                IF( NRT.GT.0 )
 837      $            CALL SLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
 838      $                         AB( KA1-L, J2+1 ), INCA, WORK( N+J2-M ),
 839      $                         WORK( J2-M ), KA1 )
 840   450       CONTINUE
 841   460    CONTINUE
 842 *
 843          IF( KB.GT.1 ) THEN
 844             DO 470 J = N - 1, I - KB + 2*KA + 1, -1
 845                WORK( N+J-M ) = WORK( N+J-KA-M )
 846                WORK( J-M ) = WORK( J-KA-M )
 847   470       CONTINUE
 848          END IF
 849 *
 850       END IF
 851 *
 852       GO TO 10
 853 *
 854   480 CONTINUE
 855 *
 856 *     **************************** Phase 2 *****************************
 857 *
 858 *     The logical structure of this phase is:
 859 *
 860 *     UPDATE = .TRUE.
 861 *     DO I = 1, M
 862 *        use S(i) to update A and create a new bulge
 863 *        apply rotations to push all bulges KA positions upward
 864 *     END DO
 865 *     UPDATE = .FALSE.
 866 *     DO I = M - KA - 1, 2, -1
 867 *        apply rotations to push all bulges KA positions upward
 868 *     END DO
 869 *
 870 *     To avoid duplicating code, the two loops are merged.
 871 *
 872       UPDATE = .TRUE.
 873       I = 0
 874   490 CONTINUE
 875       IF( UPDATE ) THEN
 876          I = I + 1
 877          KBT = MIN( KB, M-I )
 878          I0 = I + 1
 879          I1 = MAX( 1, I-KA )
 880          I2 = I + KBT - KA1
 881          IF( I.GT.M ) THEN
 882             UPDATE = .FALSE.
 883             I = I - 1
 884             I0 = M + 1
 885             IF( KA.EQ.0 )
 886      $         RETURN
 887             GO TO 490
 888          END IF
 889       ELSE
 890          I = I - KA
 891          IF( I.LT.2 )
 892      $      RETURN
 893       END IF
 894 *
 895       IF( I.LT.M-KBT ) THEN
 896          NX = M
 897       ELSE
 898          NX = N
 899       END IF
 900 *
 901       IF( UPPER ) THEN
 902 *
 903 *        Transform A, working with the upper triangle
 904 *
 905          IF( UPDATE ) THEN
 906 *
 907 *           Form  inv(S(i))**T * A * inv(S(i))
 908 *
 909             BII = BB( KB1, I )
 910             DO 500 J = I1, I
 911                AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII
 912   500       CONTINUE
 913             DO 510 J = I, MIN( N, I+KA )
 914                AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII
 915   510       CONTINUE
 916             DO 540 K = I + 1, I + KBT
 917                DO 520 J = K, I + KBT
 918                   AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
 919      $                               BB( I-J+KB1, J )*AB( I-K+KA1, K ) -
 920      $                               BB( I-K+KB1, K )*AB( I-J+KA1, J ) +
 921      $                               AB( KA1, I )*BB( I-J+KB1, J )*
 922      $                               BB( I-K+KB1, K )
 923   520          CONTINUE
 924                DO 530 J = I + KBT + 1, MIN( N, I+KA )
 925                   AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
 926      $                               BB( I-K+KB1, K )*AB( I-J+KA1, J )
 927   530          CONTINUE
 928   540       CONTINUE
 929             DO 560 J = I1, I
 930                DO 550 K = I + 1, MIN( J+KA, I+KBT )
 931                   AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
 932      $                               BB( I-K+KB1, K )*AB( J-I+KA1, I )
 933   550          CONTINUE
 934   560       CONTINUE
 935 *
 936             IF( WANTX ) THEN
 937 *
 938 *              post-multiply X by inv(S(i))
 939 *
 940                CALL SSCAL( NX, ONE / BII, X( 1, I ), 1 )
 941                IF( KBT.GT.0 )
 942      $            CALL SGER( NX, KBT, -ONE, X( 1, I ), 1, BB( KB, I+1 ),
 943      $                       LDBB-1, X( 1, I+1 ), LDX )
 944             END IF
 945 *
 946 *           store a(i1,i) in RA1 for use in next loop over K
 947 *
 948             RA1 = AB( I1-I+KA1, I )
 949          END IF
 950 *
 951 *        Generate and apply vectors of rotations to chase all the
 952 *        existing bulges KA positions up toward the top of the band
 953 *
 954          DO 610 K = 1, KB - 1
 955             IF( UPDATE ) THEN
 956 *
 957 *              Determine the rotations which would annihilate the bulge
 958 *              which has in theory just been created
 959 *
 960                IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
 961 *
 962 *                 generate rotation to annihilate a(i+k-ka-1,i)
 963 *
 964                   CALL SLARTG( AB( K+1, I ), RA1, WORK( N+I+K-KA ),
 965      $                         WORK( I+K-KA ), RA )
 966 *
 967 *                 create nonzero element a(i+k-ka-1,i+k) outside the
 968 *                 band and store it in WORK(m-kb+i+k)
 969 *
 970                   T = -BB( KB1-K, I+K )*RA1
 971                   WORK( M-KB+I+K ) = WORK( N+I+K-KA )*T -
 972      $                               WORK( I+K-KA )*AB( 1, I+K )
 973                   AB( 1, I+K ) = WORK( I+K-KA )*T +
 974      $                           WORK( N+I+K-KA )*AB( 1, I+K )
 975                   RA1 = RA
 976                END IF
 977             END IF
 978             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
 979             NR = ( J2+KA-1 ) / KA1
 980             J1 = J2 - ( NR-1 )*KA1
 981             IF( UPDATE ) THEN
 982                J2T = MIN( J2, I-2*KA+K-1 )
 983             ELSE
 984                J2T = J2
 985             END IF
 986             NRT = ( J2T+KA-1 ) / KA1
 987             DO 570 J = J1, J2T, KA1
 988 *
 989 *              create nonzero element a(j-1,j+ka) outside the band
 990 *              and store it in WORK(j)
 991 *
 992                WORK( J ) = WORK( J )*AB( 1, J+KA-1 )
 993                AB( 1, J+KA-1 ) = WORK( N+J )*AB( 1, J+KA-1 )
 994   570       CONTINUE
 995 *
 996 *           generate rotations in 1st set to annihilate elements which
 997 *           have been created outside the band
 998 *
 999             IF( NRT.GT.0 )
1000      $         CALL SLARGV( NRT, AB( 1, J1+KA ), INCA, WORK( J1 ), KA1,
1001      $                      WORK( N+J1 ), KA1 )
1002             IF( NR.GT.0 ) THEN
1003 *
1004 *              apply rotations in 1st set from the left
1005 *
1006                DO 580 L = 1, KA - 1
1007                   CALL SLARTV( NR, AB( KA1-L, J1+L ), INCA,
1008      $                         AB( KA-L, J1+L ), INCA, WORK( N+J1 ),
1009      $                         WORK( J1 ), KA1 )
1010   580          CONTINUE
1011 *
1012 *              apply rotations in 1st set from both sides to diagonal
1013 *              blocks
1014 *
1015                CALL SLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
1016      $                      AB( KA, J1 ), INCA, WORK( N+J1 ),
1017      $                      WORK( J1 ), KA1 )
1018 *
1019             END IF
1020 *
1021 *           start applying rotations in 1st set from the right
1022 *
1023             DO 590 L = KA - 1, KB - K + 1, -1
1024                NRT = ( J2+L-1 ) / KA1
1025                J1T = J2 - ( NRT-1 )*KA1
1026                IF( NRT.GT.0 )
1027      $            CALL SLARTV( NRT, AB( L, J1T ), INCA,
1028      $                         AB( L+1, J1T-1 ), INCA, WORK( N+J1T ),
1029      $                         WORK( J1T ), KA1 )
1030   590       CONTINUE
1031 *
1032             IF( WANTX ) THEN
1033 *
1034 *              post-multiply X by product of rotations in 1st set
1035 *
1036                DO 600 J = J1, J2, KA1
1037                   CALL SROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1038      $                       WORK( N+J ), WORK( J ) )
1039   600          CONTINUE
1040             END IF
1041   610    CONTINUE
1042 *
1043          IF( UPDATE ) THEN
1044             IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
1045 *
1046 *              create nonzero element a(i+kbt-ka-1,i+kbt) outside the
1047 *              band and store it in WORK(m-kb+i+kbt)
1048 *
1049                WORK( M-KB+I+KBT ) = -BB( KB1-KBT, I+KBT )*RA1
1050             END IF
1051          END IF
1052 *
1053          DO 650 K = KB, 1, -1
1054             IF( UPDATE ) THEN
1055                J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1
1056             ELSE
1057                J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1058             END IF
1059 *
1060 *           finish applying rotations in 2nd set from the right
1061 *
1062             DO 620 L = KB - K, 1, -1
1063                NRT = ( J2+KA+L-1 ) / KA1
1064                J1T = J2 - ( NRT-1 )*KA1
1065                IF( NRT.GT.0 )
1066      $            CALL SLARTV( NRT, AB( L, J1T+KA ), INCA,
1067      $                         AB( L+1, J1T+KA-1 ), INCA,
1068      $                         WORK( N+M-KB+J1T+KA ),
1069      $                         WORK( M-KB+J1T+KA ), KA1 )
1070   620       CONTINUE
1071             NR = ( J2+KA-1 ) / KA1
1072             J1 = J2 - ( NR-1 )*KA1
1073             DO 630 J = J1, J2, KA1
1074                WORK( M-KB+J ) = WORK( M-KB+J+KA )
1075                WORK( N+M-KB+J ) = WORK( N+M-KB+J+KA )
1076   630       CONTINUE
1077             DO 640 J = J1, J2, KA1
1078 *
1079 *              create nonzero element a(j-1,j+ka) outside the band
1080 *              and store it in WORK(m-kb+j)
1081 *
1082                WORK( M-KB+J ) = WORK( M-KB+J )*AB( 1, J+KA-1 )
1083                AB( 1, J+KA-1 ) = WORK( N+M-KB+J )*AB( 1, J+KA-1 )
1084   640       CONTINUE
1085             IF( UPDATE ) THEN
1086                IF( I+K.GT.KA1 .AND. K.LE.KBT )
1087      $            WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
1088             END IF
1089   650    CONTINUE
1090 *
1091          DO 690 K = KB, 1, -1
1092             J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1093             NR = ( J2+KA-1 ) / KA1
1094             J1 = J2 - ( NR-1 )*KA1
1095             IF( NR.GT.0 ) THEN
1096 *
1097 *              generate rotations in 2nd set to annihilate elements
1098 *              which have been created outside the band
1099 *
1100                CALL SLARGV( NR, AB( 1, J1+KA ), INCA, WORK( M-KB+J1 ),
1101      $                      KA1, WORK( N+M-KB+J1 ), KA1 )
1102 *
1103 *              apply rotations in 2nd set from the left
1104 *
1105                DO 660 L = 1, KA - 1
1106                   CALL SLARTV( NR, AB( KA1-L, J1+L ), INCA,
1107      $                         AB( KA-L, J1+L ), INCA,
1108      $                         WORK( N+M-KB+J1 ), WORK( M-KB+J1 ), KA1 )
1109   660          CONTINUE
1110 *
1111 *              apply rotations in 2nd set from both sides to diagonal
1112 *              blocks
1113 *
1114                CALL SLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
1115      $                      AB( KA, J1 ), INCA, WORK( N+M-KB+J1 ),
1116      $                      WORK( M-KB+J1 ), KA1 )
1117 *
1118             END IF
1119 *
1120 *           start applying rotations in 2nd set from the right
1121 *
1122             DO 670 L = KA - 1, KB - K + 1, -1
1123                NRT = ( J2+L-1 ) / KA1
1124                J1T = J2 - ( NRT-1 )*KA1
1125                IF( NRT.GT.0 )
1126      $            CALL SLARTV( NRT, AB( L, J1T ), INCA,
1127      $                         AB( L+1, J1T-1 ), INCA,
1128      $                         WORK( N+M-KB+J1T ), WORK( M-KB+J1T ),
1129      $                         KA1 )
1130   670       CONTINUE
1131 *
1132             IF( WANTX ) THEN
1133 *
1134 *              post-multiply X by product of rotations in 2nd set
1135 *
1136                DO 680 J = J1, J2, KA1
1137                   CALL SROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1138      $                       WORK( N+M-KB+J ), WORK( M-KB+J ) )
1139   680          CONTINUE
1140             END IF
1141   690    CONTINUE
1142 *
1143          DO 710 K = 1, KB - 1
1144             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
1145 *
1146 *           finish applying rotations in 1st set from the right
1147 *
1148             DO 700 L = KB - K, 1, -1
1149                NRT = ( J2+L-1 ) / KA1
1150                J1T = J2 - ( NRT-1 )*KA1
1151                IF( NRT.GT.0 )
1152      $            CALL SLARTV( NRT, AB( L, J1T ), INCA,
1153      $                         AB( L+1, J1T-1 ), INCA, WORK( N+J1T ),
1154      $                         WORK( J1T ), KA1 )
1155   700       CONTINUE
1156   710    CONTINUE
1157 *
1158          IF( KB.GT.1 ) THEN
1159             DO 720 J = 2, MIN( I+KB, M ) - 2*KA - 1
1160                WORK( N+J ) = WORK( N+J+KA )
1161                WORK( J ) = WORK( J+KA )
1162   720       CONTINUE
1163          END IF
1164 *
1165       ELSE
1166 *
1167 *        Transform A, working with the lower triangle
1168 *
1169          IF( UPDATE ) THEN
1170 *
1171 *           Form  inv(S(i))**T * A * inv(S(i))
1172 *
1173             BII = BB( 1, I )
1174             DO 730 J = I1, I
1175                AB( I-J+1, J ) = AB( I-J+1, J ) / BII
1176   730       CONTINUE
1177             DO 740 J = I, MIN( N, I+KA )
1178                AB( J-I+1, I ) = AB( J-I+1, I ) / BII
1179   740       CONTINUE
1180             DO 770 K = I + 1, I + KBT
1181                DO 750 J = K, I + KBT
1182                   AB( J-K+1, K ) = AB( J-K+1, K ) -
1183      $                             BB( J-I+1, I )*AB( K-I+1, I ) -
1184      $                             BB( K-I+1, I )*AB( J-I+1, I ) +
1185      $                             AB( 1, I )*BB( J-I+1, I )*
1186      $                             BB( K-I+1, I )
1187   750          CONTINUE
1188                DO 760 J = I + KBT + 1, MIN( N, I+KA )
1189                   AB( J-K+1, K ) = AB( J-K+1, K ) -
1190      $                             BB( K-I+1, I )*AB( J-I+1, I )
1191   760          CONTINUE
1192   770       CONTINUE
1193             DO 790 J = I1, I
1194                DO 780 K = I + 1, MIN( J+KA, I+KBT )
1195                   AB( K-J+1, J ) = AB( K-J+1, J ) -
1196      $                             BB( K-I+1, I )*AB( I-J+1, J )
1197   780          CONTINUE
1198   790       CONTINUE
1199 *
1200             IF( WANTX ) THEN
1201 *
1202 *              post-multiply X by inv(S(i))
1203 *
1204                CALL SSCAL( NX, ONE / BII, X( 1, I ), 1 )
1205                IF( KBT.GT.0 )
1206      $            CALL SGER( NX, KBT, -ONE, X( 1, I ), 1, BB( 2, I ), 1,
1207      $                       X( 1, I+1 ), LDX )
1208             END IF
1209 *
1210 *           store a(i,i1) in RA1 for use in next loop over K
1211 *
1212             RA1 = AB( I-I1+1, I1 )
1213          END IF
1214 *
1215 *        Generate and apply vectors of rotations to chase all the
1216 *        existing bulges KA positions up toward the top of the band
1217 *
1218          DO 840 K = 1, KB - 1
1219             IF( UPDATE ) THEN
1220 *
1221 *              Determine the rotations which would annihilate the bulge
1222 *              which has in theory just been created
1223 *
1224                IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
1225 *
1226 *                 generate rotation to annihilate a(i,i+k-ka-1)
1227 *
1228                   CALL SLARTG( AB( KA1-K, I+K-KA ), RA1,
1229      $                         WORK( N+I+K-KA ), WORK( I+K-KA ), RA )
1230 *
1231 *                 create nonzero element a(i+k,i+k-ka-1) outside the
1232 *                 band and store it in WORK(m-kb+i+k)
1233 *
1234                   T = -BB( K+1, I )*RA1
1235                   WORK( M-KB+I+K ) = WORK( N+I+K-KA )*T -
1236      $                               WORK( I+K-KA )*AB( KA1, I+K-KA )
1237                   AB( KA1, I+K-KA ) = WORK( I+K-KA )*T +
1238      $                                WORK( N+I+K-KA )*AB( KA1, I+K-KA )
1239                   RA1 = RA
1240                END IF
1241             END IF
1242             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
1243             NR = ( J2+KA-1 ) / KA1
1244             J1 = J2 - ( NR-1 )*KA1
1245             IF( UPDATE ) THEN
1246                J2T = MIN( J2, I-2*KA+K-1 )
1247             ELSE
1248                J2T = J2
1249             END IF
1250             NRT = ( J2T+KA-1 ) / KA1
1251             DO 800 J = J1, J2T, KA1
1252 *
1253 *              create nonzero element a(j+ka,j-1) outside the band
1254 *              and store it in WORK(j)
1255 *
1256                WORK( J ) = WORK( J )*AB( KA1, J-1 )
1257                AB( KA1, J-1 ) = WORK( N+J )*AB( KA1, J-1 )
1258   800       CONTINUE
1259 *
1260 *           generate rotations in 1st set to annihilate elements which
1261 *           have been created outside the band
1262 *
1263             IF( NRT.GT.0 )
1264      $         CALL SLARGV( NRT, AB( KA1, J1 ), INCA, WORK( J1 ), KA1,
1265      $                      WORK( N+J1 ), KA1 )
1266             IF( NR.GT.0 ) THEN
1267 *
1268 *              apply rotations in 1st set from the right
1269 *
1270                DO 810 L = 1, KA - 1
1271                   CALL SLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
1272      $                         INCA, WORK( N+J1 ), WORK( J1 ), KA1 )
1273   810          CONTINUE
1274 *
1275 *              apply rotations in 1st set from both sides to diagonal
1276 *              blocks
1277 *
1278                CALL SLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
1279      $                      AB( 2, J1-1 ), INCA, WORK( N+J1 ),
1280      $                      WORK( J1 ), KA1 )
1281 *
1282             END IF
1283 *
1284 *           start applying rotations in 1st set from the left
1285 *
1286             DO 820 L = KA - 1, KB - K + 1, -1
1287                NRT = ( J2+L-1 ) / KA1
1288                J1T = J2 - ( NRT-1 )*KA1
1289                IF( NRT.GT.0 )
1290      $            CALL SLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
1291      $                         AB( KA1-L, J1T-KA1+L ), INCA,
1292      $                         WORK( N+J1T ), WORK( J1T ), KA1 )
1293   820       CONTINUE
1294 *
1295             IF( WANTX ) THEN
1296 *
1297 *              post-multiply X by product of rotations in 1st set
1298 *
1299                DO 830 J = J1, J2, KA1
1300                   CALL SROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1301      $                       WORK( N+J ), WORK( J ) )
1302   830          CONTINUE
1303             END IF
1304   840    CONTINUE
1305 *
1306          IF( UPDATE ) THEN
1307             IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
1308 *
1309 *              create nonzero element a(i+kbt,i+kbt-ka-1) outside the
1310 *              band and store it in WORK(m-kb+i+kbt)
1311 *
1312                WORK( M-KB+I+KBT ) = -BB( KBT+1, I )*RA1
1313             END IF
1314          END IF
1315 *
1316          DO 880 K = KB, 1, -1
1317             IF( UPDATE ) THEN
1318                J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1
1319             ELSE
1320                J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1321             END IF
1322 *
1323 *           finish applying rotations in 2nd set from the left
1324 *
1325             DO 850 L = KB - K, 1, -1
1326                NRT = ( J2+KA+L-1 ) / KA1
1327                J1T = J2 - ( NRT-1 )*KA1
1328                IF( NRT.GT.0 )
1329      $            CALL SLARTV( NRT, AB( KA1-L+1, J1T+L-1 ), INCA,
1330      $                         AB( KA1-L, J1T+L-1 ), INCA,
1331      $                         WORK( N+M-KB+J1T+KA ),
1332      $                         WORK( M-KB+J1T+KA ), KA1 )
1333   850       CONTINUE
1334             NR = ( J2+KA-1 ) / KA1
1335             J1 = J2 - ( NR-1 )*KA1
1336             DO 860 J = J1, J2, KA1
1337                WORK( M-KB+J ) = WORK( M-KB+J+KA )
1338                WORK( N+M-KB+J ) = WORK( N+M-KB+J+KA )
1339   860       CONTINUE
1340             DO 870 J = J1, J2, KA1
1341 *
1342 *              create nonzero element a(j+ka,j-1) outside the band
1343 *              and store it in WORK(m-kb+j)
1344 *
1345                WORK( M-KB+J ) = WORK( M-KB+J )*AB( KA1, J-1 )
1346                AB( KA1, J-1 ) = WORK( N+M-KB+J )*AB( KA1, J-1 )
1347   870       CONTINUE
1348             IF( UPDATE ) THEN
1349                IF( I+K.GT.KA1 .AND. K.LE.KBT )
1350      $            WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
1351             END IF
1352   880    CONTINUE
1353 *
1354          DO 920 K = KB, 1, -1
1355             J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1356             NR = ( J2+KA-1 ) / KA1
1357             J1 = J2 - ( NR-1 )*KA1
1358             IF( NR.GT.0 ) THEN
1359 *
1360 *              generate rotations in 2nd set to annihilate elements
1361 *              which have been created outside the band
1362 *
1363                CALL SLARGV( NR, AB( KA1, J1 ), INCA, WORK( M-KB+J1 ),
1364      $                      KA1, WORK( N+M-KB+J1 ), KA1 )
1365 *
1366 *              apply rotations in 2nd set from the right
1367 *
1368                DO 890 L = 1, KA - 1
1369                   CALL SLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
1370      $                         INCA, WORK( N+M-KB+J1 ), WORK( M-KB+J1 ),
1371      $                         KA1 )
1372   890          CONTINUE
1373 *
1374 *              apply rotations in 2nd set from both sides to diagonal
1375 *              blocks
1376 *
1377                CALL SLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
1378      $                      AB( 2, J1-1 ), INCA, WORK( N+M-KB+J1 ),
1379      $                      WORK( M-KB+J1 ), KA1 )
1380 *
1381             END IF
1382 *
1383 *           start applying rotations in 2nd set from the left
1384 *
1385             DO 900 L = KA - 1, KB - K + 1, -1
1386                NRT = ( J2+L-1 ) / KA1
1387                J1T = J2 - ( NRT-1 )*KA1
1388                IF( NRT.GT.0 )
1389      $            CALL SLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
1390      $                         AB( KA1-L, J1T-KA1+L ), INCA,
1391      $                         WORK( N+M-KB+J1T ), WORK( M-KB+J1T ),
1392      $                         KA1 )
1393   900       CONTINUE
1394 *
1395             IF( WANTX ) THEN
1396 *
1397 *              post-multiply X by product of rotations in 2nd set
1398 *
1399                DO 910 J = J1, J2, KA1
1400                   CALL SROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1401      $                       WORK( N+M-KB+J ), WORK( M-KB+J ) )
1402   910          CONTINUE
1403             END IF
1404   920    CONTINUE
1405 *
1406          DO 940 K = 1, KB - 1
1407             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
1408 *
1409 *           finish applying rotations in 1st set from the left
1410 *
1411             DO 930 L = KB - K, 1, -1
1412                NRT = ( J2+L-1 ) / KA1
1413                J1T = J2 - ( NRT-1 )*KA1
1414                IF( NRT.GT.0 )
1415      $            CALL SLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
1416      $                         AB( KA1-L, J1T-KA1+L ), INCA,
1417      $                         WORK( N+J1T ), WORK( J1T ), KA1 )
1418   930       CONTINUE
1419   940    CONTINUE
1420 *
1421          IF( KB.GT.1 ) THEN
1422             DO 950 J = 2, MIN( I+KB, M ) - 2*KA - 1
1423                WORK( N+J ) = WORK( N+J+KA )
1424                WORK( J ) = WORK( J+KA )
1425   950       CONTINUE
1426          END IF
1427 *
1428       END IF
1429 *
1430       GO TO 490
1431 *
1432 *     End of SSBGST
1433 *
1434       END