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