SRC/zhbgst.f

   1 *> \brief \b ZHBGST
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZHBGST + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhbgst.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhbgst.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhbgst.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZHBGST( VECT, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, X,
  22 *                          LDX, WORK, RWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          UPLO, VECT
  26 *       INTEGER            INFO, KA, KB, LDAB, LDBB, LDX, N
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       DOUBLE PRECISION   RWORK( * )
  30 *       COMPLEX*16         AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
  31 *      $                   X( LDX, * )
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *> ZHBGST reduces a complex Hermitian-definite banded generalized
  41 *> eigenproblem  A*x = lambda*B*x  to standard form  C*y = lambda*y,
  42 *> such that C has the same bandwidth as A.
  43 *>
  44 *> B must have been previously factorized as S**H*S by ZPBSTF, using a
  45 *> split Cholesky factorization. A is overwritten by C = X**H*A*X, where
  46 *> X = S**(-1)*Q and Q is a unitary matrix chosen to preserve the
  47 *> bandwidth of A.
  48 *> \endverbatim
  49 *
  50 *  Arguments:
  51 *  ==========
  52 *
  53 *> \param[in] VECT
  54 *> \verbatim
  55 *>          VECT is CHARACTER*1
  56 *>          = 'N':  do not form the transformation matrix X;
  57 *>          = 'V':  form X.
  58 *> \endverbatim
  59 *>
  60 *> \param[in] UPLO
  61 *> \verbatim
  62 *>          UPLO is CHARACTER*1
  63 *>          = 'U':  Upper triangle of A is stored;
  64 *>          = 'L':  Lower triangle of A is stored.
  65 *> \endverbatim
  66 *>
  67 *> \param[in] N
  68 *> \verbatim
  69 *>          N is INTEGER
  70 *>          The order of the matrices A and B.  N >= 0.
  71 *> \endverbatim
  72 *>
  73 *> \param[in] KA
  74 *> \verbatim
  75 *>          KA is INTEGER
  76 *>          The number of superdiagonals of the matrix A if UPLO = 'U',
  77 *>          or the number of subdiagonals if UPLO = 'L'.  KA >= 0.
  78 *> \endverbatim
  79 *>
  80 *> \param[in] KB
  81 *> \verbatim
  82 *>          KB is INTEGER
  83 *>          The number of superdiagonals of the matrix B if UPLO = 'U',
  84 *>          or the number of subdiagonals if UPLO = 'L'.  KA >= KB >= 0.
  85 *> \endverbatim
  86 *>
  87 *> \param[in,out] AB
  88 *> \verbatim
  89 *>          AB is COMPLEX*16 array, dimension (LDAB,N)
  90 *>          On entry, the upper or lower triangle of the Hermitian band
  91 *>          matrix A, stored in the first ka+1 rows of the array.  The
  92 *>          j-th column of A is stored in the j-th column of the array AB
  93 *>          as follows:
  94 *>          if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;
  95 *>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+ka).
  96 *>
  97 *>          On exit, the transformed matrix X**H*A*X, stored in the same
  98 *>          format as A.
  99 *> \endverbatim
 100 *>
 101 *> \param[in] LDAB
 102 *> \verbatim
 103 *>          LDAB is INTEGER
 104 *>          The leading dimension of the array AB.  LDAB >= KA+1.
 105 *> \endverbatim
 106 *>
 107 *> \param[in] BB
 108 *> \verbatim
 109 *>          BB is COMPLEX*16 array, dimension (LDBB,N)
 110 *>          The banded factor S from the split Cholesky factorization of
 111 *>          B, as returned by ZPBSTF, stored in the first kb+1 rows of
 112 *>          the array.
 113 *> \endverbatim
 114 *>
 115 *> \param[in] LDBB
 116 *> \verbatim
 117 *>          LDBB is INTEGER
 118 *>          The leading dimension of the array BB.  LDBB >= KB+1.
 119 *> \endverbatim
 120 *>
 121 *> \param[out] X
 122 *> \verbatim
 123 *>          X is COMPLEX*16 array, dimension (LDX,N)
 124 *>          If VECT = 'V', the n-by-n matrix X.
 125 *>          If VECT = 'N', the array X is not referenced.
 126 *> \endverbatim
 127 *>
 128 *> \param[in] LDX
 129 *> \verbatim
 130 *>          LDX is INTEGER
 131 *>          The leading dimension of the array X.
 132 *>          LDX >= max(1,N) if VECT = 'V'; LDX >= 1 otherwise.
 133 *> \endverbatim
 134 *>
 135 *> \param[out] WORK
 136 *> \verbatim
 137 *>          WORK is COMPLEX*16 array, dimension (N)
 138 *> \endverbatim
 139 *>
 140 *> \param[out] RWORK
 141 *> \verbatim
 142 *>          RWORK is DOUBLE PRECISION array, dimension (N)
 143 *> \endverbatim
 144 *>
 145 *> \param[out] INFO
 146 *> \verbatim
 147 *>          INFO is INTEGER
 148 *>          = 0:  successful exit
 149 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 150 *> \endverbatim
 151 *
 152 *  Authors:
 153 *  ========
 154 *
 155 *> \author Univ. of Tennessee
 156 *> \author Univ. of California Berkeley
 157 *> \author Univ. of Colorado Denver
 158 *> \author NAG Ltd.
 159 *
 160 *> \date November 2011
 161 *
 162 *> \ingroup complex16OTHERcomputational
 163 *
 164 *  =====================================================================
 165       SUBROUTINE ZHBGST( VECT, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, X,
 166      $                   LDX, WORK, RWORK, INFO )
 167 *
 168 *  -- LAPACK computational routine (version 3.4.0) --
 169 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 170 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 171 *     November 2011
 172 *
 173 *     .. Scalar Arguments ..
 174       CHARACTER          UPLO, VECT
 175       INTEGER            INFO, KA, KB, LDAB, LDBB, LDX, N
 176 *     ..
 177 *     .. Array Arguments ..
 178       DOUBLE PRECISION   RWORK( * )
 179       COMPLEX*16         AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
 180      $                   X( LDX, * )
 181 *     ..
 182 *
 183 *  =====================================================================
 184 *
 185 *     .. Parameters ..
 186       COMPLEX*16         CZERO, CONE
 187       DOUBLE PRECISION   ONE
 188       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
 189      $                   CONE = ( 1.0D+0, 0.0D+0 ), ONE = 1.0D+0 )
 190 *     ..
 191 *     .. Local Scalars ..
 192       LOGICAL            UPDATE, UPPER, WANTX
 193       INTEGER            I, I0, I1, I2, INCA, J, J1, J1T, J2, J2T, K,
 194      $                   KA1, KB1, KBT, L, M, NR, NRT, NX
 195       DOUBLE PRECISION   BII
 196       COMPLEX*16         RA, RA1, T
 197 *     ..
 198 *     .. External Functions ..
 199       LOGICAL            LSAME
 200       EXTERNAL           LSAME
 201 *     ..
 202 *     .. External Subroutines ..
 203       EXTERNAL           XERBLA, ZDSCAL, ZGERC, ZGERU, ZLACGV, ZLAR2V,
 204      $                   ZLARGV, ZLARTG, ZLARTV, ZLASET, ZROT
 205 *     ..
 206 *     .. Intrinsic Functions ..
 207       INTRINSIC          DBLE, DCONJG, MAX, MIN
 208 *     ..
 209 *     .. Executable Statements ..
 210 *
 211 *     Test the input parameters
 212 *
 213       WANTX = LSAME( VECT, 'V' )
 214       UPPER = LSAME( UPLO, 'U' )
 215       KA1 = KA + 1
 216       KB1 = KB + 1
 217       INFO = 0
 218       IF( .NOT.WANTX .AND. .NOT.LSAME( VECT, 'N' ) ) THEN
 219          INFO = -1
 220       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 221          INFO = -2
 222       ELSE IF( N.LT.0 ) THEN
 223          INFO = -3
 224       ELSE IF( KA.LT.0 ) THEN
 225          INFO = -4
 226       ELSE IF( KB.LT.0 .OR. KB.GT.KA ) THEN
 227          INFO = -5
 228       ELSE IF( LDAB.LT.KA+1 ) THEN
 229          INFO = -7
 230       ELSE IF( LDBB.LT.KB+1 ) THEN
 231          INFO = -9
 232       ELSE IF( LDX.LT.1 .OR. WANTX .AND. LDX.LT.MAX( 1, N ) ) THEN
 233          INFO = -11
 234       END IF
 235       IF( INFO.NE.0 ) THEN
 236          CALL XERBLA( 'ZHBGST', -INFO )
 237          RETURN
 238       END IF
 239 *
 240 *     Quick return if possible
 241 *
 242       IF( N.EQ.0 )
 243      $   RETURN
 244 *
 245       INCA = LDAB*KA1
 246 *
 247 *     Initialize X to the unit matrix, if needed
 248 *
 249       IF( WANTX )
 250      $   CALL ZLASET( 'Full', N, N, CZERO, CONE, X, LDX )
 251 *
 252 *     Set M to the splitting point m. It must be the same value as is
 253 *     used in ZPBSTF. The chosen value allows the arrays WORK and RWORK
 254 *     to be of dimension (N).
 255 *
 256       M = ( N+KB ) / 2
 257 *
 258 *     The routine works in two phases, corresponding to the two halves
 259 *     of the split Cholesky factorization of B as S**H*S where
 260 *
 261 *     S = ( U    )
 262 *         ( M  L )
 263 *
 264 *     with U upper triangular of order m, and L lower triangular of
 265 *     order n-m. S has the same bandwidth as B.
 266 *
 267 *     S is treated as a product of elementary matrices:
 268 *
 269 *     S = S(m)*S(m-1)*...*S(2)*S(1)*S(m+1)*S(m+2)*...*S(n-1)*S(n)
 270 *
 271 *     where S(i) is determined by the i-th row of S.
 272 *
 273 *     In phase 1, the index i takes the values n, n-1, ... , m+1;
 274 *     in phase 2, it takes the values 1, 2, ... , m.
 275 *
 276 *     For each value of i, the current matrix A is updated by forming
 277 *     inv(S(i))**H*A*inv(S(i)). This creates a triangular bulge outside
 278 *     the band of A. The bulge is then pushed down toward the bottom of
 279 *     A in phase 1, and up toward the top of A in phase 2, by applying
 280 *     plane rotations.
 281 *
 282 *     There are kb*(kb+1)/2 elements in the bulge, but at most 2*kb-1
 283 *     of them are linearly independent, so annihilating a bulge requires
 284 *     only 2*kb-1 plane rotations. The rotations are divided into a 1st
 285 *     set of kb-1 rotations, and a 2nd set of kb rotations.
 286 *
 287 *     Wherever possible, rotations are generated and applied in vector
 288 *     operations of length NR between the indices J1 and J2 (sometimes
 289 *     replaced by modified values NRT, J1T or J2T).
 290 *
 291 *     The real cosines and complex sines of the rotations are stored in
 292 *     the arrays RWORK and WORK, those of the 1st set in elements
 293 *     2:m-kb-1, and those of the 2nd set in elements m-kb+1:n.
 294 *
 295 *     The bulges are not formed explicitly; nonzero elements outside the
 296 *     band are created only when they are required for generating new
 297 *     rotations; they are stored in the array WORK, in positions where
 298 *     they are later overwritten by the sines of the rotations which
 299 *     annihilate them.
 300 *
 301 *     **************************** Phase 1 *****************************
 302 *
 303 *     The logical structure of this phase is:
 304 *
 305 *     UPDATE = .TRUE.
 306 *     DO I = N, M + 1, -1
 307 *        use S(i) to update A and create a new bulge
 308 *        apply rotations to push all bulges KA positions downward
 309 *     END DO
 310 *     UPDATE = .FALSE.
 311 *     DO I = M + KA + 1, N - 1
 312 *        apply rotations to push all bulges KA positions downward
 313 *     END DO
 314 *
 315 *     To avoid duplicating code, the two loops are merged.
 316 *
 317       UPDATE = .TRUE.
 318       I = N + 1
 319    10 CONTINUE
 320       IF( UPDATE ) THEN
 321          I = I - 1
 322          KBT = MIN( KB, I-1 )
 323          I0 = I - 1
 324          I1 = MIN( N, I+KA )
 325          I2 = I - KBT + KA1
 326          IF( I.LT.M+1 ) THEN
 327             UPDATE = .FALSE.
 328             I = I + 1
 329             I0 = M
 330             IF( KA.EQ.0 )
 331      $         GO TO 480
 332             GO TO 10
 333          END IF
 334       ELSE
 335          I = I + KA
 336          IF( I.GT.N-1 )
 337      $      GO TO 480
 338       END IF
 339 *
 340       IF( UPPER ) THEN
 341 *
 342 *        Transform A, working with the upper triangle
 343 *
 344          IF( UPDATE ) THEN
 345 *
 346 *           Form  inv(S(i))**H * A * inv(S(i))
 347 *
 348             BII = DBLE( BB( KB1, I ) )
 349             AB( KA1, I ) = ( DBLE( AB( KA1, I ) ) / BII ) / BII
 350             DO 20 J = I + 1, I1
 351                AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII
 352    20       CONTINUE
 353             DO 30 J = MAX( 1, I-KA ), I - 1
 354                AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII
 355    30       CONTINUE
 356             DO 60 K = I - KBT, I - 1
 357                DO 40 J = I - KBT, K
 358                   AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
 359      $                               BB( J-I+KB1, I )*
 360      $                               DCONJG( AB( K-I+KA1, I ) ) -
 361      $                               DCONJG( BB( K-I+KB1, I ) )*
 362      $                               AB( J-I+KA1, I ) +
 363      $                               DBLE( AB( KA1, I ) )*
 364      $                               BB( J-I+KB1, I )*
 365      $                               DCONJG( BB( K-I+KB1, I ) )
 366    40          CONTINUE
 367                DO 50 J = MAX( 1, I-KA ), I - KBT - 1
 368                   AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
 369      $                               DCONJG( BB( K-I+KB1, I ) )*
 370      $                               AB( J-I+KA1, I )
 371    50          CONTINUE
 372    60       CONTINUE
 373             DO 80 J = I, I1
 374                DO 70 K = MAX( J-KA, I-KBT ), I - 1
 375                   AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
 376      $                               BB( K-I+KB1, I )*AB( I-J+KA1, J )
 377    70          CONTINUE
 378    80       CONTINUE
 379 *
 380             IF( WANTX ) THEN
 381 *
 382 *              post-multiply X by inv(S(i))
 383 *
 384                CALL ZDSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
 385                IF( KBT.GT.0 )
 386      $            CALL ZGERC( N-M, KBT, -CONE, X( M+1, I ), 1,
 387      $                        BB( KB1-KBT, I ), 1, X( M+1, I-KBT ),
 388      $                        LDX )
 389             END IF
 390 *
 391 *           store a(i,i1) in RA1 for use in next loop over K
 392 *
 393             RA1 = AB( I-I1+KA1, I1 )
 394          END IF
 395 *
 396 *        Generate and apply vectors of rotations to chase all the
 397 *        existing bulges KA positions down toward the bottom of the
 398 *        band
 399 *
 400          DO 130 K = 1, KB - 1
 401             IF( UPDATE ) THEN
 402 *
 403 *              Determine the rotations which would annihilate the bulge
 404 *              which has in theory just been created
 405 *
 406                IF( I-K+KA.LT.N .AND. I-K.GT.1 ) THEN
 407 *
 408 *                 generate rotation to annihilate a(i,i-k+ka+1)
 409 *
 410                   CALL ZLARTG( AB( K+1, I-K+KA ), RA1,
 411      $                         RWORK( I-K+KA-M ), WORK( I-K+KA-M ), RA )
 412 *
 413 *                 create nonzero element a(i-k,i-k+ka+1) outside the
 414 *                 band and store it in WORK(i-k)
 415 *
 416                   T = -BB( KB1-K, I )*RA1
 417                   WORK( I-K ) = RWORK( I-K+KA-M )*T -
 418      $                          DCONJG( WORK( I-K+KA-M ) )*
 419      $                          AB( 1, I-K+KA )
 420                   AB( 1, I-K+KA ) = WORK( I-K+KA-M )*T +
 421      $                              RWORK( I-K+KA-M )*AB( 1, I-K+KA )
 422                   RA1 = RA
 423                END IF
 424             END IF
 425             J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
 426             NR = ( N-J2+KA ) / KA1
 427             J1 = J2 + ( NR-1 )*KA1
 428             IF( UPDATE ) THEN
 429                J2T = MAX( J2, I+2*KA-K+1 )
 430             ELSE
 431                J2T = J2
 432             END IF
 433             NRT = ( N-J2T+KA ) / KA1
 434             DO 90 J = J2T, J1, KA1
 435 *
 436 *              create nonzero element a(j-ka,j+1) outside the band
 437 *              and store it in WORK(j-m)
 438 *
 439                WORK( J-M ) = WORK( J-M )*AB( 1, J+1 )
 440                AB( 1, J+1 ) = RWORK( J-M )*AB( 1, J+1 )
 441    90       CONTINUE
 442 *
 443 *           generate rotations in 1st set to annihilate elements which
 444 *           have been created outside the band
 445 *
 446             IF( NRT.GT.0 )
 447      $         CALL ZLARGV( NRT, AB( 1, J2T ), INCA, WORK( J2T-M ), KA1,
 448      $                      RWORK( J2T-M ), KA1 )
 449             IF( NR.GT.0 ) THEN
 450 *
 451 *              apply rotations in 1st set from the right
 452 *
 453                DO 100 L = 1, KA - 1
 454                   CALL ZLARTV( NR, AB( KA1-L, J2 ), INCA,
 455      $                         AB( KA-L, J2+1 ), INCA, RWORK( J2-M ),
 456      $                         WORK( J2-M ), KA1 )
 457   100          CONTINUE
 458 *
 459 *              apply rotations in 1st set from both sides to diagonal
 460 *              blocks
 461 *
 462                CALL ZLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ),
 463      $                      AB( KA, J2+1 ), INCA, RWORK( J2-M ),
 464      $                      WORK( J2-M ), KA1 )
 465 *
 466                CALL ZLACGV( NR, WORK( J2-M ), KA1 )
 467             END IF
 468 *
 469 *           start applying rotations in 1st set from the left
 470 *
 471             DO 110 L = KA - 1, KB - K + 1, -1
 472                NRT = ( N-J2+L ) / KA1
 473                IF( NRT.GT.0 )
 474      $            CALL ZLARTV( NRT, AB( L, J2+KA1-L ), INCA,
 475      $                         AB( L+1, J2+KA1-L ), INCA, RWORK( J2-M ),
 476      $                         WORK( J2-M ), KA1 )
 477   110       CONTINUE
 478 *
 479             IF( WANTX ) THEN
 480 *
 481 *              post-multiply X by product of rotations in 1st set
 482 *
 483                DO 120 J = J2, J1, KA1
 484                   CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
 485      $                       RWORK( J-M ), DCONJG( WORK( J-M ) ) )
 486   120          CONTINUE
 487             END IF
 488   130    CONTINUE
 489 *
 490          IF( UPDATE ) THEN
 491             IF( I2.LE.N .AND. KBT.GT.0 ) THEN
 492 *
 493 *              create nonzero element a(i-kbt,i-kbt+ka+1) outside the
 494 *              band and store it in WORK(i-kbt)
 495 *
 496                WORK( I-KBT ) = -BB( KB1-KBT, I )*RA1
 497             END IF
 498          END IF
 499 *
 500          DO 170 K = KB, 1, -1
 501             IF( UPDATE ) THEN
 502                J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1
 503             ELSE
 504                J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
 505             END IF
 506 *
 507 *           finish applying rotations in 2nd set from the left
 508 *
 509             DO 140 L = KB - K, 1, -1
 510                NRT = ( N-J2+KA+L ) / KA1
 511                IF( NRT.GT.0 )
 512      $            CALL ZLARTV( NRT, AB( L, J2-L+1 ), INCA,
 513      $                         AB( L+1, J2-L+1 ), INCA, RWORK( J2-KA ),
 514      $                         WORK( J2-KA ), KA1 )
 515   140       CONTINUE
 516             NR = ( N-J2+KA ) / KA1
 517             J1 = J2 + ( NR-1 )*KA1
 518             DO 150 J = J1, J2, -KA1
 519                WORK( J ) = WORK( J-KA )
 520                RWORK( J ) = RWORK( J-KA )
 521   150       CONTINUE
 522             DO 160 J = J2, J1, KA1
 523 *
 524 *              create nonzero element a(j-ka,j+1) outside the band
 525 *              and store it in WORK(j)
 526 *
 527                WORK( J ) = WORK( J )*AB( 1, J+1 )
 528                AB( 1, J+1 ) = RWORK( J )*AB( 1, J+1 )
 529   160       CONTINUE
 530             IF( UPDATE ) THEN
 531                IF( I-K.LT.N-KA .AND. K.LE.KBT )
 532      $            WORK( I-K+KA ) = WORK( I-K )
 533             END IF
 534   170    CONTINUE
 535 *
 536          DO 210 K = KB, 1, -1
 537             J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
 538             NR = ( N-J2+KA ) / KA1
 539             J1 = J2 + ( NR-1 )*KA1
 540             IF( NR.GT.0 ) THEN
 541 *
 542 *              generate rotations in 2nd set to annihilate elements
 543 *              which have been created outside the band
 544 *
 545                CALL ZLARGV( NR, AB( 1, J2 ), INCA, WORK( J2 ), KA1,
 546      $                      RWORK( J2 ), KA1 )
 547 *
 548 *              apply rotations in 2nd set from the right
 549 *
 550                DO 180 L = 1, KA - 1
 551                   CALL ZLARTV( NR, AB( KA1-L, J2 ), INCA,
 552      $                         AB( KA-L, J2+1 ), INCA, RWORK( J2 ),
 553      $                         WORK( J2 ), KA1 )
 554   180          CONTINUE
 555 *
 556 *              apply rotations in 2nd set from both sides to diagonal
 557 *              blocks
 558 *
 559                CALL ZLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ),
 560      $                      AB( KA, J2+1 ), INCA, RWORK( J2 ),
 561      $                      WORK( J2 ), KA1 )
 562 *
 563                CALL ZLACGV( NR, WORK( J2 ), KA1 )
 564             END IF
 565 *
 566 *           start applying rotations in 2nd set from the left
 567 *
 568             DO 190 L = KA - 1, KB - K + 1, -1
 569                NRT = ( N-J2+L ) / KA1
 570                IF( NRT.GT.0 )
 571      $            CALL ZLARTV( NRT, AB( L, J2+KA1-L ), INCA,
 572      $                         AB( L+1, J2+KA1-L ), INCA, RWORK( J2 ),
 573      $                         WORK( J2 ), KA1 )
 574   190       CONTINUE
 575 *
 576             IF( WANTX ) THEN
 577 *
 578 *              post-multiply X by product of rotations in 2nd set
 579 *
 580                DO 200 J = J2, J1, KA1
 581                   CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
 582      $                       RWORK( J ), DCONJG( WORK( J ) ) )
 583   200          CONTINUE
 584             END IF
 585   210    CONTINUE
 586 *
 587          DO 230 K = 1, KB - 1
 588             J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
 589 *
 590 *           finish applying rotations in 1st set from the left
 591 *
 592             DO 220 L = KB - K, 1, -1
 593                NRT = ( N-J2+L ) / KA1
 594                IF( NRT.GT.0 )
 595      $            CALL ZLARTV( NRT, AB( L, J2+KA1-L ), INCA,
 596      $                         AB( L+1, J2+KA1-L ), INCA, RWORK( J2-M ),
 597      $                         WORK( J2-M ), KA1 )
 598   220       CONTINUE
 599   230    CONTINUE
 600 *
 601          IF( KB.GT.1 ) THEN
 602             DO 240 J = N - 1, J2 + KA, -1
 603                RWORK( J-M ) = RWORK( J-KA-M )
 604                WORK( J-M ) = WORK( J-KA-M )
 605   240       CONTINUE
 606          END IF
 607 *
 608       ELSE
 609 *
 610 *        Transform A, working with the lower triangle
 611 *
 612          IF( UPDATE ) THEN
 613 *
 614 *           Form  inv(S(i))**H * A * inv(S(i))
 615 *
 616             BII = DBLE( BB( 1, I ) )
 617             AB( 1, I ) = ( DBLE( AB( 1, I ) ) / BII ) / BII
 618             DO 250 J = I + 1, I1
 619                AB( J-I+1, I ) = AB( J-I+1, I ) / BII
 620   250       CONTINUE
 621             DO 260 J = MAX( 1, I-KA ), I - 1
 622                AB( I-J+1, J ) = AB( I-J+1, J ) / BII
 623   260       CONTINUE
 624             DO 290 K = I - KBT, I - 1
 625                DO 270 J = I - KBT, K
 626                   AB( K-J+1, J ) = AB( K-J+1, J ) -
 627      $                             BB( I-J+1, J )*DCONJG( AB( I-K+1,
 628      $                             K ) ) - DCONJG( BB( I-K+1, K ) )*
 629      $                             AB( I-J+1, J ) + DBLE( AB( 1, I ) )*
 630      $                             BB( I-J+1, J )*DCONJG( BB( I-K+1,
 631      $                             K ) )
 632   270          CONTINUE
 633                DO 280 J = MAX( 1, I-KA ), I - KBT - 1
 634                   AB( K-J+1, J ) = AB( K-J+1, J ) -
 635      $                             DCONJG( BB( I-K+1, K ) )*
 636      $                             AB( I-J+1, J )
 637   280          CONTINUE
 638   290       CONTINUE
 639             DO 310 J = I, I1
 640                DO 300 K = MAX( J-KA, I-KBT ), I - 1
 641                   AB( J-K+1, K ) = AB( J-K+1, K ) -
 642      $                             BB( I-K+1, K )*AB( J-I+1, I )
 643   300          CONTINUE
 644   310       CONTINUE
 645 *
 646             IF( WANTX ) THEN
 647 *
 648 *              post-multiply X by inv(S(i))
 649 *
 650                CALL ZDSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
 651                IF( KBT.GT.0 )
 652      $            CALL ZGERU( N-M, KBT, -CONE, X( M+1, I ), 1,
 653      $                        BB( KBT+1, I-KBT ), LDBB-1,
 654      $                        X( M+1, I-KBT ), LDX )
 655             END IF
 656 *
 657 *           store a(i1,i) in RA1 for use in next loop over K
 658 *
 659             RA1 = AB( I1-I+1, I )
 660          END IF
 661 *
 662 *        Generate and apply vectors of rotations to chase all the
 663 *        existing bulges KA positions down toward the bottom of the
 664 *        band
 665 *
 666          DO 360 K = 1, KB - 1
 667             IF( UPDATE ) THEN
 668 *
 669 *              Determine the rotations which would annihilate the bulge
 670 *              which has in theory just been created
 671 *
 672                IF( I-K+KA.LT.N .AND. I-K.GT.1 ) THEN
 673 *
 674 *                 generate rotation to annihilate a(i-k+ka+1,i)
 675 *
 676                   CALL ZLARTG( AB( KA1-K, I ), RA1, RWORK( I-K+KA-M ),
 677      $                         WORK( I-K+KA-M ), RA )
 678 *
 679 *                 create nonzero element a(i-k+ka+1,i-k) outside the
 680 *                 band and store it in WORK(i-k)
 681 *
 682                   T = -BB( K+1, I-K )*RA1
 683                   WORK( I-K ) = RWORK( I-K+KA-M )*T -
 684      $                          DCONJG( WORK( I-K+KA-M ) )*
 685      $                          AB( KA1, I-K )
 686                   AB( KA1, I-K ) = WORK( I-K+KA-M )*T +
 687      $                             RWORK( I-K+KA-M )*AB( KA1, I-K )
 688                   RA1 = RA
 689                END IF
 690             END IF
 691             J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
 692             NR = ( N-J2+KA ) / KA1
 693             J1 = J2 + ( NR-1 )*KA1
 694             IF( UPDATE ) THEN
 695                J2T = MAX( J2, I+2*KA-K+1 )
 696             ELSE
 697                J2T = J2
 698             END IF
 699             NRT = ( N-J2T+KA ) / KA1
 700             DO 320 J = J2T, J1, KA1
 701 *
 702 *              create nonzero element a(j+1,j-ka) outside the band
 703 *              and store it in WORK(j-m)
 704 *
 705                WORK( J-M ) = WORK( J-M )*AB( KA1, J-KA+1 )
 706                AB( KA1, J-KA+1 ) = RWORK( J-M )*AB( KA1, J-KA+1 )
 707   320       CONTINUE
 708 *
 709 *           generate rotations in 1st set to annihilate elements which
 710 *           have been created outside the band
 711 *
 712             IF( NRT.GT.0 )
 713      $         CALL ZLARGV( NRT, AB( KA1, J2T-KA ), INCA, WORK( J2T-M ),
 714      $                      KA1, RWORK( J2T-M ), KA1 )
 715             IF( NR.GT.0 ) THEN
 716 *
 717 *              apply rotations in 1st set from the left
 718 *
 719                DO 330 L = 1, KA - 1
 720                   CALL ZLARTV( NR, AB( L+1, J2-L ), INCA,
 721      $                         AB( L+2, J2-L ), INCA, RWORK( J2-M ),
 722      $                         WORK( J2-M ), KA1 )
 723   330          CONTINUE
 724 *
 725 *              apply rotations in 1st set from both sides to diagonal
 726 *              blocks
 727 *
 728                CALL ZLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
 729      $                      INCA, RWORK( J2-M ), WORK( J2-M ), KA1 )
 730 *
 731                CALL ZLACGV( NR, WORK( J2-M ), KA1 )
 732             END IF
 733 *
 734 *           start applying rotations in 1st set from the right
 735 *
 736             DO 340 L = KA - 1, KB - K + 1, -1
 737                NRT = ( N-J2+L ) / KA1
 738                IF( NRT.GT.0 )
 739      $            CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
 740      $                         AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ),
 741      $                         WORK( J2-M ), KA1 )
 742   340       CONTINUE
 743 *
 744             IF( WANTX ) THEN
 745 *
 746 *              post-multiply X by product of rotations in 1st set
 747 *
 748                DO 350 J = J2, J1, KA1
 749                   CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
 750      $                       RWORK( J-M ), WORK( J-M ) )
 751   350          CONTINUE
 752             END IF
 753   360    CONTINUE
 754 *
 755          IF( UPDATE ) THEN
 756             IF( I2.LE.N .AND. KBT.GT.0 ) THEN
 757 *
 758 *              create nonzero element a(i-kbt+ka+1,i-kbt) outside the
 759 *              band and store it in WORK(i-kbt)
 760 *
 761                WORK( I-KBT ) = -BB( KBT+1, I-KBT )*RA1
 762             END IF
 763          END IF
 764 *
 765          DO 400 K = KB, 1, -1
 766             IF( UPDATE ) THEN
 767                J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1
 768             ELSE
 769                J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
 770             END IF
 771 *
 772 *           finish applying rotations in 2nd set from the right
 773 *
 774             DO 370 L = KB - K, 1, -1
 775                NRT = ( N-J2+KA+L ) / KA1
 776                IF( NRT.GT.0 )
 777      $            CALL ZLARTV( NRT, AB( KA1-L+1, J2-KA ), INCA,
 778      $                         AB( KA1-L, J2-KA+1 ), INCA,
 779      $                         RWORK( J2-KA ), WORK( J2-KA ), KA1 )
 780   370       CONTINUE
 781             NR = ( N-J2+KA ) / KA1
 782             J1 = J2 + ( NR-1 )*KA1
 783             DO 380 J = J1, J2, -KA1
 784                WORK( J ) = WORK( J-KA )
 785                RWORK( J ) = RWORK( J-KA )
 786   380       CONTINUE
 787             DO 390 J = J2, J1, KA1
 788 *
 789 *              create nonzero element a(j+1,j-ka) outside the band
 790 *              and store it in WORK(j)
 791 *
 792                WORK( J ) = WORK( J )*AB( KA1, J-KA+1 )
 793                AB( KA1, J-KA+1 ) = RWORK( J )*AB( KA1, J-KA+1 )
 794   390       CONTINUE
 795             IF( UPDATE ) THEN
 796                IF( I-K.LT.N-KA .AND. K.LE.KBT )
 797      $            WORK( I-K+KA ) = WORK( I-K )
 798             END IF
 799   400    CONTINUE
 800 *
 801          DO 440 K = KB, 1, -1
 802             J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
 803             NR = ( N-J2+KA ) / KA1
 804             J1 = J2 + ( NR-1 )*KA1
 805             IF( NR.GT.0 ) THEN
 806 *
 807 *              generate rotations in 2nd set to annihilate elements
 808 *              which have been created outside the band
 809 *
 810                CALL ZLARGV( NR, AB( KA1, J2-KA ), INCA, WORK( J2 ), KA1,
 811      $                      RWORK( J2 ), KA1 )
 812 *
 813 *              apply rotations in 2nd set from the left
 814 *
 815                DO 410 L = 1, KA - 1
 816                   CALL ZLARTV( NR, AB( L+1, J2-L ), INCA,
 817      $                         AB( L+2, J2-L ), INCA, RWORK( J2 ),
 818      $                         WORK( J2 ), KA1 )
 819   410          CONTINUE
 820 *
 821 *              apply rotations in 2nd set from both sides to diagonal
 822 *              blocks
 823 *
 824                CALL ZLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
 825      $                      INCA, RWORK( J2 ), WORK( J2 ), KA1 )
 826 *
 827                CALL ZLACGV( NR, WORK( J2 ), KA1 )
 828             END IF
 829 *
 830 *           start applying rotations in 2nd set from the right
 831 *
 832             DO 420 L = KA - 1, KB - K + 1, -1
 833                NRT = ( N-J2+L ) / KA1
 834                IF( NRT.GT.0 )
 835      $            CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
 836      $                         AB( KA1-L, J2+1 ), INCA, RWORK( J2 ),
 837      $                         WORK( J2 ), KA1 )
 838   420       CONTINUE
 839 *
 840             IF( WANTX ) THEN
 841 *
 842 *              post-multiply X by product of rotations in 2nd set
 843 *
 844                DO 430 J = J2, J1, KA1
 845                   CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
 846      $                       RWORK( J ), WORK( J ) )
 847   430          CONTINUE
 848             END IF
 849   440    CONTINUE
 850 *
 851          DO 460 K = 1, KB - 1
 852             J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
 853 *
 854 *           finish applying rotations in 1st set from the right
 855 *
 856             DO 450 L = KB - K, 1, -1
 857                NRT = ( N-J2+L ) / KA1
 858                IF( NRT.GT.0 )
 859      $            CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
 860      $                         AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ),
 861      $                         WORK( J2-M ), KA1 )
 862   450       CONTINUE
 863   460    CONTINUE
 864 *
 865          IF( KB.GT.1 ) THEN
 866             DO 470 J = N - 1, J2 + KA, -1
 867                RWORK( J-M ) = RWORK( J-KA-M )
 868                WORK( J-M ) = WORK( J-KA-M )
 869   470       CONTINUE
 870          END IF
 871 *
 872       END IF
 873 *
 874       GO TO 10
 875 *
 876   480 CONTINUE
 877 *
 878 *     **************************** Phase 2 *****************************
 879 *
 880 *     The logical structure of this phase is:
 881 *
 882 *     UPDATE = .TRUE.
 883 *     DO I = 1, M
 884 *        use S(i) to update A and create a new bulge
 885 *        apply rotations to push all bulges KA positions upward
 886 *     END DO
 887 *     UPDATE = .FALSE.
 888 *     DO I = M - KA - 1, 2, -1
 889 *        apply rotations to push all bulges KA positions upward
 890 *     END DO
 891 *
 892 *     To avoid duplicating code, the two loops are merged.
 893 *
 894       UPDATE = .TRUE.
 895       I = 0
 896   490 CONTINUE
 897       IF( UPDATE ) THEN
 898          I = I + 1
 899          KBT = MIN( KB, M-I )
 900          I0 = I + 1
 901          I1 = MAX( 1, I-KA )
 902          I2 = I + KBT - KA1
 903          IF( I.GT.M ) THEN
 904             UPDATE = .FALSE.
 905             I = I - 1
 906             I0 = M + 1
 907             IF( KA.EQ.0 )
 908      $         RETURN
 909             GO TO 490
 910          END IF
 911       ELSE
 912          I = I - KA
 913          IF( I.LT.2 )
 914      $      RETURN
 915       END IF
 916 *
 917       IF( I.LT.M-KBT ) THEN
 918          NX = M
 919       ELSE
 920          NX = N
 921       END IF
 922 *
 923       IF( UPPER ) THEN
 924 *
 925 *        Transform A, working with the upper triangle
 926 *
 927          IF( UPDATE ) THEN
 928 *
 929 *           Form  inv(S(i))**H * A * inv(S(i))
 930 *
 931             BII = DBLE( BB( KB1, I ) )
 932             AB( KA1, I ) = ( DBLE( AB( KA1, I ) ) / BII ) / BII
 933             DO 500 J = I1, I - 1
 934                AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII
 935   500       CONTINUE
 936             DO 510 J = I + 1, MIN( N, I+KA )
 937                AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII
 938   510       CONTINUE
 939             DO 540 K = I + 1, I + KBT
 940                DO 520 J = K, I + KBT
 941                   AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
 942      $                               BB( I-J+KB1, J )*
 943      $                               DCONJG( AB( I-K+KA1, K ) ) -
 944      $                               DCONJG( BB( I-K+KB1, K ) )*
 945      $                               AB( I-J+KA1, J ) +
 946      $                               DBLE( AB( KA1, I ) )*
 947      $                               BB( I-J+KB1, J )*
 948      $                               DCONJG( BB( I-K+KB1, K ) )
 949   520          CONTINUE
 950                DO 530 J = I + KBT + 1, MIN( N, I+KA )
 951                   AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
 952      $                               DCONJG( BB( I-K+KB1, K ) )*
 953      $                               AB( I-J+KA1, J )
 954   530          CONTINUE
 955   540       CONTINUE
 956             DO 560 J = I1, I
 957                DO 550 K = I + 1, MIN( J+KA, I+KBT )
 958                   AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
 959      $                               BB( I-K+KB1, K )*AB( J-I+KA1, I )
 960   550          CONTINUE
 961   560       CONTINUE
 962 *
 963             IF( WANTX ) THEN
 964 *
 965 *              post-multiply X by inv(S(i))
 966 *
 967                CALL ZDSCAL( NX, ONE / BII, X( 1, I ), 1 )
 968                IF( KBT.GT.0 )
 969      $            CALL ZGERU( NX, KBT, -CONE, X( 1, I ), 1,
 970      $                        BB( KB, I+1 ), LDBB-1, X( 1, I+1 ), LDX )
 971             END IF
 972 *
 973 *           store a(i1,i) in RA1 for use in next loop over K
 974 *
 975             RA1 = AB( I1-I+KA1, I )
 976          END IF
 977 *
 978 *        Generate and apply vectors of rotations to chase all the
 979 *        existing bulges KA positions up toward the top of the band
 980 *
 981          DO 610 K = 1, KB - 1
 982             IF( UPDATE ) THEN
 983 *
 984 *              Determine the rotations which would annihilate the bulge
 985 *              which has in theory just been created
 986 *
 987                IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
 988 *
 989 *                 generate rotation to annihilate a(i+k-ka-1,i)
 990 *
 991                   CALL ZLARTG( AB( K+1, I ), RA1, RWORK( I+K-KA ),
 992      $                         WORK( I+K-KA ), RA )
 993 *
 994 *                 create nonzero element a(i+k-ka-1,i+k) outside the
 995 *                 band and store it in WORK(m-kb+i+k)
 996 *
 997                   T = -BB( KB1-K, I+K )*RA1
 998                   WORK( M-KB+I+K ) = RWORK( I+K-KA )*T -
 999      $                               DCONJG( WORK( I+K-KA ) )*
1000      $                               AB( 1, I+K )
1001                   AB( 1, I+K ) = WORK( I+K-KA )*T +
1002      $                           RWORK( I+K-KA )*AB( 1, I+K )
1003                   RA1 = RA
1004                END IF
1005             END IF
1006             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
1007             NR = ( J2+KA-1 ) / KA1
1008             J1 = J2 - ( NR-1 )*KA1
1009             IF( UPDATE ) THEN
1010                J2T = MIN( J2, I-2*KA+K-1 )
1011             ELSE
1012                J2T = J2
1013             END IF
1014             NRT = ( J2T+KA-1 ) / KA1
1015             DO 570 J = J1, J2T, KA1
1016 *
1017 *              create nonzero element a(j-1,j+ka) outside the band
1018 *              and store it in WORK(j)
1019 *
1020                WORK( J ) = WORK( J )*AB( 1, J+KA-1 )
1021                AB( 1, J+KA-1 ) = RWORK( J )*AB( 1, J+KA-1 )
1022   570       CONTINUE
1023 *
1024 *           generate rotations in 1st set to annihilate elements which
1025 *           have been created outside the band
1026 *
1027             IF( NRT.GT.0 )
1028      $         CALL ZLARGV( NRT, AB( 1, J1+KA ), INCA, WORK( J1 ), KA1,
1029      $                      RWORK( J1 ), KA1 )
1030             IF( NR.GT.0 ) THEN
1031 *
1032 *              apply rotations in 1st set from the left
1033 *
1034                DO 580 L = 1, KA - 1
1035                   CALL ZLARTV( NR, AB( KA1-L, J1+L ), INCA,
1036      $                         AB( KA-L, J1+L ), INCA, RWORK( J1 ),
1037      $                         WORK( J1 ), KA1 )
1038   580          CONTINUE
1039 *
1040 *              apply rotations in 1st set from both sides to diagonal
1041 *              blocks
1042 *
1043                CALL ZLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
1044      $                      AB( KA, J1 ), INCA, RWORK( J1 ), WORK( J1 ),
1045      $                      KA1 )
1046 *
1047                CALL ZLACGV( NR, WORK( J1 ), KA1 )
1048             END IF
1049 *
1050 *           start applying rotations in 1st set from the right
1051 *
1052             DO 590 L = KA - 1, KB - K + 1, -1
1053                NRT = ( J2+L-1 ) / KA1
1054                J1T = J2 - ( NRT-1 )*KA1
1055                IF( NRT.GT.0 )
1056      $            CALL ZLARTV( NRT, AB( L, J1T ), INCA,
1057      $                         AB( L+1, J1T-1 ), INCA, RWORK( J1T ),
1058      $                         WORK( J1T ), KA1 )
1059   590       CONTINUE
1060 *
1061             IF( WANTX ) THEN
1062 *
1063 *              post-multiply X by product of rotations in 1st set
1064 *
1065                DO 600 J = J1, J2, KA1
1066                   CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1067      $                       RWORK( J ), WORK( J ) )
1068   600          CONTINUE
1069             END IF
1070   610    CONTINUE
1071 *
1072          IF( UPDATE ) THEN
1073             IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
1074 *
1075 *              create nonzero element a(i+kbt-ka-1,i+kbt) outside the
1076 *              band and store it in WORK(m-kb+i+kbt)
1077 *
1078                WORK( M-KB+I+KBT ) = -BB( KB1-KBT, I+KBT )*RA1
1079             END IF
1080          END IF
1081 *
1082          DO 650 K = KB, 1, -1
1083             IF( UPDATE ) THEN
1084                J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1
1085             ELSE
1086                J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1087             END IF
1088 *
1089 *           finish applying rotations in 2nd set from the right
1090 *
1091             DO 620 L = KB - K, 1, -1
1092                NRT = ( J2+KA+L-1 ) / KA1
1093                J1T = J2 - ( NRT-1 )*KA1
1094                IF( NRT.GT.0 )
1095      $            CALL ZLARTV( NRT, AB( L, J1T+KA ), INCA,
1096      $                         AB( L+1, J1T+KA-1 ), INCA,
1097      $                         RWORK( M-KB+J1T+KA ),
1098      $                         WORK( M-KB+J1T+KA ), KA1 )
1099   620       CONTINUE
1100             NR = ( J2+KA-1 ) / KA1
1101             J1 = J2 - ( NR-1 )*KA1
1102             DO 630 J = J1, J2, KA1
1103                WORK( M-KB+J ) = WORK( M-KB+J+KA )
1104                RWORK( M-KB+J ) = RWORK( M-KB+J+KA )
1105   630       CONTINUE
1106             DO 640 J = J1, J2, KA1
1107 *
1108 *              create nonzero element a(j-1,j+ka) outside the band
1109 *              and store it in WORK(m-kb+j)
1110 *
1111                WORK( M-KB+J ) = WORK( M-KB+J )*AB( 1, J+KA-1 )
1112                AB( 1, J+KA-1 ) = RWORK( M-KB+J )*AB( 1, J+KA-1 )
1113   640       CONTINUE
1114             IF( UPDATE ) THEN
1115                IF( I+K.GT.KA1 .AND. K.LE.KBT )
1116      $            WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
1117             END IF
1118   650    CONTINUE
1119 *
1120          DO 690 K = KB, 1, -1
1121             J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1122             NR = ( J2+KA-1 ) / KA1
1123             J1 = J2 - ( NR-1 )*KA1
1124             IF( NR.GT.0 ) THEN
1125 *
1126 *              generate rotations in 2nd set to annihilate elements
1127 *              which have been created outside the band
1128 *
1129                CALL ZLARGV( NR, AB( 1, J1+KA ), INCA, WORK( M-KB+J1 ),
1130      $                      KA1, RWORK( M-KB+J1 ), KA1 )
1131 *
1132 *              apply rotations in 2nd set from the left
1133 *
1134                DO 660 L = 1, KA - 1
1135                   CALL ZLARTV( NR, AB( KA1-L, J1+L ), INCA,
1136      $                         AB( KA-L, J1+L ), INCA, RWORK( M-KB+J1 ),
1137      $                         WORK( M-KB+J1 ), KA1 )
1138   660          CONTINUE
1139 *
1140 *              apply rotations in 2nd set from both sides to diagonal
1141 *              blocks
1142 *
1143                CALL ZLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
1144      $                      AB( KA, J1 ), INCA, RWORK( M-KB+J1 ),
1145      $                      WORK( M-KB+J1 ), KA1 )
1146 *
1147                CALL ZLACGV( NR, WORK( M-KB+J1 ), KA1 )
1148             END IF
1149 *
1150 *           start applying rotations in 2nd set from the right
1151 *
1152             DO 670 L = KA - 1, KB - K + 1, -1
1153                NRT = ( J2+L-1 ) / KA1
1154                J1T = J2 - ( NRT-1 )*KA1
1155                IF( NRT.GT.0 )
1156      $            CALL ZLARTV( NRT, AB( L, J1T ), INCA,
1157      $                         AB( L+1, J1T-1 ), INCA,
1158      $                         RWORK( M-KB+J1T ), WORK( M-KB+J1T ),
1159      $                         KA1 )
1160   670       CONTINUE
1161 *
1162             IF( WANTX ) THEN
1163 *
1164 *              post-multiply X by product of rotations in 2nd set
1165 *
1166                DO 680 J = J1, J2, KA1
1167                   CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1168      $                       RWORK( M-KB+J ), WORK( M-KB+J ) )
1169   680          CONTINUE
1170             END IF
1171   690    CONTINUE
1172 *
1173          DO 710 K = 1, KB - 1
1174             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
1175 *
1176 *           finish applying rotations in 1st set from the right
1177 *
1178             DO 700 L = KB - K, 1, -1
1179                NRT = ( J2+L-1 ) / KA1
1180                J1T = J2 - ( NRT-1 )*KA1
1181                IF( NRT.GT.0 )
1182      $            CALL ZLARTV( NRT, AB( L, J1T ), INCA,
1183      $                         AB( L+1, J1T-1 ), INCA, RWORK( J1T ),
1184      $                         WORK( J1T ), KA1 )
1185   700       CONTINUE
1186   710    CONTINUE
1187 *
1188          IF( KB.GT.1 ) THEN
1189             DO 720 J = 2, I2 - KA
1190                RWORK( J ) = RWORK( J+KA )
1191                WORK( J ) = WORK( J+KA )
1192   720       CONTINUE
1193          END IF
1194 *
1195       ELSE
1196 *
1197 *        Transform A, working with the lower triangle
1198 *
1199          IF( UPDATE ) THEN
1200 *
1201 *           Form  inv(S(i))**H * A * inv(S(i))
1202 *
1203             BII = DBLE( BB( 1, I ) )
1204             AB( 1, I ) = ( DBLE( AB( 1, I ) ) / BII ) / BII
1205             DO 730 J = I1, I - 1
1206                AB( I-J+1, J ) = AB( I-J+1, J ) / BII
1207   730       CONTINUE
1208             DO 740 J = I + 1, MIN( N, I+KA )
1209                AB( J-I+1, I ) = AB( J-I+1, I ) / BII
1210   740       CONTINUE
1211             DO 770 K = I + 1, I + KBT
1212                DO 750 J = K, I + KBT
1213                   AB( J-K+1, K ) = AB( J-K+1, K ) -
1214      $                             BB( J-I+1, I )*DCONJG( AB( K-I+1,
1215      $                             I ) ) - DCONJG( BB( K-I+1, I ) )*
1216      $                             AB( J-I+1, I ) + DBLE( AB( 1, I ) )*
1217      $                             BB( J-I+1, I )*DCONJG( BB( K-I+1,
1218      $                             I ) )
1219   750          CONTINUE
1220                DO 760 J = I + KBT + 1, MIN( N, I+KA )
1221                   AB( J-K+1, K ) = AB( J-K+1, K ) -
1222      $                             DCONJG( BB( K-I+1, I ) )*
1223      $                             AB( J-I+1, I )
1224   760          CONTINUE
1225   770       CONTINUE
1226             DO 790 J = I1, I
1227                DO 780 K = I + 1, MIN( J+KA, I+KBT )
1228                   AB( K-J+1, J ) = AB( K-J+1, J ) -
1229      $                             BB( K-I+1, I )*AB( I-J+1, J )
1230   780          CONTINUE
1231   790       CONTINUE
1232 *
1233             IF( WANTX ) THEN
1234 *
1235 *              post-multiply X by inv(S(i))
1236 *
1237                CALL ZDSCAL( NX, ONE / BII, X( 1, I ), 1 )
1238                IF( KBT.GT.0 )
1239      $            CALL ZGERC( NX, KBT, -CONE, X( 1, I ), 1, BB( 2, I ),
1240      $                        1, X( 1, I+1 ), LDX )
1241             END IF
1242 *
1243 *           store a(i,i1) in RA1 for use in next loop over K
1244 *
1245             RA1 = AB( I-I1+1, I1 )
1246          END IF
1247 *
1248 *        Generate and apply vectors of rotations to chase all the
1249 *        existing bulges KA positions up toward the top of the band
1250 *
1251          DO 840 K = 1, KB - 1
1252             IF( UPDATE ) THEN
1253 *
1254 *              Determine the rotations which would annihilate the bulge
1255 *              which has in theory just been created
1256 *
1257                IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
1258 *
1259 *                 generate rotation to annihilate a(i,i+k-ka-1)
1260 *
1261                   CALL ZLARTG( AB( KA1-K, I+K-KA ), RA1,
1262      $                         RWORK( I+K-KA ), WORK( I+K-KA ), RA )
1263 *
1264 *                 create nonzero element a(i+k,i+k-ka-1) outside the
1265 *                 band and store it in WORK(m-kb+i+k)
1266 *
1267                   T = -BB( K+1, I )*RA1
1268                   WORK( M-KB+I+K ) = RWORK( I+K-KA )*T -
1269      $                               DCONJG( WORK( I+K-KA ) )*
1270      $                               AB( KA1, I+K-KA )
1271                   AB( KA1, I+K-KA ) = WORK( I+K-KA )*T +
1272      $                                RWORK( I+K-KA )*AB( KA1, I+K-KA )
1273                   RA1 = RA
1274                END IF
1275             END IF
1276             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
1277             NR = ( J2+KA-1 ) / KA1
1278             J1 = J2 - ( NR-1 )*KA1
1279             IF( UPDATE ) THEN
1280                J2T = MIN( J2, I-2*KA+K-1 )
1281             ELSE
1282                J2T = J2
1283             END IF
1284             NRT = ( J2T+KA-1 ) / KA1
1285             DO 800 J = J1, J2T, KA1
1286 *
1287 *              create nonzero element a(j+ka,j-1) outside the band
1288 *              and store it in WORK(j)
1289 *
1290                WORK( J ) = WORK( J )*AB( KA1, J-1 )
1291                AB( KA1, J-1 ) = RWORK( J )*AB( KA1, J-1 )
1292   800       CONTINUE
1293 *
1294 *           generate rotations in 1st set to annihilate elements which
1295 *           have been created outside the band
1296 *
1297             IF( NRT.GT.0 )
1298      $         CALL ZLARGV( NRT, AB( KA1, J1 ), INCA, WORK( J1 ), KA1,
1299      $                      RWORK( J1 ), KA1 )
1300             IF( NR.GT.0 ) THEN
1301 *
1302 *              apply rotations in 1st set from the right
1303 *
1304                DO 810 L = 1, KA - 1
1305                   CALL ZLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
1306      $                         INCA, RWORK( J1 ), WORK( J1 ), KA1 )
1307   810          CONTINUE
1308 *
1309 *              apply rotations in 1st set from both sides to diagonal
1310 *              blocks
1311 *
1312                CALL ZLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
1313      $                      AB( 2, J1-1 ), INCA, RWORK( J1 ),
1314      $                      WORK( J1 ), KA1 )
1315 *
1316                CALL ZLACGV( NR, WORK( J1 ), KA1 )
1317             END IF
1318 *
1319 *           start applying rotations in 1st set from the left
1320 *
1321             DO 820 L = KA - 1, KB - K + 1, -1
1322                NRT = ( J2+L-1 ) / KA1
1323                J1T = J2 - ( NRT-1 )*KA1
1324                IF( NRT.GT.0 )
1325      $            CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
1326      $                         AB( KA1-L, J1T-KA1+L ), INCA,
1327      $                         RWORK( J1T ), WORK( J1T ), KA1 )
1328   820       CONTINUE
1329 *
1330             IF( WANTX ) THEN
1331 *
1332 *              post-multiply X by product of rotations in 1st set
1333 *
1334                DO 830 J = J1, J2, KA1
1335                   CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1336      $                       RWORK( J ), DCONJG( WORK( J ) ) )
1337   830          CONTINUE
1338             END IF
1339   840    CONTINUE
1340 *
1341          IF( UPDATE ) THEN
1342             IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
1343 *
1344 *              create nonzero element a(i+kbt,i+kbt-ka-1) outside the
1345 *              band and store it in WORK(m-kb+i+kbt)
1346 *
1347                WORK( M-KB+I+KBT ) = -BB( KBT+1, I )*RA1
1348             END IF
1349          END IF
1350 *
1351          DO 880 K = KB, 1, -1
1352             IF( UPDATE ) THEN
1353                J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1
1354             ELSE
1355                J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1356             END IF
1357 *
1358 *           finish applying rotations in 2nd set from the left
1359 *
1360             DO 850 L = KB - K, 1, -1
1361                NRT = ( J2+KA+L-1 ) / KA1
1362                J1T = J2 - ( NRT-1 )*KA1
1363                IF( NRT.GT.0 )
1364      $            CALL ZLARTV( NRT, AB( KA1-L+1, J1T+L-1 ), INCA,
1365      $                         AB( KA1-L, J1T+L-1 ), INCA,
1366      $                         RWORK( M-KB+J1T+KA ),
1367      $                         WORK( M-KB+J1T+KA ), KA1 )
1368   850       CONTINUE
1369             NR = ( J2+KA-1 ) / KA1
1370             J1 = J2 - ( NR-1 )*KA1
1371             DO 860 J = J1, J2, KA1
1372                WORK( M-KB+J ) = WORK( M-KB+J+KA )
1373                RWORK( M-KB+J ) = RWORK( M-KB+J+KA )
1374   860       CONTINUE
1375             DO 870 J = J1, J2, KA1
1376 *
1377 *              create nonzero element a(j+ka,j-1) outside the band
1378 *              and store it in WORK(m-kb+j)
1379 *
1380                WORK( M-KB+J ) = WORK( M-KB+J )*AB( KA1, J-1 )
1381                AB( KA1, J-1 ) = RWORK( M-KB+J )*AB( KA1, J-1 )
1382   870       CONTINUE
1383             IF( UPDATE ) THEN
1384                IF( I+K.GT.KA1 .AND. K.LE.KBT )
1385      $            WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
1386             END IF
1387   880    CONTINUE
1388 *
1389          DO 920 K = KB, 1, -1
1390             J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1391             NR = ( J2+KA-1 ) / KA1
1392             J1 = J2 - ( NR-1 )*KA1
1393             IF( NR.GT.0 ) THEN
1394 *
1395 *              generate rotations in 2nd set to annihilate elements
1396 *              which have been created outside the band
1397 *
1398                CALL ZLARGV( NR, AB( KA1, J1 ), INCA, WORK( M-KB+J1 ),
1399      $                      KA1, RWORK( M-KB+J1 ), KA1 )
1400 *
1401 *              apply rotations in 2nd set from the right
1402 *
1403                DO 890 L = 1, KA - 1
1404                   CALL ZLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
1405      $                         INCA, RWORK( M-KB+J1 ), WORK( M-KB+J1 ),
1406      $                         KA1 )
1407   890          CONTINUE
1408 *
1409 *              apply rotations in 2nd set from both sides to diagonal
1410 *              blocks
1411 *
1412                CALL ZLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
1413      $                      AB( 2, J1-1 ), INCA, RWORK( M-KB+J1 ),
1414      $                      WORK( M-KB+J1 ), KA1 )
1415 *
1416                CALL ZLACGV( NR, WORK( M-KB+J1 ), KA1 )
1417             END IF
1418 *
1419 *           start applying rotations in 2nd set from the left
1420 *
1421             DO 900 L = KA - 1, KB - K + 1, -1
1422                NRT = ( J2+L-1 ) / KA1
1423                J1T = J2 - ( NRT-1 )*KA1
1424                IF( NRT.GT.0 )
1425      $            CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
1426      $                         AB( KA1-L, J1T-KA1+L ), INCA,
1427      $                         RWORK( M-KB+J1T ), WORK( M-KB+J1T ),
1428      $                         KA1 )
1429   900       CONTINUE
1430 *
1431             IF( WANTX ) THEN
1432 *
1433 *              post-multiply X by product of rotations in 2nd set
1434 *
1435                DO 910 J = J1, J2, KA1
1436                   CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1437      $                       RWORK( M-KB+J ), DCONJG( WORK( M-KB+J ) ) )
1438   910          CONTINUE
1439             END IF
1440   920    CONTINUE
1441 *
1442          DO 940 K = 1, KB - 1
1443             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
1444 *
1445 *           finish applying rotations in 1st set from the left
1446 *
1447             DO 930 L = KB - K, 1, -1
1448                NRT = ( J2+L-1 ) / KA1
1449                J1T = J2 - ( NRT-1 )*KA1
1450                IF( NRT.GT.0 )
1451      $            CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
1452      $                         AB( KA1-L, J1T-KA1+L ), INCA,
1453      $                         RWORK( J1T ), WORK( J1T ), KA1 )
1454   930       CONTINUE
1455   940    CONTINUE
1456 *
1457          IF( KB.GT.1 ) THEN
1458             DO 950 J = 2, I2 - KA
1459                RWORK( J ) = RWORK( J+KA )
1460                WORK( J ) = WORK( J+KA )
1461   950       CONTINUE
1462          END IF
1463 *
1464       END IF
1465 *
1466       GO TO 490
1467 *
1468 *     End of ZHBGST
1469 *
1470       END