|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE ZHBGST( VECT, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, X, 00002 $ LDX, WORK, RWORK, INFO ) 00003 * 00004 * -- LAPACK routine (version 3.2.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * June 2010 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER UPLO, VECT 00011 INTEGER INFO, KA, KB, LDAB, LDBB, LDX, N 00012 * .. 00013 * .. Array Arguments .. 00014 DOUBLE PRECISION RWORK( * ) 00015 COMPLEX*16 AB( LDAB, * ), BB( LDBB, * ), WORK( * ), 00016 $ X( LDX, * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * ZHBGST reduces a complex Hermitian-definite banded generalized 00023 * eigenproblem A*x = lambda*B*x to standard form C*y = lambda*y, 00024 * such that C has the same bandwidth as A. 00025 * 00026 * B must have been previously factorized as S**H*S by ZPBSTF, using a 00027 * split Cholesky factorization. A is overwritten by C = X**H*A*X, where 00028 * X = S**(-1)*Q and Q is a unitary matrix chosen to preserve the 00029 * bandwidth of A. 00030 * 00031 * Arguments 00032 * ========= 00033 * 00034 * VECT (input) CHARACTER*1 00035 * = 'N': do not form the transformation matrix X; 00036 * = 'V': form X. 00037 * 00038 * UPLO (input) CHARACTER*1 00039 * = 'U': Upper triangle of A is stored; 00040 * = 'L': Lower triangle of A is stored. 00041 * 00042 * N (input) INTEGER 00043 * The order of the matrices A and B. N >= 0. 00044 * 00045 * KA (input) INTEGER 00046 * The number of superdiagonals of the matrix A if UPLO = 'U', 00047 * or the number of subdiagonals if UPLO = 'L'. KA >= 0. 00048 * 00049 * KB (input) INTEGER 00050 * The number of superdiagonals of the matrix B if UPLO = 'U', 00051 * or the number of subdiagonals if UPLO = 'L'. KA >= KB >= 0. 00052 * 00053 * AB (input/output) COMPLEX*16 array, dimension (LDAB,N) 00054 * On entry, the upper or lower triangle of the Hermitian band 00055 * matrix A, stored in the first ka+1 rows of the array. The 00056 * j-th column of A is stored in the j-th column of the array AB 00057 * as follows: 00058 * if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j; 00059 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+ka). 00060 * 00061 * On exit, the transformed matrix X**H*A*X, stored in the same 00062 * format as A. 00063 * 00064 * LDAB (input) INTEGER 00065 * The leading dimension of the array AB. LDAB >= KA+1. 00066 * 00067 * BB (input) COMPLEX*16 array, dimension (LDBB,N) 00068 * The banded factor S from the split Cholesky factorization of 00069 * B, as returned by ZPBSTF, stored in the first kb+1 rows of 00070 * the array. 00071 * 00072 * LDBB (input) INTEGER 00073 * The leading dimension of the array BB. LDBB >= KB+1. 00074 * 00075 * X (output) COMPLEX*16 array, dimension (LDX,N) 00076 * If VECT = 'V', the n-by-n matrix X. 00077 * If VECT = 'N', the array X is not referenced. 00078 * 00079 * LDX (input) INTEGER 00080 * The leading dimension of the array X. 00081 * LDX >= max(1,N) if VECT = 'V'; LDX >= 1 otherwise. 00082 * 00083 * WORK (workspace) COMPLEX*16 array, dimension (N) 00084 * 00085 * RWORK (workspace) DOUBLE PRECISION array, dimension (N) 00086 * 00087 * INFO (output) INTEGER 00088 * = 0: successful exit 00089 * < 0: if INFO = -i, the i-th argument had an illegal value. 00090 * 00091 * ===================================================================== 00092 * 00093 * .. Parameters .. 00094 COMPLEX*16 CZERO, CONE 00095 DOUBLE PRECISION ONE 00096 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 00097 $ CONE = ( 1.0D+0, 0.0D+0 ), ONE = 1.0D+0 ) 00098 * .. 00099 * .. Local Scalars .. 00100 LOGICAL UPDATE, UPPER, WANTX 00101 INTEGER I, I0, I1, I2, INCA, J, J1, J1T, J2, J2T, K, 00102 $ KA1, KB1, KBT, L, M, NR, NRT, NX 00103 DOUBLE PRECISION BII 00104 COMPLEX*16 RA, RA1, T 00105 * .. 00106 * .. External Functions .. 00107 LOGICAL LSAME 00108 EXTERNAL LSAME 00109 * .. 00110 * .. External Subroutines .. 00111 EXTERNAL XERBLA, ZDSCAL, ZGERC, ZGERU, ZLACGV, ZLAR2V, 00112 $ ZLARGV, ZLARTG, ZLARTV, ZLASET, ZROT 00113 * .. 00114 * .. Intrinsic Functions .. 00115 INTRINSIC DBLE, DCONJG, MAX, MIN 00116 * .. 00117 * .. Executable Statements .. 00118 * 00119 * Test the input parameters 00120 * 00121 WANTX = LSAME( VECT, 'V' ) 00122 UPPER = LSAME( UPLO, 'U' ) 00123 KA1 = KA + 1 00124 KB1 = KB + 1 00125 INFO = 0 00126 IF( .NOT.WANTX .AND. .NOT.LSAME( VECT, 'N' ) ) THEN 00127 INFO = -1 00128 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00129 INFO = -2 00130 ELSE IF( N.LT.0 ) THEN 00131 INFO = -3 00132 ELSE IF( KA.LT.0 ) THEN 00133 INFO = -4 00134 ELSE IF( KB.LT.0 .OR. KB.GT.KA ) THEN 00135 INFO = -5 00136 ELSE IF( LDAB.LT.KA+1 ) THEN 00137 INFO = -7 00138 ELSE IF( LDBB.LT.KB+1 ) THEN 00139 INFO = -9 00140 ELSE IF( LDX.LT.1 .OR. WANTX .AND. LDX.LT.MAX( 1, N ) ) THEN 00141 INFO = -11 00142 END IF 00143 IF( INFO.NE.0 ) THEN 00144 CALL XERBLA( 'ZHBGST', -INFO ) 00145 RETURN 00146 END IF 00147 * 00148 * Quick return if possible 00149 * 00150 IF( N.EQ.0 ) 00151 $ RETURN 00152 * 00153 INCA = LDAB*KA1 00154 * 00155 * Initialize X to the unit matrix, if needed 00156 * 00157 IF( WANTX ) 00158 $ CALL ZLASET( 'Full', N, N, CZERO, CONE, X, LDX ) 00159 * 00160 * Set M to the splitting point m. It must be the same value as is 00161 * used in ZPBSTF. The chosen value allows the arrays WORK and RWORK 00162 * to be of dimension (N). 00163 * 00164 M = ( N+KB ) / 2 00165 * 00166 * The routine works in two phases, corresponding to the two halves 00167 * of the split Cholesky factorization of B as S**H*S where 00168 * 00169 * S = ( U ) 00170 * ( M L ) 00171 * 00172 * with U upper triangular of order m, and L lower triangular of 00173 * order n-m. S has the same bandwidth as B. 00174 * 00175 * S is treated as a product of elementary matrices: 00176 * 00177 * S = S(m)*S(m-1)*...*S(2)*S(1)*S(m+1)*S(m+2)*...*S(n-1)*S(n) 00178 * 00179 * where S(i) is determined by the i-th row of S. 00180 * 00181 * In phase 1, the index i takes the values n, n-1, ... , m+1; 00182 * in phase 2, it takes the values 1, 2, ... , m. 00183 * 00184 * For each value of i, the current matrix A is updated by forming 00185 * inv(S(i))**H*A*inv(S(i)). This creates a triangular bulge outside 00186 * the band of A. The bulge is then pushed down toward the bottom of 00187 * A in phase 1, and up toward the top of A in phase 2, by applying 00188 * plane rotations. 00189 * 00190 * There are kb*(kb+1)/2 elements in the bulge, but at most 2*kb-1 00191 * of them are linearly independent, so annihilating a bulge requires 00192 * only 2*kb-1 plane rotations. The rotations are divided into a 1st 00193 * set of kb-1 rotations, and a 2nd set of kb rotations. 00194 * 00195 * Wherever possible, rotations are generated and applied in vector 00196 * operations of length NR between the indices J1 and J2 (sometimes 00197 * replaced by modified values NRT, J1T or J2T). 00198 * 00199 * The real cosines and complex sines of the rotations are stored in 00200 * the arrays RWORK and WORK, those of the 1st set in elements 00201 * 2:m-kb-1, and those of the 2nd set in elements m-kb+1:n. 00202 * 00203 * The bulges are not formed explicitly; nonzero elements outside the 00204 * band are created only when they are required for generating new 00205 * rotations; they are stored in the array WORK, in positions where 00206 * they are later overwritten by the sines of the rotations which 00207 * annihilate them. 00208 * 00209 * **************************** Phase 1 ***************************** 00210 * 00211 * The logical structure of this phase is: 00212 * 00213 * UPDATE = .TRUE. 00214 * DO I = N, M + 1, -1 00215 * use S(i) to update A and create a new bulge 00216 * apply rotations to push all bulges KA positions downward 00217 * END DO 00218 * UPDATE = .FALSE. 00219 * DO I = M + KA + 1, N - 1 00220 * apply rotations to push all bulges KA positions downward 00221 * END DO 00222 * 00223 * To avoid duplicating code, the two loops are merged. 00224 * 00225 UPDATE = .TRUE. 00226 I = N + 1 00227 10 CONTINUE 00228 IF( UPDATE ) THEN 00229 I = I - 1 00230 KBT = MIN( KB, I-1 ) 00231 I0 = I - 1 00232 I1 = MIN( N, I+KA ) 00233 I2 = I - KBT + KA1 00234 IF( I.LT.M+1 ) THEN 00235 UPDATE = .FALSE. 00236 I = I + 1 00237 I0 = M 00238 IF( KA.EQ.0 ) 00239 $ GO TO 480 00240 GO TO 10 00241 END IF 00242 ELSE 00243 I = I + KA 00244 IF( I.GT.N-1 ) 00245 $ GO TO 480 00246 END IF 00247 * 00248 IF( UPPER ) THEN 00249 * 00250 * Transform A, working with the upper triangle 00251 * 00252 IF( UPDATE ) THEN 00253 * 00254 * Form inv(S(i))**H * A * inv(S(i)) 00255 * 00256 BII = DBLE( BB( KB1, I ) ) 00257 AB( KA1, I ) = ( DBLE( AB( KA1, I ) ) / BII ) / BII 00258 DO 20 J = I + 1, I1 00259 AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII 00260 20 CONTINUE 00261 DO 30 J = MAX( 1, I-KA ), I - 1 00262 AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII 00263 30 CONTINUE 00264 DO 60 K = I - KBT, I - 1 00265 DO 40 J = I - KBT, K 00266 AB( J-K+KA1, K ) = AB( J-K+KA1, K ) - 00267 $ BB( J-I+KB1, I )* 00268 $ DCONJG( AB( K-I+KA1, I ) ) - 00269 $ DCONJG( BB( K-I+KB1, I ) )* 00270 $ AB( J-I+KA1, I ) + 00271 $ DBLE( AB( KA1, I ) )* 00272 $ BB( J-I+KB1, I )* 00273 $ DCONJG( BB( K-I+KB1, I ) ) 00274 40 CONTINUE 00275 DO 50 J = MAX( 1, I-KA ), I - KBT - 1 00276 AB( J-K+KA1, K ) = AB( J-K+KA1, K ) - 00277 $ DCONJG( BB( K-I+KB1, I ) )* 00278 $ AB( J-I+KA1, I ) 00279 50 CONTINUE 00280 60 CONTINUE 00281 DO 80 J = I, I1 00282 DO 70 K = MAX( J-KA, I-KBT ), I - 1 00283 AB( K-J+KA1, J ) = AB( K-J+KA1, J ) - 00284 $ BB( K-I+KB1, I )*AB( I-J+KA1, J ) 00285 70 CONTINUE 00286 80 CONTINUE 00287 * 00288 IF( WANTX ) THEN 00289 * 00290 * post-multiply X by inv(S(i)) 00291 * 00292 CALL ZDSCAL( N-M, ONE / BII, X( M+1, I ), 1 ) 00293 IF( KBT.GT.0 ) 00294 $ CALL ZGERC( N-M, KBT, -CONE, X( M+1, I ), 1, 00295 $ BB( KB1-KBT, I ), 1, X( M+1, I-KBT ), 00296 $ LDX ) 00297 END IF 00298 * 00299 * store a(i,i1) in RA1 for use in next loop over K 00300 * 00301 RA1 = AB( I-I1+KA1, I1 ) 00302 END IF 00303 * 00304 * Generate and apply vectors of rotations to chase all the 00305 * existing bulges KA positions down toward the bottom of the 00306 * band 00307 * 00308 DO 130 K = 1, KB - 1 00309 IF( UPDATE ) THEN 00310 * 00311 * Determine the rotations which would annihilate the bulge 00312 * which has in theory just been created 00313 * 00314 IF( I-K+KA.LT.N .AND. I-K.GT.1 ) THEN 00315 * 00316 * generate rotation to annihilate a(i,i-k+ka+1) 00317 * 00318 CALL ZLARTG( AB( K+1, I-K+KA ), RA1, 00319 $ RWORK( I-K+KA-M ), WORK( I-K+KA-M ), RA ) 00320 * 00321 * create nonzero element a(i-k,i-k+ka+1) outside the 00322 * band and store it in WORK(i-k) 00323 * 00324 T = -BB( KB1-K, I )*RA1 00325 WORK( I-K ) = RWORK( I-K+KA-M )*T - 00326 $ DCONJG( WORK( I-K+KA-M ) )* 00327 $ AB( 1, I-K+KA ) 00328 AB( 1, I-K+KA ) = WORK( I-K+KA-M )*T + 00329 $ RWORK( I-K+KA-M )*AB( 1, I-K+KA ) 00330 RA1 = RA 00331 END IF 00332 END IF 00333 J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1 00334 NR = ( N-J2+KA ) / KA1 00335 J1 = J2 + ( NR-1 )*KA1 00336 IF( UPDATE ) THEN 00337 J2T = MAX( J2, I+2*KA-K+1 ) 00338 ELSE 00339 J2T = J2 00340 END IF 00341 NRT = ( N-J2T+KA ) / KA1 00342 DO 90 J = J2T, J1, KA1 00343 * 00344 * create nonzero element a(j-ka,j+1) outside the band 00345 * and store it in WORK(j-m) 00346 * 00347 WORK( J-M ) = WORK( J-M )*AB( 1, J+1 ) 00348 AB( 1, J+1 ) = RWORK( J-M )*AB( 1, J+1 ) 00349 90 CONTINUE 00350 * 00351 * generate rotations in 1st set to annihilate elements which 00352 * have been created outside the band 00353 * 00354 IF( NRT.GT.0 ) 00355 $ CALL ZLARGV( NRT, AB( 1, J2T ), INCA, WORK( J2T-M ), KA1, 00356 $ RWORK( J2T-M ), KA1 ) 00357 IF( NR.GT.0 ) THEN 00358 * 00359 * apply rotations in 1st set from the right 00360 * 00361 DO 100 L = 1, KA - 1 00362 CALL ZLARTV( NR, AB( KA1-L, J2 ), INCA, 00363 $ AB( KA-L, J2+1 ), INCA, RWORK( J2-M ), 00364 $ WORK( J2-M ), KA1 ) 00365 100 CONTINUE 00366 * 00367 * apply rotations in 1st set from both sides to diagonal 00368 * blocks 00369 * 00370 CALL ZLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ), 00371 $ AB( KA, J2+1 ), INCA, RWORK( J2-M ), 00372 $ WORK( J2-M ), KA1 ) 00373 * 00374 CALL ZLACGV( NR, WORK( J2-M ), KA1 ) 00375 END IF 00376 * 00377 * start applying rotations in 1st set from the left 00378 * 00379 DO 110 L = KA - 1, KB - K + 1, -1 00380 NRT = ( N-J2+L ) / KA1 00381 IF( NRT.GT.0 ) 00382 $ CALL ZLARTV( NRT, AB( L, J2+KA1-L ), INCA, 00383 $ AB( L+1, J2+KA1-L ), INCA, RWORK( J2-M ), 00384 $ WORK( J2-M ), KA1 ) 00385 110 CONTINUE 00386 * 00387 IF( WANTX ) THEN 00388 * 00389 * post-multiply X by product of rotations in 1st set 00390 * 00391 DO 120 J = J2, J1, KA1 00392 CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1, 00393 $ RWORK( J-M ), DCONJG( WORK( J-M ) ) ) 00394 120 CONTINUE 00395 END IF 00396 130 CONTINUE 00397 * 00398 IF( UPDATE ) THEN 00399 IF( I2.LE.N .AND. KBT.GT.0 ) THEN 00400 * 00401 * create nonzero element a(i-kbt,i-kbt+ka+1) outside the 00402 * band and store it in WORK(i-kbt) 00403 * 00404 WORK( I-KBT ) = -BB( KB1-KBT, I )*RA1 00405 END IF 00406 END IF 00407 * 00408 DO 170 K = KB, 1, -1 00409 IF( UPDATE ) THEN 00410 J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1 00411 ELSE 00412 J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1 00413 END IF 00414 * 00415 * finish applying rotations in 2nd set from the left 00416 * 00417 DO 140 L = KB - K, 1, -1 00418 NRT = ( N-J2+KA+L ) / KA1 00419 IF( NRT.GT.0 ) 00420 $ CALL ZLARTV( NRT, AB( L, J2-L+1 ), INCA, 00421 $ AB( L+1, J2-L+1 ), INCA, RWORK( J2-KA ), 00422 $ WORK( J2-KA ), KA1 ) 00423 140 CONTINUE 00424 NR = ( N-J2+KA ) / KA1 00425 J1 = J2 + ( NR-1 )*KA1 00426 DO 150 J = J1, J2, -KA1 00427 WORK( J ) = WORK( J-KA ) 00428 RWORK( J ) = RWORK( J-KA ) 00429 150 CONTINUE 00430 DO 160 J = J2, J1, KA1 00431 * 00432 * create nonzero element a(j-ka,j+1) outside the band 00433 * and store it in WORK(j) 00434 * 00435 WORK( J ) = WORK( J )*AB( 1, J+1 ) 00436 AB( 1, J+1 ) = RWORK( J )*AB( 1, J+1 ) 00437 160 CONTINUE 00438 IF( UPDATE ) THEN 00439 IF( I-K.LT.N-KA .AND. K.LE.KBT ) 00440 $ WORK( I-K+KA ) = WORK( I-K ) 00441 END IF 00442 170 CONTINUE 00443 * 00444 DO 210 K = KB, 1, -1 00445 J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1 00446 NR = ( N-J2+KA ) / KA1 00447 J1 = J2 + ( NR-1 )*KA1 00448 IF( NR.GT.0 ) THEN 00449 * 00450 * generate rotations in 2nd set to annihilate elements 00451 * which have been created outside the band 00452 * 00453 CALL ZLARGV( NR, AB( 1, J2 ), INCA, WORK( J2 ), KA1, 00454 $ RWORK( J2 ), KA1 ) 00455 * 00456 * apply rotations in 2nd set from the right 00457 * 00458 DO 180 L = 1, KA - 1 00459 CALL ZLARTV( NR, AB( KA1-L, J2 ), INCA, 00460 $ AB( KA-L, J2+1 ), INCA, RWORK( J2 ), 00461 $ WORK( J2 ), KA1 ) 00462 180 CONTINUE 00463 * 00464 * apply rotations in 2nd set from both sides to diagonal 00465 * blocks 00466 * 00467 CALL ZLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ), 00468 $ AB( KA, J2+1 ), INCA, RWORK( J2 ), 00469 $ WORK( J2 ), KA1 ) 00470 * 00471 CALL ZLACGV( NR, WORK( J2 ), KA1 ) 00472 END IF 00473 * 00474 * start applying rotations in 2nd set from the left 00475 * 00476 DO 190 L = KA - 1, KB - K + 1, -1 00477 NRT = ( N-J2+L ) / KA1 00478 IF( NRT.GT.0 ) 00479 $ CALL ZLARTV( NRT, AB( L, J2+KA1-L ), INCA, 00480 $ AB( L+1, J2+KA1-L ), INCA, RWORK( J2 ), 00481 $ WORK( J2 ), KA1 ) 00482 190 CONTINUE 00483 * 00484 IF( WANTX ) THEN 00485 * 00486 * post-multiply X by product of rotations in 2nd set 00487 * 00488 DO 200 J = J2, J1, KA1 00489 CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1, 00490 $ RWORK( J ), DCONJG( WORK( J ) ) ) 00491 200 CONTINUE 00492 END IF 00493 210 CONTINUE 00494 * 00495 DO 230 K = 1, KB - 1 00496 J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1 00497 * 00498 * finish applying rotations in 1st set from the left 00499 * 00500 DO 220 L = KB - K, 1, -1 00501 NRT = ( N-J2+L ) / KA1 00502 IF( NRT.GT.0 ) 00503 $ CALL ZLARTV( NRT, AB( L, J2+KA1-L ), INCA, 00504 $ AB( L+1, J2+KA1-L ), INCA, RWORK( J2-M ), 00505 $ WORK( J2-M ), KA1 ) 00506 220 CONTINUE 00507 230 CONTINUE 00508 * 00509 IF( KB.GT.1 ) THEN 00510 DO 240 J = N - 1, J2 + KA, -1 00511 RWORK( J-M ) = RWORK( J-KA-M ) 00512 WORK( J-M ) = WORK( J-KA-M ) 00513 240 CONTINUE 00514 END IF 00515 * 00516 ELSE 00517 * 00518 * Transform A, working with the lower triangle 00519 * 00520 IF( UPDATE ) THEN 00521 * 00522 * Form inv(S(i))**H * A * inv(S(i)) 00523 * 00524 BII = DBLE( BB( 1, I ) ) 00525 AB( 1, I ) = ( DBLE( AB( 1, I ) ) / BII ) / BII 00526 DO 250 J = I + 1, I1 00527 AB( J-I+1, I ) = AB( J-I+1, I ) / BII 00528 250 CONTINUE 00529 DO 260 J = MAX( 1, I-KA ), I - 1 00530 AB( I-J+1, J ) = AB( I-J+1, J ) / BII 00531 260 CONTINUE 00532 DO 290 K = I - KBT, I - 1 00533 DO 270 J = I - KBT, K 00534 AB( K-J+1, J ) = AB( K-J+1, J ) - 00535 $ BB( I-J+1, J )*DCONJG( AB( I-K+1, 00536 $ K ) ) - DCONJG( BB( I-K+1, K ) )* 00537 $ AB( I-J+1, J ) + DBLE( AB( 1, I ) )* 00538 $ BB( I-J+1, J )*DCONJG( BB( I-K+1, 00539 $ K ) ) 00540 270 CONTINUE 00541 DO 280 J = MAX( 1, I-KA ), I - KBT - 1 00542 AB( K-J+1, J ) = AB( K-J+1, J ) - 00543 $ DCONJG( BB( I-K+1, K ) )* 00544 $ AB( I-J+1, J ) 00545 280 CONTINUE 00546 290 CONTINUE 00547 DO 310 J = I, I1 00548 DO 300 K = MAX( J-KA, I-KBT ), I - 1 00549 AB( J-K+1, K ) = AB( J-K+1, K ) - 00550 $ BB( I-K+1, K )*AB( J-I+1, I ) 00551 300 CONTINUE 00552 310 CONTINUE 00553 * 00554 IF( WANTX ) THEN 00555 * 00556 * post-multiply X by inv(S(i)) 00557 * 00558 CALL ZDSCAL( N-M, ONE / BII, X( M+1, I ), 1 ) 00559 IF( KBT.GT.0 ) 00560 $ CALL ZGERU( N-M, KBT, -CONE, X( M+1, I ), 1, 00561 $ BB( KBT+1, I-KBT ), LDBB-1, 00562 $ X( M+1, I-KBT ), LDX ) 00563 END IF 00564 * 00565 * store a(i1,i) in RA1 for use in next loop over K 00566 * 00567 RA1 = AB( I1-I+1, I ) 00568 END IF 00569 * 00570 * Generate and apply vectors of rotations to chase all the 00571 * existing bulges KA positions down toward the bottom of the 00572 * band 00573 * 00574 DO 360 K = 1, KB - 1 00575 IF( UPDATE ) THEN 00576 * 00577 * Determine the rotations which would annihilate the bulge 00578 * which has in theory just been created 00579 * 00580 IF( I-K+KA.LT.N .AND. I-K.GT.1 ) THEN 00581 * 00582 * generate rotation to annihilate a(i-k+ka+1,i) 00583 * 00584 CALL ZLARTG( AB( KA1-K, I ), RA1, RWORK( I-K+KA-M ), 00585 $ WORK( I-K+KA-M ), RA ) 00586 * 00587 * create nonzero element a(i-k+ka+1,i-k) outside the 00588 * band and store it in WORK(i-k) 00589 * 00590 T = -BB( K+1, I-K )*RA1 00591 WORK( I-K ) = RWORK( I-K+KA-M )*T - 00592 $ DCONJG( WORK( I-K+KA-M ) )* 00593 $ AB( KA1, I-K ) 00594 AB( KA1, I-K ) = WORK( I-K+KA-M )*T + 00595 $ RWORK( I-K+KA-M )*AB( KA1, I-K ) 00596 RA1 = RA 00597 END IF 00598 END IF 00599 J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1 00600 NR = ( N-J2+KA ) / KA1 00601 J1 = J2 + ( NR-1 )*KA1 00602 IF( UPDATE ) THEN 00603 J2T = MAX( J2, I+2*KA-K+1 ) 00604 ELSE 00605 J2T = J2 00606 END IF 00607 NRT = ( N-J2T+KA ) / KA1 00608 DO 320 J = J2T, J1, KA1 00609 * 00610 * create nonzero element a(j+1,j-ka) outside the band 00611 * and store it in WORK(j-m) 00612 * 00613 WORK( J-M ) = WORK( J-M )*AB( KA1, J-KA+1 ) 00614 AB( KA1, J-KA+1 ) = RWORK( J-M )*AB( KA1, J-KA+1 ) 00615 320 CONTINUE 00616 * 00617 * generate rotations in 1st set to annihilate elements which 00618 * have been created outside the band 00619 * 00620 IF( NRT.GT.0 ) 00621 $ CALL ZLARGV( NRT, AB( KA1, J2T-KA ), INCA, WORK( J2T-M ), 00622 $ KA1, RWORK( J2T-M ), KA1 ) 00623 IF( NR.GT.0 ) THEN 00624 * 00625 * apply rotations in 1st set from the left 00626 * 00627 DO 330 L = 1, KA - 1 00628 CALL ZLARTV( NR, AB( L+1, J2-L ), INCA, 00629 $ AB( L+2, J2-L ), INCA, RWORK( J2-M ), 00630 $ WORK( J2-M ), KA1 ) 00631 330 CONTINUE 00632 * 00633 * apply rotations in 1st set from both sides to diagonal 00634 * blocks 00635 * 00636 CALL ZLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ), 00637 $ INCA, RWORK( J2-M ), WORK( J2-M ), KA1 ) 00638 * 00639 CALL ZLACGV( NR, WORK( J2-M ), KA1 ) 00640 END IF 00641 * 00642 * start applying rotations in 1st set from the right 00643 * 00644 DO 340 L = KA - 1, KB - K + 1, -1 00645 NRT = ( N-J2+L ) / KA1 00646 IF( NRT.GT.0 ) 00647 $ CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA, 00648 $ AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ), 00649 $ WORK( J2-M ), KA1 ) 00650 340 CONTINUE 00651 * 00652 IF( WANTX ) THEN 00653 * 00654 * post-multiply X by product of rotations in 1st set 00655 * 00656 DO 350 J = J2, J1, KA1 00657 CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1, 00658 $ RWORK( J-M ), WORK( J-M ) ) 00659 350 CONTINUE 00660 END IF 00661 360 CONTINUE 00662 * 00663 IF( UPDATE ) THEN 00664 IF( I2.LE.N .AND. KBT.GT.0 ) THEN 00665 * 00666 * create nonzero element a(i-kbt+ka+1,i-kbt) outside the 00667 * band and store it in WORK(i-kbt) 00668 * 00669 WORK( I-KBT ) = -BB( KBT+1, I-KBT )*RA1 00670 END IF 00671 END IF 00672 * 00673 DO 400 K = KB, 1, -1 00674 IF( UPDATE ) THEN 00675 J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1 00676 ELSE 00677 J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1 00678 END IF 00679 * 00680 * finish applying rotations in 2nd set from the right 00681 * 00682 DO 370 L = KB - K, 1, -1 00683 NRT = ( N-J2+KA+L ) / KA1 00684 IF( NRT.GT.0 ) 00685 $ CALL ZLARTV( NRT, AB( KA1-L+1, J2-KA ), INCA, 00686 $ AB( KA1-L, J2-KA+1 ), INCA, 00687 $ RWORK( J2-KA ), WORK( J2-KA ), KA1 ) 00688 370 CONTINUE 00689 NR = ( N-J2+KA ) / KA1 00690 J1 = J2 + ( NR-1 )*KA1 00691 DO 380 J = J1, J2, -KA1 00692 WORK( J ) = WORK( J-KA ) 00693 RWORK( J ) = RWORK( J-KA ) 00694 380 CONTINUE 00695 DO 390 J = J2, J1, KA1 00696 * 00697 * create nonzero element a(j+1,j-ka) outside the band 00698 * and store it in WORK(j) 00699 * 00700 WORK( J ) = WORK( J )*AB( KA1, J-KA+1 ) 00701 AB( KA1, J-KA+1 ) = RWORK( J )*AB( KA1, J-KA+1 ) 00702 390 CONTINUE 00703 IF( UPDATE ) THEN 00704 IF( I-K.LT.N-KA .AND. K.LE.KBT ) 00705 $ WORK( I-K+KA ) = WORK( I-K ) 00706 END IF 00707 400 CONTINUE 00708 * 00709 DO 440 K = KB, 1, -1 00710 J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1 00711 NR = ( N-J2+KA ) / KA1 00712 J1 = J2 + ( NR-1 )*KA1 00713 IF( NR.GT.0 ) THEN 00714 * 00715 * generate rotations in 2nd set to annihilate elements 00716 * which have been created outside the band 00717 * 00718 CALL ZLARGV( NR, AB( KA1, J2-KA ), INCA, WORK( J2 ), KA1, 00719 $ RWORK( J2 ), KA1 ) 00720 * 00721 * apply rotations in 2nd set from the left 00722 * 00723 DO 410 L = 1, KA - 1 00724 CALL ZLARTV( NR, AB( L+1, J2-L ), INCA, 00725 $ AB( L+2, J2-L ), INCA, RWORK( J2 ), 00726 $ WORK( J2 ), KA1 ) 00727 410 CONTINUE 00728 * 00729 * apply rotations in 2nd set from both sides to diagonal 00730 * blocks 00731 * 00732 CALL ZLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ), 00733 $ INCA, RWORK( J2 ), WORK( J2 ), KA1 ) 00734 * 00735 CALL ZLACGV( NR, WORK( J2 ), KA1 ) 00736 END IF 00737 * 00738 * start applying rotations in 2nd set from the right 00739 * 00740 DO 420 L = KA - 1, KB - K + 1, -1 00741 NRT = ( N-J2+L ) / KA1 00742 IF( NRT.GT.0 ) 00743 $ CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA, 00744 $ AB( KA1-L, J2+1 ), INCA, RWORK( J2 ), 00745 $ WORK( J2 ), KA1 ) 00746 420 CONTINUE 00747 * 00748 IF( WANTX ) THEN 00749 * 00750 * post-multiply X by product of rotations in 2nd set 00751 * 00752 DO 430 J = J2, J1, KA1 00753 CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1, 00754 $ RWORK( J ), WORK( J ) ) 00755 430 CONTINUE 00756 END IF 00757 440 CONTINUE 00758 * 00759 DO 460 K = 1, KB - 1 00760 J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1 00761 * 00762 * finish applying rotations in 1st set from the right 00763 * 00764 DO 450 L = KB - K, 1, -1 00765 NRT = ( N-J2+L ) / KA1 00766 IF( NRT.GT.0 ) 00767 $ CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA, 00768 $ AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ), 00769 $ WORK( J2-M ), KA1 ) 00770 450 CONTINUE 00771 460 CONTINUE 00772 * 00773 IF( KB.GT.1 ) THEN 00774 DO 470 J = N - 1, J2 + KA, -1 00775 RWORK( J-M ) = RWORK( J-KA-M ) 00776 WORK( J-M ) = WORK( J-KA-M ) 00777 470 CONTINUE 00778 END IF 00779 * 00780 END IF 00781 * 00782 GO TO 10 00783 * 00784 480 CONTINUE 00785 * 00786 * **************************** Phase 2 ***************************** 00787 * 00788 * The logical structure of this phase is: 00789 * 00790 * UPDATE = .TRUE. 00791 * DO I = 1, M 00792 * use S(i) to update A and create a new bulge 00793 * apply rotations to push all bulges KA positions upward 00794 * END DO 00795 * UPDATE = .FALSE. 00796 * DO I = M - KA - 1, 2, -1 00797 * apply rotations to push all bulges KA positions upward 00798 * END DO 00799 * 00800 * To avoid duplicating code, the two loops are merged. 00801 * 00802 UPDATE = .TRUE. 00803 I = 0 00804 490 CONTINUE 00805 IF( UPDATE ) THEN 00806 I = I + 1 00807 KBT = MIN( KB, M-I ) 00808 I0 = I + 1 00809 I1 = MAX( 1, I-KA ) 00810 I2 = I + KBT - KA1 00811 IF( I.GT.M ) THEN 00812 UPDATE = .FALSE. 00813 I = I - 1 00814 I0 = M + 1 00815 IF( KA.EQ.0 ) 00816 $ RETURN 00817 GO TO 490 00818 END IF 00819 ELSE 00820 I = I - KA 00821 IF( I.LT.2 ) 00822 $ RETURN 00823 END IF 00824 * 00825 IF( I.LT.M-KBT ) THEN 00826 NX = M 00827 ELSE 00828 NX = N 00829 END IF 00830 * 00831 IF( UPPER ) THEN 00832 * 00833 * Transform A, working with the upper triangle 00834 * 00835 IF( UPDATE ) THEN 00836 * 00837 * Form inv(S(i))**H * A * inv(S(i)) 00838 * 00839 BII = DBLE( BB( KB1, I ) ) 00840 AB( KA1, I ) = ( DBLE( AB( KA1, I ) ) / BII ) / BII 00841 DO 500 J = I1, I - 1 00842 AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII 00843 500 CONTINUE 00844 DO 510 J = I + 1, MIN( N, I+KA ) 00845 AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII 00846 510 CONTINUE 00847 DO 540 K = I + 1, I + KBT 00848 DO 520 J = K, I + KBT 00849 AB( K-J+KA1, J ) = AB( K-J+KA1, J ) - 00850 $ BB( I-J+KB1, J )* 00851 $ DCONJG( AB( I-K+KA1, K ) ) - 00852 $ DCONJG( BB( I-K+KB1, K ) )* 00853 $ AB( I-J+KA1, J ) + 00854 $ DBLE( AB( KA1, I ) )* 00855 $ BB( I-J+KB1, J )* 00856 $ DCONJG( BB( I-K+KB1, K ) ) 00857 520 CONTINUE 00858 DO 530 J = I + KBT + 1, MIN( N, I+KA ) 00859 AB( K-J+KA1, J ) = AB( K-J+KA1, J ) - 00860 $ DCONJG( BB( I-K+KB1, K ) )* 00861 $ AB( I-J+KA1, J ) 00862 530 CONTINUE 00863 540 CONTINUE 00864 DO 560 J = I1, I 00865 DO 550 K = I + 1, MIN( J+KA, I+KBT ) 00866 AB( J-K+KA1, K ) = AB( J-K+KA1, K ) - 00867 $ BB( I-K+KB1, K )*AB( J-I+KA1, I ) 00868 550 CONTINUE 00869 560 CONTINUE 00870 * 00871 IF( WANTX ) THEN 00872 * 00873 * post-multiply X by inv(S(i)) 00874 * 00875 CALL ZDSCAL( NX, ONE / BII, X( 1, I ), 1 ) 00876 IF( KBT.GT.0 ) 00877 $ CALL ZGERU( NX, KBT, -CONE, X( 1, I ), 1, 00878 $ BB( KB, I+1 ), LDBB-1, X( 1, I+1 ), LDX ) 00879 END IF 00880 * 00881 * store a(i1,i) in RA1 for use in next loop over K 00882 * 00883 RA1 = AB( I1-I+KA1, I ) 00884 END IF 00885 * 00886 * Generate and apply vectors of rotations to chase all the 00887 * existing bulges KA positions up toward the top of the band 00888 * 00889 DO 610 K = 1, KB - 1 00890 IF( UPDATE ) THEN 00891 * 00892 * Determine the rotations which would annihilate the bulge 00893 * which has in theory just been created 00894 * 00895 IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN 00896 * 00897 * generate rotation to annihilate a(i+k-ka-1,i) 00898 * 00899 CALL ZLARTG( AB( K+1, I ), RA1, RWORK( I+K-KA ), 00900 $ WORK( I+K-KA ), RA ) 00901 * 00902 * create nonzero element a(i+k-ka-1,i+k) outside the 00903 * band and store it in WORK(m-kb+i+k) 00904 * 00905 T = -BB( KB1-K, I+K )*RA1 00906 WORK( M-KB+I+K ) = RWORK( I+K-KA )*T - 00907 $ DCONJG( WORK( I+K-KA ) )* 00908 $ AB( 1, I+K ) 00909 AB( 1, I+K ) = WORK( I+K-KA )*T + 00910 $ RWORK( I+K-KA )*AB( 1, I+K ) 00911 RA1 = RA 00912 END IF 00913 END IF 00914 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1 00915 NR = ( J2+KA-1 ) / KA1 00916 J1 = J2 - ( NR-1 )*KA1 00917 IF( UPDATE ) THEN 00918 J2T = MIN( J2, I-2*KA+K-1 ) 00919 ELSE 00920 J2T = J2 00921 END IF 00922 NRT = ( J2T+KA-1 ) / KA1 00923 DO 570 J = J1, J2T, KA1 00924 * 00925 * create nonzero element a(j-1,j+ka) outside the band 00926 * and store it in WORK(j) 00927 * 00928 WORK( J ) = WORK( J )*AB( 1, J+KA-1 ) 00929 AB( 1, J+KA-1 ) = RWORK( J )*AB( 1, J+KA-1 ) 00930 570 CONTINUE 00931 * 00932 * generate rotations in 1st set to annihilate elements which 00933 * have been created outside the band 00934 * 00935 IF( NRT.GT.0 ) 00936 $ CALL ZLARGV( NRT, AB( 1, J1+KA ), INCA, WORK( J1 ), KA1, 00937 $ RWORK( J1 ), KA1 ) 00938 IF( NR.GT.0 ) THEN 00939 * 00940 * apply rotations in 1st set from the left 00941 * 00942 DO 580 L = 1, KA - 1 00943 CALL ZLARTV( NR, AB( KA1-L, J1+L ), INCA, 00944 $ AB( KA-L, J1+L ), INCA, RWORK( J1 ), 00945 $ WORK( J1 ), KA1 ) 00946 580 CONTINUE 00947 * 00948 * apply rotations in 1st set from both sides to diagonal 00949 * blocks 00950 * 00951 CALL ZLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ), 00952 $ AB( KA, J1 ), INCA, RWORK( J1 ), WORK( J1 ), 00953 $ KA1 ) 00954 * 00955 CALL ZLACGV( NR, WORK( J1 ), KA1 ) 00956 END IF 00957 * 00958 * start applying rotations in 1st set from the right 00959 * 00960 DO 590 L = KA - 1, KB - K + 1, -1 00961 NRT = ( J2+L-1 ) / KA1 00962 J1T = J2 - ( NRT-1 )*KA1 00963 IF( NRT.GT.0 ) 00964 $ CALL ZLARTV( NRT, AB( L, J1T ), INCA, 00965 $ AB( L+1, J1T-1 ), INCA, RWORK( J1T ), 00966 $ WORK( J1T ), KA1 ) 00967 590 CONTINUE 00968 * 00969 IF( WANTX ) THEN 00970 * 00971 * post-multiply X by product of rotations in 1st set 00972 * 00973 DO 600 J = J1, J2, KA1 00974 CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1, 00975 $ RWORK( J ), WORK( J ) ) 00976 600 CONTINUE 00977 END IF 00978 610 CONTINUE 00979 * 00980 IF( UPDATE ) THEN 00981 IF( I2.GT.0 .AND. KBT.GT.0 ) THEN 00982 * 00983 * create nonzero element a(i+kbt-ka-1,i+kbt) outside the 00984 * band and store it in WORK(m-kb+i+kbt) 00985 * 00986 WORK( M-KB+I+KBT ) = -BB( KB1-KBT, I+KBT )*RA1 00987 END IF 00988 END IF 00989 * 00990 DO 650 K = KB, 1, -1 00991 IF( UPDATE ) THEN 00992 J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1 00993 ELSE 00994 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1 00995 END IF 00996 * 00997 * finish applying rotations in 2nd set from the right 00998 * 00999 DO 620 L = KB - K, 1, -1 01000 NRT = ( J2+KA+L-1 ) / KA1 01001 J1T = J2 - ( NRT-1 )*KA1 01002 IF( NRT.GT.0 ) 01003 $ CALL ZLARTV( NRT, AB( L, J1T+KA ), INCA, 01004 $ AB( L+1, J1T+KA-1 ), INCA, 01005 $ RWORK( M-KB+J1T+KA ), 01006 $ WORK( M-KB+J1T+KA ), KA1 ) 01007 620 CONTINUE 01008 NR = ( J2+KA-1 ) / KA1 01009 J1 = J2 - ( NR-1 )*KA1 01010 DO 630 J = J1, J2, KA1 01011 WORK( M-KB+J ) = WORK( M-KB+J+KA ) 01012 RWORK( M-KB+J ) = RWORK( M-KB+J+KA ) 01013 630 CONTINUE 01014 DO 640 J = J1, J2, KA1 01015 * 01016 * create nonzero element a(j-1,j+ka) outside the band 01017 * and store it in WORK(m-kb+j) 01018 * 01019 WORK( M-KB+J ) = WORK( M-KB+J )*AB( 1, J+KA-1 ) 01020 AB( 1, J+KA-1 ) = RWORK( M-KB+J )*AB( 1, J+KA-1 ) 01021 640 CONTINUE 01022 IF( UPDATE ) THEN 01023 IF( I+K.GT.KA1 .AND. K.LE.KBT ) 01024 $ WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K ) 01025 END IF 01026 650 CONTINUE 01027 * 01028 DO 690 K = KB, 1, -1 01029 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1 01030 NR = ( J2+KA-1 ) / KA1 01031 J1 = J2 - ( NR-1 )*KA1 01032 IF( NR.GT.0 ) THEN 01033 * 01034 * generate rotations in 2nd set to annihilate elements 01035 * which have been created outside the band 01036 * 01037 CALL ZLARGV( NR, AB( 1, J1+KA ), INCA, WORK( M-KB+J1 ), 01038 $ KA1, RWORK( M-KB+J1 ), KA1 ) 01039 * 01040 * apply rotations in 2nd set from the left 01041 * 01042 DO 660 L = 1, KA - 1 01043 CALL ZLARTV( NR, AB( KA1-L, J1+L ), INCA, 01044 $ AB( KA-L, J1+L ), INCA, RWORK( M-KB+J1 ), 01045 $ WORK( M-KB+J1 ), KA1 ) 01046 660 CONTINUE 01047 * 01048 * apply rotations in 2nd set from both sides to diagonal 01049 * blocks 01050 * 01051 CALL ZLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ), 01052 $ AB( KA, J1 ), INCA, RWORK( M-KB+J1 ), 01053 $ WORK( M-KB+J1 ), KA1 ) 01054 * 01055 CALL ZLACGV( NR, WORK( M-KB+J1 ), KA1 ) 01056 END IF 01057 * 01058 * start applying rotations in 2nd set from the right 01059 * 01060 DO 670 L = KA - 1, KB - K + 1, -1 01061 NRT = ( J2+L-1 ) / KA1 01062 J1T = J2 - ( NRT-1 )*KA1 01063 IF( NRT.GT.0 ) 01064 $ CALL ZLARTV( NRT, AB( L, J1T ), INCA, 01065 $ AB( L+1, J1T-1 ), INCA, 01066 $ RWORK( M-KB+J1T ), WORK( M-KB+J1T ), 01067 $ KA1 ) 01068 670 CONTINUE 01069 * 01070 IF( WANTX ) THEN 01071 * 01072 * post-multiply X by product of rotations in 2nd set 01073 * 01074 DO 680 J = J1, J2, KA1 01075 CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1, 01076 $ RWORK( M-KB+J ), WORK( M-KB+J ) ) 01077 680 CONTINUE 01078 END IF 01079 690 CONTINUE 01080 * 01081 DO 710 K = 1, KB - 1 01082 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1 01083 * 01084 * finish applying rotations in 1st set from the right 01085 * 01086 DO 700 L = KB - K, 1, -1 01087 NRT = ( J2+L-1 ) / KA1 01088 J1T = J2 - ( NRT-1 )*KA1 01089 IF( NRT.GT.0 ) 01090 $ CALL ZLARTV( NRT, AB( L, J1T ), INCA, 01091 $ AB( L+1, J1T-1 ), INCA, RWORK( J1T ), 01092 $ WORK( J1T ), KA1 ) 01093 700 CONTINUE 01094 710 CONTINUE 01095 * 01096 IF( KB.GT.1 ) THEN 01097 DO 720 J = 2, I2 - KA 01098 RWORK( J ) = RWORK( J+KA ) 01099 WORK( J ) = WORK( J+KA ) 01100 720 CONTINUE 01101 END IF 01102 * 01103 ELSE 01104 * 01105 * Transform A, working with the lower triangle 01106 * 01107 IF( UPDATE ) THEN 01108 * 01109 * Form inv(S(i))**H * A * inv(S(i)) 01110 * 01111 BII = DBLE( BB( 1, I ) ) 01112 AB( 1, I ) = ( DBLE( AB( 1, I ) ) / BII ) / BII 01113 DO 730 J = I1, I - 1 01114 AB( I-J+1, J ) = AB( I-J+1, J ) / BII 01115 730 CONTINUE 01116 DO 740 J = I + 1, MIN( N, I+KA ) 01117 AB( J-I+1, I ) = AB( J-I+1, I ) / BII 01118 740 CONTINUE 01119 DO 770 K = I + 1, I + KBT 01120 DO 750 J = K, I + KBT 01121 AB( J-K+1, K ) = AB( J-K+1, K ) - 01122 $ BB( J-I+1, I )*DCONJG( AB( K-I+1, 01123 $ I ) ) - DCONJG( BB( K-I+1, I ) )* 01124 $ AB( J-I+1, I ) + DBLE( AB( 1, I ) )* 01125 $ BB( J-I+1, I )*DCONJG( BB( K-I+1, 01126 $ I ) ) 01127 750 CONTINUE 01128 DO 760 J = I + KBT + 1, MIN( N, I+KA ) 01129 AB( J-K+1, K ) = AB( J-K+1, K ) - 01130 $ DCONJG( BB( K-I+1, I ) )* 01131 $ AB( J-I+1, I ) 01132 760 CONTINUE 01133 770 CONTINUE 01134 DO 790 J = I1, I 01135 DO 780 K = I + 1, MIN( J+KA, I+KBT ) 01136 AB( K-J+1, J ) = AB( K-J+1, J ) - 01137 $ BB( K-I+1, I )*AB( I-J+1, J ) 01138 780 CONTINUE 01139 790 CONTINUE 01140 * 01141 IF( WANTX ) THEN 01142 * 01143 * post-multiply X by inv(S(i)) 01144 * 01145 CALL ZDSCAL( NX, ONE / BII, X( 1, I ), 1 ) 01146 IF( KBT.GT.0 ) 01147 $ CALL ZGERC( NX, KBT, -CONE, X( 1, I ), 1, BB( 2, I ), 01148 $ 1, X( 1, I+1 ), LDX ) 01149 END IF 01150 * 01151 * store a(i,i1) in RA1 for use in next loop over K 01152 * 01153 RA1 = AB( I-I1+1, I1 ) 01154 END IF 01155 * 01156 * Generate and apply vectors of rotations to chase all the 01157 * existing bulges KA positions up toward the top of the band 01158 * 01159 DO 840 K = 1, KB - 1 01160 IF( UPDATE ) THEN 01161 * 01162 * Determine the rotations which would annihilate the bulge 01163 * which has in theory just been created 01164 * 01165 IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN 01166 * 01167 * generate rotation to annihilate a(i,i+k-ka-1) 01168 * 01169 CALL ZLARTG( AB( KA1-K, I+K-KA ), RA1, 01170 $ RWORK( I+K-KA ), WORK( I+K-KA ), RA ) 01171 * 01172 * create nonzero element a(i+k,i+k-ka-1) outside the 01173 * band and store it in WORK(m-kb+i+k) 01174 * 01175 T = -BB( K+1, I )*RA1 01176 WORK( M-KB+I+K ) = RWORK( I+K-KA )*T - 01177 $ DCONJG( WORK( I+K-KA ) )* 01178 $ AB( KA1, I+K-KA ) 01179 AB( KA1, I+K-KA ) = WORK( I+K-KA )*T + 01180 $ RWORK( I+K-KA )*AB( KA1, I+K-KA ) 01181 RA1 = RA 01182 END IF 01183 END IF 01184 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1 01185 NR = ( J2+KA-1 ) / KA1 01186 J1 = J2 - ( NR-1 )*KA1 01187 IF( UPDATE ) THEN 01188 J2T = MIN( J2, I-2*KA+K-1 ) 01189 ELSE 01190 J2T = J2 01191 END IF 01192 NRT = ( J2T+KA-1 ) / KA1 01193 DO 800 J = J1, J2T, KA1 01194 * 01195 * create nonzero element a(j+ka,j-1) outside the band 01196 * and store it in WORK(j) 01197 * 01198 WORK( J ) = WORK( J )*AB( KA1, J-1 ) 01199 AB( KA1, J-1 ) = RWORK( J )*AB( KA1, J-1 ) 01200 800 CONTINUE 01201 * 01202 * generate rotations in 1st set to annihilate elements which 01203 * have been created outside the band 01204 * 01205 IF( NRT.GT.0 ) 01206 $ CALL ZLARGV( NRT, AB( KA1, J1 ), INCA, WORK( J1 ), KA1, 01207 $ RWORK( J1 ), KA1 ) 01208 IF( NR.GT.0 ) THEN 01209 * 01210 * apply rotations in 1st set from the right 01211 * 01212 DO 810 L = 1, KA - 1 01213 CALL ZLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ), 01214 $ INCA, RWORK( J1 ), WORK( J1 ), KA1 ) 01215 810 CONTINUE 01216 * 01217 * apply rotations in 1st set from both sides to diagonal 01218 * blocks 01219 * 01220 CALL ZLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ), 01221 $ AB( 2, J1-1 ), INCA, RWORK( J1 ), 01222 $ WORK( J1 ), KA1 ) 01223 * 01224 CALL ZLACGV( NR, WORK( J1 ), KA1 ) 01225 END IF 01226 * 01227 * start applying rotations in 1st set from the left 01228 * 01229 DO 820 L = KA - 1, KB - K + 1, -1 01230 NRT = ( J2+L-1 ) / KA1 01231 J1T = J2 - ( NRT-1 )*KA1 01232 IF( NRT.GT.0 ) 01233 $ CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA, 01234 $ AB( KA1-L, J1T-KA1+L ), INCA, 01235 $ RWORK( J1T ), WORK( J1T ), KA1 ) 01236 820 CONTINUE 01237 * 01238 IF( WANTX ) THEN 01239 * 01240 * post-multiply X by product of rotations in 1st set 01241 * 01242 DO 830 J = J1, J2, KA1 01243 CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1, 01244 $ RWORK( J ), DCONJG( WORK( J ) ) ) 01245 830 CONTINUE 01246 END IF 01247 840 CONTINUE 01248 * 01249 IF( UPDATE ) THEN 01250 IF( I2.GT.0 .AND. KBT.GT.0 ) THEN 01251 * 01252 * create nonzero element a(i+kbt,i+kbt-ka-1) outside the 01253 * band and store it in WORK(m-kb+i+kbt) 01254 * 01255 WORK( M-KB+I+KBT ) = -BB( KBT+1, I )*RA1 01256 END IF 01257 END IF 01258 * 01259 DO 880 K = KB, 1, -1 01260 IF( UPDATE ) THEN 01261 J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1 01262 ELSE 01263 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1 01264 END IF 01265 * 01266 * finish applying rotations in 2nd set from the left 01267 * 01268 DO 850 L = KB - K, 1, -1 01269 NRT = ( J2+KA+L-1 ) / KA1 01270 J1T = J2 - ( NRT-1 )*KA1 01271 IF( NRT.GT.0 ) 01272 $ CALL ZLARTV( NRT, AB( KA1-L+1, J1T+L-1 ), INCA, 01273 $ AB( KA1-L, J1T+L-1 ), INCA, 01274 $ RWORK( M-KB+J1T+KA ), 01275 $ WORK( M-KB+J1T+KA ), KA1 ) 01276 850 CONTINUE 01277 NR = ( J2+KA-1 ) / KA1 01278 J1 = J2 - ( NR-1 )*KA1 01279 DO 860 J = J1, J2, KA1 01280 WORK( M-KB+J ) = WORK( M-KB+J+KA ) 01281 RWORK( M-KB+J ) = RWORK( M-KB+J+KA ) 01282 860 CONTINUE 01283 DO 870 J = J1, J2, KA1 01284 * 01285 * create nonzero element a(j+ka,j-1) outside the band 01286 * and store it in WORK(m-kb+j) 01287 * 01288 WORK( M-KB+J ) = WORK( M-KB+J )*AB( KA1, J-1 ) 01289 AB( KA1, J-1 ) = RWORK( M-KB+J )*AB( KA1, J-1 ) 01290 870 CONTINUE 01291 IF( UPDATE ) THEN 01292 IF( I+K.GT.KA1 .AND. K.LE.KBT ) 01293 $ WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K ) 01294 END IF 01295 880 CONTINUE 01296 * 01297 DO 920 K = KB, 1, -1 01298 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1 01299 NR = ( J2+KA-1 ) / KA1 01300 J1 = J2 - ( NR-1 )*KA1 01301 IF( NR.GT.0 ) THEN 01302 * 01303 * generate rotations in 2nd set to annihilate elements 01304 * which have been created outside the band 01305 * 01306 CALL ZLARGV( NR, AB( KA1, J1 ), INCA, WORK( M-KB+J1 ), 01307 $ KA1, RWORK( M-KB+J1 ), KA1 ) 01308 * 01309 * apply rotations in 2nd set from the right 01310 * 01311 DO 890 L = 1, KA - 1 01312 CALL ZLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ), 01313 $ INCA, RWORK( M-KB+J1 ), WORK( M-KB+J1 ), 01314 $ KA1 ) 01315 890 CONTINUE 01316 * 01317 * apply rotations in 2nd set from both sides to diagonal 01318 * blocks 01319 * 01320 CALL ZLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ), 01321 $ AB( 2, J1-1 ), INCA, RWORK( M-KB+J1 ), 01322 $ WORK( M-KB+J1 ), KA1 ) 01323 * 01324 CALL ZLACGV( NR, WORK( M-KB+J1 ), KA1 ) 01325 END IF 01326 * 01327 * start applying rotations in 2nd set from the left 01328 * 01329 DO 900 L = KA - 1, KB - K + 1, -1 01330 NRT = ( J2+L-1 ) / KA1 01331 J1T = J2 - ( NRT-1 )*KA1 01332 IF( NRT.GT.0 ) 01333 $ CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA, 01334 $ AB( KA1-L, J1T-KA1+L ), INCA, 01335 $ RWORK( M-KB+J1T ), WORK( M-KB+J1T ), 01336 $ KA1 ) 01337 900 CONTINUE 01338 * 01339 IF( WANTX ) THEN 01340 * 01341 * post-multiply X by product of rotations in 2nd set 01342 * 01343 DO 910 J = J1, J2, KA1 01344 CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1, 01345 $ RWORK( M-KB+J ), DCONJG( WORK( M-KB+J ) ) ) 01346 910 CONTINUE 01347 END IF 01348 920 CONTINUE 01349 * 01350 DO 940 K = 1, KB - 1 01351 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1 01352 * 01353 * finish applying rotations in 1st set from the left 01354 * 01355 DO 930 L = KB - K, 1, -1 01356 NRT = ( J2+L-1 ) / KA1 01357 J1T = J2 - ( NRT-1 )*KA1 01358 IF( NRT.GT.0 ) 01359 $ CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA, 01360 $ AB( KA1-L, J1T-KA1+L ), INCA, 01361 $ RWORK( J1T ), WORK( J1T ), KA1 ) 01362 930 CONTINUE 01363 940 CONTINUE 01364 * 01365 IF( KB.GT.1 ) THEN 01366 DO 950 J = 2, I2 - KA 01367 RWORK( J ) = RWORK( J+KA ) 01368 WORK( J ) = WORK( J+KA ) 01369 950 CONTINUE 01370 END IF 01371 * 01372 END IF 01373 * 01374 GO TO 490 01375 * 01376 * End of ZHBGST 01377 * 01378 END
1.7.3 
