LAPACK 3.3.0

dsbgst.f

Go to the documentation of this file.
00001       SUBROUTINE DSBGST( VECT, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, X,
00002      $                   LDX, WORK, INFO )
00003 *
00004 *  -- LAPACK routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          UPLO, VECT
00011       INTEGER            INFO, KA, KB, LDAB, LDBB, LDX, N
00012 *     ..
00013 *     .. Array Arguments ..
00014       DOUBLE PRECISION   AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
00015      $                   X( LDX, * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  DSBGST reduces a real symmetric-definite banded generalized
00022 *  eigenproblem  A*x = lambda*B*x  to standard form  C*y = lambda*y,
00023 *  such that C has the same bandwidth as A.
00024 *
00025 *  B must have been previously factorized as S**T*S by DPBSTF, using a
00026 *  split Cholesky factorization. A is overwritten by C = X**T*A*X, where
00027 *  X = S**(-1)*Q and Q is an orthogonal matrix chosen to preserve the
00028 *  bandwidth of A.
00029 *
00030 *  Arguments
00031 *  =========
00032 *
00033 *  VECT    (input) CHARACTER*1
00034 *          = 'N':  do not form the transformation matrix X;
00035 *          = 'V':  form X.
00036 *
00037 *  UPLO    (input) CHARACTER*1
00038 *          = 'U':  Upper triangle of A is stored;
00039 *          = 'L':  Lower triangle of A is stored.
00040 *
00041 *  N       (input) INTEGER
00042 *          The order of the matrices A and B.  N >= 0.
00043 *
00044 *  KA      (input) INTEGER
00045 *          The number of superdiagonals of the matrix A if UPLO = 'U',
00046 *          or the number of subdiagonals if UPLO = 'L'.  KA >= 0.
00047 *
00048 *  KB      (input) INTEGER
00049 *          The number of superdiagonals of the matrix B if UPLO = 'U',
00050 *          or the number of subdiagonals if UPLO = 'L'.  KA >= KB >= 0.
00051 *
00052 *  AB      (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
00053 *          On entry, the upper or lower triangle of the symmetric band
00054 *          matrix A, stored in the first ka+1 rows of the array.  The
00055 *          j-th column of A is stored in the j-th column of the array AB
00056 *          as follows:
00057 *          if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;
00058 *          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+ka).
00059 *
00060 *          On exit, the transformed matrix X**T*A*X, stored in the same
00061 *          format as A.
00062 *
00063 *  LDAB    (input) INTEGER
00064 *          The leading dimension of the array AB.  LDAB >= KA+1.
00065 *
00066 *  BB      (input) DOUBLE PRECISION array, dimension (LDBB,N)
00067 *          The banded factor S from the split Cholesky factorization of
00068 *          B, as returned by DPBSTF, stored in the first KB+1 rows of
00069 *          the array.
00070 *
00071 *  LDBB    (input) INTEGER
00072 *          The leading dimension of the array BB.  LDBB >= KB+1.
00073 *
00074 *  X       (output) DOUBLE PRECISION array, dimension (LDX,N)
00075 *          If VECT = 'V', the n-by-n matrix X.
00076 *          If VECT = 'N', the array X is not referenced.
00077 *
00078 *  LDX     (input) INTEGER
00079 *          The leading dimension of the array X.
00080 *          LDX >= max(1,N) if VECT = 'V'; LDX >= 1 otherwise.
00081 *
00082 *  WORK    (workspace) DOUBLE PRECISION array, dimension (2*N)
00083 *
00084 *  INFO    (output) INTEGER
00085 *          = 0:  successful exit
00086 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00087 *
00088 *  =====================================================================
00089 *
00090 *     .. Parameters ..
00091       DOUBLE PRECISION   ZERO, ONE
00092       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00093 *     ..
00094 *     .. Local Scalars ..
00095       LOGICAL            UPDATE, UPPER, WANTX
00096       INTEGER            I, I0, I1, I2, INCA, J, J1, J1T, J2, J2T, K,
00097      $                   KA1, KB1, KBT, L, M, NR, NRT, NX
00098       DOUBLE PRECISION   BII, RA, RA1, T
00099 *     ..
00100 *     .. External Functions ..
00101       LOGICAL            LSAME
00102       EXTERNAL           LSAME
00103 *     ..
00104 *     .. External Subroutines ..
00105       EXTERNAL           DGER, DLAR2V, DLARGV, DLARTG, DLARTV, DLASET,
00106      $                   DROT, DSCAL, XERBLA
00107 *     ..
00108 *     .. Intrinsic Functions ..
00109       INTRINSIC          MAX, MIN
00110 *     ..
00111 *     .. Executable Statements ..
00112 *
00113 *     Test the input parameters
00114 *
00115       WANTX = LSAME( VECT, 'V' )
00116       UPPER = LSAME( UPLO, 'U' )
00117       KA1 = KA + 1
00118       KB1 = KB + 1
00119       INFO = 0
00120       IF( .NOT.WANTX .AND. .NOT.LSAME( VECT, 'N' ) ) THEN
00121          INFO = -1
00122       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00123          INFO = -2
00124       ELSE IF( N.LT.0 ) THEN
00125          INFO = -3
00126       ELSE IF( KA.LT.0 ) THEN
00127          INFO = -4
00128       ELSE IF( KB.LT.0 .OR. KB.GT.KA ) THEN
00129          INFO = -5
00130       ELSE IF( LDAB.LT.KA+1 ) THEN
00131          INFO = -7
00132       ELSE IF( LDBB.LT.KB+1 ) THEN
00133          INFO = -9
00134       ELSE IF( LDX.LT.1 .OR. WANTX .AND. LDX.LT.MAX( 1, N ) ) THEN
00135          INFO = -11
00136       END IF
00137       IF( INFO.NE.0 ) THEN
00138          CALL XERBLA( 'DSBGST', -INFO )
00139          RETURN
00140       END IF
00141 *
00142 *     Quick return if possible
00143 *
00144       IF( N.EQ.0 )
00145      $   RETURN
00146 *
00147       INCA = LDAB*KA1
00148 *
00149 *     Initialize X to the unit matrix, if needed
00150 *
00151       IF( WANTX )
00152      $   CALL DLASET( 'Full', N, N, ZERO, ONE, X, LDX )
00153 *
00154 *     Set M to the splitting point m. It must be the same value as is
00155 *     used in DPBSTF. The chosen value allows the arrays WORK and RWORK
00156 *     to be of dimension (N).
00157 *
00158       M = ( N+KB ) / 2
00159 *
00160 *     The routine works in two phases, corresponding to the two halves
00161 *     of the split Cholesky factorization of B as S**T*S where
00162 *
00163 *     S = ( U    )
00164 *         ( M  L )
00165 *
00166 *     with U upper triangular of order m, and L lower triangular of
00167 *     order n-m. S has the same bandwidth as B.
00168 *
00169 *     S is treated as a product of elementary matrices:
00170 *
00171 *     S = S(m)*S(m-1)*...*S(2)*S(1)*S(m+1)*S(m+2)*...*S(n-1)*S(n)
00172 *
00173 *     where S(i) is determined by the i-th row of S.
00174 *
00175 *     In phase 1, the index i takes the values n, n-1, ... , m+1;
00176 *     in phase 2, it takes the values 1, 2, ... , m.
00177 *
00178 *     For each value of i, the current matrix A is updated by forming
00179 *     inv(S(i))**T*A*inv(S(i)). This creates a triangular bulge outside
00180 *     the band of A. The bulge is then pushed down toward the bottom of
00181 *     A in phase 1, and up toward the top of A in phase 2, by applying
00182 *     plane rotations.
00183 *
00184 *     There are kb*(kb+1)/2 elements in the bulge, but at most 2*kb-1
00185 *     of them are linearly independent, so annihilating a bulge requires
00186 *     only 2*kb-1 plane rotations. The rotations are divided into a 1st
00187 *     set of kb-1 rotations, and a 2nd set of kb rotations.
00188 *
00189 *     Wherever possible, rotations are generated and applied in vector
00190 *     operations of length NR between the indices J1 and J2 (sometimes
00191 *     replaced by modified values NRT, J1T or J2T).
00192 *
00193 *     The cosines and sines of the rotations are stored in the array
00194 *     WORK. The cosines of the 1st set of rotations are stored in
00195 *     elements n+2:n+m-kb-1 and the sines of the 1st set in elements
00196 *     2:m-kb-1; the cosines of the 2nd set are stored in elements
00197 *     n+m-kb+1:2*n and the sines of the second set in elements m-kb+1:n.
00198 *
00199 *     The bulges are not formed explicitly; nonzero elements outside the
00200 *     band are created only when they are required for generating new
00201 *     rotations; they are stored in the array WORK, in positions where
00202 *     they are later overwritten by the sines of the rotations which
00203 *     annihilate them.
00204 *
00205 *     **************************** Phase 1 *****************************
00206 *
00207 *     The logical structure of this phase is:
00208 *
00209 *     UPDATE = .TRUE.
00210 *     DO I = N, M + 1, -1
00211 *        use S(i) to update A and create a new bulge
00212 *        apply rotations to push all bulges KA positions downward
00213 *     END DO
00214 *     UPDATE = .FALSE.
00215 *     DO I = M + KA + 1, N - 1
00216 *        apply rotations to push all bulges KA positions downward
00217 *     END DO
00218 *
00219 *     To avoid duplicating code, the two loops are merged.
00220 *
00221       UPDATE = .TRUE.
00222       I = N + 1
00223    10 CONTINUE
00224       IF( UPDATE ) THEN
00225          I = I - 1
00226          KBT = MIN( KB, I-1 )
00227          I0 = I - 1
00228          I1 = MIN( N, I+KA )
00229          I2 = I - KBT + KA1
00230          IF( I.LT.M+1 ) THEN
00231             UPDATE = .FALSE.
00232             I = I + 1
00233             I0 = M
00234             IF( KA.EQ.0 )
00235      $         GO TO 480
00236             GO TO 10
00237          END IF
00238       ELSE
00239          I = I + KA
00240          IF( I.GT.N-1 )
00241      $      GO TO 480
00242       END IF
00243 *
00244       IF( UPPER ) THEN
00245 *
00246 *        Transform A, working with the upper triangle
00247 *
00248          IF( UPDATE ) THEN
00249 *
00250 *           Form  inv(S(i))**T * A * inv(S(i))
00251 *
00252             BII = BB( KB1, I )
00253             DO 20 J = I, I1
00254                AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII
00255    20       CONTINUE
00256             DO 30 J = MAX( 1, I-KA ), I
00257                AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII
00258    30       CONTINUE
00259             DO 60 K = I - KBT, I - 1
00260                DO 40 J = I - KBT, K
00261                   AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
00262      $                               BB( J-I+KB1, I )*AB( K-I+KA1, I ) -
00263      $                               BB( K-I+KB1, I )*AB( J-I+KA1, I ) +
00264      $                               AB( KA1, I )*BB( J-I+KB1, I )*
00265      $                               BB( K-I+KB1, I )
00266    40          CONTINUE
00267                DO 50 J = MAX( 1, I-KA ), I - KBT - 1
00268                   AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
00269      $                               BB( K-I+KB1, I )*AB( J-I+KA1, I )
00270    50          CONTINUE
00271    60       CONTINUE
00272             DO 80 J = I, I1
00273                DO 70 K = MAX( J-KA, I-KBT ), I - 1
00274                   AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
00275      $                               BB( K-I+KB1, I )*AB( I-J+KA1, J )
00276    70          CONTINUE
00277    80       CONTINUE
00278 *
00279             IF( WANTX ) THEN
00280 *
00281 *              post-multiply X by inv(S(i))
00282 *
00283                CALL DSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
00284                IF( KBT.GT.0 )
00285      $            CALL DGER( N-M, KBT, -ONE, X( M+1, I ), 1,
00286      $                       BB( KB1-KBT, I ), 1, X( M+1, I-KBT ), LDX )
00287             END IF
00288 *
00289 *           store a(i,i1) in RA1 for use in next loop over K
00290 *
00291             RA1 = AB( I-I1+KA1, I1 )
00292          END IF
00293 *
00294 *        Generate and apply vectors of rotations to chase all the
00295 *        existing bulges KA positions down toward the bottom of the
00296 *        band
00297 *
00298          DO 130 K = 1, KB - 1
00299             IF( UPDATE ) THEN
00300 *
00301 *              Determine the rotations which would annihilate the bulge
00302 *              which has in theory just been created
00303 *
00304                IF( I-K+KA.LT.N .AND. I-K.GT.1 ) THEN
00305 *
00306 *                 generate rotation to annihilate a(i,i-k+ka+1)
00307 *
00308                   CALL DLARTG( AB( K+1, I-K+KA ), RA1,
00309      $                         WORK( N+I-K+KA-M ), WORK( I-K+KA-M ),
00310      $                         RA )
00311 *
00312 *                 create nonzero element a(i-k,i-k+ka+1) outside the
00313 *                 band and store it in WORK(i-k)
00314 *
00315                   T = -BB( KB1-K, I )*RA1
00316                   WORK( I-K ) = WORK( N+I-K+KA-M )*T -
00317      $                          WORK( I-K+KA-M )*AB( 1, I-K+KA )
00318                   AB( 1, I-K+KA ) = WORK( I-K+KA-M )*T +
00319      $                              WORK( N+I-K+KA-M )*AB( 1, I-K+KA )
00320                   RA1 = RA
00321                END IF
00322             END IF
00323             J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
00324             NR = ( N-J2+KA ) / KA1
00325             J1 = J2 + ( NR-1 )*KA1
00326             IF( UPDATE ) THEN
00327                J2T = MAX( J2, I+2*KA-K+1 )
00328             ELSE
00329                J2T = J2
00330             END IF
00331             NRT = ( N-J2T+KA ) / KA1
00332             DO 90 J = J2T, J1, KA1
00333 *
00334 *              create nonzero element a(j-ka,j+1) outside the band
00335 *              and store it in WORK(j-m)
00336 *
00337                WORK( J-M ) = WORK( J-M )*AB( 1, J+1 )
00338                AB( 1, J+1 ) = WORK( N+J-M )*AB( 1, J+1 )
00339    90       CONTINUE
00340 *
00341 *           generate rotations in 1st set to annihilate elements which
00342 *           have been created outside the band
00343 *
00344             IF( NRT.GT.0 )
00345      $         CALL DLARGV( NRT, AB( 1, J2T ), INCA, WORK( J2T-M ), KA1,
00346      $                      WORK( N+J2T-M ), KA1 )
00347             IF( NR.GT.0 ) THEN
00348 *
00349 *              apply rotations in 1st set from the right
00350 *
00351                DO 100 L = 1, KA - 1
00352                   CALL DLARTV( NR, AB( KA1-L, J2 ), INCA,
00353      $                         AB( KA-L, J2+1 ), INCA, WORK( N+J2-M ),
00354      $                         WORK( J2-M ), KA1 )
00355   100          CONTINUE
00356 *
00357 *              apply rotations in 1st set from both sides to diagonal
00358 *              blocks
00359 *
00360                CALL DLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ),
00361      $                      AB( KA, J2+1 ), INCA, WORK( N+J2-M ),
00362      $                      WORK( J2-M ), KA1 )
00363 *
00364             END IF
00365 *
00366 *           start applying rotations in 1st set from the left
00367 *
00368             DO 110 L = KA - 1, KB - K + 1, -1
00369                NRT = ( N-J2+L ) / KA1
00370                IF( NRT.GT.0 )
00371      $            CALL DLARTV( NRT, AB( L, J2+KA1-L ), INCA,
00372      $                         AB( L+1, J2+KA1-L ), INCA,
00373      $                         WORK( N+J2-M ), WORK( J2-M ), KA1 )
00374   110       CONTINUE
00375 *
00376             IF( WANTX ) THEN
00377 *
00378 *              post-multiply X by product of rotations in 1st set
00379 *
00380                DO 120 J = J2, J1, KA1
00381                   CALL DROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
00382      $                       WORK( N+J-M ), WORK( J-M ) )
00383   120          CONTINUE
00384             END IF
00385   130    CONTINUE
00386 *
00387          IF( UPDATE ) THEN
00388             IF( I2.LE.N .AND. KBT.GT.0 ) THEN
00389 *
00390 *              create nonzero element a(i-kbt,i-kbt+ka+1) outside the
00391 *              band and store it in WORK(i-kbt)
00392 *
00393                WORK( I-KBT ) = -BB( KB1-KBT, I )*RA1
00394             END IF
00395          END IF
00396 *
00397          DO 170 K = KB, 1, -1
00398             IF( UPDATE ) THEN
00399                J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1
00400             ELSE
00401                J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
00402             END IF
00403 *
00404 *           finish applying rotations in 2nd set from the left
00405 *
00406             DO 140 L = KB - K, 1, -1
00407                NRT = ( N-J2+KA+L ) / KA1
00408                IF( NRT.GT.0 )
00409      $            CALL DLARTV( NRT, AB( L, J2-L+1 ), INCA,
00410      $                         AB( L+1, J2-L+1 ), INCA, WORK( N+J2-KA ),
00411      $                         WORK( J2-KA ), KA1 )
00412   140       CONTINUE
00413             NR = ( N-J2+KA ) / KA1
00414             J1 = J2 + ( NR-1 )*KA1
00415             DO 150 J = J1, J2, -KA1
00416                WORK( J ) = WORK( J-KA )
00417                WORK( N+J ) = WORK( N+J-KA )
00418   150       CONTINUE
00419             DO 160 J = J2, J1, KA1
00420 *
00421 *              create nonzero element a(j-ka,j+1) outside the band
00422 *              and store it in WORK(j)
00423 *
00424                WORK( J ) = WORK( J )*AB( 1, J+1 )
00425                AB( 1, J+1 ) = WORK( N+J )*AB( 1, J+1 )
00426   160       CONTINUE
00427             IF( UPDATE ) THEN
00428                IF( I-K.LT.N-KA .AND. K.LE.KBT )
00429      $            WORK( I-K+KA ) = WORK( I-K )
00430             END IF
00431   170    CONTINUE
00432 *
00433          DO 210 K = KB, 1, -1
00434             J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
00435             NR = ( N-J2+KA ) / KA1
00436             J1 = J2 + ( NR-1 )*KA1
00437             IF( NR.GT.0 ) THEN
00438 *
00439 *              generate rotations in 2nd set to annihilate elements
00440 *              which have been created outside the band
00441 *
00442                CALL DLARGV( NR, AB( 1, J2 ), INCA, WORK( J2 ), KA1,
00443      $                      WORK( N+J2 ), KA1 )
00444 *
00445 *              apply rotations in 2nd set from the right
00446 *
00447                DO 180 L = 1, KA - 1
00448                   CALL DLARTV( NR, AB( KA1-L, J2 ), INCA,
00449      $                         AB( KA-L, J2+1 ), INCA, WORK( N+J2 ),
00450      $                         WORK( J2 ), KA1 )
00451   180          CONTINUE
00452 *
00453 *              apply rotations in 2nd set from both sides to diagonal
00454 *              blocks
00455 *
00456                CALL DLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ),
00457      $                      AB( KA, J2+1 ), INCA, WORK( N+J2 ),
00458      $                      WORK( J2 ), KA1 )
00459 *
00460             END IF
00461 *
00462 *           start applying rotations in 2nd set from the left
00463 *
00464             DO 190 L = KA - 1, KB - K + 1, -1
00465                NRT = ( N-J2+L ) / KA1
00466                IF( NRT.GT.0 )
00467      $            CALL DLARTV( NRT, AB( L, J2+KA1-L ), INCA,
00468      $                         AB( L+1, J2+KA1-L ), INCA, WORK( N+J2 ),
00469      $                         WORK( J2 ), KA1 )
00470   190       CONTINUE
00471 *
00472             IF( WANTX ) THEN
00473 *
00474 *              post-multiply X by product of rotations in 2nd set
00475 *
00476                DO 200 J = J2, J1, KA1
00477                   CALL DROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
00478      $                       WORK( N+J ), WORK( J ) )
00479   200          CONTINUE
00480             END IF
00481   210    CONTINUE
00482 *
00483          DO 230 K = 1, KB - 1
00484             J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
00485 *
00486 *           finish applying rotations in 1st set from the left
00487 *
00488             DO 220 L = KB - K, 1, -1
00489                NRT = ( N-J2+L ) / KA1
00490                IF( NRT.GT.0 )
00491      $            CALL DLARTV( NRT, AB( L, J2+KA1-L ), INCA,
00492      $                         AB( L+1, J2+KA1-L ), INCA,
00493      $                         WORK( N+J2-M ), WORK( J2-M ), KA1 )
00494   220       CONTINUE
00495   230    CONTINUE
00496 *
00497          IF( KB.GT.1 ) THEN
00498             DO 240 J = N - 1, I - KB + 2*KA + 1, -1
00499                WORK( N+J-M ) = WORK( N+J-KA-M )
00500                WORK( J-M ) = WORK( J-KA-M )
00501   240       CONTINUE
00502          END IF
00503 *
00504       ELSE
00505 *
00506 *        Transform A, working with the lower triangle
00507 *
00508          IF( UPDATE ) THEN
00509 *
00510 *           Form  inv(S(i))**T * A * inv(S(i))
00511 *
00512             BII = BB( 1, I )
00513             DO 250 J = I, I1
00514                AB( J-I+1, I ) = AB( J-I+1, I ) / BII
00515   250       CONTINUE
00516             DO 260 J = MAX( 1, I-KA ), I
00517                AB( I-J+1, J ) = AB( I-J+1, J ) / BII
00518   260       CONTINUE
00519             DO 290 K = I - KBT, I - 1
00520                DO 270 J = I - KBT, K
00521                   AB( K-J+1, J ) = AB( K-J+1, J ) -
00522      $                             BB( I-J+1, J )*AB( I-K+1, K ) -
00523      $                             BB( I-K+1, K )*AB( I-J+1, J ) +
00524      $                             AB( 1, I )*BB( I-J+1, J )*
00525      $                             BB( I-K+1, K )
00526   270          CONTINUE
00527                DO 280 J = MAX( 1, I-KA ), I - KBT - 1
00528                   AB( K-J+1, J ) = AB( K-J+1, J ) -
00529      $                             BB( I-K+1, K )*AB( I-J+1, J )
00530   280          CONTINUE
00531   290       CONTINUE
00532             DO 310 J = I, I1
00533                DO 300 K = MAX( J-KA, I-KBT ), I - 1
00534                   AB( J-K+1, K ) = AB( J-K+1, K ) -
00535      $                             BB( I-K+1, K )*AB( J-I+1, I )
00536   300          CONTINUE
00537   310       CONTINUE
00538 *
00539             IF( WANTX ) THEN
00540 *
00541 *              post-multiply X by inv(S(i))
00542 *
00543                CALL DSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
00544                IF( KBT.GT.0 )
00545      $            CALL DGER( N-M, KBT, -ONE, X( M+1, I ), 1,
00546      $                       BB( KBT+1, I-KBT ), LDBB-1,
00547      $                       X( M+1, I-KBT ), LDX )
00548             END IF
00549 *
00550 *           store a(i1,i) in RA1 for use in next loop over K
00551 *
00552             RA1 = AB( I1-I+1, I )
00553          END IF
00554 *
00555 *        Generate and apply vectors of rotations to chase all the
00556 *        existing bulges KA positions down toward the bottom of the
00557 *        band
00558 *
00559          DO 360 K = 1, KB - 1
00560             IF( UPDATE ) THEN
00561 *
00562 *              Determine the rotations which would annihilate the bulge
00563 *              which has in theory just been created
00564 *
00565                IF( I-K+KA.LT.N .AND. I-K.GT.1 ) THEN
00566 *
00567 *                 generate rotation to annihilate a(i-k+ka+1,i)
00568 *
00569                   CALL DLARTG( AB( KA1-K, I ), RA1, WORK( N+I-K+KA-M ),
00570      $                         WORK( I-K+KA-M ), RA )
00571 *
00572 *                 create nonzero element a(i-k+ka+1,i-k) outside the
00573 *                 band and store it in WORK(i-k)
00574 *
00575                   T = -BB( K+1, I-K )*RA1
00576                   WORK( I-K ) = WORK( N+I-K+KA-M )*T -
00577      $                          WORK( I-K+KA-M )*AB( KA1, I-K )
00578                   AB( KA1, I-K ) = WORK( I-K+KA-M )*T +
00579      $                             WORK( N+I-K+KA-M )*AB( KA1, I-K )
00580                   RA1 = RA
00581                END IF
00582             END IF
00583             J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
00584             NR = ( N-J2+KA ) / KA1
00585             J1 = J2 + ( NR-1 )*KA1
00586             IF( UPDATE ) THEN
00587                J2T = MAX( J2, I+2*KA-K+1 )
00588             ELSE
00589                J2T = J2
00590             END IF
00591             NRT = ( N-J2T+KA ) / KA1
00592             DO 320 J = J2T, J1, KA1
00593 *
00594 *              create nonzero element a(j+1,j-ka) outside the band
00595 *              and store it in WORK(j-m)
00596 *
00597                WORK( J-M ) = WORK( J-M )*AB( KA1, J-KA+1 )
00598                AB( KA1, J-KA+1 ) = WORK( N+J-M )*AB( KA1, J-KA+1 )
00599   320       CONTINUE
00600 *
00601 *           generate rotations in 1st set to annihilate elements which
00602 *           have been created outside the band
00603 *
00604             IF( NRT.GT.0 )
00605      $         CALL DLARGV( NRT, AB( KA1, J2T-KA ), INCA, WORK( J2T-M ),
00606      $                      KA1, WORK( N+J2T-M ), KA1 )
00607             IF( NR.GT.0 ) THEN
00608 *
00609 *              apply rotations in 1st set from the left
00610 *
00611                DO 330 L = 1, KA - 1
00612                   CALL DLARTV( NR, AB( L+1, J2-L ), INCA,
00613      $                         AB( L+2, J2-L ), INCA, WORK( N+J2-M ),
00614      $                         WORK( J2-M ), KA1 )
00615   330          CONTINUE
00616 *
00617 *              apply rotations in 1st set from both sides to diagonal
00618 *              blocks
00619 *
00620                CALL DLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
00621      $                      INCA, WORK( N+J2-M ), WORK( J2-M ), KA1 )
00622 *
00623             END IF
00624 *
00625 *           start applying rotations in 1st set from the right
00626 *
00627             DO 340 L = KA - 1, KB - K + 1, -1
00628                NRT = ( N-J2+L ) / KA1
00629                IF( NRT.GT.0 )
00630      $            CALL DLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
00631      $                         AB( KA1-L, J2+1 ), INCA, WORK( N+J2-M ),
00632      $                         WORK( J2-M ), KA1 )
00633   340       CONTINUE
00634 *
00635             IF( WANTX ) THEN
00636 *
00637 *              post-multiply X by product of rotations in 1st set
00638 *
00639                DO 350 J = J2, J1, KA1
00640                   CALL DROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
00641      $                       WORK( N+J-M ), WORK( J-M ) )
00642   350          CONTINUE
00643             END IF
00644   360    CONTINUE
00645 *
00646          IF( UPDATE ) THEN
00647             IF( I2.LE.N .AND. KBT.GT.0 ) THEN
00648 *
00649 *              create nonzero element a(i-kbt+ka+1,i-kbt) outside the
00650 *              band and store it in WORK(i-kbt)
00651 *
00652                WORK( I-KBT ) = -BB( KBT+1, I-KBT )*RA1
00653             END IF
00654          END IF
00655 *
00656          DO 400 K = KB, 1, -1
00657             IF( UPDATE ) THEN
00658                J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1
00659             ELSE
00660                J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
00661             END IF
00662 *
00663 *           finish applying rotations in 2nd set from the right
00664 *
00665             DO 370 L = KB - K, 1, -1
00666                NRT = ( N-J2+KA+L ) / KA1
00667                IF( NRT.GT.0 )
00668      $            CALL DLARTV( NRT, AB( KA1-L+1, J2-KA ), INCA,
00669      $                         AB( KA1-L, J2-KA+1 ), INCA,
00670      $                         WORK( N+J2-KA ), WORK( J2-KA ), KA1 )
00671   370       CONTINUE
00672             NR = ( N-J2+KA ) / KA1
00673             J1 = J2 + ( NR-1 )*KA1
00674             DO 380 J = J1, J2, -KA1
00675                WORK( J ) = WORK( J-KA )
00676                WORK( N+J ) = WORK( N+J-KA )
00677   380       CONTINUE
00678             DO 390 J = J2, J1, KA1
00679 *
00680 *              create nonzero element a(j+1,j-ka) outside the band
00681 *              and store it in WORK(j)
00682 *
00683                WORK( J ) = WORK( J )*AB( KA1, J-KA+1 )
00684                AB( KA1, J-KA+1 ) = WORK( N+J )*AB( KA1, J-KA+1 )
00685   390       CONTINUE
00686             IF( UPDATE ) THEN
00687                IF( I-K.LT.N-KA .AND. K.LE.KBT )
00688      $            WORK( I-K+KA ) = WORK( I-K )
00689             END IF
00690   400    CONTINUE
00691 *
00692          DO 440 K = KB, 1, -1
00693             J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
00694             NR = ( N-J2+KA ) / KA1
00695             J1 = J2 + ( NR-1 )*KA1
00696             IF( NR.GT.0 ) THEN
00697 *
00698 *              generate rotations in 2nd set to annihilate elements
00699 *              which have been created outside the band
00700 *
00701                CALL DLARGV( NR, AB( KA1, J2-KA ), INCA, WORK( J2 ), KA1,
00702      $                      WORK( N+J2 ), KA1 )
00703 *
00704 *              apply rotations in 2nd set from the left
00705 *
00706                DO 410 L = 1, KA - 1
00707                   CALL DLARTV( NR, AB( L+1, J2-L ), INCA,
00708      $                         AB( L+2, J2-L ), INCA, WORK( N+J2 ),
00709      $                         WORK( J2 ), KA1 )
00710   410          CONTINUE
00711 *
00712 *              apply rotations in 2nd set from both sides to diagonal
00713 *              blocks
00714 *
00715                CALL DLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
00716      $                      INCA, WORK( N+J2 ), WORK( J2 ), KA1 )
00717 *
00718             END IF
00719 *
00720 *           start applying rotations in 2nd set from the right
00721 *
00722             DO 420 L = KA - 1, KB - K + 1, -1
00723                NRT = ( N-J2+L ) / KA1
00724                IF( NRT.GT.0 )
00725      $            CALL DLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
00726      $                         AB( KA1-L, J2+1 ), INCA, WORK( N+J2 ),
00727      $                         WORK( J2 ), KA1 )
00728   420       CONTINUE
00729 *
00730             IF( WANTX ) THEN
00731 *
00732 *              post-multiply X by product of rotations in 2nd set
00733 *
00734                DO 430 J = J2, J1, KA1
00735                   CALL DROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
00736      $                       WORK( N+J ), WORK( J ) )
00737   430          CONTINUE
00738             END IF
00739   440    CONTINUE
00740 *
00741          DO 460 K = 1, KB - 1
00742             J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
00743 *
00744 *           finish applying rotations in 1st set from the right
00745 *
00746             DO 450 L = KB - K, 1, -1
00747                NRT = ( N-J2+L ) / KA1
00748                IF( NRT.GT.0 )
00749      $            CALL DLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
00750      $                         AB( KA1-L, J2+1 ), INCA, WORK( N+J2-M ),
00751      $                         WORK( J2-M ), KA1 )
00752   450       CONTINUE
00753   460    CONTINUE
00754 *
00755          IF( KB.GT.1 ) THEN
00756             DO 470 J = N - 1, I - KB + 2*KA + 1, -1
00757                WORK( N+J-M ) = WORK( N+J-KA-M )
00758                WORK( J-M ) = WORK( J-KA-M )
00759   470       CONTINUE
00760          END IF
00761 *
00762       END IF
00763 *
00764       GO TO 10
00765 *
00766   480 CONTINUE
00767 *
00768 *     **************************** Phase 2 *****************************
00769 *
00770 *     The logical structure of this phase is:
00771 *
00772 *     UPDATE = .TRUE.
00773 *     DO I = 1, M
00774 *        use S(i) to update A and create a new bulge
00775 *        apply rotations to push all bulges KA positions upward
00776 *     END DO
00777 *     UPDATE = .FALSE.
00778 *     DO I = M - KA - 1, 2, -1
00779 *        apply rotations to push all bulges KA positions upward
00780 *     END DO
00781 *
00782 *     To avoid duplicating code, the two loops are merged.
00783 *
00784       UPDATE = .TRUE.
00785       I = 0
00786   490 CONTINUE
00787       IF( UPDATE ) THEN
00788          I = I + 1
00789          KBT = MIN( KB, M-I )
00790          I0 = I + 1
00791          I1 = MAX( 1, I-KA )
00792          I2 = I + KBT - KA1
00793          IF( I.GT.M ) THEN
00794             UPDATE = .FALSE.
00795             I = I - 1
00796             I0 = M + 1
00797             IF( KA.EQ.0 )
00798      $         RETURN
00799             GO TO 490
00800          END IF
00801       ELSE
00802          I = I - KA
00803          IF( I.LT.2 )
00804      $      RETURN
00805       END IF
00806 *
00807       IF( I.LT.M-KBT ) THEN
00808          NX = M
00809       ELSE
00810          NX = N
00811       END IF
00812 *
00813       IF( UPPER ) THEN
00814 *
00815 *        Transform A, working with the upper triangle
00816 *
00817          IF( UPDATE ) THEN
00818 *
00819 *           Form  inv(S(i))**T * A * inv(S(i))
00820 *
00821             BII = BB( KB1, I )
00822             DO 500 J = I1, I
00823                AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII
00824   500       CONTINUE
00825             DO 510 J = I, MIN( N, I+KA )
00826                AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII
00827   510       CONTINUE
00828             DO 540 K = I + 1, I + KBT
00829                DO 520 J = K, I + KBT
00830                   AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
00831      $                               BB( I-J+KB1, J )*AB( I-K+KA1, K ) -
00832      $                               BB( I-K+KB1, K )*AB( I-J+KA1, J ) +
00833      $                               AB( KA1, I )*BB( I-J+KB1, J )*
00834      $                               BB( I-K+KB1, K )
00835   520          CONTINUE
00836                DO 530 J = I + KBT + 1, MIN( N, I+KA )
00837                   AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
00838      $                               BB( I-K+KB1, K )*AB( I-J+KA1, J )
00839   530          CONTINUE
00840   540       CONTINUE
00841             DO 560 J = I1, I
00842                DO 550 K = I + 1, MIN( J+KA, I+KBT )
00843                   AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
00844      $                               BB( I-K+KB1, K )*AB( J-I+KA1, I )
00845   550          CONTINUE
00846   560       CONTINUE
00847 *
00848             IF( WANTX ) THEN
00849 *
00850 *              post-multiply X by inv(S(i))
00851 *
00852                CALL DSCAL( NX, ONE / BII, X( 1, I ), 1 )
00853                IF( KBT.GT.0 )
00854      $            CALL DGER( NX, KBT, -ONE, X( 1, I ), 1, BB( KB, I+1 ),
00855      $                       LDBB-1, X( 1, I+1 ), LDX )
00856             END IF
00857 *
00858 *           store a(i1,i) in RA1 for use in next loop over K
00859 *
00860             RA1 = AB( I1-I+KA1, I )
00861          END IF
00862 *
00863 *        Generate and apply vectors of rotations to chase all the
00864 *        existing bulges KA positions up toward the top of the band
00865 *
00866          DO 610 K = 1, KB - 1
00867             IF( UPDATE ) THEN
00868 *
00869 *              Determine the rotations which would annihilate the bulge
00870 *              which has in theory just been created
00871 *
00872                IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
00873 *
00874 *                 generate rotation to annihilate a(i+k-ka-1,i)
00875 *
00876                   CALL DLARTG( AB( K+1, I ), RA1, WORK( N+I+K-KA ),
00877      $                         WORK( I+K-KA ), RA )
00878 *
00879 *                 create nonzero element a(i+k-ka-1,i+k) outside the
00880 *                 band and store it in WORK(m-kb+i+k)
00881 *
00882                   T = -BB( KB1-K, I+K )*RA1
00883                   WORK( M-KB+I+K ) = WORK( N+I+K-KA )*T -
00884      $                               WORK( I+K-KA )*AB( 1, I+K )
00885                   AB( 1, I+K ) = WORK( I+K-KA )*T +
00886      $                           WORK( N+I+K-KA )*AB( 1, I+K )
00887                   RA1 = RA
00888                END IF
00889             END IF
00890             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
00891             NR = ( J2+KA-1 ) / KA1
00892             J1 = J2 - ( NR-1 )*KA1
00893             IF( UPDATE ) THEN
00894                J2T = MIN( J2, I-2*KA+K-1 )
00895             ELSE
00896                J2T = J2
00897             END IF
00898             NRT = ( J2T+KA-1 ) / KA1
00899             DO 570 J = J1, J2T, KA1
00900 *
00901 *              create nonzero element a(j-1,j+ka) outside the band
00902 *              and store it in WORK(j)
00903 *
00904                WORK( J ) = WORK( J )*AB( 1, J+KA-1 )
00905                AB( 1, J+KA-1 ) = WORK( N+J )*AB( 1, J+KA-1 )
00906   570       CONTINUE
00907 *
00908 *           generate rotations in 1st set to annihilate elements which
00909 *           have been created outside the band
00910 *
00911             IF( NRT.GT.0 )
00912      $         CALL DLARGV( NRT, AB( 1, J1+KA ), INCA, WORK( J1 ), KA1,
00913      $                      WORK( N+J1 ), KA1 )
00914             IF( NR.GT.0 ) THEN
00915 *
00916 *              apply rotations in 1st set from the left
00917 *
00918                DO 580 L = 1, KA - 1
00919                   CALL DLARTV( NR, AB( KA1-L, J1+L ), INCA,
00920      $                         AB( KA-L, J1+L ), INCA, WORK( N+J1 ),
00921      $                         WORK( J1 ), KA1 )
00922   580          CONTINUE
00923 *
00924 *              apply rotations in 1st set from both sides to diagonal
00925 *              blocks
00926 *
00927                CALL DLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
00928      $                      AB( KA, J1 ), INCA, WORK( N+J1 ),
00929      $                      WORK( J1 ), KA1 )
00930 *
00931             END IF
00932 *
00933 *           start applying rotations in 1st set from the right
00934 *
00935             DO 590 L = KA - 1, KB - K + 1, -1
00936                NRT = ( J2+L-1 ) / KA1
00937                J1T = J2 - ( NRT-1 )*KA1
00938                IF( NRT.GT.0 )
00939      $            CALL DLARTV( NRT, AB( L, J1T ), INCA,
00940      $                         AB( L+1, J1T-1 ), INCA, WORK( N+J1T ),
00941      $                         WORK( J1T ), KA1 )
00942   590       CONTINUE
00943 *
00944             IF( WANTX ) THEN
00945 *
00946 *              post-multiply X by product of rotations in 1st set
00947 *
00948                DO 600 J = J1, J2, KA1
00949                   CALL DROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
00950      $                       WORK( N+J ), WORK( J ) )
00951   600          CONTINUE
00952             END IF
00953   610    CONTINUE
00954 *
00955          IF( UPDATE ) THEN
00956             IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
00957 *
00958 *              create nonzero element a(i+kbt-ka-1,i+kbt) outside the
00959 *              band and store it in WORK(m-kb+i+kbt)
00960 *
00961                WORK( M-KB+I+KBT ) = -BB( KB1-KBT, I+KBT )*RA1
00962             END IF
00963          END IF
00964 *
00965          DO 650 K = KB, 1, -1
00966             IF( UPDATE ) THEN
00967                J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1
00968             ELSE
00969                J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
00970             END IF
00971 *
00972 *           finish applying rotations in 2nd set from the right
00973 *
00974             DO 620 L = KB - K, 1, -1
00975                NRT = ( J2+KA+L-1 ) / KA1
00976                J1T = J2 - ( NRT-1 )*KA1
00977                IF( NRT.GT.0 )
00978      $            CALL DLARTV( NRT, AB( L, J1T+KA ), INCA,
00979      $                         AB( L+1, J1T+KA-1 ), INCA,
00980      $                         WORK( N+M-KB+J1T+KA ),
00981      $                         WORK( M-KB+J1T+KA ), KA1 )
00982   620       CONTINUE
00983             NR = ( J2+KA-1 ) / KA1
00984             J1 = J2 - ( NR-1 )*KA1
00985             DO 630 J = J1, J2, KA1
00986                WORK( M-KB+J ) = WORK( M-KB+J+KA )
00987                WORK( N+M-KB+J ) = WORK( N+M-KB+J+KA )
00988   630       CONTINUE
00989             DO 640 J = J1, J2, KA1
00990 *
00991 *              create nonzero element a(j-1,j+ka) outside the band
00992 *              and store it in WORK(m-kb+j)
00993 *
00994                WORK( M-KB+J ) = WORK( M-KB+J )*AB( 1, J+KA-1 )
00995                AB( 1, J+KA-1 ) = WORK( N+M-KB+J )*AB( 1, J+KA-1 )
00996   640       CONTINUE
00997             IF( UPDATE ) THEN
00998                IF( I+K.GT.KA1 .AND. K.LE.KBT )
00999      $            WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
01000             END IF
01001   650    CONTINUE
01002 *
01003          DO 690 K = KB, 1, -1
01004             J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
01005             NR = ( J2+KA-1 ) / KA1
01006             J1 = J2 - ( NR-1 )*KA1
01007             IF( NR.GT.0 ) THEN
01008 *
01009 *              generate rotations in 2nd set to annihilate elements
01010 *              which have been created outside the band
01011 *
01012                CALL DLARGV( NR, AB( 1, J1+KA ), INCA, WORK( M-KB+J1 ),
01013      $                      KA1, WORK( N+M-KB+J1 ), KA1 )
01014 *
01015 *              apply rotations in 2nd set from the left
01016 *
01017                DO 660 L = 1, KA - 1
01018                   CALL DLARTV( NR, AB( KA1-L, J1+L ), INCA,
01019      $                         AB( KA-L, J1+L ), INCA,
01020      $                         WORK( N+M-KB+J1 ), WORK( M-KB+J1 ), KA1 )
01021   660          CONTINUE
01022 *
01023 *              apply rotations in 2nd set from both sides to diagonal
01024 *              blocks
01025 *
01026                CALL DLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
01027      $                      AB( KA, J1 ), INCA, WORK( N+M-KB+J1 ),
01028      $                      WORK( M-KB+J1 ), KA1 )
01029 *
01030             END IF
01031 *
01032 *           start applying rotations in 2nd set from the right
01033 *
01034             DO 670 L = KA - 1, KB - K + 1, -1
01035                NRT = ( J2+L-1 ) / KA1
01036                J1T = J2 - ( NRT-1 )*KA1
01037                IF( NRT.GT.0 )
01038      $            CALL DLARTV( NRT, AB( L, J1T ), INCA,
01039      $                         AB( L+1, J1T-1 ), INCA,
01040      $                         WORK( N+M-KB+J1T ), WORK( M-KB+J1T ),
01041      $                         KA1 )
01042   670       CONTINUE
01043 *
01044             IF( WANTX ) THEN
01045 *
01046 *              post-multiply X by product of rotations in 2nd set
01047 *
01048                DO 680 J = J1, J2, KA1
01049                   CALL DROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
01050      $                       WORK( N+M-KB+J ), WORK( M-KB+J ) )
01051   680          CONTINUE
01052             END IF
01053   690    CONTINUE
01054 *
01055          DO 710 K = 1, KB - 1
01056             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
01057 *
01058 *           finish applying rotations in 1st set from the right
01059 *
01060             DO 700 L = KB - K, 1, -1
01061                NRT = ( J2+L-1 ) / KA1
01062                J1T = J2 - ( NRT-1 )*KA1
01063                IF( NRT.GT.0 )
01064      $            CALL DLARTV( NRT, AB( L, J1T ), INCA,
01065      $                         AB( L+1, J1T-1 ), INCA, WORK( N+J1T ),
01066      $                         WORK( J1T ), KA1 )
01067   700       CONTINUE
01068   710    CONTINUE
01069 *
01070          IF( KB.GT.1 ) THEN
01071             DO 720 J = 2, MIN( I+KB, M ) - 2*KA - 1
01072                WORK( N+J ) = WORK( N+J+KA )
01073                WORK( J ) = WORK( J+KA )
01074   720       CONTINUE
01075          END IF
01076 *
01077       ELSE
01078 *
01079 *        Transform A, working with the lower triangle
01080 *
01081          IF( UPDATE ) THEN
01082 *
01083 *           Form  inv(S(i))**T * A * inv(S(i))
01084 *
01085             BII = BB( 1, I )
01086             DO 730 J = I1, I
01087                AB( I-J+1, J ) = AB( I-J+1, J ) / BII
01088   730       CONTINUE
01089             DO 740 J = I, MIN( N, I+KA )
01090                AB( J-I+1, I ) = AB( J-I+1, I ) / BII
01091   740       CONTINUE
01092             DO 770 K = I + 1, I + KBT
01093                DO 750 J = K, I + KBT
01094                   AB( J-K+1, K ) = AB( J-K+1, K ) -
01095      $                             BB( J-I+1, I )*AB( K-I+1, I ) -
01096      $                             BB( K-I+1, I )*AB( J-I+1, I ) +
01097      $                             AB( 1, I )*BB( J-I+1, I )*
01098      $                             BB( K-I+1, I )
01099   750          CONTINUE
01100                DO 760 J = I + KBT + 1, MIN( N, I+KA )
01101                   AB( J-K+1, K ) = AB( J-K+1, K ) -
01102      $                             BB( K-I+1, I )*AB( J-I+1, I )
01103   760          CONTINUE
01104   770       CONTINUE
01105             DO 790 J = I1, I
01106                DO 780 K = I + 1, MIN( J+KA, I+KBT )
01107                   AB( K-J+1, J ) = AB( K-J+1, J ) -
01108      $                             BB( K-I+1, I )*AB( I-J+1, J )
01109   780          CONTINUE
01110   790       CONTINUE
01111 *
01112             IF( WANTX ) THEN
01113 *
01114 *              post-multiply X by inv(S(i))
01115 *
01116                CALL DSCAL( NX, ONE / BII, X( 1, I ), 1 )
01117                IF( KBT.GT.0 )
01118      $            CALL DGER( NX, KBT, -ONE, X( 1, I ), 1, BB( 2, I ), 1,
01119      $                       X( 1, I+1 ), LDX )
01120             END IF
01121 *
01122 *           store a(i,i1) in RA1 for use in next loop over K
01123 *
01124             RA1 = AB( I-I1+1, I1 )
01125          END IF
01126 *
01127 *        Generate and apply vectors of rotations to chase all the
01128 *        existing bulges KA positions up toward the top of the band
01129 *
01130          DO 840 K = 1, KB - 1
01131             IF( UPDATE ) THEN
01132 *
01133 *              Determine the rotations which would annihilate the bulge
01134 *              which has in theory just been created
01135 *
01136                IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
01137 *
01138 *                 generate rotation to annihilate a(i,i+k-ka-1)
01139 *
01140                   CALL DLARTG( AB( KA1-K, I+K-KA ), RA1,
01141      $                         WORK( N+I+K-KA ), WORK( I+K-KA ), RA )
01142 *
01143 *                 create nonzero element a(i+k,i+k-ka-1) outside the
01144 *                 band and store it in WORK(m-kb+i+k)
01145 *
01146                   T = -BB( K+1, I )*RA1
01147                   WORK( M-KB+I+K ) = WORK( N+I+K-KA )*T -
01148      $                               WORK( I+K-KA )*AB( KA1, I+K-KA )
01149                   AB( KA1, I+K-KA ) = WORK( I+K-KA )*T +
01150      $                                WORK( N+I+K-KA )*AB( KA1, I+K-KA )
01151                   RA1 = RA
01152                END IF
01153             END IF
01154             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
01155             NR = ( J2+KA-1 ) / KA1
01156             J1 = J2 - ( NR-1 )*KA1
01157             IF( UPDATE ) THEN
01158                J2T = MIN( J2, I-2*KA+K-1 )
01159             ELSE
01160                J2T = J2
01161             END IF
01162             NRT = ( J2T+KA-1 ) / KA1
01163             DO 800 J = J1, J2T, KA1
01164 *
01165 *              create nonzero element a(j+ka,j-1) outside the band
01166 *              and store it in WORK(j)
01167 *
01168                WORK( J ) = WORK( J )*AB( KA1, J-1 )
01169                AB( KA1, J-1 ) = WORK( N+J )*AB( KA1, J-1 )
01170   800       CONTINUE
01171 *
01172 *           generate rotations in 1st set to annihilate elements which
01173 *           have been created outside the band
01174 *
01175             IF( NRT.GT.0 )
01176      $         CALL DLARGV( NRT, AB( KA1, J1 ), INCA, WORK( J1 ), KA1,
01177      $                      WORK( N+J1 ), KA1 )
01178             IF( NR.GT.0 ) THEN
01179 *
01180 *              apply rotations in 1st set from the right
01181 *
01182                DO 810 L = 1, KA - 1
01183                   CALL DLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
01184      $                         INCA, WORK( N+J1 ), WORK( J1 ), KA1 )
01185   810          CONTINUE
01186 *
01187 *              apply rotations in 1st set from both sides to diagonal
01188 *              blocks
01189 *
01190                CALL DLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
01191      $                      AB( 2, J1-1 ), INCA, WORK( N+J1 ),
01192      $                      WORK( J1 ), KA1 )
01193 *
01194             END IF
01195 *
01196 *           start applying rotations in 1st set from the left
01197 *
01198             DO 820 L = KA - 1, KB - K + 1, -1
01199                NRT = ( J2+L-1 ) / KA1
01200                J1T = J2 - ( NRT-1 )*KA1
01201                IF( NRT.GT.0 )
01202      $            CALL DLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
01203      $                         AB( KA1-L, J1T-KA1+L ), INCA,
01204      $                         WORK( N+J1T ), WORK( J1T ), KA1 )
01205   820       CONTINUE
01206 *
01207             IF( WANTX ) THEN
01208 *
01209 *              post-multiply X by product of rotations in 1st set
01210 *
01211                DO 830 J = J1, J2, KA1
01212                   CALL DROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
01213      $                       WORK( N+J ), WORK( J ) )
01214   830          CONTINUE
01215             END IF
01216   840    CONTINUE
01217 *
01218          IF( UPDATE ) THEN
01219             IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
01220 *
01221 *              create nonzero element a(i+kbt,i+kbt-ka-1) outside the
01222 *              band and store it in WORK(m-kb+i+kbt)
01223 *
01224                WORK( M-KB+I+KBT ) = -BB( KBT+1, I )*RA1
01225             END IF
01226          END IF
01227 *
01228          DO 880 K = KB, 1, -1
01229             IF( UPDATE ) THEN
01230                J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1
01231             ELSE
01232                J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
01233             END IF
01234 *
01235 *           finish applying rotations in 2nd set from the left
01236 *
01237             DO 850 L = KB - K, 1, -1
01238                NRT = ( J2+KA+L-1 ) / KA1
01239                J1T = J2 - ( NRT-1 )*KA1
01240                IF( NRT.GT.0 )
01241      $            CALL DLARTV( NRT, AB( KA1-L+1, J1T+L-1 ), INCA,
01242      $                         AB( KA1-L, J1T+L-1 ), INCA,
01243      $                         WORK( N+M-KB+J1T+KA ),
01244      $                         WORK( M-KB+J1T+KA ), KA1 )
01245   850       CONTINUE
01246             NR = ( J2+KA-1 ) / KA1
01247             J1 = J2 - ( NR-1 )*KA1
01248             DO 860 J = J1, J2, KA1
01249                WORK( M-KB+J ) = WORK( M-KB+J+KA )
01250                WORK( N+M-KB+J ) = WORK( N+M-KB+J+KA )
01251   860       CONTINUE
01252             DO 870 J = J1, J2, KA1
01253 *
01254 *              create nonzero element a(j+ka,j-1) outside the band
01255 *              and store it in WORK(m-kb+j)
01256 *
01257                WORK( M-KB+J ) = WORK( M-KB+J )*AB( KA1, J-1 )
01258                AB( KA1, J-1 ) = WORK( N+M-KB+J )*AB( KA1, J-1 )
01259   870       CONTINUE
01260             IF( UPDATE ) THEN
01261                IF( I+K.GT.KA1 .AND. K.LE.KBT )
01262      $            WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
01263             END IF
01264   880    CONTINUE
01265 *
01266          DO 920 K = KB, 1, -1
01267             J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
01268             NR = ( J2+KA-1 ) / KA1
01269             J1 = J2 - ( NR-1 )*KA1
01270             IF( NR.GT.0 ) THEN
01271 *
01272 *              generate rotations in 2nd set to annihilate elements
01273 *              which have been created outside the band
01274 *
01275                CALL DLARGV( NR, AB( KA1, J1 ), INCA, WORK( M-KB+J1 ),
01276      $                      KA1, WORK( N+M-KB+J1 ), KA1 )
01277 *
01278 *              apply rotations in 2nd set from the right
01279 *
01280                DO 890 L = 1, KA - 1
01281                   CALL DLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
01282      $                         INCA, WORK( N+M-KB+J1 ), WORK( M-KB+J1 ),
01283      $                         KA1 )
01284   890          CONTINUE
01285 *
01286 *              apply rotations in 2nd set from both sides to diagonal
01287 *              blocks
01288 *
01289                CALL DLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
01290      $                      AB( 2, J1-1 ), INCA, WORK( N+M-KB+J1 ),
01291      $                      WORK( M-KB+J1 ), KA1 )
01292 *
01293             END IF
01294 *
01295 *           start applying rotations in 2nd set from the left
01296 *
01297             DO 900 L = KA - 1, KB - K + 1, -1
01298                NRT = ( J2+L-1 ) / KA1
01299                J1T = J2 - ( NRT-1 )*KA1
01300                IF( NRT.GT.0 )
01301      $            CALL DLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
01302      $                         AB( KA1-L, J1T-KA1+L ), INCA,
01303      $                         WORK( N+M-KB+J1T ), WORK( M-KB+J1T ),
01304      $                         KA1 )
01305   900       CONTINUE
01306 *
01307             IF( WANTX ) THEN
01308 *
01309 *              post-multiply X by product of rotations in 2nd set
01310 *
01311                DO 910 J = J1, J2, KA1
01312                   CALL DROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
01313      $                       WORK( N+M-KB+J ), WORK( M-KB+J ) )
01314   910          CONTINUE
01315             END IF
01316   920    CONTINUE
01317 *
01318          DO 940 K = 1, KB - 1
01319             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
01320 *
01321 *           finish applying rotations in 1st set from the left
01322 *
01323             DO 930 L = KB - K, 1, -1
01324                NRT = ( J2+L-1 ) / KA1
01325                J1T = J2 - ( NRT-1 )*KA1
01326                IF( NRT.GT.0 )
01327      $            CALL DLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
01328      $                         AB( KA1-L, J1T-KA1+L ), INCA,
01329      $                         WORK( N+J1T ), WORK( J1T ), KA1 )
01330   930       CONTINUE
01331   940    CONTINUE
01332 *
01333          IF( KB.GT.1 ) THEN
01334             DO 950 J = 2, MIN( I+KB, M ) - 2*KA - 1
01335                WORK( N+J ) = WORK( N+J+KA )
01336                WORK( J ) = WORK( J+KA )
01337   950       CONTINUE
01338          END IF
01339 *
01340       END IF
01341 *
01342       GO TO 490
01343 *
01344 *     End of DSBGST
01345 *
01346       END
 All Files Functions