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*> \ingroup complex16OTHERcomputational 161* 162* ===================================================================== 163 SUBROUTINE ZHBGST( 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 DOUBLE PRECISION RWORK( * ) 176 COMPLEX*16 AB( LDAB, * ), BB( LDBB, * ), WORK( * ), 177 $ X( LDX, * ) 178* .. 179* 180* ===================================================================== 181* 182* .. Parameters .. 183 COMPLEX*16 CZERO, CONE 184 DOUBLE PRECISION ONE 185 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 186 $ CONE = ( 1.0D+0, 0.0D+0 ), ONE = 1.0D+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 DOUBLE PRECISION BII 193 COMPLEX*16 RA, RA1, T 194* .. 195* .. External Functions .. 196 LOGICAL LSAME 197 EXTERNAL LSAME 198* .. 199* .. External Subroutines .. 200 EXTERNAL XERBLA, ZDSCAL, ZGERC, ZGERU, ZLACGV, ZLAR2V, 201 $ ZLARGV, ZLARTG, ZLARTV, ZLASET, ZROT 202* .. 203* .. Intrinsic Functions .. 204 INTRINSIC DBLE, DCONJG, MAX, MIN 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( 'ZHBGST', -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 ZLASET( '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 ZPBSTF. 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 = DBLE( BB( KB1, I ) ) 346 AB( KA1, I ) = ( DBLE( 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 $ DCONJG( AB( K-I+KA1, I ) ) - 358 $ DCONJG( BB( K-I+KB1, I ) )* 359 $ AB( J-I+KA1, I ) + 360 $ DBLE( AB( KA1, I ) )* 361 $ BB( J-I+KB1, I )* 362 $ DCONJG( 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 $ DCONJG( 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 ZDSCAL( N-M, ONE / BII, X( M+1, I ), 1 ) 382 IF( KBT.GT.0 ) 383 $ CALL ZGERC( 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 ZLARTG( 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 $ DCONJG( 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 ZLARGV( 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 ZLARTV( 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 ZLAR2V( 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 ZLACGV( 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 ZLARTV( 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 ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1, 482 $ RWORK( J-M ), DCONJG( 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 ZLARTV( 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 ZLARGV( 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 ZLARTV( 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 ZLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ), 557 $ AB( KA, J2+1 ), INCA, RWORK( J2 ), 558 $ WORK( J2 ), KA1 ) 559* 560 CALL ZLACGV( 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 ZLARTV( 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 ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1, 579 $ RWORK( J ), DCONJG( 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 ZLARTV( 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 = DBLE( BB( 1, I ) ) 614 AB( 1, I ) = ( DBLE( 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 )*DCONJG( AB( I-K+1, 625 $ K ) ) - DCONJG( BB( I-K+1, K ) )* 626 $ AB( I-J+1, J ) + DBLE( AB( 1, I ) )* 627 $ BB( I-J+1, J )*DCONJG( 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 $ DCONJG( 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 ZDSCAL( N-M, ONE / BII, X( M+1, I ), 1 ) 648 IF( KBT.GT.0 ) 649 $ CALL ZGERU( 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 ZLARTG( 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 $ DCONJG( WORK( I-K+KA-M ) )* 682 $ AB( KA1, I-K ) 683 AB( KA1, I-K ) = WORK( I-K+KA-M )*T + 684 $ RWORK( I-K+KA-M )*AB( KA1, I-K ) 685 RA1 = RA 686 END IF 687 END IF 688 J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1 689 NR = ( N-J2+KA ) / KA1 690 J1 = J2 + ( NR-1 )*KA1 691 IF( UPDATE ) THEN 692 J2T = MAX( J2, I+2*KA-K+1 ) 693 ELSE 694 J2T = J2 695 END IF 696 NRT = ( N-J2T+KA ) / KA1 697 DO 320 J = J2T, J1, KA1 698* 699* create nonzero element a(j+1,j-ka) outside the band 700* and store it in WORK(j-m) 701* 702 WORK( J-M ) = WORK( J-M )*AB( KA1, J-KA+1 ) 703 AB( KA1, J-KA+1 ) = RWORK( J-M )*AB( KA1, J-KA+1 ) 704 320 CONTINUE 705* 706* generate rotations in 1st set to annihilate elements which 707* have been created outside the band 708* 709 IF( NRT.GT.0 ) 710 $ CALL ZLARGV( NRT, AB( KA1, J2T-KA ), INCA, WORK( J2T-M ), 711 $ KA1, RWORK( J2T-M ), KA1 ) 712 IF( NR.GT.0 ) THEN 713* 714* apply rotations in 1st set from the left 715* 716 DO 330 L = 1, KA - 1 717 CALL ZLARTV( NR, AB( L+1, J2-L ), INCA, 718 $ AB( L+2, J2-L ), INCA, RWORK( J2-M ), 719 $ WORK( J2-M ), KA1 ) 720 330 CONTINUE 721* 722* apply rotations in 1st set from both sides to diagonal 723* blocks 724* 725 CALL ZLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ), 726 $ INCA, RWORK( J2-M ), WORK( J2-M ), KA1 ) 727* 728 CALL ZLACGV( NR, WORK( J2-M ), KA1 ) 729 END IF 730* 731* start applying rotations in 1st set from the right 732* 733 DO 340 L = KA - 1, KB - K + 1, -1 734 NRT = ( N-J2+L ) / KA1 735 IF( NRT.GT.0 ) 736 $ CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA, 737 $ AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ), 738 $ WORK( J2-M ), KA1 ) 739 340 CONTINUE 740* 741 IF( WANTX ) THEN 742* 743* post-multiply X by product of rotations in 1st set 744* 745 DO 350 J = J2, J1, KA1 746 CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1, 747 $ RWORK( J-M ), WORK( J-M ) ) 748 350 CONTINUE 749 END IF 750 360 CONTINUE 751* 752 IF( UPDATE ) THEN 753 IF( I2.LE.N .AND. KBT.GT.0 ) THEN 754* 755* create nonzero element a(i-kbt+ka+1,i-kbt) outside the 756* band and store it in WORK(i-kbt) 757* 758 WORK( I-KBT ) = -BB( KBT+1, I-KBT )*RA1 759 END IF 760 END IF 761* 762 DO 400 K = KB, 1, -1 763 IF( UPDATE ) THEN 764 J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1 765 ELSE 766 J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1 767 END IF 768* 769* finish applying rotations in 2nd set from the right 770* 771 DO 370 L = KB - K, 1, -1 772 NRT = ( N-J2+KA+L ) / KA1 773 IF( NRT.GT.0 ) 774 $ CALL ZLARTV( NRT, AB( KA1-L+1, J2-KA ), INCA, 775 $ AB( KA1-L, J2-KA+1 ), INCA, 776 $ RWORK( J2-KA ), WORK( J2-KA ), KA1 ) 777 370 CONTINUE 778 NR = ( N-J2+KA ) / KA1 779 J1 = J2 + ( NR-1 )*KA1 780 DO 380 J = J1, J2, -KA1 781 WORK( J ) = WORK( J-KA ) 782 RWORK( J ) = RWORK( J-KA ) 783 380 CONTINUE 784 DO 390 J = J2, J1, KA1 785* 786* create nonzero element a(j+1,j-ka) outside the band 787* and store it in WORK(j) 788* 789 WORK( J ) = WORK( J )*AB( KA1, J-KA+1 ) 790 AB( KA1, J-KA+1 ) = RWORK( J )*AB( KA1, J-KA+1 ) 791 390 CONTINUE 792 IF( UPDATE ) THEN 793 IF( I-K.LT.N-KA .AND. K.LE.KBT ) 794 $ WORK( I-K+KA ) = WORK( I-K ) 795 END IF 796 400 CONTINUE 797* 798 DO 440 K = KB, 1, -1 799 J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1 800 NR = ( N-J2+KA ) / KA1 801 J1 = J2 + ( NR-1 )*KA1 802 IF( NR.GT.0 ) THEN 803* 804* generate rotations in 2nd set to annihilate elements 805* which have been created outside the band 806* 807 CALL ZLARGV( NR, AB( KA1, J2-KA ), INCA, WORK( J2 ), KA1, 808 $ RWORK( J2 ), KA1 ) 809* 810* apply rotations in 2nd set from the left 811* 812 DO 410 L = 1, KA - 1 813 CALL ZLARTV( NR, AB( L+1, J2-L ), INCA, 814 $ AB( L+2, J2-L ), INCA, RWORK( J2 ), 815 $ WORK( J2 ), KA1 ) 816 410 CONTINUE 817* 818* apply rotations in 2nd set from both sides to diagonal 819* blocks 820* 821 CALL ZLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ), 822 $ INCA, RWORK( J2 ), WORK( J2 ), KA1 ) 823* 824 CALL ZLACGV( NR, WORK( J2 ), KA1 ) 825 END IF 826* 827* start applying rotations in 2nd set from the right 828* 829 DO 420 L = KA - 1, KB - K + 1, -1 830 NRT = ( N-J2+L ) / KA1 831 IF( NRT.GT.0 ) 832 $ CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA, 833 $ AB( KA1-L, J2+1 ), INCA, RWORK( J2 ), 834 $ WORK( J2 ), KA1 ) 835 420 CONTINUE 836* 837 IF( WANTX ) THEN 838* 839* post-multiply X by product of rotations in 2nd set 840* 841 DO 430 J = J2, J1, KA1 842 CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1, 843 $ RWORK( J ), WORK( J ) ) 844 430 CONTINUE 845 END IF 846 440 CONTINUE 847* 848 DO 460 K = 1, KB - 1 849 J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1 850* 851* finish applying rotations in 1st set from the right 852* 853 DO 450 L = KB - K, 1, -1 854 NRT = ( N-J2+L ) / KA1 855 IF( NRT.GT.0 ) 856 $ CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA, 857 $ AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ), 858 $ WORK( J2-M ), KA1 ) 859 450 CONTINUE 860 460 CONTINUE 861* 862 IF( KB.GT.1 ) THEN 863 DO 470 J = N - 1, J2 + KA, -1 864 RWORK( J-M ) = RWORK( J-KA-M ) 865 WORK( J-M ) = WORK( J-KA-M ) 866 470 CONTINUE 867 END IF 868* 869 END IF 870* 871 GO TO 10 872* 873 480 CONTINUE 874* 875* **************************** Phase 2 ***************************** 876* 877* The logical structure of this phase is: 878* 879* UPDATE = .TRUE. 880* DO I = 1, M 881* use S(i) to update A and create a new bulge 882* apply rotations to push all bulges KA positions upward 883* END DO 884* UPDATE = .FALSE. 885* DO I = M - KA - 1, 2, -1 886* apply rotations to push all bulges KA positions upward 887* END DO 888* 889* To avoid duplicating code, the two loops are merged. 890* 891 UPDATE = .TRUE. 892 I = 0 893 490 CONTINUE 894 IF( UPDATE ) THEN 895 I = I + 1 896 KBT = MIN( KB, M-I ) 897 I0 = I + 1 898 I1 = MAX( 1, I-KA ) 899 I2 = I + KBT - KA1 900 IF( I.GT.M ) THEN 901 UPDATE = .FALSE. 902 I = I - 1 903 I0 = M + 1 904 IF( KA.EQ.0 ) 905 $ RETURN 906 GO TO 490 907 END IF 908 ELSE 909 I = I - KA 910 IF( I.LT.2 ) 911 $ RETURN 912 END IF 913* 914 IF( I.LT.M-KBT ) THEN 915 NX = M 916 ELSE 917 NX = N 918 END IF 919* 920 IF( UPPER ) THEN 921* 922* Transform A, working with the upper triangle 923* 924 IF( UPDATE ) THEN 925* 926* Form inv(S(i))**H * A * inv(S(i)) 927* 928 BII = DBLE( BB( KB1, I ) ) 929 AB( KA1, I ) = ( DBLE( AB( KA1, I ) ) / BII ) / BII 930 DO 500 J = I1, I - 1 931 AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII 932 500 CONTINUE 933 DO 510 J = I + 1, MIN( N, I+KA ) 934 AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII 935 510 CONTINUE 936 DO 540 K = I + 1, I + KBT 937 DO 520 J = K, I + KBT 938 AB( K-J+KA1, J ) = AB( K-J+KA1, J ) - 939 $ BB( I-J+KB1, J )* 940 $ DCONJG( AB( I-K+KA1, K ) ) - 941 $ DCONJG( BB( I-K+KB1, K ) )* 942 $ AB( I-J+KA1, J ) + 943 $ DBLE( AB( KA1, I ) )* 944 $ BB( I-J+KB1, J )* 945 $ DCONJG( BB( I-K+KB1, K ) ) 946 520 CONTINUE 947 DO 530 J = I + KBT + 1, MIN( N, I+KA ) 948 AB( K-J+KA1, J ) = AB( K-J+KA1, J ) - 949 $ DCONJG( BB( I-K+KB1, K ) )* 950 $ AB( I-J+KA1, J ) 951 530 CONTINUE 952 540 CONTINUE 953 DO 560 J = I1, I 954 DO 550 K = I + 1, MIN( J+KA, I+KBT ) 955 AB( J-K+KA1, K ) = AB( J-K+KA1, K ) - 956 $ BB( I-K+KB1, K )*AB( J-I+KA1, I ) 957 550 CONTINUE 958 560 CONTINUE 959* 960 IF( WANTX ) THEN 961* 962* post-multiply X by inv(S(i)) 963* 964 CALL ZDSCAL( NX, ONE / BII, X( 1, I ), 1 ) 965 IF( KBT.GT.0 ) 966 $ CALL ZGERU( NX, KBT, -CONE, X( 1, I ), 1, 967 $ BB( KB, I+1 ), LDBB-1, X( 1, I+1 ), LDX ) 968 END IF 969* 970* store a(i1,i) in RA1 for use in next loop over K 971* 972 RA1 = AB( I1-I+KA1, I ) 973 END IF 974* 975* Generate and apply vectors of rotations to chase all the 976* existing bulges KA positions up toward the top of the band 977* 978 DO 610 K = 1, KB - 1 979 IF( UPDATE ) THEN 980* 981* Determine the rotations which would annihilate the bulge 982* which has in theory just been created 983* 984 IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN 985* 986* generate rotation to annihilate a(i+k-ka-1,i) 987* 988 CALL ZLARTG( AB( K+1, I ), RA1, RWORK( I+K-KA ), 989 $ WORK( I+K-KA ), RA ) 990* 991* create nonzero element a(i+k-ka-1,i+k) outside the 992* band and store it in WORK(m-kb+i+k) 993* 994 T = -BB( KB1-K, I+K )*RA1 995 WORK( M-KB+I+K ) = RWORK( I+K-KA )*T - 996 $ DCONJG( WORK( I+K-KA ) )* 997 $ AB( 1, I+K ) 998 AB( 1, I+K ) = WORK( I+K-KA )*T + 999 $ RWORK( I+K-KA )*AB( 1, I+K ) 1000 RA1 = RA 1001 END IF 1002 END IF 1003 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1 1004 NR = ( J2+KA-1 ) / KA1 1005 J1 = J2 - ( NR-1 )*KA1 1006 IF( UPDATE ) THEN 1007 J2T = MIN( J2, I-2*KA+K-1 ) 1008 ELSE 1009 J2T = J2 1010 END IF 1011 NRT = ( J2T+KA-1 ) / KA1 1012 DO 570 J = J1, J2T, KA1 1013* 1014* create nonzero element a(j-1,j+ka) outside the band 1015* and store it in WORK(j) 1016* 1017 WORK( J ) = WORK( J )*AB( 1, J+KA-1 ) 1018 AB( 1, J+KA-1 ) = RWORK( J )*AB( 1, J+KA-1 ) 1019 570 CONTINUE 1020* 1021* generate rotations in 1st set to annihilate elements which 1022* have been created outside the band 1023* 1024 IF( NRT.GT.0 ) 1025 $ CALL ZLARGV( NRT, AB( 1, J1+KA ), INCA, WORK( J1 ), KA1, 1026 $ RWORK( J1 ), KA1 ) 1027 IF( NR.GT.0 ) THEN 1028* 1029* apply rotations in 1st set from the left 1030* 1031 DO 580 L = 1, KA - 1 1032 CALL ZLARTV( NR, AB( KA1-L, J1+L ), INCA, 1033 $ AB( KA-L, J1+L ), INCA, RWORK( J1 ), 1034 $ WORK( J1 ), KA1 ) 1035 580 CONTINUE 1036* 1037* apply rotations in 1st set from both sides to diagonal 1038* blocks 1039* 1040 CALL ZLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ), 1041 $ AB( KA, J1 ), INCA, RWORK( J1 ), WORK( J1 ), 1042 $ KA1 ) 1043* 1044 CALL ZLACGV( NR, WORK( J1 ), KA1 ) 1045 END IF 1046* 1047* start applying rotations in 1st set from the right 1048* 1049 DO 590 L = KA - 1, KB - K + 1, -1 1050 NRT = ( J2+L-1 ) / KA1 1051 J1T = J2 - ( NRT-1 )*KA1 1052 IF( NRT.GT.0 ) 1053 $ CALL ZLARTV( NRT, AB( L, J1T ), INCA, 1054 $ AB( L+1, J1T-1 ), INCA, RWORK( J1T ), 1055 $ WORK( J1T ), KA1 ) 1056 590 CONTINUE 1057* 1058 IF( WANTX ) THEN 1059* 1060* post-multiply X by product of rotations in 1st set 1061* 1062 DO 600 J = J1, J2, KA1 1063 CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1, 1064 $ RWORK( J ), WORK( J ) ) 1065 600 CONTINUE 1066 END IF 1067 610 CONTINUE 1068* 1069 IF( UPDATE ) THEN 1070 IF( I2.GT.0 .AND. KBT.GT.0 ) THEN 1071* 1072* create nonzero element a(i+kbt-ka-1,i+kbt) outside the 1073* band and store it in WORK(m-kb+i+kbt) 1074* 1075 WORK( M-KB+I+KBT ) = -BB( KB1-KBT, I+KBT )*RA1 1076 END IF 1077 END IF 1078* 1079 DO 650 K = KB, 1, -1 1080 IF( UPDATE ) THEN 1081 J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1 1082 ELSE 1083 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1 1084 END IF 1085* 1086* finish applying rotations in 2nd set from the right 1087* 1088 DO 620 L = KB - K, 1, -1 1089 NRT = ( J2+KA+L-1 ) / KA1 1090 J1T = J2 - ( NRT-1 )*KA1 1091 IF( NRT.GT.0 ) 1092 $ CALL ZLARTV( NRT, AB( L, J1T+KA ), INCA, 1093 $ AB( L+1, J1T+KA-1 ), INCA, 1094 $ RWORK( M-KB+J1T+KA ), 1095 $ WORK( M-KB+J1T+KA ), KA1 ) 1096 620 CONTINUE 1097 NR = ( J2+KA-1 ) / KA1 1098 J1 = J2 - ( NR-1 )*KA1 1099 DO 630 J = J1, J2, KA1 1100 WORK( M-KB+J ) = WORK( M-KB+J+KA ) 1101 RWORK( M-KB+J ) = RWORK( M-KB+J+KA ) 1102 630 CONTINUE 1103 DO 640 J = J1, J2, KA1 1104* 1105* create nonzero element a(j-1,j+ka) outside the band 1106* and store it in WORK(m-kb+j) 1107* 1108 WORK( M-KB+J ) = WORK( M-KB+J )*AB( 1, J+KA-1 ) 1109 AB( 1, J+KA-1 ) = RWORK( M-KB+J )*AB( 1, J+KA-1 ) 1110 640 CONTINUE 1111 IF( UPDATE ) THEN 1112 IF( I+K.GT.KA1 .AND. K.LE.KBT ) 1113 $ WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K ) 1114 END IF 1115 650 CONTINUE 1116* 1117 DO 690 K = KB, 1, -1 1118 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1 1119 NR = ( J2+KA-1 ) / KA1 1120 J1 = J2 - ( NR-1 )*KA1 1121 IF( NR.GT.0 ) THEN 1122* 1123* generate rotations in 2nd set to annihilate elements 1124* which have been created outside the band 1125* 1126 CALL ZLARGV( NR, AB( 1, J1+KA ), INCA, WORK( M-KB+J1 ), 1127 $ KA1, RWORK( M-KB+J1 ), KA1 ) 1128* 1129* apply rotations in 2nd set from the left 1130* 1131 DO 660 L = 1, KA - 1 1132 CALL ZLARTV( NR, AB( KA1-L, J1+L ), INCA, 1133 $ AB( KA-L, J1+L ), INCA, RWORK( M-KB+J1 ), 1134 $ WORK( M-KB+J1 ), KA1 ) 1135 660 CONTINUE 1136* 1137* apply rotations in 2nd set from both sides to diagonal 1138* blocks 1139* 1140 CALL ZLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ), 1141 $ AB( KA, J1 ), INCA, RWORK( M-KB+J1 ), 1142 $ WORK( M-KB+J1 ), KA1 ) 1143* 1144 CALL ZLACGV( NR, WORK( M-KB+J1 ), KA1 ) 1145 END IF 1146* 1147* start applying rotations in 2nd set from the right 1148* 1149 DO 670 L = KA - 1, KB - K + 1, -1 1150 NRT = ( J2+L-1 ) / KA1 1151 J1T = J2 - ( NRT-1 )*KA1 1152 IF( NRT.GT.0 ) 1153 $ CALL ZLARTV( NRT, AB( L, J1T ), INCA, 1154 $ AB( L+1, J1T-1 ), INCA, 1155 $ RWORK( M-KB+J1T ), WORK( M-KB+J1T ), 1156 $ KA1 ) 1157 670 CONTINUE 1158* 1159 IF( WANTX ) THEN 1160* 1161* post-multiply X by product of rotations in 2nd set 1162* 1163 DO 680 J = J1, J2, KA1 1164 CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1, 1165 $ RWORK( M-KB+J ), WORK( M-KB+J ) ) 1166 680 CONTINUE 1167 END IF 1168 690 CONTINUE 1169* 1170 DO 710 K = 1, KB - 1 1171 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1 1172* 1173* finish applying rotations in 1st set from the right 1174* 1175 DO 700 L = KB - K, 1, -1 1176 NRT = ( J2+L-1 ) / KA1 1177 J1T = J2 - ( NRT-1 )*KA1 1178 IF( NRT.GT.0 ) 1179 $ CALL ZLARTV( NRT, AB( L, J1T ), INCA, 1180 $ AB( L+1, J1T-1 ), INCA, RWORK( J1T ), 1181 $ WORK( J1T ), KA1 ) 1182 700 CONTINUE 1183 710 CONTINUE 1184* 1185 IF( KB.GT.1 ) THEN 1186 DO 720 J = 2, I2 - KA 1187 RWORK( J ) = RWORK( J+KA ) 1188 WORK( J ) = WORK( J+KA ) 1189 720 CONTINUE 1190 END IF 1191* 1192 ELSE 1193* 1194* Transform A, working with the lower triangle 1195* 1196 IF( UPDATE ) THEN 1197* 1198* Form inv(S(i))**H * A * inv(S(i)) 1199* 1200 BII = DBLE( BB( 1, I ) ) 1201 AB( 1, I ) = ( DBLE( AB( 1, I ) ) / BII ) / BII 1202 DO 730 J = I1, I - 1 1203 AB( I-J+1, J ) = AB( I-J+1, J ) / BII 1204 730 CONTINUE 1205 DO 740 J = I + 1, MIN( N, I+KA ) 1206 AB( J-I+1, I ) = AB( J-I+1, I ) / BII 1207 740 CONTINUE 1208 DO 770 K = I + 1, I + KBT 1209 DO 750 J = K, I + KBT 1210 AB( J-K+1, K ) = AB( J-K+1, K ) - 1211 $ BB( J-I+1, I )*DCONJG( AB( K-I+1, 1212 $ I ) ) - DCONJG( BB( K-I+1, I ) )* 1213 $ AB( J-I+1, I ) + DBLE( AB( 1, I ) )* 1214 $ BB( J-I+1, I )*DCONJG( BB( K-I+1, 1215 $ I ) ) 1216 750 CONTINUE 1217 DO 760 J = I + KBT + 1, MIN( N, I+KA ) 1218 AB( J-K+1, K ) = AB( J-K+1, K ) - 1219 $ DCONJG( BB( K-I+1, I ) )* 1220 $ AB( J-I+1, I ) 1221 760 CONTINUE 1222 770 CONTINUE 1223 DO 790 J = I1, I 1224 DO 780 K = I + 1, MIN( J+KA, I+KBT ) 1225 AB( K-J+1, J ) = AB( K-J+1, J ) - 1226 $ BB( K-I+1, I )*AB( I-J+1, J ) 1227 780 CONTINUE 1228 790 CONTINUE 1229* 1230 IF( WANTX ) THEN 1231* 1232* post-multiply X by inv(S(i)) 1233* 1234 CALL ZDSCAL( NX, ONE / BII, X( 1, I ), 1 ) 1235 IF( KBT.GT.0 ) 1236 $ CALL ZGERC( NX, KBT, -CONE, X( 1, I ), 1, BB( 2, I ), 1237 $ 1, X( 1, I+1 ), LDX ) 1238 END IF 1239* 1240* store a(i,i1) in RA1 for use in next loop over K 1241* 1242 RA1 = AB( I-I1+1, I1 ) 1243 END IF 1244* 1245* Generate and apply vectors of rotations to chase all the 1246* existing bulges KA positions up toward the top of the band 1247* 1248 DO 840 K = 1, KB - 1 1249 IF( UPDATE ) THEN 1250* 1251* Determine the rotations which would annihilate the bulge 1252* which has in theory just been created 1253* 1254 IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN 1255* 1256* generate rotation to annihilate a(i,i+k-ka-1) 1257* 1258 CALL ZLARTG( AB( KA1-K, I+K-KA ), RA1, 1259 $ RWORK( I+K-KA ), WORK( I+K-KA ), RA ) 1260* 1261* create nonzero element a(i+k,i+k-ka-1) outside the 1262* band and store it in WORK(m-kb+i+k) 1263* 1264 T = -BB( K+1, I )*RA1 1265 WORK( M-KB+I+K ) = RWORK( I+K-KA )*T - 1266 $ DCONJG( WORK( I+K-KA ) )* 1267 $ AB( KA1, I+K-KA ) 1268 AB( KA1, I+K-KA ) = WORK( I+K-KA )*T + 1269 $ RWORK( I+K-KA )*AB( KA1, I+K-KA ) 1270 RA1 = RA 1271 END IF 1272 END IF 1273 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1 1274 NR = ( J2+KA-1 ) / KA1 1275 J1 = J2 - ( NR-1 )*KA1 1276 IF( UPDATE ) THEN 1277 J2T = MIN( J2, I-2*KA+K-1 ) 1278 ELSE 1279 J2T = J2 1280 END IF 1281 NRT = ( J2T+KA-1 ) / KA1 1282 DO 800 J = J1, J2T, KA1 1283* 1284* create nonzero element a(j+ka,j-1) outside the band 1285* and store it in WORK(j) 1286* 1287 WORK( J ) = WORK( J )*AB( KA1, J-1 ) 1288 AB( KA1, J-1 ) = RWORK( J )*AB( KA1, J-1 ) 1289 800 CONTINUE 1290* 1291* generate rotations in 1st set to annihilate elements which 1292* have been created outside the band 1293* 1294 IF( NRT.GT.0 ) 1295 $ CALL ZLARGV( NRT, AB( KA1, J1 ), INCA, WORK( J1 ), KA1, 1296 $ RWORK( J1 ), KA1 ) 1297 IF( NR.GT.0 ) THEN 1298* 1299* apply rotations in 1st set from the right 1300* 1301 DO 810 L = 1, KA - 1 1302 CALL ZLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ), 1303 $ INCA, RWORK( J1 ), WORK( J1 ), KA1 ) 1304 810 CONTINUE 1305* 1306* apply rotations in 1st set from both sides to diagonal 1307* blocks 1308* 1309 CALL ZLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ), 1310 $ AB( 2, J1-1 ), INCA, RWORK( J1 ), 1311 $ WORK( J1 ), KA1 ) 1312* 1313 CALL ZLACGV( NR, WORK( J1 ), KA1 ) 1314 END IF 1315* 1316* start applying rotations in 1st set from the left 1317* 1318 DO 820 L = KA - 1, KB - K + 1, -1 1319 NRT = ( J2+L-1 ) / KA1 1320 J1T = J2 - ( NRT-1 )*KA1 1321 IF( NRT.GT.0 ) 1322 $ CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA, 1323 $ AB( KA1-L, J1T-KA1+L ), INCA, 1324 $ RWORK( J1T ), WORK( J1T ), KA1 ) 1325 820 CONTINUE 1326* 1327 IF( WANTX ) THEN 1328* 1329* post-multiply X by product of rotations in 1st set 1330* 1331 DO 830 J = J1, J2, KA1 1332 CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1, 1333 $ RWORK( J ), DCONJG( WORK( J ) ) ) 1334 830 CONTINUE 1335 END IF 1336 840 CONTINUE 1337* 1338 IF( UPDATE ) THEN 1339 IF( I2.GT.0 .AND. KBT.GT.0 ) THEN 1340* 1341* create nonzero element a(i+kbt,i+kbt-ka-1) outside the 1342* band and store it in WORK(m-kb+i+kbt) 1343* 1344 WORK( M-KB+I+KBT ) = -BB( KBT+1, I )*RA1 1345 END IF 1346 END IF 1347* 1348 DO 880 K = KB, 1, -1 1349 IF( UPDATE ) THEN 1350 J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1 1351 ELSE 1352 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1 1353 END IF 1354* 1355* finish applying rotations in 2nd set from the left 1356* 1357 DO 850 L = KB - K, 1, -1 1358 NRT = ( J2+KA+L-1 ) / KA1 1359 J1T = J2 - ( NRT-1 )*KA1 1360 IF( NRT.GT.0 ) 1361 $ CALL ZLARTV( NRT, AB( KA1-L+1, J1T+L-1 ), INCA, 1362 $ AB( KA1-L, J1T+L-1 ), INCA, 1363 $ RWORK( M-KB+J1T+KA ), 1364 $ WORK( M-KB+J1T+KA ), KA1 ) 1365 850 CONTINUE 1366 NR = ( J2+KA-1 ) / KA1 1367 J1 = J2 - ( NR-1 )*KA1 1368 DO 860 J = J1, J2, KA1 1369 WORK( M-KB+J ) = WORK( M-KB+J+KA ) 1370 RWORK( M-KB+J ) = RWORK( M-KB+J+KA ) 1371 860 CONTINUE 1372 DO 870 J = J1, J2, KA1 1373* 1374* create nonzero element a(j+ka,j-1) outside the band 1375* and store it in WORK(m-kb+j) 1376* 1377 WORK( M-KB+J ) = WORK( M-KB+J )*AB( KA1, J-1 ) 1378 AB( KA1, J-1 ) = RWORK( M-KB+J )*AB( KA1, J-1 ) 1379 870 CONTINUE 1380 IF( UPDATE ) THEN 1381 IF( I+K.GT.KA1 .AND. K.LE.KBT ) 1382 $ WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K ) 1383 END IF 1384 880 CONTINUE 1385* 1386 DO 920 K = KB, 1, -1 1387 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1 1388 NR = ( J2+KA-1 ) / KA1 1389 J1 = J2 - ( NR-1 )*KA1 1390 IF( NR.GT.0 ) THEN 1391* 1392* generate rotations in 2nd set to annihilate elements 1393* which have been created outside the band 1394* 1395 CALL ZLARGV( NR, AB( KA1, J1 ), INCA, WORK( M-KB+J1 ), 1396 $ KA1, RWORK( M-KB+J1 ), KA1 ) 1397* 1398* apply rotations in 2nd set from the right 1399* 1400 DO 890 L = 1, KA - 1 1401 CALL ZLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ), 1402 $ INCA, RWORK( M-KB+J1 ), WORK( M-KB+J1 ), 1403 $ KA1 ) 1404 890 CONTINUE 1405* 1406* apply rotations in 2nd set from both sides to diagonal 1407* blocks 1408* 1409 CALL ZLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ), 1410 $ AB( 2, J1-1 ), INCA, RWORK( M-KB+J1 ), 1411 $ WORK( M-KB+J1 ), KA1 ) 1412* 1413 CALL ZLACGV( NR, WORK( M-KB+J1 ), KA1 ) 1414 END IF 1415* 1416* start applying rotations in 2nd set from the left 1417* 1418 DO 900 L = KA - 1, KB - K + 1, -1 1419 NRT = ( J2+L-1 ) / KA1 1420 J1T = J2 - ( NRT-1 )*KA1 1421 IF( NRT.GT.0 ) 1422 $ CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA, 1423 $ AB( KA1-L, J1T-KA1+L ), INCA, 1424 $ RWORK( M-KB+J1T ), WORK( M-KB+J1T ), 1425 $ KA1 ) 1426 900 CONTINUE 1427* 1428 IF( WANTX ) THEN 1429* 1430* post-multiply X by product of rotations in 2nd set 1431* 1432 DO 910 J = J1, J2, KA1 1433 CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1, 1434 $ RWORK( M-KB+J ), DCONJG( WORK( M-KB+J ) ) ) 1435 910 CONTINUE 1436 END IF 1437 920 CONTINUE 1438* 1439 DO 940 K = 1, KB - 1 1440 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1 1441* 1442* finish applying rotations in 1st set from the left 1443* 1444 DO 930 L = KB - K, 1, -1 1445 NRT = ( J2+L-1 ) / KA1 1446 J1T = J2 - ( NRT-1 )*KA1 1447 IF( NRT.GT.0 ) 1448 $ CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA, 1449 $ AB( KA1-L, J1T-KA1+L ), INCA, 1450 $ RWORK( J1T ), WORK( J1T ), KA1 ) 1451 930 CONTINUE 1452 940 CONTINUE 1453* 1454 IF( KB.GT.1 ) THEN 1455 DO 950 J = 2, I2 - KA 1456 RWORK( J ) = RWORK( J+KA ) 1457 WORK( J ) = WORK( J+KA ) 1458 950 CONTINUE 1459 END IF 1460* 1461 END IF 1462* 1463 GO TO 490 1464* 1465* End of ZHBGST 1466* 1467 END 1468