1 SUBROUTINE ZHBGST( 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 DOUBLE PRECISION RWORK( * ) 15 COMPLEX*16 AB( LDAB, * ), BB( LDBB, * ), WORK( * ), 16 $ X( LDX, * ) 17* .. 18* 19* Purpose 20* ======= 21* 22* ZHBGST 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 ZPBSTF, 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*16 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*16 array, dimension (LDBB,N) 68* The banded factor S from the split Cholesky factorization of 69* B, as returned by ZPBSTF, 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*16 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*16 array, dimension (N) 84* 85* RWORK (workspace) DOUBLE PRECISION 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*16 CZERO, CONE 95 DOUBLE PRECISION ONE 96 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 97 $ CONE = ( 1.0D+0, 0.0D+0 ), ONE = 1.0D+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 DOUBLE PRECISION BII 104 COMPLEX*16 RA, RA1, T 105* .. 106* .. External Functions .. 107 LOGICAL LSAME 108 EXTERNAL LSAME 109* .. 110* .. External Subroutines .. 111 EXTERNAL XERBLA, ZDSCAL, ZGERC, ZGERU, ZLACGV, ZLAR2V, 112 $ ZLARGV, ZLARTG, ZLARTV, ZLASET, ZROT 113* .. 114* .. Intrinsic Functions .. 115 INTRINSIC DBLE, DCONJG, MAX, MIN 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( 'ZHBGST', -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 ZLASET( '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 ZPBSTF. 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 = DBLE( BB( KB1, I ) ) 257 AB( KA1, I ) = ( DBLE( 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 $ DCONJG( AB( K-I+KA1, I ) ) - 269 $ DCONJG( BB( K-I+KB1, I ) )* 270 $ AB( J-I+KA1, I ) + 271 $ DBLE( AB( KA1, I ) )* 272 $ BB( J-I+KB1, I )* 273 $ DCONJG( 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 $ DCONJG( 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 ZDSCAL( N-M, ONE / BII, X( M+1, I ), 1 ) 293 IF( KBT.GT.0 ) 294 $ CALL ZGERC( 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 ZLARTG( 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 $ DCONJG( 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 ZLARGV( 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 ZLARTV( 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 ZLAR2V( 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 ZLACGV( 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 ZLARTV( 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 ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1, 393 $ RWORK( J-M ), DCONJG( 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 ZLARTV( 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 ZLARGV( 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 ZLARTV( 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 ZLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ), 468 $ AB( KA, J2+1 ), INCA, RWORK( J2 ), 469 $ WORK( J2 ), KA1 ) 470* 471 CALL ZLACGV( 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 ZLARTV( 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 ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1, 490 $ RWORK( J ), DCONJG( 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 ZLARTV( 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 = DBLE( BB( 1, I ) ) 525 AB( 1, I ) = ( DBLE( 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 )*DCONJG( AB( I-K+1, 536 $ K ) ) - DCONJG( BB( I-K+1, K ) )* 537 $ AB( I-J+1, J ) + DBLE( AB( 1, I ) )* 538 $ BB( I-J+1, J )*DCONJG( 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 $ DCONJG( 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 ZDSCAL( N-M, ONE / BII, X( M+1, I ), 1 ) 559 IF( KBT.GT.0 ) 560 $ CALL ZGERU( 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 ZLARTG( 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 $ DCONJG( WORK( I-K+KA-M ) )* 593 $ AB( KA1, I-K ) 594 AB( KA1, I-K ) = WORK( I-K+KA-M )*T + 595 $ RWORK( I-K+KA-M )*AB( KA1, I-K ) 596 RA1 = RA 597 END IF 598 END IF 599 J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1 600 NR = ( N-J2+KA ) / KA1 601 J1 = J2 + ( NR-1 )*KA1 602 IF( UPDATE ) THEN 603 J2T = MAX( J2, I+2*KA-K+1 ) 604 ELSE 605 J2T = J2 606 END IF 607 NRT = ( N-J2T+KA ) / KA1 608 DO 320 J = J2T, J1, KA1 609* 610* create nonzero element a(j+1,j-ka) outside the band 611* and store it in WORK(j-m) 612* 613 WORK( J-M ) = WORK( J-M )*AB( KA1, J-KA+1 ) 614 AB( KA1, J-KA+1 ) = RWORK( J-M )*AB( KA1, J-KA+1 ) 615 320 CONTINUE 616* 617* generate rotations in 1st set to annihilate elements which 618* have been created outside the band 619* 620 IF( NRT.GT.0 ) 621 $ CALL ZLARGV( NRT, AB( KA1, J2T-KA ), INCA, WORK( J2T-M ), 622 $ KA1, RWORK( J2T-M ), KA1 ) 623 IF( NR.GT.0 ) THEN 624* 625* apply rotations in 1st set from the left 626* 627 DO 330 L = 1, KA - 1 628 CALL ZLARTV( NR, AB( L+1, J2-L ), INCA, 629 $ AB( L+2, J2-L ), INCA, RWORK( J2-M ), 630 $ WORK( J2-M ), KA1 ) 631 330 CONTINUE 632* 633* apply rotations in 1st set from both sides to diagonal 634* blocks 635* 636 CALL ZLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ), 637 $ INCA, RWORK( J2-M ), WORK( J2-M ), KA1 ) 638* 639 CALL ZLACGV( NR, WORK( J2-M ), KA1 ) 640 END IF 641* 642* start applying rotations in 1st set from the right 643* 644 DO 340 L = KA - 1, KB - K + 1, -1 645 NRT = ( N-J2+L ) / KA1 646 IF( NRT.GT.0 ) 647 $ CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA, 648 $ AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ), 649 $ WORK( J2-M ), KA1 ) 650 340 CONTINUE 651* 652 IF( WANTX ) THEN 653* 654* post-multiply X by product of rotations in 1st set 655* 656 DO 350 J = J2, J1, KA1 657 CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1, 658 $ RWORK( J-M ), WORK( J-M ) ) 659 350 CONTINUE 660 END IF 661 360 CONTINUE 662* 663 IF( UPDATE ) THEN 664 IF( I2.LE.N .AND. KBT.GT.0 ) THEN 665* 666* create nonzero element a(i-kbt+ka+1,i-kbt) outside the 667* band and store it in WORK(i-kbt) 668* 669 WORK( I-KBT ) = -BB( KBT+1, I-KBT )*RA1 670 END IF 671 END IF 672* 673 DO 400 K = KB, 1, -1 674 IF( UPDATE ) THEN 675 J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1 676 ELSE 677 J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1 678 END IF 679* 680* finish applying rotations in 2nd set from the right 681* 682 DO 370 L = KB - K, 1, -1 683 NRT = ( N-J2+KA+L ) / KA1 684 IF( NRT.GT.0 ) 685 $ CALL ZLARTV( NRT, AB( KA1-L+1, J2-KA ), INCA, 686 $ AB( KA1-L, J2-KA+1 ), INCA, 687 $ RWORK( J2-KA ), WORK( J2-KA ), KA1 ) 688 370 CONTINUE 689 NR = ( N-J2+KA ) / KA1 690 J1 = J2 + ( NR-1 )*KA1 691 DO 380 J = J1, J2, -KA1 692 WORK( J ) = WORK( J-KA ) 693 RWORK( J ) = RWORK( J-KA ) 694 380 CONTINUE 695 DO 390 J = J2, J1, KA1 696* 697* create nonzero element a(j+1,j-ka) outside the band 698* and store it in WORK(j) 699* 700 WORK( J ) = WORK( J )*AB( KA1, J-KA+1 ) 701 AB( KA1, J-KA+1 ) = RWORK( J )*AB( KA1, J-KA+1 ) 702 390 CONTINUE 703 IF( UPDATE ) THEN 704 IF( I-K.LT.N-KA .AND. K.LE.KBT ) 705 $ WORK( I-K+KA ) = WORK( I-K ) 706 END IF 707 400 CONTINUE 708* 709 DO 440 K = KB, 1, -1 710 J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1 711 NR = ( N-J2+KA ) / KA1 712 J1 = J2 + ( NR-1 )*KA1 713 IF( NR.GT.0 ) THEN 714* 715* generate rotations in 2nd set to annihilate elements 716* which have been created outside the band 717* 718 CALL ZLARGV( NR, AB( KA1, J2-KA ), INCA, WORK( J2 ), KA1, 719 $ RWORK( J2 ), KA1 ) 720* 721* apply rotations in 2nd set from the left 722* 723 DO 410 L = 1, KA - 1 724 CALL ZLARTV( NR, AB( L+1, J2-L ), INCA, 725 $ AB( L+2, J2-L ), INCA, RWORK( J2 ), 726 $ WORK( J2 ), KA1 ) 727 410 CONTINUE 728* 729* apply rotations in 2nd set from both sides to diagonal 730* blocks 731* 732 CALL ZLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ), 733 $ INCA, RWORK( J2 ), WORK( J2 ), KA1 ) 734* 735 CALL ZLACGV( NR, WORK( J2 ), KA1 ) 736 END IF 737* 738* start applying rotations in 2nd set from the right 739* 740 DO 420 L = KA - 1, KB - K + 1, -1 741 NRT = ( N-J2+L ) / KA1 742 IF( NRT.GT.0 ) 743 $ CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA, 744 $ AB( KA1-L, J2+1 ), INCA, RWORK( J2 ), 745 $ WORK( J2 ), KA1 ) 746 420 CONTINUE 747* 748 IF( WANTX ) THEN 749* 750* post-multiply X by product of rotations in 2nd set 751* 752 DO 430 J = J2, J1, KA1 753 CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1, 754 $ RWORK( J ), WORK( J ) ) 755 430 CONTINUE 756 END IF 757 440 CONTINUE 758* 759 DO 460 K = 1, KB - 1 760 J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1 761* 762* finish applying rotations in 1st set from the right 763* 764 DO 450 L = KB - K, 1, -1 765 NRT = ( N-J2+L ) / KA1 766 IF( NRT.GT.0 ) 767 $ CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA, 768 $ AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ), 769 $ WORK( J2-M ), KA1 ) 770 450 CONTINUE 771 460 CONTINUE 772* 773 IF( KB.GT.1 ) THEN 774 DO 470 J = N - 1, I2 + KA, -1 775 RWORK( J-M ) = RWORK( J-KA-M ) 776 WORK( J-M ) = WORK( J-KA-M ) 777 470 CONTINUE 778 END IF 779* 780 END IF 781* 782 GO TO 10 783* 784 480 CONTINUE 785* 786* **************************** Phase 2 ***************************** 787* 788* The logical structure of this phase is: 789* 790* UPDATE = .TRUE. 791* DO I = 1, M 792* use S(i) to update A and create a new bulge 793* apply rotations to push all bulges KA positions upward 794* END DO 795* UPDATE = .FALSE. 796* DO I = M - KA - 1, 2, -1 797* apply rotations to push all bulges KA positions upward 798* END DO 799* 800* To avoid duplicating code, the two loops are merged. 801* 802 UPDATE = .TRUE. 803 I = 0 804 490 CONTINUE 805 IF( UPDATE ) THEN 806 I = I + 1 807 KBT = MIN( KB, M-I ) 808 I0 = I + 1 809 I1 = MAX( 1, I-KA ) 810 I2 = I + KBT - KA1 811 IF( I.GT.M ) THEN 812 UPDATE = .FALSE. 813 I = I - 1 814 I0 = M + 1 815 IF( KA.EQ.0 ) 816 $ RETURN 817 GO TO 490 818 END IF 819 ELSE 820 I = I - KA 821 IF( I.LT.2 ) 822 $ RETURN 823 END IF 824* 825 IF( I.LT.M-KBT ) THEN 826 NX = M 827 ELSE 828 NX = N 829 END IF 830* 831 IF( UPPER ) THEN 832* 833* Transform A, working with the upper triangle 834* 835 IF( UPDATE ) THEN 836* 837* Form inv(S(i))**H * A * inv(S(i)) 838* 839 BII = DBLE( BB( KB1, I ) ) 840 AB( KA1, I ) = ( DBLE( AB( KA1, I ) ) / BII ) / BII 841 DO 500 J = I1, I - 1 842 AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII 843 500 CONTINUE 844 DO 510 J = I + 1, MIN( N, I+KA ) 845 AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII 846 510 CONTINUE 847 DO 540 K = I + 1, I + KBT 848 DO 520 J = K, I + KBT 849 AB( K-J+KA1, J ) = AB( K-J+KA1, J ) - 850 $ BB( I-J+KB1, J )* 851 $ DCONJG( AB( I-K+KA1, K ) ) - 852 $ DCONJG( BB( I-K+KB1, K ) )* 853 $ AB( I-J+KA1, J ) + 854 $ DBLE( AB( KA1, I ) )* 855 $ BB( I-J+KB1, J )* 856 $ DCONJG( BB( I-K+KB1, K ) ) 857 520 CONTINUE 858 DO 530 J = I + KBT + 1, MIN( N, I+KA ) 859 AB( K-J+KA1, J ) = AB( K-J+KA1, J ) - 860 $ DCONJG( BB( I-K+KB1, K ) )* 861 $ AB( I-J+KA1, J ) 862 530 CONTINUE 863 540 CONTINUE 864 DO 560 J = I1, I 865 DO 550 K = I + 1, MIN( J+KA, I+KBT ) 866 AB( J-K+KA1, K ) = AB( J-K+KA1, K ) - 867 $ BB( I-K+KB1, K )*AB( J-I+KA1, I ) 868 550 CONTINUE 869 560 CONTINUE 870* 871 IF( WANTX ) THEN 872* 873* post-multiply X by inv(S(i)) 874* 875 CALL ZDSCAL( NX, ONE / BII, X( 1, I ), 1 ) 876 IF( KBT.GT.0 ) 877 $ CALL ZGERU( NX, KBT, -CONE, X( 1, I ), 1, 878 $ BB( KB, I+1 ), LDBB-1, X( 1, I+1 ), LDX ) 879 END IF 880* 881* store a(i1,i) in RA1 for use in next loop over K 882* 883 RA1 = AB( I1-I+KA1, I ) 884 END IF 885* 886* Generate and apply vectors of rotations to chase all the 887* existing bulges KA positions up toward the top of the band 888* 889 DO 610 K = 1, KB - 1 890 IF( UPDATE ) THEN 891* 892* Determine the rotations which would annihilate the bulge 893* which has in theory just been created 894* 895 IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN 896* 897* generate rotation to annihilate a(i+k-ka-1,i) 898* 899 CALL ZLARTG( AB( K+1, I ), RA1, RWORK( I+K-KA ), 900 $ WORK( I+K-KA ), RA ) 901* 902* create nonzero element a(i+k-ka-1,i+k) outside the 903* band and store it in WORK(m-kb+i+k) 904* 905 T = -BB( KB1-K, I+K )*RA1 906 WORK( M-KB+I+K ) = RWORK( I+K-KA )*T - 907 $ DCONJG( WORK( I+K-KA ) )* 908 $ AB( 1, I+K ) 909 AB( 1, I+K ) = WORK( I+K-KA )*T + 910 $ RWORK( I+K-KA )*AB( 1, I+K ) 911 RA1 = RA 912 END IF 913 END IF 914 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1 915 NR = ( J2+KA-1 ) / KA1 916 J1 = J2 - ( NR-1 )*KA1 917 IF( UPDATE ) THEN 918 J2T = MIN( J2, I-2*KA+K-1 ) 919 ELSE 920 J2T = J2 921 END IF 922 NRT = ( J2T+KA-1 ) / KA1 923 DO 570 J = J1, J2T, KA1 924* 925* create nonzero element a(j-1,j+ka) outside the band 926* and store it in WORK(j) 927* 928 WORK( J ) = WORK( J )*AB( 1, J+KA-1 ) 929 AB( 1, J+KA-1 ) = RWORK( J )*AB( 1, J+KA-1 ) 930 570 CONTINUE 931* 932* generate rotations in 1st set to annihilate elements which 933* have been created outside the band 934* 935 IF( NRT.GT.0 ) 936 $ CALL ZLARGV( NRT, AB( 1, J1+KA ), INCA, WORK( J1 ), KA1, 937 $ RWORK( J1 ), KA1 ) 938 IF( NR.GT.0 ) THEN 939* 940* apply rotations in 1st set from the left 941* 942 DO 580 L = 1, KA - 1 943 CALL ZLARTV( NR, AB( KA1-L, J1+L ), INCA, 944 $ AB( KA-L, J1+L ), INCA, RWORK( J1 ), 945 $ WORK( J1 ), KA1 ) 946 580 CONTINUE 947* 948* apply rotations in 1st set from both sides to diagonal 949* blocks 950* 951 CALL ZLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ), 952 $ AB( KA, J1 ), INCA, RWORK( J1 ), WORK( J1 ), 953 $ KA1 ) 954* 955 CALL ZLACGV( NR, WORK( J1 ), KA1 ) 956 END IF 957* 958* start applying rotations in 1st set from the right 959* 960 DO 590 L = KA - 1, KB - K + 1, -1 961 NRT = ( J2+L-1 ) / KA1 962 J1T = J2 - ( NRT-1 )*KA1 963 IF( NRT.GT.0 ) 964 $ CALL ZLARTV( NRT, AB( L, J1T ), INCA, 965 $ AB( L+1, J1T-1 ), INCA, RWORK( J1T ), 966 $ WORK( J1T ), KA1 ) 967 590 CONTINUE 968* 969 IF( WANTX ) THEN 970* 971* post-multiply X by product of rotations in 1st set 972* 973 DO 600 J = J1, J2, KA1 974 CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1, 975 $ RWORK( J ), WORK( J ) ) 976 600 CONTINUE 977 END IF 978 610 CONTINUE 979* 980 IF( UPDATE ) THEN 981 IF( I2.GT.0 .AND. KBT.GT.0 ) THEN 982* 983* create nonzero element a(i+kbt-ka-1,i+kbt) outside the 984* band and store it in WORK(m-kb+i+kbt) 985* 986 WORK( M-KB+I+KBT ) = -BB( KB1-KBT, I+KBT )*RA1 987 END IF 988 END IF 989* 990 DO 650 K = KB, 1, -1 991 IF( UPDATE ) THEN 992 J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1 993 ELSE 994 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1 995 END IF 996* 997* finish applying rotations in 2nd set from the right 998* 999 DO 620 L = KB - K, 1, -1 1000 NRT = ( J2+KA+L-1 ) / KA1 1001 J1T = J2 - ( NRT-1 )*KA1 1002 IF( NRT.GT.0 ) 1003 $ CALL ZLARTV( NRT, AB( L, J1T+KA ), INCA, 1004 $ AB( L+1, J1T+KA-1 ), INCA, 1005 $ RWORK( M-KB+J1T+KA ), 1006 $ WORK( M-KB+J1T+KA ), KA1 ) 1007 620 CONTINUE 1008 NR = ( J2+KA-1 ) / KA1 1009 J1 = J2 - ( NR-1 )*KA1 1010 DO 630 J = J1, J2, KA1 1011 WORK( M-KB+J ) = WORK( M-KB+J+KA ) 1012 RWORK( M-KB+J ) = RWORK( M-KB+J+KA ) 1013 630 CONTINUE 1014 DO 640 J = J1, J2, KA1 1015* 1016* create nonzero element a(j-1,j+ka) outside the band 1017* and store it in WORK(m-kb+j) 1018* 1019 WORK( M-KB+J ) = WORK( M-KB+J )*AB( 1, J+KA-1 ) 1020 AB( 1, J+KA-1 ) = RWORK( M-KB+J )*AB( 1, J+KA-1 ) 1021 640 CONTINUE 1022 IF( UPDATE ) THEN 1023 IF( I+K.GT.KA1 .AND. K.LE.KBT ) 1024 $ WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K ) 1025 END IF 1026 650 CONTINUE 1027* 1028 DO 690 K = KB, 1, -1 1029 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1 1030 NR = ( J2+KA-1 ) / KA1 1031 J1 = J2 - ( NR-1 )*KA1 1032 IF( NR.GT.0 ) THEN 1033* 1034* generate rotations in 2nd set to annihilate elements 1035* which have been created outside the band 1036* 1037 CALL ZLARGV( NR, AB( 1, J1+KA ), INCA, WORK( M-KB+J1 ), 1038 $ KA1, RWORK( M-KB+J1 ), KA1 ) 1039* 1040* apply rotations in 2nd set from the left 1041* 1042 DO 660 L = 1, KA - 1 1043 CALL ZLARTV( NR, AB( KA1-L, J1+L ), INCA, 1044 $ AB( KA-L, J1+L ), INCA, RWORK( M-KB+J1 ), 1045 $ WORK( M-KB+J1 ), KA1 ) 1046 660 CONTINUE 1047* 1048* apply rotations in 2nd set from both sides to diagonal 1049* blocks 1050* 1051 CALL ZLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ), 1052 $ AB( KA, J1 ), INCA, RWORK( M-KB+J1 ), 1053 $ WORK( M-KB+J1 ), KA1 ) 1054* 1055 CALL ZLACGV( NR, WORK( M-KB+J1 ), KA1 ) 1056 END IF 1057* 1058* start applying rotations in 2nd set from the right 1059* 1060 DO 670 L = KA - 1, KB - K + 1, -1 1061 NRT = ( J2+L-1 ) / KA1 1062 J1T = J2 - ( NRT-1 )*KA1 1063 IF( NRT.GT.0 ) 1064 $ CALL ZLARTV( NRT, AB( L, J1T ), INCA, 1065 $ AB( L+1, J1T-1 ), INCA, 1066 $ RWORK( M-KB+J1T ), WORK( M-KB+J1T ), 1067 $ KA1 ) 1068 670 CONTINUE 1069* 1070 IF( WANTX ) THEN 1071* 1072* post-multiply X by product of rotations in 2nd set 1073* 1074 DO 680 J = J1, J2, KA1 1075 CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1, 1076 $ RWORK( M-KB+J ), WORK( M-KB+J ) ) 1077 680 CONTINUE 1078 END IF 1079 690 CONTINUE 1080* 1081 DO 710 K = 1, KB - 1 1082 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1 1083* 1084* finish applying rotations in 1st set from the right 1085* 1086 DO 700 L = KB - K, 1, -1 1087 NRT = ( J2+L-1 ) / KA1 1088 J1T = J2 - ( NRT-1 )*KA1 1089 IF( NRT.GT.0 ) 1090 $ CALL ZLARTV( NRT, AB( L, J1T ), INCA, 1091 $ AB( L+1, J1T-1 ), INCA, RWORK( J1T ), 1092 $ WORK( J1T ), KA1 ) 1093 700 CONTINUE 1094 710 CONTINUE 1095* 1096 IF( KB.GT.1 ) THEN 1097 DO 720 J = 2, I2 - KA 1098 RWORK( J ) = RWORK( J+KA ) 1099 WORK( J ) = WORK( J+KA ) 1100 720 CONTINUE 1101 END IF 1102* 1103 ELSE 1104* 1105* Transform A, working with the lower triangle 1106* 1107 IF( UPDATE ) THEN 1108* 1109* Form inv(S(i))**H * A * inv(S(i)) 1110* 1111 BII = DBLE( BB( 1, I ) ) 1112 AB( 1, I ) = ( DBLE( AB( 1, I ) ) / BII ) / BII 1113 DO 730 J = I1, I - 1 1114 AB( I-J+1, J ) = AB( I-J+1, J ) / BII 1115 730 CONTINUE 1116 DO 740 J = I + 1, MIN( N, I+KA ) 1117 AB( J-I+1, I ) = AB( J-I+1, I ) / BII 1118 740 CONTINUE 1119 DO 770 K = I + 1, I + KBT 1120 DO 750 J = K, I + KBT 1121 AB( J-K+1, K ) = AB( J-K+1, K ) - 1122 $ BB( J-I+1, I )*DCONJG( AB( K-I+1, 1123 $ I ) ) - DCONJG( BB( K-I+1, I ) )* 1124 $ AB( J-I+1, I ) + DBLE( AB( 1, I ) )* 1125 $ BB( J-I+1, I )*DCONJG( BB( K-I+1, 1126 $ I ) ) 1127 750 CONTINUE 1128 DO 760 J = I + KBT + 1, MIN( N, I+KA ) 1129 AB( J-K+1, K ) = AB( J-K+1, K ) - 1130 $ DCONJG( BB( K-I+1, I ) )* 1131 $ AB( J-I+1, I ) 1132 760 CONTINUE 1133 770 CONTINUE 1134 DO 790 J = I1, I 1135 DO 780 K = I + 1, MIN( J+KA, I+KBT ) 1136 AB( K-J+1, J ) = AB( K-J+1, J ) - 1137 $ BB( K-I+1, I )*AB( I-J+1, J ) 1138 780 CONTINUE 1139 790 CONTINUE 1140* 1141 IF( WANTX ) THEN 1142* 1143* post-multiply X by inv(S(i)) 1144* 1145 CALL ZDSCAL( NX, ONE / BII, X( 1, I ), 1 ) 1146 IF( KBT.GT.0 ) 1147 $ CALL ZGERC( NX, KBT, -CONE, X( 1, I ), 1, BB( 2, I ), 1148 $ 1, X( 1, I+1 ), LDX ) 1149 END IF 1150* 1151* store a(i,i1) in RA1 for use in next loop over K 1152* 1153 RA1 = AB( I-I1+1, I1 ) 1154 END IF 1155* 1156* Generate and apply vectors of rotations to chase all the 1157* existing bulges KA positions up toward the top of the band 1158* 1159 DO 840 K = 1, KB - 1 1160 IF( UPDATE ) THEN 1161* 1162* Determine the rotations which would annihilate the bulge 1163* which has in theory just been created 1164* 1165 IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN 1166* 1167* generate rotation to annihilate a(i,i+k-ka-1) 1168* 1169 CALL ZLARTG( AB( KA1-K, I+K-KA ), RA1, 1170 $ RWORK( I+K-KA ), WORK( I+K-KA ), RA ) 1171* 1172* create nonzero element a(i+k,i+k-ka-1) outside the 1173* band and store it in WORK(m-kb+i+k) 1174* 1175 T = -BB( K+1, I )*RA1 1176 WORK( M-KB+I+K ) = RWORK( I+K-KA )*T - 1177 $ DCONJG( WORK( I+K-KA ) )* 1178 $ AB( KA1, I+K-KA ) 1179 AB( KA1, I+K-KA ) = WORK( I+K-KA )*T + 1180 $ RWORK( I+K-KA )*AB( KA1, I+K-KA ) 1181 RA1 = RA 1182 END IF 1183 END IF 1184 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1 1185 NR = ( J2+KA-1 ) / KA1 1186 J1 = J2 - ( NR-1 )*KA1 1187 IF( UPDATE ) THEN 1188 J2T = MIN( J2, I-2*KA+K-1 ) 1189 ELSE 1190 J2T = J2 1191 END IF 1192 NRT = ( J2T+KA-1 ) / KA1 1193 DO 800 J = J1, J2T, KA1 1194* 1195* create nonzero element a(j+ka,j-1) outside the band 1196* and store it in WORK(j) 1197* 1198 WORK( J ) = WORK( J )*AB( KA1, J-1 ) 1199 AB( KA1, J-1 ) = RWORK( J )*AB( KA1, J-1 ) 1200 800 CONTINUE 1201* 1202* generate rotations in 1st set to annihilate elements which 1203* have been created outside the band 1204* 1205 IF( NRT.GT.0 ) 1206 $ CALL ZLARGV( NRT, AB( KA1, J1 ), INCA, WORK( J1 ), KA1, 1207 $ RWORK( J1 ), KA1 ) 1208 IF( NR.GT.0 ) THEN 1209* 1210* apply rotations in 1st set from the right 1211* 1212 DO 810 L = 1, KA - 1 1213 CALL ZLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ), 1214 $ INCA, RWORK( J1 ), WORK( J1 ), KA1 ) 1215 810 CONTINUE 1216* 1217* apply rotations in 1st set from both sides to diagonal 1218* blocks 1219* 1220 CALL ZLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ), 1221 $ AB( 2, J1-1 ), INCA, RWORK( J1 ), 1222 $ WORK( J1 ), KA1 ) 1223* 1224 CALL ZLACGV( NR, WORK( J1 ), KA1 ) 1225 END IF 1226* 1227* start applying rotations in 1st set from the left 1228* 1229 DO 820 L = KA - 1, KB - K + 1, -1 1230 NRT = ( J2+L-1 ) / KA1 1231 J1T = J2 - ( NRT-1 )*KA1 1232 IF( NRT.GT.0 ) 1233 $ CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA, 1234 $ AB( KA1-L, J1T-KA1+L ), INCA, 1235 $ RWORK( J1T ), WORK( J1T ), KA1 ) 1236 820 CONTINUE 1237* 1238 IF( WANTX ) THEN 1239* 1240* post-multiply X by product of rotations in 1st set 1241* 1242 DO 830 J = J1, J2, KA1 1243 CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1, 1244 $ RWORK( J ), DCONJG( WORK( J ) ) ) 1245 830 CONTINUE 1246 END IF 1247 840 CONTINUE 1248* 1249 IF( UPDATE ) THEN 1250 IF( I2.GT.0 .AND. KBT.GT.0 ) THEN 1251* 1252* create nonzero element a(i+kbt,i+kbt-ka-1) outside the 1253* band and store it in WORK(m-kb+i+kbt) 1254* 1255 WORK( M-KB+I+KBT ) = -BB( KBT+1, I )*RA1 1256 END IF 1257 END IF 1258* 1259 DO 880 K = KB, 1, -1 1260 IF( UPDATE ) THEN 1261 J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1 1262 ELSE 1263 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1 1264 END IF 1265* 1266* finish applying rotations in 2nd set from the left 1267* 1268 DO 850 L = KB - K, 1, -1 1269 NRT = ( J2+KA+L-1 ) / KA1 1270 J1T = J2 - ( NRT-1 )*KA1 1271 IF( NRT.GT.0 ) 1272 $ CALL ZLARTV( NRT, AB( KA1-L+1, J1T+L-1 ), INCA, 1273 $ AB( KA1-L, J1T+L-1 ), INCA, 1274 $ RWORK( M-KB+J1T+KA ), 1275 $ WORK( M-KB+J1T+KA ), KA1 ) 1276 850 CONTINUE 1277 NR = ( J2+KA-1 ) / KA1 1278 J1 = J2 - ( NR-1 )*KA1 1279 DO 860 J = J1, J2, KA1 1280 WORK( M-KB+J ) = WORK( M-KB+J+KA ) 1281 RWORK( M-KB+J ) = RWORK( M-KB+J+KA ) 1282 860 CONTINUE 1283 DO 870 J = J1, J2, KA1 1284* 1285* create nonzero element a(j+ka,j-1) outside the band 1286* and store it in WORK(m-kb+j) 1287* 1288 WORK( M-KB+J ) = WORK( M-KB+J )*AB( KA1, J-1 ) 1289 AB( KA1, J-1 ) = RWORK( M-KB+J )*AB( KA1, J-1 ) 1290 870 CONTINUE 1291 IF( UPDATE ) THEN 1292 IF( I+K.GT.KA1 .AND. K.LE.KBT ) 1293 $ WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K ) 1294 END IF 1295 880 CONTINUE 1296* 1297 DO 920 K = KB, 1, -1 1298 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1 1299 NR = ( J2+KA-1 ) / KA1 1300 J1 = J2 - ( NR-1 )*KA1 1301 IF( NR.GT.0 ) THEN 1302* 1303* generate rotations in 2nd set to annihilate elements 1304* which have been created outside the band 1305* 1306 CALL ZLARGV( NR, AB( KA1, J1 ), INCA, WORK( M-KB+J1 ), 1307 $ KA1, RWORK( M-KB+J1 ), KA1 ) 1308* 1309* apply rotations in 2nd set from the right 1310* 1311 DO 890 L = 1, KA - 1 1312 CALL ZLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ), 1313 $ INCA, RWORK( M-KB+J1 ), WORK( M-KB+J1 ), 1314 $ KA1 ) 1315 890 CONTINUE 1316* 1317* apply rotations in 2nd set from both sides to diagonal 1318* blocks 1319* 1320 CALL ZLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ), 1321 $ AB( 2, J1-1 ), INCA, RWORK( M-KB+J1 ), 1322 $ WORK( M-KB+J1 ), KA1 ) 1323* 1324 CALL ZLACGV( NR, WORK( M-KB+J1 ), KA1 ) 1325 END IF 1326* 1327* start applying rotations in 2nd set from the left 1328* 1329 DO 900 L = KA - 1, KB - K + 1, -1 1330 NRT = ( J2+L-1 ) / KA1 1331 J1T = J2 - ( NRT-1 )*KA1 1332 IF( NRT.GT.0 ) 1333 $ CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA, 1334 $ AB( KA1-L, J1T-KA1+L ), INCA, 1335 $ RWORK( M-KB+J1T ), WORK( M-KB+J1T ), 1336 $ KA1 ) 1337 900 CONTINUE 1338* 1339 IF( WANTX ) THEN 1340* 1341* post-multiply X by product of rotations in 2nd set 1342* 1343 DO 910 J = J1, J2, KA1 1344 CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1, 1345 $ RWORK( M-KB+J ), DCONJG( WORK( M-KB+J ) ) ) 1346 910 CONTINUE 1347 END IF 1348 920 CONTINUE 1349* 1350 DO 940 K = 1, KB - 1 1351 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1 1352* 1353* finish applying rotations in 1st set from the left 1354* 1355 DO 930 L = KB - K, 1, -1 1356 NRT = ( J2+L-1 ) / KA1 1357 J1T = J2 - ( NRT-1 )*KA1 1358 IF( NRT.GT.0 ) 1359 $ CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA, 1360 $ AB( KA1-L, J1T-KA1+L ), INCA, 1361 $ RWORK( J1T ), WORK( J1T ), KA1 ) 1362 930 CONTINUE 1363 940 CONTINUE 1364* 1365 IF( KB.GT.1 ) THEN 1366 DO 950 J = 2, I2 - KA 1367 RWORK( J ) = RWORK( J+KA ) 1368 WORK( J ) = WORK( J+KA ) 1369 950 CONTINUE 1370 END IF 1371* 1372 END IF 1373* 1374 GO TO 490 1375* 1376* End of ZHBGST 1377* 1378 END 1379