SRC/chbgst.f

   1 *> \brief \b CHBGST
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CHBGST + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chbgst.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chbgst.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chbgst.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CHBGST( 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 *       REAL               RWORK( * )
  30 *       COMPLEX            AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
  31 *      $                   X( LDX, * )
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *> CHBGST 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 CPBSTF, 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 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 array, dimension (LDBB,N)
 110 *>          The banded factor S from the split Cholesky factorization of
 111 *>          B, as returned by CPBSTF, 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 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 array, dimension (N)
 138 *> \endverbatim
 139 *>
 140 *> \param[out] RWORK
 141 *> \verbatim
 142 *>          RWORK is REAL 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 complexOTHERcomputational
 163 *
 164 *  =====================================================================
 165       SUBROUTINE CHBGST( 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       REAL               RWORK( * )
 179       COMPLEX            AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
 180      $                   X( LDX, * )
 181 *     ..
 182 *
 183 *  =====================================================================
 184 *
 185 *     .. Parameters ..
 186       COMPLEX            CZERO, CONE
 187       REAL               ONE
 188       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
 189      $                   CONE = ( 1.0E+0, 0.0E+0 ), ONE = 1.0E+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       REAL               BII
 196       COMPLEX            RA, RA1, T
 197 *     ..
 198 *     .. External Functions ..
 199       LOGICAL            LSAME
 200       EXTERNAL           LSAME
 201 *     ..
 202 *     .. External Subroutines ..
 203       EXTERNAL           CGERC, CGERU, CLACGV, CLAR2V, CLARGV, CLARTG,
 204      $                   CLARTV, CLASET, CROT, CSSCAL, XERBLA
 205 *     ..
 206 *     .. Intrinsic Functions ..
 207       INTRINSIC          CONJG, MAX, MIN, REAL
 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( 'CHBGST', -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 CLASET( '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 CPBSTF. 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 = REAL( BB( KB1, I ) )
 349             AB( KA1, I ) = ( REAL( 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      $                               CONJG( AB( K-I+KA1, I ) ) -
 361      $                               CONJG( BB( K-I+KB1, I ) )*
 362      $                               AB( J-I+KA1, I ) +
 363      $                               REAL( AB( KA1, I ) )*
 364      $                               BB( J-I+KB1, I )*
 365      $                               CONJG( 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      $                               CONJG( 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 CSSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
 385                IF( KBT.GT.0 )
 386      $            CALL CGERC( 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 CLARTG( 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      $                          CONJG( 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 CLARGV( 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 CLARTV( 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 CLAR2V( 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 CLACGV( 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 CLARTV( 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 CROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
 485      $                       RWORK( J-M ), CONJG( 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 CLARTV( 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 CLARGV( 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 CLARTV( 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 CLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ),
 560      $                      AB( KA, J2+1 ), INCA, RWORK( J2 ),
 561      $                      WORK( J2 ), KA1 )
 562 *
 563                CALL CLACGV( 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 CLARTV( 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 CROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
 582      $                       RWORK( J ), CONJG( 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 CLARTV( 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 = REAL( BB( 1, I ) )
 617             AB( 1, I ) = ( REAL( 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 )*CONJG( AB( I-K+1,
 628      $                             K ) ) - CONJG( BB( I-K+1, K ) )*
 629      $                             AB( I-J+1, J ) + REAL( AB( 1, I ) )*
 630      $                             BB( I-J+1, J )*CONJG( 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      $                             CONJG( 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 CSSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
 651                IF( KBT.GT.0 )
 652      $            CALL CGERU( 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 CLARTG( 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      $                          CONJG( WORK( I-K+KA-M ) )*AB( KA1, I-K )
 685                   AB( KA1, I-K ) = WORK( I-K+KA-M )*T +
 686      $                             RWORK( I-K+KA-M )*AB( KA1, I-K )
 687                   RA1 = RA
 688                END IF
 689             END IF
 690             J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
 691             NR = ( N-J2+KA ) / KA1
 692             J1 = J2 + ( NR-1 )*KA1
 693             IF( UPDATE ) THEN
 694                J2T = MAX( J2, I+2*KA-K+1 )
 695             ELSE
 696                J2T = J2
 697             END IF
 698             NRT = ( N-J2T+KA ) / KA1
 699             DO 320 J = J2T, J1, KA1
 700 *
 701 *              create nonzero element a(j+1,j-ka) outside the band
 702 *              and store it in WORK(j-m)
 703 *
 704                WORK( J-M ) = WORK( J-M )*AB( KA1, J-KA+1 )
 705                AB( KA1, J-KA+1 ) = RWORK( J-M )*AB( KA1, J-KA+1 )
 706   320       CONTINUE
 707 *
 708 *           generate rotations in 1st set to annihilate elements which
 709 *           have been created outside the band
 710 *
 711             IF( NRT.GT.0 )
 712      $         CALL CLARGV( NRT, AB( KA1, J2T-KA ), INCA, WORK( J2T-M ),
 713      $                      KA1, RWORK( J2T-M ), KA1 )
 714             IF( NR.GT.0 ) THEN
 715 *
 716 *              apply rotations in 1st set from the left
 717 *
 718                DO 330 L = 1, KA - 1
 719                   CALL CLARTV( NR, AB( L+1, J2-L ), INCA,
 720      $                         AB( L+2, J2-L ), INCA, RWORK( J2-M ),
 721      $                         WORK( J2-M ), KA1 )
 722   330          CONTINUE
 723 *
 724 *              apply rotations in 1st set from both sides to diagonal
 725 *              blocks
 726 *
 727                CALL CLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
 728      $                      INCA, RWORK( J2-M ), WORK( J2-M ), KA1 )
 729 *
 730                CALL CLACGV( NR, WORK( J2-M ), KA1 )
 731             END IF
 732 *
 733 *           start applying rotations in 1st set from the right
 734 *
 735             DO 340 L = KA - 1, KB - K + 1, -1
 736                NRT = ( N-J2+L ) / KA1
 737                IF( NRT.GT.0 )
 738      $            CALL CLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
 739      $                         AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ),
 740      $                         WORK( J2-M ), KA1 )
 741   340       CONTINUE
 742 *
 743             IF( WANTX ) THEN
 744 *
 745 *              post-multiply X by product of rotations in 1st set
 746 *
 747                DO 350 J = J2, J1, KA1
 748                   CALL CROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
 749      $                       RWORK( J-M ), WORK( J-M ) )
 750   350          CONTINUE
 751             END IF
 752   360    CONTINUE
 753 *
 754          IF( UPDATE ) THEN
 755             IF( I2.LE.N .AND. KBT.GT.0 ) THEN
 756 *
 757 *              create nonzero element a(i-kbt+ka+1,i-kbt) outside the
 758 *              band and store it in WORK(i-kbt)
 759 *
 760                WORK( I-KBT ) = -BB( KBT+1, I-KBT )*RA1
 761             END IF
 762          END IF
 763 *
 764          DO 400 K = KB, 1, -1
 765             IF( UPDATE ) THEN
 766                J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1
 767             ELSE
 768                J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
 769             END IF
 770 *
 771 *           finish applying rotations in 2nd set from the right
 772 *
 773             DO 370 L = KB - K, 1, -1
 774                NRT = ( N-J2+KA+L ) / KA1
 775                IF( NRT.GT.0 )
 776      $            CALL CLARTV( NRT, AB( KA1-L+1, J2-KA ), INCA,
 777      $                         AB( KA1-L, J2-KA+1 ), INCA,
 778      $                         RWORK( J2-KA ), WORK( J2-KA ), KA1 )
 779   370       CONTINUE
 780             NR = ( N-J2+KA ) / KA1
 781             J1 = J2 + ( NR-1 )*KA1
 782             DO 380 J = J1, J2, -KA1
 783                WORK( J ) = WORK( J-KA )
 784                RWORK( J ) = RWORK( J-KA )
 785   380       CONTINUE
 786             DO 390 J = J2, J1, KA1
 787 *
 788 *              create nonzero element a(j+1,j-ka) outside the band
 789 *              and store it in WORK(j)
 790 *
 791                WORK( J ) = WORK( J )*AB( KA1, J-KA+1 )
 792                AB( KA1, J-KA+1 ) = RWORK( J )*AB( KA1, J-KA+1 )
 793   390       CONTINUE
 794             IF( UPDATE ) THEN
 795                IF( I-K.LT.N-KA .AND. K.LE.KBT )
 796      $            WORK( I-K+KA ) = WORK( I-K )
 797             END IF
 798   400    CONTINUE
 799 *
 800          DO 440 K = KB, 1, -1
 801             J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
 802             NR = ( N-J2+KA ) / KA1
 803             J1 = J2 + ( NR-1 )*KA1
 804             IF( NR.GT.0 ) THEN
 805 *
 806 *              generate rotations in 2nd set to annihilate elements
 807 *              which have been created outside the band
 808 *
 809                CALL CLARGV( NR, AB( KA1, J2-KA ), INCA, WORK( J2 ), KA1,
 810      $                      RWORK( J2 ), KA1 )
 811 *
 812 *              apply rotations in 2nd set from the left
 813 *
 814                DO 410 L = 1, KA - 1
 815                   CALL CLARTV( NR, AB( L+1, J2-L ), INCA,
 816      $                         AB( L+2, J2-L ), INCA, RWORK( J2 ),
 817      $                         WORK( J2 ), KA1 )
 818   410          CONTINUE
 819 *
 820 *              apply rotations in 2nd set from both sides to diagonal
 821 *              blocks
 822 *
 823                CALL CLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
 824      $                      INCA, RWORK( J2 ), WORK( J2 ), KA1 )
 825 *
 826                CALL CLACGV( NR, WORK( J2 ), KA1 )
 827             END IF
 828 *
 829 *           start applying rotations in 2nd set from the right
 830 *
 831             DO 420 L = KA - 1, KB - K + 1, -1
 832                NRT = ( N-J2+L ) / KA1
 833                IF( NRT.GT.0 )
 834      $            CALL CLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
 835      $                         AB( KA1-L, J2+1 ), INCA, RWORK( J2 ),
 836      $                         WORK( J2 ), KA1 )
 837   420       CONTINUE
 838 *
 839             IF( WANTX ) THEN
 840 *
 841 *              post-multiply X by product of rotations in 2nd set
 842 *
 843                DO 430 J = J2, J1, KA1
 844                   CALL CROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
 845      $                       RWORK( J ), WORK( J ) )
 846   430          CONTINUE
 847             END IF
 848   440    CONTINUE
 849 *
 850          DO 460 K = 1, KB - 1
 851             J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
 852 *
 853 *           finish applying rotations in 1st set from the right
 854 *
 855             DO 450 L = KB - K, 1, -1
 856                NRT = ( N-J2+L ) / KA1
 857                IF( NRT.GT.0 )
 858      $            CALL CLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
 859      $                         AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ),
 860      $                         WORK( J2-M ), KA1 )
 861   450       CONTINUE
 862   460    CONTINUE
 863 *
 864          IF( KB.GT.1 ) THEN
 865             DO 470 J = N - 1, J2 + KA, -1
 866                RWORK( J-M ) = RWORK( J-KA-M )
 867                WORK( J-M ) = WORK( J-KA-M )
 868   470       CONTINUE
 869          END IF
 870 *
 871       END IF
 872 *
 873       GO TO 10
 874 *
 875   480 CONTINUE
 876 *
 877 *     **************************** Phase 2 *****************************
 878 *
 879 *     The logical structure of this phase is:
 880 *
 881 *     UPDATE = .TRUE.
 882 *     DO I = 1, M
 883 *        use S(i) to update A and create a new bulge
 884 *        apply rotations to push all bulges KA positions upward
 885 *     END DO
 886 *     UPDATE = .FALSE.
 887 *     DO I = M - KA - 1, 2, -1
 888 *        apply rotations to push all bulges KA positions upward
 889 *     END DO
 890 *
 891 *     To avoid duplicating code, the two loops are merged.
 892 *
 893       UPDATE = .TRUE.
 894       I = 0
 895   490 CONTINUE
 896       IF( UPDATE ) THEN
 897          I = I + 1
 898          KBT = MIN( KB, M-I )
 899          I0 = I + 1
 900          I1 = MAX( 1, I-KA )
 901          I2 = I + KBT - KA1
 902          IF( I.GT.M ) THEN
 903             UPDATE = .FALSE.
 904             I = I - 1
 905             I0 = M + 1
 906             IF( KA.EQ.0 )
 907      $         RETURN
 908             GO TO 490
 909          END IF
 910       ELSE
 911          I = I - KA
 912          IF( I.LT.2 )
 913      $      RETURN
 914       END IF
 915 *
 916       IF( I.LT.M-KBT ) THEN
 917          NX = M
 918       ELSE
 919          NX = N
 920       END IF
 921 *
 922       IF( UPPER ) THEN
 923 *
 924 *        Transform A, working with the upper triangle
 925 *
 926          IF( UPDATE ) THEN
 927 *
 928 *           Form  inv(S(i))**H * A * inv(S(i))
 929 *
 930             BII = REAL( BB( KB1, I ) )
 931             AB( KA1, I ) = ( REAL( AB( KA1, I ) ) / BII ) / BII
 932             DO 500 J = I1, I - 1
 933                AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII
 934   500       CONTINUE
 935             DO 510 J = I + 1, MIN( N, I+KA )
 936                AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII
 937   510       CONTINUE
 938             DO 540 K = I + 1, I + KBT
 939                DO 520 J = K, I + KBT
 940                   AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
 941      $                               BB( I-J+KB1, J )*
 942      $                               CONJG( AB( I-K+KA1, K ) ) -
 943      $                               CONJG( BB( I-K+KB1, K ) )*
 944      $                               AB( I-J+KA1, J ) +
 945      $                               REAL( AB( KA1, I ) )*
 946      $                               BB( I-J+KB1, J )*
 947      $                               CONJG( BB( I-K+KB1, K ) )
 948   520          CONTINUE
 949                DO 530 J = I + KBT + 1, MIN( N, I+KA )
 950                   AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
 951      $                               CONJG( BB( I-K+KB1, K ) )*
 952      $                               AB( I-J+KA1, J )
 953   530          CONTINUE
 954   540       CONTINUE
 955             DO 560 J = I1, I
 956                DO 550 K = I + 1, MIN( J+KA, I+KBT )
 957                   AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
 958      $                               BB( I-K+KB1, K )*AB( J-I+KA1, I )
 959   550          CONTINUE
 960   560       CONTINUE
 961 *
 962             IF( WANTX ) THEN
 963 *
 964 *              post-multiply X by inv(S(i))
 965 *
 966                CALL CSSCAL( NX, ONE / BII, X( 1, I ), 1 )
 967                IF( KBT.GT.0 )
 968      $            CALL CGERU( NX, KBT, -CONE, X( 1, I ), 1,
 969      $                        BB( KB, I+1 ), LDBB-1, X( 1, I+1 ), LDX )
 970             END IF
 971 *
 972 *           store a(i1,i) in RA1 for use in next loop over K
 973 *
 974             RA1 = AB( I1-I+KA1, I )
 975          END IF
 976 *
 977 *        Generate and apply vectors of rotations to chase all the
 978 *        existing bulges KA positions up toward the top of the band
 979 *
 980          DO 610 K = 1, KB - 1
 981             IF( UPDATE ) THEN
 982 *
 983 *              Determine the rotations which would annihilate the bulge
 984 *              which has in theory just been created
 985 *
 986                IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
 987 *
 988 *                 generate rotation to annihilate a(i+k-ka-1,i)
 989 *
 990                   CALL CLARTG( AB( K+1, I ), RA1, RWORK( I+K-KA ),
 991      $                         WORK( I+K-KA ), RA )
 992 *
 993 *                 create nonzero element a(i+k-ka-1,i+k) outside the
 994 *                 band and store it in WORK(m-kb+i+k)
 995 *
 996                   T = -BB( KB1-K, I+K )*RA1
 997                   WORK( M-KB+I+K ) = RWORK( I+K-KA )*T -
 998      $                               CONJG( WORK( I+K-KA ) )*
 999      $                               AB( 1, I+K )
1000                   AB( 1, I+K ) = WORK( I+K-KA )*T +
1001      $                           RWORK( I+K-KA )*AB( 1, I+K )
1002                   RA1 = RA
1003                END IF
1004             END IF
1005             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
1006             NR = ( J2+KA-1 ) / KA1
1007             J1 = J2 - ( NR-1 )*KA1
1008             IF( UPDATE ) THEN
1009                J2T = MIN( J2, I-2*KA+K-1 )
1010             ELSE
1011                J2T = J2
1012             END IF
1013             NRT = ( J2T+KA-1 ) / KA1
1014             DO 570 J = J1, J2T, KA1
1015 *
1016 *              create nonzero element a(j-1,j+ka) outside the band
1017 *              and store it in WORK(j)
1018 *
1019                WORK( J ) = WORK( J )*AB( 1, J+KA-1 )
1020                AB( 1, J+KA-1 ) = RWORK( J )*AB( 1, J+KA-1 )
1021   570       CONTINUE
1022 *
1023 *           generate rotations in 1st set to annihilate elements which
1024 *           have been created outside the band
1025 *
1026             IF( NRT.GT.0 )
1027      $         CALL CLARGV( NRT, AB( 1, J1+KA ), INCA, WORK( J1 ), KA1,
1028      $                      RWORK( J1 ), KA1 )
1029             IF( NR.GT.0 ) THEN
1030 *
1031 *              apply rotations in 1st set from the left
1032 *
1033                DO 580 L = 1, KA - 1
1034                   CALL CLARTV( NR, AB( KA1-L, J1+L ), INCA,
1035      $                         AB( KA-L, J1+L ), INCA, RWORK( J1 ),
1036      $                         WORK( J1 ), KA1 )
1037   580          CONTINUE
1038 *
1039 *              apply rotations in 1st set from both sides to diagonal
1040 *              blocks
1041 *
1042                CALL CLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
1043      $                      AB( KA, J1 ), INCA, RWORK( J1 ), WORK( J1 ),
1044      $                      KA1 )
1045 *
1046                CALL CLACGV( NR, WORK( J1 ), KA1 )
1047             END IF
1048 *
1049 *           start applying rotations in 1st set from the right
1050 *
1051             DO 590 L = KA - 1, KB - K + 1, -1
1052                NRT = ( J2+L-1 ) / KA1
1053                J1T = J2 - ( NRT-1 )*KA1
1054                IF( NRT.GT.0 )
1055      $            CALL CLARTV( NRT, AB( L, J1T ), INCA,
1056      $                         AB( L+1, J1T-1 ), INCA, RWORK( J1T ),
1057      $                         WORK( J1T ), KA1 )
1058   590       CONTINUE
1059 *
1060             IF( WANTX ) THEN
1061 *
1062 *              post-multiply X by product of rotations in 1st set
1063 *
1064                DO 600 J = J1, J2, KA1
1065                   CALL CROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1066      $                       RWORK( J ), WORK( J ) )
1067   600          CONTINUE
1068             END IF
1069   610    CONTINUE
1070 *
1071          IF( UPDATE ) THEN
1072             IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
1073 *
1074 *              create nonzero element a(i+kbt-ka-1,i+kbt) outside the
1075 *              band and store it in WORK(m-kb+i+kbt)
1076 *
1077                WORK( M-KB+I+KBT ) = -BB( KB1-KBT, I+KBT )*RA1
1078             END IF
1079          END IF
1080 *
1081          DO 650 K = KB, 1, -1
1082             IF( UPDATE ) THEN
1083                J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1
1084             ELSE
1085                J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1086             END IF
1087 *
1088 *           finish applying rotations in 2nd set from the right
1089 *
1090             DO 620 L = KB - K, 1, -1
1091                NRT = ( J2+KA+L-1 ) / KA1
1092                J1T = J2 - ( NRT-1 )*KA1
1093                IF( NRT.GT.0 )
1094      $            CALL CLARTV( NRT, AB( L, J1T+KA ), INCA,
1095      $                         AB( L+1, J1T+KA-1 ), INCA,
1096      $                         RWORK( M-KB+J1T+KA ),
1097      $                         WORK( M-KB+J1T+KA ), KA1 )
1098   620       CONTINUE
1099             NR = ( J2+KA-1 ) / KA1
1100             J1 = J2 - ( NR-1 )*KA1
1101             DO 630 J = J1, J2, KA1
1102                WORK( M-KB+J ) = WORK( M-KB+J+KA )
1103                RWORK( M-KB+J ) = RWORK( M-KB+J+KA )
1104   630       CONTINUE
1105             DO 640 J = J1, J2, KA1
1106 *
1107 *              create nonzero element a(j-1,j+ka) outside the band
1108 *              and store it in WORK(m-kb+j)
1109 *
1110                WORK( M-KB+J ) = WORK( M-KB+J )*AB( 1, J+KA-1 )
1111                AB( 1, J+KA-1 ) = RWORK( M-KB+J )*AB( 1, J+KA-1 )
1112   640       CONTINUE
1113             IF( UPDATE ) THEN
1114                IF( I+K.GT.KA1 .AND. K.LE.KBT )
1115      $            WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
1116             END IF
1117   650    CONTINUE
1118 *
1119          DO 690 K = KB, 1, -1
1120             J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1121             NR = ( J2+KA-1 ) / KA1
1122             J1 = J2 - ( NR-1 )*KA1
1123             IF( NR.GT.0 ) THEN
1124 *
1125 *              generate rotations in 2nd set to annihilate elements
1126 *              which have been created outside the band
1127 *
1128                CALL CLARGV( NR, AB( 1, J1+KA ), INCA, WORK( M-KB+J1 ),
1129      $                      KA1, RWORK( M-KB+J1 ), KA1 )
1130 *
1131 *              apply rotations in 2nd set from the left
1132 *
1133                DO 660 L = 1, KA - 1
1134                   CALL CLARTV( NR, AB( KA1-L, J1+L ), INCA,
1135      $                         AB( KA-L, J1+L ), INCA, RWORK( M-KB+J1 ),
1136      $                         WORK( M-KB+J1 ), KA1 )
1137   660          CONTINUE
1138 *
1139 *              apply rotations in 2nd set from both sides to diagonal
1140 *              blocks
1141 *
1142                CALL CLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
1143      $                      AB( KA, J1 ), INCA, RWORK( M-KB+J1 ),
1144      $                      WORK( M-KB+J1 ), KA1 )
1145 *
1146                CALL CLACGV( NR, WORK( M-KB+J1 ), KA1 )
1147             END IF
1148 *
1149 *           start applying rotations in 2nd set from the right
1150 *
1151             DO 670 L = KA - 1, KB - K + 1, -1
1152                NRT = ( J2+L-1 ) / KA1
1153                J1T = J2 - ( NRT-1 )*KA1
1154                IF( NRT.GT.0 )
1155      $            CALL CLARTV( NRT, AB( L, J1T ), INCA,
1156      $                         AB( L+1, J1T-1 ), INCA,
1157      $                         RWORK( M-KB+J1T ), WORK( M-KB+J1T ),
1158      $                         KA1 )
1159   670       CONTINUE
1160 *
1161             IF( WANTX ) THEN
1162 *
1163 *              post-multiply X by product of rotations in 2nd set
1164 *
1165                DO 680 J = J1, J2, KA1
1166                   CALL CROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1167      $                       RWORK( M-KB+J ), WORK( M-KB+J ) )
1168   680          CONTINUE
1169             END IF
1170   690    CONTINUE
1171 *
1172          DO 710 K = 1, KB - 1
1173             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
1174 *
1175 *           finish applying rotations in 1st set from the right
1176 *
1177             DO 700 L = KB - K, 1, -1
1178                NRT = ( J2+L-1 ) / KA1
1179                J1T = J2 - ( NRT-1 )*KA1
1180                IF( NRT.GT.0 )
1181      $            CALL CLARTV( NRT, AB( L, J1T ), INCA,
1182      $                         AB( L+1, J1T-1 ), INCA, RWORK( J1T ),
1183      $                         WORK( J1T ), KA1 )
1184   700       CONTINUE
1185   710    CONTINUE
1186 *
1187          IF( KB.GT.1 ) THEN
1188             DO 720 J = 2, I2 - KA
1189                RWORK( J ) = RWORK( J+KA )
1190                WORK( J ) = WORK( J+KA )
1191   720       CONTINUE
1192          END IF
1193 *
1194       ELSE
1195 *
1196 *        Transform A, working with the lower triangle
1197 *
1198          IF( UPDATE ) THEN
1199 *
1200 *           Form  inv(S(i))**H * A * inv(S(i))
1201 *
1202             BII = REAL( BB( 1, I ) )
1203             AB( 1, I ) = ( REAL( AB( 1, I ) ) / BII ) / BII
1204             DO 730 J = I1, I - 1
1205                AB( I-J+1, J ) = AB( I-J+1, J ) / BII
1206   730       CONTINUE
1207             DO 740 J = I + 1, MIN( N, I+KA )
1208                AB( J-I+1, I ) = AB( J-I+1, I ) / BII
1209   740       CONTINUE
1210             DO 770 K = I + 1, I + KBT
1211                DO 750 J = K, I + KBT
1212                   AB( J-K+1, K ) = AB( J-K+1, K ) -
1213      $                             BB( J-I+1, I )*CONJG( AB( K-I+1,
1214      $                             I ) ) - CONJG( BB( K-I+1, I ) )*
1215      $                             AB( J-I+1, I ) + REAL( AB( 1, I ) )*
1216      $                             BB( J-I+1, I )*CONJG( BB( K-I+1,
1217      $                             I ) )
1218   750          CONTINUE
1219                DO 760 J = I + KBT + 1, MIN( N, I+KA )
1220                   AB( J-K+1, K ) = AB( J-K+1, K ) -
1221      $                             CONJG( BB( K-I+1, I ) )*
1222      $                             AB( J-I+1, I )
1223   760          CONTINUE
1224   770       CONTINUE
1225             DO 790 J = I1, I
1226                DO 780 K = I + 1, MIN( J+KA, I+KBT )
1227                   AB( K-J+1, J ) = AB( K-J+1, J ) -
1228      $                             BB( K-I+1, I )*AB( I-J+1, J )
1229   780          CONTINUE
1230   790       CONTINUE
1231 *
1232             IF( WANTX ) THEN
1233 *
1234 *              post-multiply X by inv(S(i))
1235 *
1236                CALL CSSCAL( NX, ONE / BII, X( 1, I ), 1 )
1237                IF( KBT.GT.0 )
1238      $            CALL CGERC( NX, KBT, -CONE, X( 1, I ), 1, BB( 2, I ),
1239      $                        1, X( 1, I+1 ), LDX )
1240             END IF
1241 *
1242 *           store a(i,i1) in RA1 for use in next loop over K
1243 *
1244             RA1 = AB( I-I1+1, I1 )
1245          END IF
1246 *
1247 *        Generate and apply vectors of rotations to chase all the
1248 *        existing bulges KA positions up toward the top of the band
1249 *
1250          DO 840 K = 1, KB - 1
1251             IF( UPDATE ) THEN
1252 *
1253 *              Determine the rotations which would annihilate the bulge
1254 *              which has in theory just been created
1255 *
1256                IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
1257 *
1258 *                 generate rotation to annihilate a(i,i+k-ka-1)
1259 *
1260                   CALL CLARTG( AB( KA1-K, I+K-KA ), RA1,
1261      $                         RWORK( I+K-KA ), WORK( I+K-KA ), RA )
1262 *
1263 *                 create nonzero element a(i+k,i+k-ka-1) outside the
1264 *                 band and store it in WORK(m-kb+i+k)
1265 *
1266                   T = -BB( K+1, I )*RA1
1267                   WORK( M-KB+I+K ) = RWORK( I+K-KA )*T -
1268      $                               CONJG( WORK( I+K-KA ) )*
1269      $                               AB( KA1, I+K-KA )
1270                   AB( KA1, I+K-KA ) = WORK( I+K-KA )*T +
1271      $                                RWORK( I+K-KA )*AB( KA1, I+K-KA )
1272                   RA1 = RA
1273                END IF
1274             END IF
1275             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
1276             NR = ( J2+KA-1 ) / KA1
1277             J1 = J2 - ( NR-1 )*KA1
1278             IF( UPDATE ) THEN
1279                J2T = MIN( J2, I-2*KA+K-1 )
1280             ELSE
1281                J2T = J2
1282             END IF
1283             NRT = ( J2T+KA-1 ) / KA1
1284             DO 800 J = J1, J2T, KA1
1285 *
1286 *              create nonzero element a(j+ka,j-1) outside the band
1287 *              and store it in WORK(j)
1288 *
1289                WORK( J ) = WORK( J )*AB( KA1, J-1 )
1290                AB( KA1, J-1 ) = RWORK( J )*AB( KA1, J-1 )
1291   800       CONTINUE
1292 *
1293 *           generate rotations in 1st set to annihilate elements which
1294 *           have been created outside the band
1295 *
1296             IF( NRT.GT.0 )
1297      $         CALL CLARGV( NRT, AB( KA1, J1 ), INCA, WORK( J1 ), KA1,
1298      $                      RWORK( J1 ), KA1 )
1299             IF( NR.GT.0 ) THEN
1300 *
1301 *              apply rotations in 1st set from the right
1302 *
1303                DO 810 L = 1, KA - 1
1304                   CALL CLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
1305      $                         INCA, RWORK( J1 ), WORK( J1 ), KA1 )
1306   810          CONTINUE
1307 *
1308 *              apply rotations in 1st set from both sides to diagonal
1309 *              blocks
1310 *
1311                CALL CLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
1312      $                      AB( 2, J1-1 ), INCA, RWORK( J1 ),
1313      $                      WORK( J1 ), KA1 )
1314 *
1315                CALL CLACGV( NR, WORK( J1 ), KA1 )
1316             END IF
1317 *
1318 *           start applying rotations in 1st set from the left
1319 *
1320             DO 820 L = KA - 1, KB - K + 1, -1
1321                NRT = ( J2+L-1 ) / KA1
1322                J1T = J2 - ( NRT-1 )*KA1
1323                IF( NRT.GT.0 )
1324      $            CALL CLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
1325      $                         AB( KA1-L, J1T-KA1+L ), INCA,
1326      $                         RWORK( J1T ), WORK( J1T ), KA1 )
1327   820       CONTINUE
1328 *
1329             IF( WANTX ) THEN
1330 *
1331 *              post-multiply X by product of rotations in 1st set
1332 *
1333                DO 830 J = J1, J2, KA1
1334                   CALL CROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1335      $                       RWORK( J ), CONJG( WORK( J ) ) )
1336   830          CONTINUE
1337             END IF
1338   840    CONTINUE
1339 *
1340          IF( UPDATE ) THEN
1341             IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
1342 *
1343 *              create nonzero element a(i+kbt,i+kbt-ka-1) outside the
1344 *              band and store it in WORK(m-kb+i+kbt)
1345 *
1346                WORK( M-KB+I+KBT ) = -BB( KBT+1, I )*RA1
1347             END IF
1348          END IF
1349 *
1350          DO 880 K = KB, 1, -1
1351             IF( UPDATE ) THEN
1352                J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1
1353             ELSE
1354                J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1355             END IF
1356 *
1357 *           finish applying rotations in 2nd set from the left
1358 *
1359             DO 850 L = KB - K, 1, -1
1360                NRT = ( J2+KA+L-1 ) / KA1
1361                J1T = J2 - ( NRT-1 )*KA1
1362                IF( NRT.GT.0 )
1363      $            CALL CLARTV( NRT, AB( KA1-L+1, J1T+L-1 ), INCA,
1364      $                         AB( KA1-L, J1T+L-1 ), INCA,
1365      $                         RWORK( M-KB+J1T+KA ),
1366      $                         WORK( M-KB+J1T+KA ), KA1 )
1367   850       CONTINUE
1368             NR = ( J2+KA-1 ) / KA1
1369             J1 = J2 - ( NR-1 )*KA1
1370             DO 860 J = J1, J2, KA1
1371                WORK( M-KB+J ) = WORK( M-KB+J+KA )
1372                RWORK( M-KB+J ) = RWORK( M-KB+J+KA )
1373   860       CONTINUE
1374             DO 870 J = J1, J2, KA1
1375 *
1376 *              create nonzero element a(j+ka,j-1) outside the band
1377 *              and store it in WORK(m-kb+j)
1378 *
1379                WORK( M-KB+J ) = WORK( M-KB+J )*AB( KA1, J-1 )
1380                AB( KA1, J-1 ) = RWORK( M-KB+J )*AB( KA1, J-1 )
1381   870       CONTINUE
1382             IF( UPDATE ) THEN
1383                IF( I+K.GT.KA1 .AND. K.LE.KBT )
1384      $            WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
1385             END IF
1386   880    CONTINUE
1387 *
1388          DO 920 K = KB, 1, -1
1389             J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1390             NR = ( J2+KA-1 ) / KA1
1391             J1 = J2 - ( NR-1 )*KA1
1392             IF( NR.GT.0 ) THEN
1393 *
1394 *              generate rotations in 2nd set to annihilate elements
1395 *              which have been created outside the band
1396 *
1397                CALL CLARGV( NR, AB( KA1, J1 ), INCA, WORK( M-KB+J1 ),
1398      $                      KA1, RWORK( M-KB+J1 ), KA1 )
1399 *
1400 *              apply rotations in 2nd set from the right
1401 *
1402                DO 890 L = 1, KA - 1
1403                   CALL CLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
1404      $                         INCA, RWORK( M-KB+J1 ), WORK( M-KB+J1 ),
1405      $                         KA1 )
1406   890          CONTINUE
1407 *
1408 *              apply rotations in 2nd set from both sides to diagonal
1409 *              blocks
1410 *
1411                CALL CLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
1412      $                      AB( 2, J1-1 ), INCA, RWORK( M-KB+J1 ),
1413      $                      WORK( M-KB+J1 ), KA1 )
1414 *
1415                CALL CLACGV( NR, WORK( M-KB+J1 ), KA1 )
1416             END IF
1417 *
1418 *           start applying rotations in 2nd set from the left
1419 *
1420             DO 900 L = KA - 1, KB - K + 1, -1
1421                NRT = ( J2+L-1 ) / KA1
1422                J1T = J2 - ( NRT-1 )*KA1
1423                IF( NRT.GT.0 )
1424      $            CALL CLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
1425      $                         AB( KA1-L, J1T-KA1+L ), INCA,
1426      $                         RWORK( M-KB+J1T ), WORK( M-KB+J1T ),
1427      $                         KA1 )
1428   900       CONTINUE
1429 *
1430             IF( WANTX ) THEN
1431 *
1432 *              post-multiply X by product of rotations in 2nd set
1433 *
1434                DO 910 J = J1, J2, KA1
1435                   CALL CROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1436      $                       RWORK( M-KB+J ), CONJG( WORK( M-KB+J ) ) )
1437   910          CONTINUE
1438             END IF
1439   920    CONTINUE
1440 *
1441          DO 940 K = 1, KB - 1
1442             J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
1443 *
1444 *           finish applying rotations in 1st set from the left
1445 *
1446             DO 930 L = KB - K, 1, -1
1447                NRT = ( J2+L-1 ) / KA1
1448                J1T = J2 - ( NRT-1 )*KA1
1449                IF( NRT.GT.0 )
1450      $            CALL CLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
1451      $                         AB( KA1-L, J1T-KA1+L ), INCA,
1452      $                         RWORK( J1T ), WORK( J1T ), KA1 )
1453   930       CONTINUE
1454   940    CONTINUE
1455 *
1456          IF( KB.GT.1 ) THEN
1457             DO 950 J = 2, I2 - KA
1458                RWORK( J ) = RWORK( J+KA )
1459                WORK( J ) = WORK( J+KA )
1460   950       CONTINUE
1461          END IF
1462 *
1463       END IF
1464 *
1465       GO TO 490
1466 *
1467 *     End of CHBGST
1468 *
1469       END