LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zhbgst.f
Go to the documentation of this file.
00001 *> \brief \b ZHBGST
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download ZHBGST + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhbgst.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhbgst.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhbgst.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZHBGST( 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 *       DOUBLE PRECISION   RWORK( * )
00030 *       COMPLEX*16         AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
00031 *      $                   X( LDX, * )
00032 *       ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *> ZHBGST 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 ZPBSTF, 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*16 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*16 array, dimension (LDBB,N)
00110 *>          The banded factor S from the split Cholesky factorization of
00111 *>          B, as returned by ZPBSTF, 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*16 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*16 array, dimension (N)
00138 *> \endverbatim
00139 *>
00140 *> \param[out] RWORK
00141 *> \verbatim
00142 *>          RWORK is DOUBLE PRECISION 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 complex16OTHERcomputational
00163 *
00164 *  =====================================================================
00165       SUBROUTINE ZHBGST( 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       DOUBLE PRECISION   RWORK( * )
00179       COMPLEX*16         AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
00180      $                   X( LDX, * )
00181 *     ..
00182 *
00183 *  =====================================================================
00184 *
00185 *     .. Parameters ..
00186       COMPLEX*16         CZERO, CONE
00187       DOUBLE PRECISION   ONE
00188       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
00189      $                   CONE = ( 1.0D+0, 0.0D+0 ), ONE = 1.0D+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       DOUBLE PRECISION   BII
00196       COMPLEX*16         RA, RA1, T
00197 *     ..
00198 *     .. External Functions ..
00199       LOGICAL            LSAME
00200       EXTERNAL           LSAME
00201 *     ..
00202 *     .. External Subroutines ..
00203       EXTERNAL           XERBLA, ZDSCAL, ZGERC, ZGERU, ZLACGV, ZLAR2V,
00204      $                   ZLARGV, ZLARTG, ZLARTV, ZLASET, ZROT
00205 *     ..
00206 *     .. Intrinsic Functions ..
00207       INTRINSIC          DBLE, DCONJG, MAX, MIN
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( 'ZHBGST', -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 ZLASET( '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 ZPBSTF. 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 = DBLE( BB( KB1, I ) )
00349             AB( KA1, I ) = ( DBLE( 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      $                               DCONJG( AB( K-I+KA1, I ) ) -
00361      $                               DCONJG( BB( K-I+KB1, I ) )*
00362      $                               AB( J-I+KA1, I ) +
00363      $                               DBLE( AB( KA1, I ) )*
00364      $                               BB( J-I+KB1, I )*
00365      $                               DCONJG( 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      $                               DCONJG( 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 ZDSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
00385                IF( KBT.GT.0 )
00386      $            CALL ZGERC( 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 ZLARTG( 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      $                          DCONJG( 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 ZLARGV( 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 ZLARTV( 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 ZLAR2V( 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 ZLACGV( 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 ZLARTV( 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 ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
00485      $                       RWORK( J-M ), DCONJG( 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 ZLARTV( 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 ZLARGV( 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 ZLARTV( 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 ZLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ),
00560      $                      AB( KA, J2+1 ), INCA, RWORK( J2 ),
00561      $                      WORK( J2 ), KA1 )
00562 *
00563                CALL ZLACGV( 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 ZLARTV( 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 ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
00582      $                       RWORK( J ), DCONJG( 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 ZLARTV( 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 = DBLE( BB( 1, I ) )
00617             AB( 1, I ) = ( DBLE( 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 )*DCONJG( AB( I-K+1,
00628      $                             K ) ) - DCONJG( BB( I-K+1, K ) )*
00629      $                             AB( I-J+1, J ) + DBLE( AB( 1, I ) )*
00630      $                             BB( I-J+1, J )*DCONJG( 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      $                             DCONJG( 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 ZDSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
00651                IF( KBT.GT.0 )
00652      $            CALL ZGERU( 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 ZLARTG( 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      $                          DCONJG( WORK( I-K+KA-M ) )*
00685      $                          AB( KA1, I-K )
00686                   AB( KA1, I-K ) = WORK( I-K+KA-M )*T +
00687      $                             RWORK( I-K+KA-M )*AB( KA1, I-K )
00688                   RA1 = RA
00689                END IF
00690             END IF
00691             J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
00692             NR = ( N-J2+KA ) / KA1
00693             J1 = J2 + ( NR-1 )*KA1
00694             IF( UPDATE ) THEN
00695                J2T = MAX( J2, I+2*KA-K+1 )
00696             ELSE
00697                J2T = J2
00698             END IF
00699             NRT = ( N-J2T+KA ) / KA1
00700             DO 320 J = J2T, J1, KA1
00701 *
00702 *              create nonzero element a(j+1,j-ka) outside the band
00703 *              and store it in WORK(j-m)
00704 *
00705                WORK( J-M ) = WORK( J-M )*AB( KA1, J-KA+1 )
00706                AB( KA1, J-KA+1 ) = RWORK( J-M )*AB( KA1, J-KA+1 )
00707   320       CONTINUE
00708 *
00709 *           generate rotations in 1st set to annihilate elements which
00710 *           have been created outside the band
00711 *
00712             IF( NRT.GT.0 )
00713      $         CALL ZLARGV( NRT, AB( KA1, J2T-KA ), INCA, WORK( J2T-M ),
00714      $                      KA1, RWORK( J2T-M ), KA1 )
00715             IF( NR.GT.0 ) THEN
00716 *
00717 *              apply rotations in 1st set from the left
00718 *
00719                DO 330 L = 1, KA - 1
00720                   CALL ZLARTV( NR, AB( L+1, J2-L ), INCA,
00721      $                         AB( L+2, J2-L ), INCA, RWORK( J2-M ),
00722      $                         WORK( J2-M ), KA1 )
00723   330          CONTINUE
00724 *
00725 *              apply rotations in 1st set from both sides to diagonal
00726 *              blocks
00727 *
00728                CALL ZLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
00729      $                      INCA, RWORK( J2-M ), WORK( J2-M ), KA1 )
00730 *
00731                CALL ZLACGV( NR, WORK( J2-M ), KA1 )
00732             END IF
00733 *
00734 *           start applying rotations in 1st set from the right
00735 *
00736             DO 340 L = KA - 1, KB - K + 1, -1
00737                NRT = ( N-J2+L ) / KA1
00738                IF( NRT.GT.0 )
00739      $            CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
00740      $                         AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ),
00741      $                         WORK( J2-M ), KA1 )
00742   340       CONTINUE
00743 *
00744             IF( WANTX ) THEN
00745 *
00746 *              post-multiply X by product of rotations in 1st set
00747 *
00748                DO 350 J = J2, J1, KA1
00749                   CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
00750      $                       RWORK( J-M ), WORK( J-M ) )
00751   350          CONTINUE
00752             END IF
00753   360    CONTINUE
00754 *
00755          IF( UPDATE ) THEN
00756             IF( I2.LE.N .AND. KBT.GT.0 ) THEN
00757 *
00758 *              create nonzero element a(i-kbt+ka+1,i-kbt) outside the
00759 *              band and store it in WORK(i-kbt)
00760 *
00761                WORK( I-KBT ) = -BB( KBT+1, I-KBT )*RA1
00762             END IF
00763          END IF
00764 *
00765          DO 400 K = KB, 1, -1
00766             IF( UPDATE ) THEN
00767                J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1
00768             ELSE
00769                J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
00770             END IF
00771 *
00772 *           finish applying rotations in 2nd set from the right
00773 *
00774             DO 370 L = KB - K, 1, -1
00775                NRT = ( N-J2+KA+L ) / KA1
00776                IF( NRT.GT.0 )
00777      $            CALL ZLARTV( NRT, AB( KA1-L+1, J2-KA ), INCA,
00778      $                         AB( KA1-L, J2-KA+1 ), INCA,
00779      $                         RWORK( J2-KA ), WORK( J2-KA ), KA1 )
00780   370       CONTINUE
00781             NR = ( N-J2+KA ) / KA1
00782             J1 = J2 + ( NR-1 )*KA1
00783             DO 380 J = J1, J2, -KA1
00784                WORK( J ) = WORK( J-KA )
00785                RWORK( J ) = RWORK( J-KA )
00786   380       CONTINUE
00787             DO 390 J = J2, J1, KA1
00788 *
00789 *              create nonzero element a(j+1,j-ka) outside the band
00790 *              and store it in WORK(j)
00791 *
00792                WORK( J ) = WORK( J )*AB( KA1, J-KA+1 )
00793                AB( KA1, J-KA+1 ) = RWORK( J )*AB( KA1, J-KA+1 )
00794   390       CONTINUE
00795             IF( UPDATE ) THEN
00796                IF( I-K.LT.N-KA .AND. K.LE.KBT )
00797      $            WORK( I-K+KA ) = WORK( I-K )
00798             END IF
00799   400    CONTINUE
00800 *
00801          DO 440 K = KB, 1, -1
00802             J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
00803             NR = ( N-J2+KA ) / KA1
00804             J1 = J2 + ( NR-1 )*KA1
00805             IF( NR.GT.0 ) THEN
00806 *
00807 *              generate rotations in 2nd set to annihilate elements
00808 *              which have been created outside the band
00809 *
00810                CALL ZLARGV( NR, AB( KA1, J2-KA ), INCA, WORK( J2 ), KA1,
00811      $                      RWORK( J2 ), KA1 )
00812 *
00813 *              apply rotations in 2nd set from the left
00814 *
00815                DO 410 L = 1, KA - 1
00816                   CALL ZLARTV( NR, AB( L+1, J2-L ), INCA,
00817      $                         AB( L+2, J2-L ), INCA, RWORK( J2 ),
00818      $                         WORK( J2 ), KA1 )
00819   410          CONTINUE
00820 *
00821 *              apply rotations in 2nd set from both sides to diagonal
00822 *              blocks
00823 *
00824                CALL ZLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
00825      $                      INCA, RWORK( J2 ), WORK( J2 ), KA1 )
00826 *
00827                CALL ZLACGV( NR, WORK( J2 ), KA1 )
00828             END IF
00829 *
00830 *           start applying rotations in 2nd set from the right
00831 *
00832             DO 420 L = KA - 1, KB - K + 1, -1
00833                NRT = ( N-J2+L ) / KA1
00834                IF( NRT.GT.0 )
00835      $            CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
00836      $                         AB( KA1-L, J2+1 ), INCA, RWORK( J2 ),
00837      $                         WORK( J2 ), KA1 )
00838   420       CONTINUE
00839 *
00840             IF( WANTX ) THEN
00841 *
00842 *              post-multiply X by product of rotations in 2nd set
00843 *
00844                DO 430 J = J2, J1, KA1
00845                   CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
00846      $                       RWORK( J ), WORK( J ) )
00847   430          CONTINUE
00848             END IF
00849   440    CONTINUE
00850 *
00851          DO 460 K = 1, KB - 1
00852             J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
00853 *
00854 *           finish applying rotations in 1st set from the right
00855 *
00856             DO 450 L = KB - K, 1, -1
00857                NRT = ( N-J2+L ) / KA1
00858                IF( NRT.GT.0 )
00859      $            CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
00860      $                         AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ),
00861      $                         WORK( J2-M ), KA1 )
00862   450       CONTINUE
00863   460    CONTINUE
00864 *
00865          IF( KB.GT.1 ) THEN
00866             DO 470 J = N - 1, J2 + KA, -1
00867                RWORK( J-M ) = RWORK( J-KA-M )
00868                WORK( J-M ) = WORK( J-KA-M )
00869   470       CONTINUE
00870          END IF
00871 *
00872       END IF
00873 *
00874       GO TO 10
00875 *
00876   480 CONTINUE
00877 *
00878 *     **************************** Phase 2 *****************************
00879 *
00880 *     The logical structure of this phase is:
00881 *
00882 *     UPDATE = .TRUE.
00883 *     DO I = 1, M
00884 *        use S(i) to update A and create a new bulge
00885 *        apply rotations to push all bulges KA positions upward
00886 *     END DO
00887 *     UPDATE = .FALSE.
00888 *     DO I = M - KA - 1, 2, -1
00889 *        apply rotations to push all bulges KA positions upward
00890 *     END DO
00891 *
00892 *     To avoid duplicating code, the two loops are merged.
00893 *
00894       UPDATE = .TRUE.
00895       I = 0
00896   490 CONTINUE
00897       IF( UPDATE ) THEN
00898          I = I + 1
00899          KBT = MIN( KB, M-I )
00900          I0 = I + 1
00901          I1 = MAX( 1, I-KA )
00902          I2 = I + KBT - KA1
00903          IF( I.GT.M ) THEN
00904             UPDATE = .FALSE.
00905             I = I - 1
00906             I0 = M + 1
00907             IF( KA.EQ.0 )
00908      $         RETURN
00909             GO TO 490
00910          END IF
00911       ELSE
00912          I = I - KA
00913          IF( I.LT.2 )
00914      $      RETURN
00915       END IF
00916 *
00917       IF( I.LT.M-KBT ) THEN
00918          NX = M
00919       ELSE
00920          NX = N
00921       END IF
00922 *
00923       IF( UPPER ) THEN
00924 *
00925 *        Transform A, working with the upper triangle
00926 *
00927          IF( UPDATE ) THEN
00928 *
00929 *           Form  inv(S(i))**H * A * inv(S(i))
00930 *
00931             BII = DBLE( BB( KB1, I ) )
00932             AB( KA1, I ) = ( DBLE( AB( KA1, I ) ) / BII ) / BII
00933             DO 500 J = I1, I - 1
00934                AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII
00935   500       CONTINUE
00936             DO 510 J = I + 1, MIN( N, I+KA )
00937                AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII
00938   510       CONTINUE
00939             DO 540 K = I + 1, I + KBT
00940                DO 520 J = K, I + KBT
00941                   AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
00942      $                               BB( I-J+KB1, J )*
00943      $                               DCONJG( AB( I-K+KA1, K ) ) -
00944      $                               DCONJG( BB( I-K+KB1, K ) )*
00945      $                               AB( I-J+KA1, J ) +
00946      $                               DBLE( AB( KA1, I ) )*
00947      $                               BB( I-J+KB1, J )*
00948      $                               DCONJG( BB( I-K+KB1, K ) )
00949   520          CONTINUE
00950                DO 530 J = I + KBT + 1, MIN( N, I+KA )
00951                   AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
00952      $                               DCONJG( BB( I-K+KB1, K ) )*
00953      $                               AB( I-J+KA1, J )
00954   530          CONTINUE
00955   540       CONTINUE
00956             DO 560 J = I1, I
00957                DO 550 K = I + 1, MIN( J+KA, I+KBT )
00958                   AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
00959      $                               BB( I-K+KB1, K )*AB( J-I+KA1, I )
00960   550          CONTINUE
00961   560       CONTINUE
00962 *
00963             IF( WANTX ) THEN
00964 *
00965 *              post-multiply X by inv(S(i))
00966 *
00967                CALL ZDSCAL( NX, ONE / BII, X( 1, I ), 1 )
00968                IF( KBT.GT.0 )
00969      $            CALL ZGERU( NX, KBT, -CONE, X( 1, I ), 1,
00970      $                        BB( KB, I+1 ), LDBB-1, X( 1, I+1 ), LDX )
00971             END IF
00972 *
00973 *           store a(i1,i) in RA1 for use in next loop over K
00974 *
00975             RA1 = AB( I1-I+KA1, I )
00976          END IF
00977 *
00978 *        Generate and apply vectors of rotations to chase all the
00979 *        existing bulges KA positions up toward the top of the band
00980 *
00981          DO 610 K = 1, KB - 1
00982             IF( UPDATE ) THEN
00983 *
00984 *              Determine the rotations which would annihilate the bulge
00985 *              which has in theory just been created
00986 *
00987                IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
00988 *
00989 *                 generate rotation to annihilate a(i+k-ka-1,i)
00990 *
00991                   CALL ZLARTG( AB( K+1, I ), RA1, RWORK( I+K-KA ),
00992      $                         WORK( I+K-KA ), RA )
00993 *
00994 *                 create nonzero element a(i+k-ka-1,i+k) outside the
00995 *                 band and store it in WORK(m-kb+i+k)
00996 *
00997                   T = -BB( KB1-K, I+K )*RA1
00998                   WORK( M-KB+I+K ) = RWORK( I+K-KA )*T -
00999      $                               DCONJG( WORK( I+K-KA ) )*
01000      $                               AB( 1, I+K )
01001                   AB( 1, I+K ) = WORK( I+K-KA )*T +
01002      $                           RWORK( I+K-KA )*AB( 1, I+K )
01003                   RA1 = RA
01004                END IF
01005             END IF
01006             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
01007             NR = ( J2+KA-1 ) / KA1
01008             J1 = J2 - ( NR-1 )*KA1
01009             IF( UPDATE ) THEN
01010                J2T = MIN( J2, I-2*KA+K-1 )
01011             ELSE
01012                J2T = J2
01013             END IF
01014             NRT = ( J2T+KA-1 ) / KA1
01015             DO 570 J = J1, J2T, KA1
01016 *
01017 *              create nonzero element a(j-1,j+ka) outside the band
01018 *              and store it in WORK(j)
01019 *
01020                WORK( J ) = WORK( J )*AB( 1, J+KA-1 )
01021                AB( 1, J+KA-1 ) = RWORK( J )*AB( 1, J+KA-1 )
01022   570       CONTINUE
01023 *
01024 *           generate rotations in 1st set to annihilate elements which
01025 *           have been created outside the band
01026 *
01027             IF( NRT.GT.0 )
01028      $         CALL ZLARGV( NRT, AB( 1, J1+KA ), INCA, WORK( J1 ), KA1,
01029      $                      RWORK( J1 ), KA1 )
01030             IF( NR.GT.0 ) THEN
01031 *
01032 *              apply rotations in 1st set from the left
01033 *
01034                DO 580 L = 1, KA - 1
01035                   CALL ZLARTV( NR, AB( KA1-L, J1+L ), INCA,
01036      $                         AB( KA-L, J1+L ), INCA, RWORK( J1 ),
01037      $                         WORK( J1 ), KA1 )
01038   580          CONTINUE
01039 *
01040 *              apply rotations in 1st set from both sides to diagonal
01041 *              blocks
01042 *
01043                CALL ZLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
01044      $                      AB( KA, J1 ), INCA, RWORK( J1 ), WORK( J1 ),
01045      $                      KA1 )
01046 *
01047                CALL ZLACGV( NR, WORK( J1 ), KA1 )
01048             END IF
01049 *
01050 *           start applying rotations in 1st set from the right
01051 *
01052             DO 590 L = KA - 1, KB - K + 1, -1
01053                NRT = ( J2+L-1 ) / KA1
01054                J1T = J2 - ( NRT-1 )*KA1
01055                IF( NRT.GT.0 )
01056      $            CALL ZLARTV( NRT, AB( L, J1T ), INCA,
01057      $                         AB( L+1, J1T-1 ), INCA, RWORK( J1T ),
01058      $                         WORK( J1T ), KA1 )
01059   590       CONTINUE
01060 *
01061             IF( WANTX ) THEN
01062 *
01063 *              post-multiply X by product of rotations in 1st set
01064 *
01065                DO 600 J = J1, J2, KA1
01066                   CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
01067      $                       RWORK( J ), WORK( J ) )
01068   600          CONTINUE
01069             END IF
01070   610    CONTINUE
01071 *
01072          IF( UPDATE ) THEN
01073             IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
01074 *
01075 *              create nonzero element a(i+kbt-ka-1,i+kbt) outside the
01076 *              band and store it in WORK(m-kb+i+kbt)
01077 *
01078                WORK( M-KB+I+KBT ) = -BB( KB1-KBT, I+KBT )*RA1
01079             END IF
01080          END IF
01081 *
01082          DO 650 K = KB, 1, -1
01083             IF( UPDATE ) THEN
01084                J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1
01085             ELSE
01086                J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
01087             END IF
01088 *
01089 *           finish applying rotations in 2nd set from the right
01090 *
01091             DO 620 L = KB - K, 1, -1
01092                NRT = ( J2+KA+L-1 ) / KA1
01093                J1T = J2 - ( NRT-1 )*KA1
01094                IF( NRT.GT.0 )
01095      $            CALL ZLARTV( NRT, AB( L, J1T+KA ), INCA,
01096      $                         AB( L+1, J1T+KA-1 ), INCA,
01097      $                         RWORK( M-KB+J1T+KA ),
01098      $                         WORK( M-KB+J1T+KA ), KA1 )
01099   620       CONTINUE
01100             NR = ( J2+KA-1 ) / KA1
01101             J1 = J2 - ( NR-1 )*KA1
01102             DO 630 J = J1, J2, KA1
01103                WORK( M-KB+J ) = WORK( M-KB+J+KA )
01104                RWORK( M-KB+J ) = RWORK( M-KB+J+KA )
01105   630       CONTINUE
01106             DO 640 J = J1, J2, KA1
01107 *
01108 *              create nonzero element a(j-1,j+ka) outside the band
01109 *              and store it in WORK(m-kb+j)
01110 *
01111                WORK( M-KB+J ) = WORK( M-KB+J )*AB( 1, J+KA-1 )
01112                AB( 1, J+KA-1 ) = RWORK( M-KB+J )*AB( 1, J+KA-1 )
01113   640       CONTINUE
01114             IF( UPDATE ) THEN
01115                IF( I+K.GT.KA1 .AND. K.LE.KBT )
01116      $            WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
01117             END IF
01118   650    CONTINUE
01119 *
01120          DO 690 K = KB, 1, -1
01121             J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
01122             NR = ( J2+KA-1 ) / KA1
01123             J1 = J2 - ( NR-1 )*KA1
01124             IF( NR.GT.0 ) THEN
01125 *
01126 *              generate rotations in 2nd set to annihilate elements
01127 *              which have been created outside the band
01128 *
01129                CALL ZLARGV( NR, AB( 1, J1+KA ), INCA, WORK( M-KB+J1 ),
01130      $                      KA1, RWORK( M-KB+J1 ), KA1 )
01131 *
01132 *              apply rotations in 2nd set from the left
01133 *
01134                DO 660 L = 1, KA - 1
01135                   CALL ZLARTV( NR, AB( KA1-L, J1+L ), INCA,
01136      $                         AB( KA-L, J1+L ), INCA, RWORK( M-KB+J1 ),
01137      $                         WORK( M-KB+J1 ), KA1 )
01138   660          CONTINUE
01139 *
01140 *              apply rotations in 2nd set from both sides to diagonal
01141 *              blocks
01142 *
01143                CALL ZLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
01144      $                      AB( KA, J1 ), INCA, RWORK( M-KB+J1 ),
01145      $                      WORK( M-KB+J1 ), KA1 )
01146 *
01147                CALL ZLACGV( NR, WORK( M-KB+J1 ), KA1 )
01148             END IF
01149 *
01150 *           start applying rotations in 2nd set from the right
01151 *
01152             DO 670 L = KA - 1, KB - K + 1, -1
01153                NRT = ( J2+L-1 ) / KA1
01154                J1T = J2 - ( NRT-1 )*KA1
01155                IF( NRT.GT.0 )
01156      $            CALL ZLARTV( NRT, AB( L, J1T ), INCA,
01157      $                         AB( L+1, J1T-1 ), INCA,
01158      $                         RWORK( M-KB+J1T ), WORK( M-KB+J1T ),
01159      $                         KA1 )
01160   670       CONTINUE
01161 *
01162             IF( WANTX ) THEN
01163 *
01164 *              post-multiply X by product of rotations in 2nd set
01165 *
01166                DO 680 J = J1, J2, KA1
01167                   CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
01168      $                       RWORK( M-KB+J ), WORK( M-KB+J ) )
01169   680          CONTINUE
01170             END IF
01171   690    CONTINUE
01172 *
01173          DO 710 K = 1, KB - 1
01174             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
01175 *
01176 *           finish applying rotations in 1st set from the right
01177 *
01178             DO 700 L = KB - K, 1, -1
01179                NRT = ( J2+L-1 ) / KA1
01180                J1T = J2 - ( NRT-1 )*KA1
01181                IF( NRT.GT.0 )
01182      $            CALL ZLARTV( NRT, AB( L, J1T ), INCA,
01183      $                         AB( L+1, J1T-1 ), INCA, RWORK( J1T ),
01184      $                         WORK( J1T ), KA1 )
01185   700       CONTINUE
01186   710    CONTINUE
01187 *
01188          IF( KB.GT.1 ) THEN
01189             DO 720 J = 2, I2 - KA
01190                RWORK( J ) = RWORK( J+KA )
01191                WORK( J ) = WORK( J+KA )
01192   720       CONTINUE
01193          END IF
01194 *
01195       ELSE
01196 *
01197 *        Transform A, working with the lower triangle
01198 *
01199          IF( UPDATE ) THEN
01200 *
01201 *           Form  inv(S(i))**H * A * inv(S(i))
01202 *
01203             BII = DBLE( BB( 1, I ) )
01204             AB( 1, I ) = ( DBLE( AB( 1, I ) ) / BII ) / BII
01205             DO 730 J = I1, I - 1
01206                AB( I-J+1, J ) = AB( I-J+1, J ) / BII
01207   730       CONTINUE
01208             DO 740 J = I + 1, MIN( N, I+KA )
01209                AB( J-I+1, I ) = AB( J-I+1, I ) / BII
01210   740       CONTINUE
01211             DO 770 K = I + 1, I + KBT
01212                DO 750 J = K, I + KBT
01213                   AB( J-K+1, K ) = AB( J-K+1, K ) -
01214      $                             BB( J-I+1, I )*DCONJG( AB( K-I+1,
01215      $                             I ) ) - DCONJG( BB( K-I+1, I ) )*
01216      $                             AB( J-I+1, I ) + DBLE( AB( 1, I ) )*
01217      $                             BB( J-I+1, I )*DCONJG( BB( K-I+1,
01218      $                             I ) )
01219   750          CONTINUE
01220                DO 760 J = I + KBT + 1, MIN( N, I+KA )
01221                   AB( J-K+1, K ) = AB( J-K+1, K ) -
01222      $                             DCONJG( BB( K-I+1, I ) )*
01223      $                             AB( J-I+1, I )
01224   760          CONTINUE
01225   770       CONTINUE
01226             DO 790 J = I1, I
01227                DO 780 K = I + 1, MIN( J+KA, I+KBT )
01228                   AB( K-J+1, J ) = AB( K-J+1, J ) -
01229      $                             BB( K-I+1, I )*AB( I-J+1, J )
01230   780          CONTINUE
01231   790       CONTINUE
01232 *
01233             IF( WANTX ) THEN
01234 *
01235 *              post-multiply X by inv(S(i))
01236 *
01237                CALL ZDSCAL( NX, ONE / BII, X( 1, I ), 1 )
01238                IF( KBT.GT.0 )
01239      $            CALL ZGERC( NX, KBT, -CONE, X( 1, I ), 1, BB( 2, I ),
01240      $                        1, X( 1, I+1 ), LDX )
01241             END IF
01242 *
01243 *           store a(i,i1) in RA1 for use in next loop over K
01244 *
01245             RA1 = AB( I-I1+1, I1 )
01246          END IF
01247 *
01248 *        Generate and apply vectors of rotations to chase all the
01249 *        existing bulges KA positions up toward the top of the band
01250 *
01251          DO 840 K = 1, KB - 1
01252             IF( UPDATE ) THEN
01253 *
01254 *              Determine the rotations which would annihilate the bulge
01255 *              which has in theory just been created
01256 *
01257                IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
01258 *
01259 *                 generate rotation to annihilate a(i,i+k-ka-1)
01260 *
01261                   CALL ZLARTG( AB( KA1-K, I+K-KA ), RA1,
01262      $                         RWORK( I+K-KA ), WORK( I+K-KA ), RA )
01263 *
01264 *                 create nonzero element a(i+k,i+k-ka-1) outside the
01265 *                 band and store it in WORK(m-kb+i+k)
01266 *
01267                   T = -BB( K+1, I )*RA1
01268                   WORK( M-KB+I+K ) = RWORK( I+K-KA )*T -
01269      $                               DCONJG( WORK( I+K-KA ) )*
01270      $                               AB( KA1, I+K-KA )
01271                   AB( KA1, I+K-KA ) = WORK( I+K-KA )*T +
01272      $                                RWORK( I+K-KA )*AB( KA1, I+K-KA )
01273                   RA1 = RA
01274                END IF
01275             END IF
01276             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
01277             NR = ( J2+KA-1 ) / KA1
01278             J1 = J2 - ( NR-1 )*KA1
01279             IF( UPDATE ) THEN
01280                J2T = MIN( J2, I-2*KA+K-1 )
01281             ELSE
01282                J2T = J2
01283             END IF
01284             NRT = ( J2T+KA-1 ) / KA1
01285             DO 800 J = J1, J2T, KA1
01286 *
01287 *              create nonzero element a(j+ka,j-1) outside the band
01288 *              and store it in WORK(j)
01289 *
01290                WORK( J ) = WORK( J )*AB( KA1, J-1 )
01291                AB( KA1, J-1 ) = RWORK( J )*AB( KA1, J-1 )
01292   800       CONTINUE
01293 *
01294 *           generate rotations in 1st set to annihilate elements which
01295 *           have been created outside the band
01296 *
01297             IF( NRT.GT.0 )
01298      $         CALL ZLARGV( NRT, AB( KA1, J1 ), INCA, WORK( J1 ), KA1,
01299      $                      RWORK( J1 ), KA1 )
01300             IF( NR.GT.0 ) THEN
01301 *
01302 *              apply rotations in 1st set from the right
01303 *
01304                DO 810 L = 1, KA - 1
01305                   CALL ZLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
01306      $                         INCA, RWORK( J1 ), WORK( J1 ), KA1 )
01307   810          CONTINUE
01308 *
01309 *              apply rotations in 1st set from both sides to diagonal
01310 *              blocks
01311 *
01312                CALL ZLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
01313      $                      AB( 2, J1-1 ), INCA, RWORK( J1 ),
01314      $                      WORK( J1 ), KA1 )
01315 *
01316                CALL ZLACGV( NR, WORK( J1 ), KA1 )
01317             END IF
01318 *
01319 *           start applying rotations in 1st set from the left
01320 *
01321             DO 820 L = KA - 1, KB - K + 1, -1
01322                NRT = ( J2+L-1 ) / KA1
01323                J1T = J2 - ( NRT-1 )*KA1
01324                IF( NRT.GT.0 )
01325      $            CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
01326      $                         AB( KA1-L, J1T-KA1+L ), INCA,
01327      $                         RWORK( J1T ), WORK( J1T ), KA1 )
01328   820       CONTINUE
01329 *
01330             IF( WANTX ) THEN
01331 *
01332 *              post-multiply X by product of rotations in 1st set
01333 *
01334                DO 830 J = J1, J2, KA1
01335                   CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
01336      $                       RWORK( J ), DCONJG( WORK( J ) ) )
01337   830          CONTINUE
01338             END IF
01339   840    CONTINUE
01340 *
01341          IF( UPDATE ) THEN
01342             IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
01343 *
01344 *              create nonzero element a(i+kbt,i+kbt-ka-1) outside the
01345 *              band and store it in WORK(m-kb+i+kbt)
01346 *
01347                WORK( M-KB+I+KBT ) = -BB( KBT+1, I )*RA1
01348             END IF
01349          END IF
01350 *
01351          DO 880 K = KB, 1, -1
01352             IF( UPDATE ) THEN
01353                J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1
01354             ELSE
01355                J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
01356             END IF
01357 *
01358 *           finish applying rotations in 2nd set from the left
01359 *
01360             DO 850 L = KB - K, 1, -1
01361                NRT = ( J2+KA+L-1 ) / KA1
01362                J1T = J2 - ( NRT-1 )*KA1
01363                IF( NRT.GT.0 )
01364      $            CALL ZLARTV( NRT, AB( KA1-L+1, J1T+L-1 ), INCA,
01365      $                         AB( KA1-L, J1T+L-1 ), INCA,
01366      $                         RWORK( M-KB+J1T+KA ),
01367      $                         WORK( M-KB+J1T+KA ), KA1 )
01368   850       CONTINUE
01369             NR = ( J2+KA-1 ) / KA1
01370             J1 = J2 - ( NR-1 )*KA1
01371             DO 860 J = J1, J2, KA1
01372                WORK( M-KB+J ) = WORK( M-KB+J+KA )
01373                RWORK( M-KB+J ) = RWORK( M-KB+J+KA )
01374   860       CONTINUE
01375             DO 870 J = J1, J2, KA1
01376 *
01377 *              create nonzero element a(j+ka,j-1) outside the band
01378 *              and store it in WORK(m-kb+j)
01379 *
01380                WORK( M-KB+J ) = WORK( M-KB+J )*AB( KA1, J-1 )
01381                AB( KA1, J-1 ) = RWORK( M-KB+J )*AB( KA1, J-1 )
01382   870       CONTINUE
01383             IF( UPDATE ) THEN
01384                IF( I+K.GT.KA1 .AND. K.LE.KBT )
01385      $            WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
01386             END IF
01387   880    CONTINUE
01388 *
01389          DO 920 K = KB, 1, -1
01390             J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
01391             NR = ( J2+KA-1 ) / KA1
01392             J1 = J2 - ( NR-1 )*KA1
01393             IF( NR.GT.0 ) THEN
01394 *
01395 *              generate rotations in 2nd set to annihilate elements
01396 *              which have been created outside the band
01397 *
01398                CALL ZLARGV( NR, AB( KA1, J1 ), INCA, WORK( M-KB+J1 ),
01399      $                      KA1, RWORK( M-KB+J1 ), KA1 )
01400 *
01401 *              apply rotations in 2nd set from the right
01402 *
01403                DO 890 L = 1, KA - 1
01404                   CALL ZLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
01405      $                         INCA, RWORK( M-KB+J1 ), WORK( M-KB+J1 ),
01406      $                         KA1 )
01407   890          CONTINUE
01408 *
01409 *              apply rotations in 2nd set from both sides to diagonal
01410 *              blocks
01411 *
01412                CALL ZLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
01413      $                      AB( 2, J1-1 ), INCA, RWORK( M-KB+J1 ),
01414      $                      WORK( M-KB+J1 ), KA1 )
01415 *
01416                CALL ZLACGV( NR, WORK( M-KB+J1 ), KA1 )
01417             END IF
01418 *
01419 *           start applying rotations in 2nd set from the left
01420 *
01421             DO 900 L = KA - 1, KB - K + 1, -1
01422                NRT = ( J2+L-1 ) / KA1
01423                J1T = J2 - ( NRT-1 )*KA1
01424                IF( NRT.GT.0 )
01425      $            CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
01426      $                         AB( KA1-L, J1T-KA1+L ), INCA,
01427      $                         RWORK( M-KB+J1T ), WORK( M-KB+J1T ),
01428      $                         KA1 )
01429   900       CONTINUE
01430 *
01431             IF( WANTX ) THEN
01432 *
01433 *              post-multiply X by product of rotations in 2nd set
01434 *
01435                DO 910 J = J1, J2, KA1
01436                   CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
01437      $                       RWORK( M-KB+J ), DCONJG( WORK( M-KB+J ) ) )
01438   910          CONTINUE
01439             END IF
01440   920    CONTINUE
01441 *
01442          DO 940 K = 1, KB - 1
01443             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
01444 *
01445 *           finish applying rotations in 1st set from the left
01446 *
01447             DO 930 L = KB - K, 1, -1
01448                NRT = ( J2+L-1 ) / KA1
01449                J1T = J2 - ( NRT-1 )*KA1
01450                IF( NRT.GT.0 )
01451      $            CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
01452      $                         AB( KA1-L, J1T-KA1+L ), INCA,
01453      $                         RWORK( J1T ), WORK( J1T ), KA1 )
01454   930       CONTINUE
01455   940    CONTINUE
01456 *
01457          IF( KB.GT.1 ) THEN
01458             DO 950 J = 2, I2 - KA
01459                RWORK( J ) = RWORK( J+KA )
01460                WORK( J ) = WORK( J+KA )
01461   950       CONTINUE
01462          END IF
01463 *
01464       END IF
01465 *
01466       GO TO 490
01467 *
01468 *     End of ZHBGST
01469 *
01470       END
 All Files Functions