LAPACK 3.3.1
Linear Algebra PACKage

chbgst.f

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