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