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 1470