LAPACK 3.3.0

zhbgst.f

Go to the documentation of this file.
00001       SUBROUTINE ZHBGST( 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       DOUBLE PRECISION   RWORK( * )
00015       COMPLEX*16         AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
00016      $                   X( LDX, * )
00017 *     ..
00018 *
00019 *  Purpose
00020 *  =======
00021 *
00022 *  ZHBGST 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 ZPBSTF, 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*16 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*16 array, dimension (LDBB,N)
00068 *          The banded factor S from the split Cholesky factorization of
00069 *          B, as returned by ZPBSTF, 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*16 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*16 array, dimension (N)
00084 *
00085 *  RWORK   (workspace) DOUBLE PRECISION 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*16         CZERO, CONE
00095       DOUBLE PRECISION   ONE
00096       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
00097      $                   CONE = ( 1.0D+0, 0.0D+0 ), ONE = 1.0D+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       DOUBLE PRECISION   BII
00104       COMPLEX*16         RA, RA1, T
00105 *     ..
00106 *     .. External Functions ..
00107       LOGICAL            LSAME
00108       EXTERNAL           LSAME
00109 *     ..
00110 *     .. External Subroutines ..
00111       EXTERNAL           XERBLA, ZDSCAL, ZGERC, ZGERU, ZLACGV, ZLAR2V,
00112      $                   ZLARGV, ZLARTG, ZLARTV, ZLASET, ZROT
00113 *     ..
00114 *     .. Intrinsic Functions ..
00115       INTRINSIC          DBLE, DCONJG, MAX, MIN
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( 'ZHBGST', -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 ZLASET( '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 ZPBSTF. 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 = DBLE( BB( KB1, I ) )
00257             AB( KA1, I ) = ( DBLE( 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      $                               DCONJG( AB( K-I+KA1, I ) ) -
00269      $                               DCONJG( BB( K-I+KB1, I ) )*
00270      $                               AB( J-I+KA1, I ) +
00271      $                               DBLE( AB( KA1, I ) )*
00272      $                               BB( J-I+KB1, I )*
00273      $                               DCONJG( 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      $                               DCONJG( 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 ZDSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
00293                IF( KBT.GT.0 )
00294      $            CALL ZGERC( 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 ZLARTG( 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      $                          DCONJG( 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 ZLARGV( 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 ZLARTV( 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 ZLAR2V( 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 ZLACGV( 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 ZLARTV( 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 ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
00393      $                       RWORK( J-M ), DCONJG( 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 ZLARTV( 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 ZLARGV( 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 ZLARTV( 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 ZLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ),
00468      $                      AB( KA, J2+1 ), INCA, RWORK( J2 ),
00469      $                      WORK( J2 ), KA1 )
00470 *
00471                CALL ZLACGV( 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 ZLARTV( 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 ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
00490      $                       RWORK( J ), DCONJG( 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 ZLARTV( 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 = DBLE( BB( 1, I ) )
00525             AB( 1, I ) = ( DBLE( 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 )*DCONJG( AB( I-K+1,
00536      $                             K ) ) - DCONJG( BB( I-K+1, K ) )*
00537      $                             AB( I-J+1, J ) + DBLE( AB( 1, I ) )*
00538      $                             BB( I-J+1, J )*DCONJG( 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      $                             DCONJG( 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 ZDSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
00559                IF( KBT.GT.0 )
00560      $            CALL ZGERU( 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 ZLARTG( 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      $                          DCONJG( WORK( I-K+KA-M ) )*
00593      $                          AB( KA1, I-K )
00594                   AB( KA1, I-K ) = WORK( I-K+KA-M )*T +
00595      $                             RWORK( I-K+KA-M )*AB( KA1, I-K )
00596                   RA1 = RA
00597                END IF
00598             END IF
00599             J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
00600             NR = ( N-J2+KA ) / KA1
00601             J1 = J2 + ( NR-1 )*KA1
00602             IF( UPDATE ) THEN
00603                J2T = MAX( J2, I+2*KA-K+1 )
00604             ELSE
00605                J2T = J2
00606             END IF
00607             NRT = ( N-J2T+KA ) / KA1
00608             DO 320 J = J2T, J1, KA1
00609 *
00610 *              create nonzero element a(j+1,j-ka) outside the band
00611 *              and store it in WORK(j-m)
00612 *
00613                WORK( J-M ) = WORK( J-M )*AB( KA1, J-KA+1 )
00614                AB( KA1, J-KA+1 ) = RWORK( J-M )*AB( KA1, J-KA+1 )
00615   320       CONTINUE
00616 *
00617 *           generate rotations in 1st set to annihilate elements which
00618 *           have been created outside the band
00619 *
00620             IF( NRT.GT.0 )
00621      $         CALL ZLARGV( NRT, AB( KA1, J2T-KA ), INCA, WORK( J2T-M ),
00622      $                      KA1, RWORK( J2T-M ), KA1 )
00623             IF( NR.GT.0 ) THEN
00624 *
00625 *              apply rotations in 1st set from the left
00626 *
00627                DO 330 L = 1, KA - 1
00628                   CALL ZLARTV( NR, AB( L+1, J2-L ), INCA,
00629      $                         AB( L+2, J2-L ), INCA, RWORK( J2-M ),
00630      $                         WORK( J2-M ), KA1 )
00631   330          CONTINUE
00632 *
00633 *              apply rotations in 1st set from both sides to diagonal
00634 *              blocks
00635 *
00636                CALL ZLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
00637      $                      INCA, RWORK( J2-M ), WORK( J2-M ), KA1 )
00638 *
00639                CALL ZLACGV( NR, WORK( J2-M ), KA1 )
00640             END IF
00641 *
00642 *           start applying rotations in 1st set from the right
00643 *
00644             DO 340 L = KA - 1, KB - K + 1, -1
00645                NRT = ( N-J2+L ) / KA1
00646                IF( NRT.GT.0 )
00647      $            CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
00648      $                         AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ),
00649      $                         WORK( J2-M ), KA1 )
00650   340       CONTINUE
00651 *
00652             IF( WANTX ) THEN
00653 *
00654 *              post-multiply X by product of rotations in 1st set
00655 *
00656                DO 350 J = J2, J1, KA1
00657                   CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
00658      $                       RWORK( J-M ), WORK( J-M ) )
00659   350          CONTINUE
00660             END IF
00661   360    CONTINUE
00662 *
00663          IF( UPDATE ) THEN
00664             IF( I2.LE.N .AND. KBT.GT.0 ) THEN
00665 *
00666 *              create nonzero element a(i-kbt+ka+1,i-kbt) outside the
00667 *              band and store it in WORK(i-kbt)
00668 *
00669                WORK( I-KBT ) = -BB( KBT+1, I-KBT )*RA1
00670             END IF
00671          END IF
00672 *
00673          DO 400 K = KB, 1, -1
00674             IF( UPDATE ) THEN
00675                J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1
00676             ELSE
00677                J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
00678             END IF
00679 *
00680 *           finish applying rotations in 2nd set from the right
00681 *
00682             DO 370 L = KB - K, 1, -1
00683                NRT = ( N-J2+KA+L ) / KA1
00684                IF( NRT.GT.0 )
00685      $            CALL ZLARTV( NRT, AB( KA1-L+1, J2-KA ), INCA,
00686      $                         AB( KA1-L, J2-KA+1 ), INCA,
00687      $                         RWORK( J2-KA ), WORK( J2-KA ), KA1 )
00688   370       CONTINUE
00689             NR = ( N-J2+KA ) / KA1
00690             J1 = J2 + ( NR-1 )*KA1
00691             DO 380 J = J1, J2, -KA1
00692                WORK( J ) = WORK( J-KA )
00693                RWORK( J ) = RWORK( J-KA )
00694   380       CONTINUE
00695             DO 390 J = J2, J1, KA1
00696 *
00697 *              create nonzero element a(j+1,j-ka) outside the band
00698 *              and store it in WORK(j)
00699 *
00700                WORK( J ) = WORK( J )*AB( KA1, J-KA+1 )
00701                AB( KA1, J-KA+1 ) = RWORK( J )*AB( KA1, J-KA+1 )
00702   390       CONTINUE
00703             IF( UPDATE ) THEN
00704                IF( I-K.LT.N-KA .AND. K.LE.KBT )
00705      $            WORK( I-K+KA ) = WORK( I-K )
00706             END IF
00707   400    CONTINUE
00708 *
00709          DO 440 K = KB, 1, -1
00710             J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
00711             NR = ( N-J2+KA ) / KA1
00712             J1 = J2 + ( NR-1 )*KA1
00713             IF( NR.GT.0 ) THEN
00714 *
00715 *              generate rotations in 2nd set to annihilate elements
00716 *              which have been created outside the band
00717 *
00718                CALL ZLARGV( NR, AB( KA1, J2-KA ), INCA, WORK( J2 ), KA1,
00719      $                      RWORK( J2 ), KA1 )
00720 *
00721 *              apply rotations in 2nd set from the left
00722 *
00723                DO 410 L = 1, KA - 1
00724                   CALL ZLARTV( NR, AB( L+1, J2-L ), INCA,
00725      $                         AB( L+2, J2-L ), INCA, RWORK( J2 ),
00726      $                         WORK( J2 ), KA1 )
00727   410          CONTINUE
00728 *
00729 *              apply rotations in 2nd set from both sides to diagonal
00730 *              blocks
00731 *
00732                CALL ZLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
00733      $                      INCA, RWORK( J2 ), WORK( J2 ), KA1 )
00734 *
00735                CALL ZLACGV( NR, WORK( J2 ), KA1 )
00736             END IF
00737 *
00738 *           start applying rotations in 2nd set from the right
00739 *
00740             DO 420 L = KA - 1, KB - K + 1, -1
00741                NRT = ( N-J2+L ) / KA1
00742                IF( NRT.GT.0 )
00743      $            CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
00744      $                         AB( KA1-L, J2+1 ), INCA, RWORK( J2 ),
00745      $                         WORK( J2 ), KA1 )
00746   420       CONTINUE
00747 *
00748             IF( WANTX ) THEN
00749 *
00750 *              post-multiply X by product of rotations in 2nd set
00751 *
00752                DO 430 J = J2, J1, KA1
00753                   CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
00754      $                       RWORK( J ), WORK( J ) )
00755   430          CONTINUE
00756             END IF
00757   440    CONTINUE
00758 *
00759          DO 460 K = 1, KB - 1
00760             J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
00761 *
00762 *           finish applying rotations in 1st set from the right
00763 *
00764             DO 450 L = KB - K, 1, -1
00765                NRT = ( N-J2+L ) / KA1
00766                IF( NRT.GT.0 )
00767      $            CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
00768      $                         AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ),
00769      $                         WORK( J2-M ), KA1 )
00770   450       CONTINUE
00771   460    CONTINUE
00772 *
00773          IF( KB.GT.1 ) THEN
00774             DO 470 J = N - 1, J2 + KA, -1
00775                RWORK( J-M ) = RWORK( J-KA-M )
00776                WORK( J-M ) = WORK( J-KA-M )
00777   470       CONTINUE
00778          END IF
00779 *
00780       END IF
00781 *
00782       GO TO 10
00783 *
00784   480 CONTINUE
00785 *
00786 *     **************************** Phase 2 *****************************
00787 *
00788 *     The logical structure of this phase is:
00789 *
00790 *     UPDATE = .TRUE.
00791 *     DO I = 1, M
00792 *        use S(i) to update A and create a new bulge
00793 *        apply rotations to push all bulges KA positions upward
00794 *     END DO
00795 *     UPDATE = .FALSE.
00796 *     DO I = M - KA - 1, 2, -1
00797 *        apply rotations to push all bulges KA positions upward
00798 *     END DO
00799 *
00800 *     To avoid duplicating code, the two loops are merged.
00801 *
00802       UPDATE = .TRUE.
00803       I = 0
00804   490 CONTINUE
00805       IF( UPDATE ) THEN
00806          I = I + 1
00807          KBT = MIN( KB, M-I )
00808          I0 = I + 1
00809          I1 = MAX( 1, I-KA )
00810          I2 = I + KBT - KA1
00811          IF( I.GT.M ) THEN
00812             UPDATE = .FALSE.
00813             I = I - 1
00814             I0 = M + 1
00815             IF( KA.EQ.0 )
00816      $         RETURN
00817             GO TO 490
00818          END IF
00819       ELSE
00820          I = I - KA
00821          IF( I.LT.2 )
00822      $      RETURN
00823       END IF
00824 *
00825       IF( I.LT.M-KBT ) THEN
00826          NX = M
00827       ELSE
00828          NX = N
00829       END IF
00830 *
00831       IF( UPPER ) THEN
00832 *
00833 *        Transform A, working with the upper triangle
00834 *
00835          IF( UPDATE ) THEN
00836 *
00837 *           Form  inv(S(i))**H * A * inv(S(i))
00838 *
00839             BII = DBLE( BB( KB1, I ) )
00840             AB( KA1, I ) = ( DBLE( AB( KA1, I ) ) / BII ) / BII
00841             DO 500 J = I1, I - 1
00842                AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII
00843   500       CONTINUE
00844             DO 510 J = I + 1, MIN( N, I+KA )
00845                AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII
00846   510       CONTINUE
00847             DO 540 K = I + 1, I + KBT
00848                DO 520 J = K, I + KBT
00849                   AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
00850      $                               BB( I-J+KB1, J )*
00851      $                               DCONJG( AB( I-K+KA1, K ) ) -
00852      $                               DCONJG( BB( I-K+KB1, K ) )*
00853      $                               AB( I-J+KA1, J ) +
00854      $                               DBLE( AB( KA1, I ) )*
00855      $                               BB( I-J+KB1, J )*
00856      $                               DCONJG( BB( I-K+KB1, K ) )
00857   520          CONTINUE
00858                DO 530 J = I + KBT + 1, MIN( N, I+KA )
00859                   AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
00860      $                               DCONJG( BB( I-K+KB1, K ) )*
00861      $                               AB( I-J+KA1, J )
00862   530          CONTINUE
00863   540       CONTINUE
00864             DO 560 J = I1, I
00865                DO 550 K = I + 1, MIN( J+KA, I+KBT )
00866                   AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
00867      $                               BB( I-K+KB1, K )*AB( J-I+KA1, I )
00868   550          CONTINUE
00869   560       CONTINUE
00870 *
00871             IF( WANTX ) THEN
00872 *
00873 *              post-multiply X by inv(S(i))
00874 *
00875                CALL ZDSCAL( NX, ONE / BII, X( 1, I ), 1 )
00876                IF( KBT.GT.0 )
00877      $            CALL ZGERU( NX, KBT, -CONE, X( 1, I ), 1,
00878      $                        BB( KB, I+1 ), LDBB-1, X( 1, I+1 ), LDX )
00879             END IF
00880 *
00881 *           store a(i1,i) in RA1 for use in next loop over K
00882 *
00883             RA1 = AB( I1-I+KA1, I )
00884          END IF
00885 *
00886 *        Generate and apply vectors of rotations to chase all the
00887 *        existing bulges KA positions up toward the top of the band
00888 *
00889          DO 610 K = 1, KB - 1
00890             IF( UPDATE ) THEN
00891 *
00892 *              Determine the rotations which would annihilate the bulge
00893 *              which has in theory just been created
00894 *
00895                IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
00896 *
00897 *                 generate rotation to annihilate a(i+k-ka-1,i)
00898 *
00899                   CALL ZLARTG( AB( K+1, I ), RA1, RWORK( I+K-KA ),
00900      $                         WORK( I+K-KA ), RA )
00901 *
00902 *                 create nonzero element a(i+k-ka-1,i+k) outside the
00903 *                 band and store it in WORK(m-kb+i+k)
00904 *
00905                   T = -BB( KB1-K, I+K )*RA1
00906                   WORK( M-KB+I+K ) = RWORK( I+K-KA )*T -
00907      $                               DCONJG( WORK( I+K-KA ) )*
00908      $                               AB( 1, I+K )
00909                   AB( 1, I+K ) = WORK( I+K-KA )*T +
00910      $                           RWORK( I+K-KA )*AB( 1, I+K )
00911                   RA1 = RA
00912                END IF
00913             END IF
00914             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
00915             NR = ( J2+KA-1 ) / KA1
00916             J1 = J2 - ( NR-1 )*KA1
00917             IF( UPDATE ) THEN
00918                J2T = MIN( J2, I-2*KA+K-1 )
00919             ELSE
00920                J2T = J2
00921             END IF
00922             NRT = ( J2T+KA-1 ) / KA1
00923             DO 570 J = J1, J2T, KA1
00924 *
00925 *              create nonzero element a(j-1,j+ka) outside the band
00926 *              and store it in WORK(j)
00927 *
00928                WORK( J ) = WORK( J )*AB( 1, J+KA-1 )
00929                AB( 1, J+KA-1 ) = RWORK( J )*AB( 1, J+KA-1 )
00930   570       CONTINUE
00931 *
00932 *           generate rotations in 1st set to annihilate elements which
00933 *           have been created outside the band
00934 *
00935             IF( NRT.GT.0 )
00936      $         CALL ZLARGV( NRT, AB( 1, J1+KA ), INCA, WORK( J1 ), KA1,
00937      $                      RWORK( J1 ), KA1 )
00938             IF( NR.GT.0 ) THEN
00939 *
00940 *              apply rotations in 1st set from the left
00941 *
00942                DO 580 L = 1, KA - 1
00943                   CALL ZLARTV( NR, AB( KA1-L, J1+L ), INCA,
00944      $                         AB( KA-L, J1+L ), INCA, RWORK( J1 ),
00945      $                         WORK( J1 ), KA1 )
00946   580          CONTINUE
00947 *
00948 *              apply rotations in 1st set from both sides to diagonal
00949 *              blocks
00950 *
00951                CALL ZLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
00952      $                      AB( KA, J1 ), INCA, RWORK( J1 ), WORK( J1 ),
00953      $                      KA1 )
00954 *
00955                CALL ZLACGV( NR, WORK( J1 ), KA1 )
00956             END IF
00957 *
00958 *           start applying rotations in 1st set from the right
00959 *
00960             DO 590 L = KA - 1, KB - K + 1, -1
00961                NRT = ( J2+L-1 ) / KA1
00962                J1T = J2 - ( NRT-1 )*KA1
00963                IF( NRT.GT.0 )
00964      $            CALL ZLARTV( NRT, AB( L, J1T ), INCA,
00965      $                         AB( L+1, J1T-1 ), INCA, RWORK( J1T ),
00966      $                         WORK( J1T ), KA1 )
00967   590       CONTINUE
00968 *
00969             IF( WANTX ) THEN
00970 *
00971 *              post-multiply X by product of rotations in 1st set
00972 *
00973                DO 600 J = J1, J2, KA1
00974                   CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
00975      $                       RWORK( J ), WORK( J ) )
00976   600          CONTINUE
00977             END IF
00978   610    CONTINUE
00979 *
00980          IF( UPDATE ) THEN
00981             IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
00982 *
00983 *              create nonzero element a(i+kbt-ka-1,i+kbt) outside the
00984 *              band and store it in WORK(m-kb+i+kbt)
00985 *
00986                WORK( M-KB+I+KBT ) = -BB( KB1-KBT, I+KBT )*RA1
00987             END IF
00988          END IF
00989 *
00990          DO 650 K = KB, 1, -1
00991             IF( UPDATE ) THEN
00992                J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1
00993             ELSE
00994                J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
00995             END IF
00996 *
00997 *           finish applying rotations in 2nd set from the right
00998 *
00999             DO 620 L = KB - K, 1, -1
01000                NRT = ( J2+KA+L-1 ) / KA1
01001                J1T = J2 - ( NRT-1 )*KA1
01002                IF( NRT.GT.0 )
01003      $            CALL ZLARTV( NRT, AB( L, J1T+KA ), INCA,
01004      $                         AB( L+1, J1T+KA-1 ), INCA,
01005      $                         RWORK( M-KB+J1T+KA ),
01006      $                         WORK( M-KB+J1T+KA ), KA1 )
01007   620       CONTINUE
01008             NR = ( J2+KA-1 ) / KA1
01009             J1 = J2 - ( NR-1 )*KA1
01010             DO 630 J = J1, J2, KA1
01011                WORK( M-KB+J ) = WORK( M-KB+J+KA )
01012                RWORK( M-KB+J ) = RWORK( M-KB+J+KA )
01013   630       CONTINUE
01014             DO 640 J = J1, J2, KA1
01015 *
01016 *              create nonzero element a(j-1,j+ka) outside the band
01017 *              and store it in WORK(m-kb+j)
01018 *
01019                WORK( M-KB+J ) = WORK( M-KB+J )*AB( 1, J+KA-1 )
01020                AB( 1, J+KA-1 ) = RWORK( M-KB+J )*AB( 1, J+KA-1 )
01021   640       CONTINUE
01022             IF( UPDATE ) THEN
01023                IF( I+K.GT.KA1 .AND. K.LE.KBT )
01024      $            WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
01025             END IF
01026   650    CONTINUE
01027 *
01028          DO 690 K = KB, 1, -1
01029             J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
01030             NR = ( J2+KA-1 ) / KA1
01031             J1 = J2 - ( NR-1 )*KA1
01032             IF( NR.GT.0 ) THEN
01033 *
01034 *              generate rotations in 2nd set to annihilate elements
01035 *              which have been created outside the band
01036 *
01037                CALL ZLARGV( NR, AB( 1, J1+KA ), INCA, WORK( M-KB+J1 ),
01038      $                      KA1, RWORK( M-KB+J1 ), KA1 )
01039 *
01040 *              apply rotations in 2nd set from the left
01041 *
01042                DO 660 L = 1, KA - 1
01043                   CALL ZLARTV( NR, AB( KA1-L, J1+L ), INCA,
01044      $                         AB( KA-L, J1+L ), INCA, RWORK( M-KB+J1 ),
01045      $                         WORK( M-KB+J1 ), KA1 )
01046   660          CONTINUE
01047 *
01048 *              apply rotations in 2nd set from both sides to diagonal
01049 *              blocks
01050 *
01051                CALL ZLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
01052      $                      AB( KA, J1 ), INCA, RWORK( M-KB+J1 ),
01053      $                      WORK( M-KB+J1 ), KA1 )
01054 *
01055                CALL ZLACGV( NR, WORK( M-KB+J1 ), KA1 )
01056             END IF
01057 *
01058 *           start applying rotations in 2nd set from the right
01059 *
01060             DO 670 L = KA - 1, KB - K + 1, -1
01061                NRT = ( J2+L-1 ) / KA1
01062                J1T = J2 - ( NRT-1 )*KA1
01063                IF( NRT.GT.0 )
01064      $            CALL ZLARTV( NRT, AB( L, J1T ), INCA,
01065      $                         AB( L+1, J1T-1 ), INCA,
01066      $                         RWORK( M-KB+J1T ), WORK( M-KB+J1T ),
01067      $                         KA1 )
01068   670       CONTINUE
01069 *
01070             IF( WANTX ) THEN
01071 *
01072 *              post-multiply X by product of rotations in 2nd set
01073 *
01074                DO 680 J = J1, J2, KA1
01075                   CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
01076      $                       RWORK( M-KB+J ), WORK( M-KB+J ) )
01077   680          CONTINUE
01078             END IF
01079   690    CONTINUE
01080 *
01081          DO 710 K = 1, KB - 1
01082             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
01083 *
01084 *           finish applying rotations in 1st set from the right
01085 *
01086             DO 700 L = KB - K, 1, -1
01087                NRT = ( J2+L-1 ) / KA1
01088                J1T = J2 - ( NRT-1 )*KA1
01089                IF( NRT.GT.0 )
01090      $            CALL ZLARTV( NRT, AB( L, J1T ), INCA,
01091      $                         AB( L+1, J1T-1 ), INCA, RWORK( J1T ),
01092      $                         WORK( J1T ), KA1 )
01093   700       CONTINUE
01094   710    CONTINUE
01095 *
01096          IF( KB.GT.1 ) THEN
01097             DO 720 J = 2, I2 - KA
01098                RWORK( J ) = RWORK( J+KA )
01099                WORK( J ) = WORK( J+KA )
01100   720       CONTINUE
01101          END IF
01102 *
01103       ELSE
01104 *
01105 *        Transform A, working with the lower triangle
01106 *
01107          IF( UPDATE ) THEN
01108 *
01109 *           Form  inv(S(i))**H * A * inv(S(i))
01110 *
01111             BII = DBLE( BB( 1, I ) )
01112             AB( 1, I ) = ( DBLE( AB( 1, I ) ) / BII ) / BII
01113             DO 730 J = I1, I - 1
01114                AB( I-J+1, J ) = AB( I-J+1, J ) / BII
01115   730       CONTINUE
01116             DO 740 J = I + 1, MIN( N, I+KA )
01117                AB( J-I+1, I ) = AB( J-I+1, I ) / BII
01118   740       CONTINUE
01119             DO 770 K = I + 1, I + KBT
01120                DO 750 J = K, I + KBT
01121                   AB( J-K+1, K ) = AB( J-K+1, K ) -
01122      $                             BB( J-I+1, I )*DCONJG( AB( K-I+1,
01123      $                             I ) ) - DCONJG( BB( K-I+1, I ) )*
01124      $                             AB( J-I+1, I ) + DBLE( AB( 1, I ) )*
01125      $                             BB( J-I+1, I )*DCONJG( BB( K-I+1,
01126      $                             I ) )
01127   750          CONTINUE
01128                DO 760 J = I + KBT + 1, MIN( N, I+KA )
01129                   AB( J-K+1, K ) = AB( J-K+1, K ) -
01130      $                             DCONJG( BB( K-I+1, I ) )*
01131      $                             AB( J-I+1, I )
01132   760          CONTINUE
01133   770       CONTINUE
01134             DO 790 J = I1, I
01135                DO 780 K = I + 1, MIN( J+KA, I+KBT )
01136                   AB( K-J+1, J ) = AB( K-J+1, J ) -
01137      $                             BB( K-I+1, I )*AB( I-J+1, J )
01138   780          CONTINUE
01139   790       CONTINUE
01140 *
01141             IF( WANTX ) THEN
01142 *
01143 *              post-multiply X by inv(S(i))
01144 *
01145                CALL ZDSCAL( NX, ONE / BII, X( 1, I ), 1 )
01146                IF( KBT.GT.0 )
01147      $            CALL ZGERC( NX, KBT, -CONE, X( 1, I ), 1, BB( 2, I ),
01148      $                        1, X( 1, I+1 ), LDX )
01149             END IF
01150 *
01151 *           store a(i,i1) in RA1 for use in next loop over K
01152 *
01153             RA1 = AB( I-I1+1, I1 )
01154          END IF
01155 *
01156 *        Generate and apply vectors of rotations to chase all the
01157 *        existing bulges KA positions up toward the top of the band
01158 *
01159          DO 840 K = 1, KB - 1
01160             IF( UPDATE ) THEN
01161 *
01162 *              Determine the rotations which would annihilate the bulge
01163 *              which has in theory just been created
01164 *
01165                IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
01166 *
01167 *                 generate rotation to annihilate a(i,i+k-ka-1)
01168 *
01169                   CALL ZLARTG( AB( KA1-K, I+K-KA ), RA1,
01170      $                         RWORK( I+K-KA ), WORK( I+K-KA ), RA )
01171 *
01172 *                 create nonzero element a(i+k,i+k-ka-1) outside the
01173 *                 band and store it in WORK(m-kb+i+k)
01174 *
01175                   T = -BB( K+1, I )*RA1
01176                   WORK( M-KB+I+K ) = RWORK( I+K-KA )*T -
01177      $                               DCONJG( WORK( I+K-KA ) )*
01178      $                               AB( KA1, I+K-KA )
01179                   AB( KA1, I+K-KA ) = WORK( I+K-KA )*T +
01180      $                                RWORK( I+K-KA )*AB( KA1, I+K-KA )
01181                   RA1 = RA
01182                END IF
01183             END IF
01184             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
01185             NR = ( J2+KA-1 ) / KA1
01186             J1 = J2 - ( NR-1 )*KA1
01187             IF( UPDATE ) THEN
01188                J2T = MIN( J2, I-2*KA+K-1 )
01189             ELSE
01190                J2T = J2
01191             END IF
01192             NRT = ( J2T+KA-1 ) / KA1
01193             DO 800 J = J1, J2T, KA1
01194 *
01195 *              create nonzero element a(j+ka,j-1) outside the band
01196 *              and store it in WORK(j)
01197 *
01198                WORK( J ) = WORK( J )*AB( KA1, J-1 )
01199                AB( KA1, J-1 ) = RWORK( J )*AB( KA1, J-1 )
01200   800       CONTINUE
01201 *
01202 *           generate rotations in 1st set to annihilate elements which
01203 *           have been created outside the band
01204 *
01205             IF( NRT.GT.0 )
01206      $         CALL ZLARGV( NRT, AB( KA1, J1 ), INCA, WORK( J1 ), KA1,
01207      $                      RWORK( J1 ), KA1 )
01208             IF( NR.GT.0 ) THEN
01209 *
01210 *              apply rotations in 1st set from the right
01211 *
01212                DO 810 L = 1, KA - 1
01213                   CALL ZLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
01214      $                         INCA, RWORK( J1 ), WORK( J1 ), KA1 )
01215   810          CONTINUE
01216 *
01217 *              apply rotations in 1st set from both sides to diagonal
01218 *              blocks
01219 *
01220                CALL ZLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
01221      $                      AB( 2, J1-1 ), INCA, RWORK( J1 ),
01222      $                      WORK( J1 ), KA1 )
01223 *
01224                CALL ZLACGV( NR, WORK( J1 ), KA1 )
01225             END IF
01226 *
01227 *           start applying rotations in 1st set from the left
01228 *
01229             DO 820 L = KA - 1, KB - K + 1, -1
01230                NRT = ( J2+L-1 ) / KA1
01231                J1T = J2 - ( NRT-1 )*KA1
01232                IF( NRT.GT.0 )
01233      $            CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
01234      $                         AB( KA1-L, J1T-KA1+L ), INCA,
01235      $                         RWORK( J1T ), WORK( J1T ), KA1 )
01236   820       CONTINUE
01237 *
01238             IF( WANTX ) THEN
01239 *
01240 *              post-multiply X by product of rotations in 1st set
01241 *
01242                DO 830 J = J1, J2, KA1
01243                   CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
01244      $                       RWORK( J ), DCONJG( WORK( J ) ) )
01245   830          CONTINUE
01246             END IF
01247   840    CONTINUE
01248 *
01249          IF( UPDATE ) THEN
01250             IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
01251 *
01252 *              create nonzero element a(i+kbt,i+kbt-ka-1) outside the
01253 *              band and store it in WORK(m-kb+i+kbt)
01254 *
01255                WORK( M-KB+I+KBT ) = -BB( KBT+1, I )*RA1
01256             END IF
01257          END IF
01258 *
01259          DO 880 K = KB, 1, -1
01260             IF( UPDATE ) THEN
01261                J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1
01262             ELSE
01263                J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
01264             END IF
01265 *
01266 *           finish applying rotations in 2nd set from the left
01267 *
01268             DO 850 L = KB - K, 1, -1
01269                NRT = ( J2+KA+L-1 ) / KA1
01270                J1T = J2 - ( NRT-1 )*KA1
01271                IF( NRT.GT.0 )
01272      $            CALL ZLARTV( NRT, AB( KA1-L+1, J1T+L-1 ), INCA,
01273      $                         AB( KA1-L, J1T+L-1 ), INCA,
01274      $                         RWORK( M-KB+J1T+KA ),
01275      $                         WORK( M-KB+J1T+KA ), KA1 )
01276   850       CONTINUE
01277             NR = ( J2+KA-1 ) / KA1
01278             J1 = J2 - ( NR-1 )*KA1
01279             DO 860 J = J1, J2, KA1
01280                WORK( M-KB+J ) = WORK( M-KB+J+KA )
01281                RWORK( M-KB+J ) = RWORK( M-KB+J+KA )
01282   860       CONTINUE
01283             DO 870 J = J1, J2, KA1
01284 *
01285 *              create nonzero element a(j+ka,j-1) outside the band
01286 *              and store it in WORK(m-kb+j)
01287 *
01288                WORK( M-KB+J ) = WORK( M-KB+J )*AB( KA1, J-1 )
01289                AB( KA1, J-1 ) = RWORK( M-KB+J )*AB( KA1, J-1 )
01290   870       CONTINUE
01291             IF( UPDATE ) THEN
01292                IF( I+K.GT.KA1 .AND. K.LE.KBT )
01293      $            WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
01294             END IF
01295   880    CONTINUE
01296 *
01297          DO 920 K = KB, 1, -1
01298             J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
01299             NR = ( J2+KA-1 ) / KA1
01300             J1 = J2 - ( NR-1 )*KA1
01301             IF( NR.GT.0 ) THEN
01302 *
01303 *              generate rotations in 2nd set to annihilate elements
01304 *              which have been created outside the band
01305 *
01306                CALL ZLARGV( NR, AB( KA1, J1 ), INCA, WORK( M-KB+J1 ),
01307      $                      KA1, RWORK( M-KB+J1 ), KA1 )
01308 *
01309 *              apply rotations in 2nd set from the right
01310 *
01311                DO 890 L = 1, KA - 1
01312                   CALL ZLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
01313      $                         INCA, RWORK( M-KB+J1 ), WORK( M-KB+J1 ),
01314      $                         KA1 )
01315   890          CONTINUE
01316 *
01317 *              apply rotations in 2nd set from both sides to diagonal
01318 *              blocks
01319 *
01320                CALL ZLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
01321      $                      AB( 2, J1-1 ), INCA, RWORK( M-KB+J1 ),
01322      $                      WORK( M-KB+J1 ), KA1 )
01323 *
01324                CALL ZLACGV( NR, WORK( M-KB+J1 ), KA1 )
01325             END IF
01326 *
01327 *           start applying rotations in 2nd set from the left
01328 *
01329             DO 900 L = KA - 1, KB - K + 1, -1
01330                NRT = ( J2+L-1 ) / KA1
01331                J1T = J2 - ( NRT-1 )*KA1
01332                IF( NRT.GT.0 )
01333      $            CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
01334      $                         AB( KA1-L, J1T-KA1+L ), INCA,
01335      $                         RWORK( M-KB+J1T ), WORK( M-KB+J1T ),
01336      $                         KA1 )
01337   900       CONTINUE
01338 *
01339             IF( WANTX ) THEN
01340 *
01341 *              post-multiply X by product of rotations in 2nd set
01342 *
01343                DO 910 J = J1, J2, KA1
01344                   CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
01345      $                       RWORK( M-KB+J ), DCONJG( WORK( M-KB+J ) ) )
01346   910          CONTINUE
01347             END IF
01348   920    CONTINUE
01349 *
01350          DO 940 K = 1, KB - 1
01351             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
01352 *
01353 *           finish applying rotations in 1st set from the left
01354 *
01355             DO 930 L = KB - K, 1, -1
01356                NRT = ( J2+L-1 ) / KA1
01357                J1T = J2 - ( NRT-1 )*KA1
01358                IF( NRT.GT.0 )
01359      $            CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
01360      $                         AB( KA1-L, J1T-KA1+L ), INCA,
01361      $                         RWORK( J1T ), WORK( J1T ), KA1 )
01362   930       CONTINUE
01363   940    CONTINUE
01364 *
01365          IF( KB.GT.1 ) THEN
01366             DO 950 J = 2, I2 - KA
01367                RWORK( J ) = RWORK( J+KA )
01368                WORK( J ) = WORK( J+KA )
01369   950       CONTINUE
01370          END IF
01371 *
01372       END IF
01373 *
01374       GO TO 490
01375 *
01376 *     End of ZHBGST
01377 *
01378       END
 All Files Functions