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