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