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