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