|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE CHBGST( 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 REAL RWORK( * ) 00015 COMPLEX AB( LDAB, * ), BB( LDBB, * ), WORK( * ), 00016 $ X( LDX, * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * CHBGST 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 CPBSTF, 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 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 array, dimension (LDBB,N) 00068 * The banded factor S from the split Cholesky factorization of 00069 * B, as returned by CPBSTF, 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 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 array, dimension (N) 00084 * 00085 * RWORK (workspace) REAL 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 CZERO, CONE 00095 REAL ONE 00096 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 00097 $ CONE = ( 1.0E+0, 0.0E+0 ), ONE = 1.0E+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 REAL BII 00104 COMPLEX RA, RA1, T 00105 * .. 00106 * .. External Functions .. 00107 LOGICAL LSAME 00108 EXTERNAL LSAME 00109 * .. 00110 * .. External Subroutines .. 00111 EXTERNAL CGERC, CGERU, CLACGV, CLAR2V, CLARGV, CLARTG, 00112 $ CLARTV, CLASET, CROT, CSSCAL, XERBLA 00113 * .. 00114 * .. Intrinsic Functions .. 00115 INTRINSIC CONJG, MAX, MIN, REAL 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( 'CHBGST', -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 CLASET( '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 CPBSTF. 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 = REAL( BB( KB1, I ) ) 00257 AB( KA1, I ) = ( REAL( 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 $ CONJG( AB( K-I+KA1, I ) ) - 00269 $ CONJG( BB( K-I+KB1, I ) )* 00270 $ AB( J-I+KA1, I ) + 00271 $ REAL( AB( KA1, I ) )* 00272 $ BB( J-I+KB1, I )* 00273 $ CONJG( 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 $ CONJG( 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 CSSCAL( N-M, ONE / BII, X( M+1, I ), 1 ) 00293 IF( KBT.GT.0 ) 00294 $ CALL CGERC( 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 CLARTG( 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 $ CONJG( 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 CLARGV( 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 CLARTV( 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 CLAR2V( 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 CLACGV( 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 CLARTV( 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 CROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1, 00393 $ RWORK( J-M ), CONJG( 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 CLARTV( 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 CLARGV( 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 CLARTV( 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 CLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ), 00468 $ AB( KA, J2+1 ), INCA, RWORK( J2 ), 00469 $ WORK( J2 ), KA1 ) 00470 * 00471 CALL CLACGV( 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 CLARTV( 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 CROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1, 00490 $ RWORK( J ), CONJG( 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 CLARTV( 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 = REAL( BB( 1, I ) ) 00525 AB( 1, I ) = ( REAL( 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 )*CONJG( AB( I-K+1, 00536 $ K ) ) - CONJG( BB( I-K+1, K ) )* 00537 $ AB( I-J+1, J ) + REAL( AB( 1, I ) )* 00538 $ BB( I-J+1, J )*CONJG( 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 $ CONJG( 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 CSSCAL( N-M, ONE / BII, X( M+1, I ), 1 ) 00559 IF( KBT.GT.0 ) 00560 $ CALL CGERU( 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 CLARTG( 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 $ CONJG( WORK( I-K+KA-M ) )*AB( KA1, I-K ) 00593 AB( KA1, I-K ) = WORK( I-K+KA-M )*T + 00594 $ RWORK( I-K+KA-M )*AB( KA1, I-K ) 00595 RA1 = RA 00596 END IF 00597 END IF 00598 J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1 00599 NR = ( N-J2+KA ) / KA1 00600 J1 = J2 + ( NR-1 )*KA1 00601 IF( UPDATE ) THEN 00602 J2T = MAX( J2, I+2*KA-K+1 ) 00603 ELSE 00604 J2T = J2 00605 END IF 00606 NRT = ( N-J2T+KA ) / KA1 00607 DO 320 J = J2T, J1, KA1 00608 * 00609 * create nonzero element a(j+1,j-ka) outside the band 00610 * and store it in WORK(j-m) 00611 * 00612 WORK( J-M ) = WORK( J-M )*AB( KA1, J-KA+1 ) 00613 AB( KA1, J-KA+1 ) = RWORK( J-M )*AB( KA1, J-KA+1 ) 00614 320 CONTINUE 00615 * 00616 * generate rotations in 1st set to annihilate elements which 00617 * have been created outside the band 00618 * 00619 IF( NRT.GT.0 ) 00620 $ CALL CLARGV( NRT, AB( KA1, J2T-KA ), INCA, WORK( J2T-M ), 00621 $ KA1, RWORK( J2T-M ), KA1 ) 00622 IF( NR.GT.0 ) THEN 00623 * 00624 * apply rotations in 1st set from the left 00625 * 00626 DO 330 L = 1, KA - 1 00627 CALL CLARTV( NR, AB( L+1, J2-L ), INCA, 00628 $ AB( L+2, J2-L ), INCA, RWORK( J2-M ), 00629 $ WORK( J2-M ), KA1 ) 00630 330 CONTINUE 00631 * 00632 * apply rotations in 1st set from both sides to diagonal 00633 * blocks 00634 * 00635 CALL CLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ), 00636 $ INCA, RWORK( J2-M ), WORK( J2-M ), KA1 ) 00637 * 00638 CALL CLACGV( NR, WORK( J2-M ), KA1 ) 00639 END IF 00640 * 00641 * start applying rotations in 1st set from the right 00642 * 00643 DO 340 L = KA - 1, KB - K + 1, -1 00644 NRT = ( N-J2+L ) / KA1 00645 IF( NRT.GT.0 ) 00646 $ CALL CLARTV( NRT, AB( KA1-L+1, J2 ), INCA, 00647 $ AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ), 00648 $ WORK( J2-M ), KA1 ) 00649 340 CONTINUE 00650 * 00651 IF( WANTX ) THEN 00652 * 00653 * post-multiply X by product of rotations in 1st set 00654 * 00655 DO 350 J = J2, J1, KA1 00656 CALL CROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1, 00657 $ RWORK( J-M ), WORK( J-M ) ) 00658 350 CONTINUE 00659 END IF 00660 360 CONTINUE 00661 * 00662 IF( UPDATE ) THEN 00663 IF( I2.LE.N .AND. KBT.GT.0 ) THEN 00664 * 00665 * create nonzero element a(i-kbt+ka+1,i-kbt) outside the 00666 * band and store it in WORK(i-kbt) 00667 * 00668 WORK( I-KBT ) = -BB( KBT+1, I-KBT )*RA1 00669 END IF 00670 END IF 00671 * 00672 DO 400 K = KB, 1, -1 00673 IF( UPDATE ) THEN 00674 J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1 00675 ELSE 00676 J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1 00677 END IF 00678 * 00679 * finish applying rotations in 2nd set from the right 00680 * 00681 DO 370 L = KB - K, 1, -1 00682 NRT = ( N-J2+KA+L ) / KA1 00683 IF( NRT.GT.0 ) 00684 $ CALL CLARTV( NRT, AB( KA1-L+1, J2-KA ), INCA, 00685 $ AB( KA1-L, J2-KA+1 ), INCA, 00686 $ RWORK( J2-KA ), WORK( J2-KA ), KA1 ) 00687 370 CONTINUE 00688 NR = ( N-J2+KA ) / KA1 00689 J1 = J2 + ( NR-1 )*KA1 00690 DO 380 J = J1, J2, -KA1 00691 WORK( J ) = WORK( J-KA ) 00692 RWORK( J ) = RWORK( J-KA ) 00693 380 CONTINUE 00694 DO 390 J = J2, J1, KA1 00695 * 00696 * create nonzero element a(j+1,j-ka) outside the band 00697 * and store it in WORK(j) 00698 * 00699 WORK( J ) = WORK( J )*AB( KA1, J-KA+1 ) 00700 AB( KA1, J-KA+1 ) = RWORK( J )*AB( KA1, J-KA+1 ) 00701 390 CONTINUE 00702 IF( UPDATE ) THEN 00703 IF( I-K.LT.N-KA .AND. K.LE.KBT ) 00704 $ WORK( I-K+KA ) = WORK( I-K ) 00705 END IF 00706 400 CONTINUE 00707 * 00708 DO 440 K = KB, 1, -1 00709 J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1 00710 NR = ( N-J2+KA ) / KA1 00711 J1 = J2 + ( NR-1 )*KA1 00712 IF( NR.GT.0 ) THEN 00713 * 00714 * generate rotations in 2nd set to annihilate elements 00715 * which have been created outside the band 00716 * 00717 CALL CLARGV( NR, AB( KA1, J2-KA ), INCA, WORK( J2 ), KA1, 00718 $ RWORK( J2 ), KA1 ) 00719 * 00720 * apply rotations in 2nd set from the left 00721 * 00722 DO 410 L = 1, KA - 1 00723 CALL CLARTV( NR, AB( L+1, J2-L ), INCA, 00724 $ AB( L+2, J2-L ), INCA, RWORK( J2 ), 00725 $ WORK( J2 ), KA1 ) 00726 410 CONTINUE 00727 * 00728 * apply rotations in 2nd set from both sides to diagonal 00729 * blocks 00730 * 00731 CALL CLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ), 00732 $ INCA, RWORK( J2 ), WORK( J2 ), KA1 ) 00733 * 00734 CALL CLACGV( NR, WORK( J2 ), KA1 ) 00735 END IF 00736 * 00737 * start applying rotations in 2nd set from the right 00738 * 00739 DO 420 L = KA - 1, KB - K + 1, -1 00740 NRT = ( N-J2+L ) / KA1 00741 IF( NRT.GT.0 ) 00742 $ CALL CLARTV( NRT, AB( KA1-L+1, J2 ), INCA, 00743 $ AB( KA1-L, J2+1 ), INCA, RWORK( J2 ), 00744 $ WORK( J2 ), KA1 ) 00745 420 CONTINUE 00746 * 00747 IF( WANTX ) THEN 00748 * 00749 * post-multiply X by product of rotations in 2nd set 00750 * 00751 DO 430 J = J2, J1, KA1 00752 CALL CROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1, 00753 $ RWORK( J ), WORK( J ) ) 00754 430 CONTINUE 00755 END IF 00756 440 CONTINUE 00757 * 00758 DO 460 K = 1, KB - 1 00759 J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1 00760 * 00761 * finish applying rotations in 1st set from the right 00762 * 00763 DO 450 L = KB - K, 1, -1 00764 NRT = ( N-J2+L ) / KA1 00765 IF( NRT.GT.0 ) 00766 $ CALL CLARTV( NRT, AB( KA1-L+1, J2 ), INCA, 00767 $ AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ), 00768 $ WORK( J2-M ), KA1 ) 00769 450 CONTINUE 00770 460 CONTINUE 00771 * 00772 IF( KB.GT.1 ) THEN 00773 DO 470 J = N - 1, J2 + KA, -1 00774 RWORK( J-M ) = RWORK( J-KA-M ) 00775 WORK( J-M ) = WORK( J-KA-M ) 00776 470 CONTINUE 00777 END IF 00778 * 00779 END IF 00780 * 00781 GO TO 10 00782 * 00783 480 CONTINUE 00784 * 00785 * **************************** Phase 2 ***************************** 00786 * 00787 * The logical structure of this phase is: 00788 * 00789 * UPDATE = .TRUE. 00790 * DO I = 1, M 00791 * use S(i) to update A and create a new bulge 00792 * apply rotations to push all bulges KA positions upward 00793 * END DO 00794 * UPDATE = .FALSE. 00795 * DO I = M - KA - 1, 2, -1 00796 * apply rotations to push all bulges KA positions upward 00797 * END DO 00798 * 00799 * To avoid duplicating code, the two loops are merged. 00800 * 00801 UPDATE = .TRUE. 00802 I = 0 00803 490 CONTINUE 00804 IF( UPDATE ) THEN 00805 I = I + 1 00806 KBT = MIN( KB, M-I ) 00807 I0 = I + 1 00808 I1 = MAX( 1, I-KA ) 00809 I2 = I + KBT - KA1 00810 IF( I.GT.M ) THEN 00811 UPDATE = .FALSE. 00812 I = I - 1 00813 I0 = M + 1 00814 IF( KA.EQ.0 ) 00815 $ RETURN 00816 GO TO 490 00817 END IF 00818 ELSE 00819 I = I - KA 00820 IF( I.LT.2 ) 00821 $ RETURN 00822 END IF 00823 * 00824 IF( I.LT.M-KBT ) THEN 00825 NX = M 00826 ELSE 00827 NX = N 00828 END IF 00829 * 00830 IF( UPPER ) THEN 00831 * 00832 * Transform A, working with the upper triangle 00833 * 00834 IF( UPDATE ) THEN 00835 * 00836 * Form inv(S(i))**H * A * inv(S(i)) 00837 * 00838 BII = REAL( BB( KB1, I ) ) 00839 AB( KA1, I ) = ( REAL( AB( KA1, I ) ) / BII ) / BII 00840 DO 500 J = I1, I - 1 00841 AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII 00842 500 CONTINUE 00843 DO 510 J = I + 1, MIN( N, I+KA ) 00844 AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII 00845 510 CONTINUE 00846 DO 540 K = I + 1, I + KBT 00847 DO 520 J = K, I + KBT 00848 AB( K-J+KA1, J ) = AB( K-J+KA1, J ) - 00849 $ BB( I-J+KB1, J )* 00850 $ CONJG( AB( I-K+KA1, K ) ) - 00851 $ CONJG( BB( I-K+KB1, K ) )* 00852 $ AB( I-J+KA1, J ) + 00853 $ REAL( AB( KA1, I ) )* 00854 $ BB( I-J+KB1, J )* 00855 $ CONJG( BB( I-K+KB1, K ) ) 00856 520 CONTINUE 00857 DO 530 J = I + KBT + 1, MIN( N, I+KA ) 00858 AB( K-J+KA1, J ) = AB( K-J+KA1, J ) - 00859 $ CONJG( BB( I-K+KB1, K ) )* 00860 $ AB( I-J+KA1, J ) 00861 530 CONTINUE 00862 540 CONTINUE 00863 DO 560 J = I1, I 00864 DO 550 K = I + 1, MIN( J+KA, I+KBT ) 00865 AB( J-K+KA1, K ) = AB( J-K+KA1, K ) - 00866 $ BB( I-K+KB1, K )*AB( J-I+KA1, I ) 00867 550 CONTINUE 00868 560 CONTINUE 00869 * 00870 IF( WANTX ) THEN 00871 * 00872 * post-multiply X by inv(S(i)) 00873 * 00874 CALL CSSCAL( NX, ONE / BII, X( 1, I ), 1 ) 00875 IF( KBT.GT.0 ) 00876 $ CALL CGERU( NX, KBT, -CONE, X( 1, I ), 1, 00877 $ BB( KB, I+1 ), LDBB-1, X( 1, I+1 ), LDX ) 00878 END IF 00879 * 00880 * store a(i1,i) in RA1 for use in next loop over K 00881 * 00882 RA1 = AB( I1-I+KA1, I ) 00883 END IF 00884 * 00885 * Generate and apply vectors of rotations to chase all the 00886 * existing bulges KA positions up toward the top of the band 00887 * 00888 DO 610 K = 1, KB - 1 00889 IF( UPDATE ) THEN 00890 * 00891 * Determine the rotations which would annihilate the bulge 00892 * which has in theory just been created 00893 * 00894 IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN 00895 * 00896 * generate rotation to annihilate a(i+k-ka-1,i) 00897 * 00898 CALL CLARTG( AB( K+1, I ), RA1, RWORK( I+K-KA ), 00899 $ WORK( I+K-KA ), RA ) 00900 * 00901 * create nonzero element a(i+k-ka-1,i+k) outside the 00902 * band and store it in WORK(m-kb+i+k) 00903 * 00904 T = -BB( KB1-K, I+K )*RA1 00905 WORK( M-KB+I+K ) = RWORK( I+K-KA )*T - 00906 $ CONJG( WORK( I+K-KA ) )* 00907 $ AB( 1, I+K ) 00908 AB( 1, I+K ) = WORK( I+K-KA )*T + 00909 $ RWORK( I+K-KA )*AB( 1, I+K ) 00910 RA1 = RA 00911 END IF 00912 END IF 00913 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1 00914 NR = ( J2+KA-1 ) / KA1 00915 J1 = J2 - ( NR-1 )*KA1 00916 IF( UPDATE ) THEN 00917 J2T = MIN( J2, I-2*KA+K-1 ) 00918 ELSE 00919 J2T = J2 00920 END IF 00921 NRT = ( J2T+KA-1 ) / KA1 00922 DO 570 J = J1, J2T, KA1 00923 * 00924 * create nonzero element a(j-1,j+ka) outside the band 00925 * and store it in WORK(j) 00926 * 00927 WORK( J ) = WORK( J )*AB( 1, J+KA-1 ) 00928 AB( 1, J+KA-1 ) = RWORK( J )*AB( 1, J+KA-1 ) 00929 570 CONTINUE 00930 * 00931 * generate rotations in 1st set to annihilate elements which 00932 * have been created outside the band 00933 * 00934 IF( NRT.GT.0 ) 00935 $ CALL CLARGV( NRT, AB( 1, J1+KA ), INCA, WORK( J1 ), KA1, 00936 $ RWORK( J1 ), KA1 ) 00937 IF( NR.GT.0 ) THEN 00938 * 00939 * apply rotations in 1st set from the left 00940 * 00941 DO 580 L = 1, KA - 1 00942 CALL CLARTV( NR, AB( KA1-L, J1+L ), INCA, 00943 $ AB( KA-L, J1+L ), INCA, RWORK( J1 ), 00944 $ WORK( J1 ), KA1 ) 00945 580 CONTINUE 00946 * 00947 * apply rotations in 1st set from both sides to diagonal 00948 * blocks 00949 * 00950 CALL CLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ), 00951 $ AB( KA, J1 ), INCA, RWORK( J1 ), WORK( J1 ), 00952 $ KA1 ) 00953 * 00954 CALL CLACGV( NR, WORK( J1 ), KA1 ) 00955 END IF 00956 * 00957 * start applying rotations in 1st set from the right 00958 * 00959 DO 590 L = KA - 1, KB - K + 1, -1 00960 NRT = ( J2+L-1 ) / KA1 00961 J1T = J2 - ( NRT-1 )*KA1 00962 IF( NRT.GT.0 ) 00963 $ CALL CLARTV( NRT, AB( L, J1T ), INCA, 00964 $ AB( L+1, J1T-1 ), INCA, RWORK( J1T ), 00965 $ WORK( J1T ), KA1 ) 00966 590 CONTINUE 00967 * 00968 IF( WANTX ) THEN 00969 * 00970 * post-multiply X by product of rotations in 1st set 00971 * 00972 DO 600 J = J1, J2, KA1 00973 CALL CROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1, 00974 $ RWORK( J ), WORK( J ) ) 00975 600 CONTINUE 00976 END IF 00977 610 CONTINUE 00978 * 00979 IF( UPDATE ) THEN 00980 IF( I2.GT.0 .AND. KBT.GT.0 ) THEN 00981 * 00982 * create nonzero element a(i+kbt-ka-1,i+kbt) outside the 00983 * band and store it in WORK(m-kb+i+kbt) 00984 * 00985 WORK( M-KB+I+KBT ) = -BB( KB1-KBT, I+KBT )*RA1 00986 END IF 00987 END IF 00988 * 00989 DO 650 K = KB, 1, -1 00990 IF( UPDATE ) THEN 00991 J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1 00992 ELSE 00993 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1 00994 END IF 00995 * 00996 * finish applying rotations in 2nd set from the right 00997 * 00998 DO 620 L = KB - K, 1, -1 00999 NRT = ( J2+KA+L-1 ) / KA1 01000 J1T = J2 - ( NRT-1 )*KA1 01001 IF( NRT.GT.0 ) 01002 $ CALL CLARTV( NRT, AB( L, J1T+KA ), INCA, 01003 $ AB( L+1, J1T+KA-1 ), INCA, 01004 $ RWORK( M-KB+J1T+KA ), 01005 $ WORK( M-KB+J1T+KA ), KA1 ) 01006 620 CONTINUE 01007 NR = ( J2+KA-1 ) / KA1 01008 J1 = J2 - ( NR-1 )*KA1 01009 DO 630 J = J1, J2, KA1 01010 WORK( M-KB+J ) = WORK( M-KB+J+KA ) 01011 RWORK( M-KB+J ) = RWORK( M-KB+J+KA ) 01012 630 CONTINUE 01013 DO 640 J = J1, J2, KA1 01014 * 01015 * create nonzero element a(j-1,j+ka) outside the band 01016 * and store it in WORK(m-kb+j) 01017 * 01018 WORK( M-KB+J ) = WORK( M-KB+J )*AB( 1, J+KA-1 ) 01019 AB( 1, J+KA-1 ) = RWORK( M-KB+J )*AB( 1, J+KA-1 ) 01020 640 CONTINUE 01021 IF( UPDATE ) THEN 01022 IF( I+K.GT.KA1 .AND. K.LE.KBT ) 01023 $ WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K ) 01024 END IF 01025 650 CONTINUE 01026 * 01027 DO 690 K = KB, 1, -1 01028 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1 01029 NR = ( J2+KA-1 ) / KA1 01030 J1 = J2 - ( NR-1 )*KA1 01031 IF( NR.GT.0 ) THEN 01032 * 01033 * generate rotations in 2nd set to annihilate elements 01034 * which have been created outside the band 01035 * 01036 CALL CLARGV( NR, AB( 1, J1+KA ), INCA, WORK( M-KB+J1 ), 01037 $ KA1, RWORK( M-KB+J1 ), KA1 ) 01038 * 01039 * apply rotations in 2nd set from the left 01040 * 01041 DO 660 L = 1, KA - 1 01042 CALL CLARTV( NR, AB( KA1-L, J1+L ), INCA, 01043 $ AB( KA-L, J1+L ), INCA, RWORK( M-KB+J1 ), 01044 $ WORK( M-KB+J1 ), KA1 ) 01045 660 CONTINUE 01046 * 01047 * apply rotations in 2nd set from both sides to diagonal 01048 * blocks 01049 * 01050 CALL CLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ), 01051 $ AB( KA, J1 ), INCA, RWORK( M-KB+J1 ), 01052 $ WORK( M-KB+J1 ), KA1 ) 01053 * 01054 CALL CLACGV( NR, WORK( M-KB+J1 ), KA1 ) 01055 END IF 01056 * 01057 * start applying rotations in 2nd set from the right 01058 * 01059 DO 670 L = KA - 1, KB - K + 1, -1 01060 NRT = ( J2+L-1 ) / KA1 01061 J1T = J2 - ( NRT-1 )*KA1 01062 IF( NRT.GT.0 ) 01063 $ CALL CLARTV( NRT, AB( L, J1T ), INCA, 01064 $ AB( L+1, J1T-1 ), INCA, 01065 $ RWORK( M-KB+J1T ), WORK( M-KB+J1T ), 01066 $ KA1 ) 01067 670 CONTINUE 01068 * 01069 IF( WANTX ) THEN 01070 * 01071 * post-multiply X by product of rotations in 2nd set 01072 * 01073 DO 680 J = J1, J2, KA1 01074 CALL CROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1, 01075 $ RWORK( M-KB+J ), WORK( M-KB+J ) ) 01076 680 CONTINUE 01077 END IF 01078 690 CONTINUE 01079 * 01080 DO 710 K = 1, KB - 1 01081 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1 01082 * 01083 * finish applying rotations in 1st set from the right 01084 * 01085 DO 700 L = KB - K, 1, -1 01086 NRT = ( J2+L-1 ) / KA1 01087 J1T = J2 - ( NRT-1 )*KA1 01088 IF( NRT.GT.0 ) 01089 $ CALL CLARTV( NRT, AB( L, J1T ), INCA, 01090 $ AB( L+1, J1T-1 ), INCA, RWORK( J1T ), 01091 $ WORK( J1T ), KA1 ) 01092 700 CONTINUE 01093 710 CONTINUE 01094 * 01095 IF( KB.GT.1 ) THEN 01096 DO 720 J = 2, I2 - KA 01097 RWORK( J ) = RWORK( J+KA ) 01098 WORK( J ) = WORK( J+KA ) 01099 720 CONTINUE 01100 END IF 01101 * 01102 ELSE 01103 * 01104 * Transform A, working with the lower triangle 01105 * 01106 IF( UPDATE ) THEN 01107 * 01108 * Form inv(S(i))**H * A * inv(S(i)) 01109 * 01110 BII = REAL( BB( 1, I ) ) 01111 AB( 1, I ) = ( REAL( AB( 1, I ) ) / BII ) / BII 01112 DO 730 J = I1, I - 1 01113 AB( I-J+1, J ) = AB( I-J+1, J ) / BII 01114 730 CONTINUE 01115 DO 740 J = I + 1, MIN( N, I+KA ) 01116 AB( J-I+1, I ) = AB( J-I+1, I ) / BII 01117 740 CONTINUE 01118 DO 770 K = I + 1, I + KBT 01119 DO 750 J = K, I + KBT 01120 AB( J-K+1, K ) = AB( J-K+1, K ) - 01121 $ BB( J-I+1, I )*CONJG( AB( K-I+1, 01122 $ I ) ) - CONJG( BB( K-I+1, I ) )* 01123 $ AB( J-I+1, I ) + REAL( AB( 1, I ) )* 01124 $ BB( J-I+1, I )*CONJG( BB( K-I+1, 01125 $ I ) ) 01126 750 CONTINUE 01127 DO 760 J = I + KBT + 1, MIN( N, I+KA ) 01128 AB( J-K+1, K ) = AB( J-K+1, K ) - 01129 $ CONJG( BB( K-I+1, I ) )* 01130 $ AB( J-I+1, I ) 01131 760 CONTINUE 01132 770 CONTINUE 01133 DO 790 J = I1, I 01134 DO 780 K = I + 1, MIN( J+KA, I+KBT ) 01135 AB( K-J+1, J ) = AB( K-J+1, J ) - 01136 $ BB( K-I+1, I )*AB( I-J+1, J ) 01137 780 CONTINUE 01138 790 CONTINUE 01139 * 01140 IF( WANTX ) THEN 01141 * 01142 * post-multiply X by inv(S(i)) 01143 * 01144 CALL CSSCAL( NX, ONE / BII, X( 1, I ), 1 ) 01145 IF( KBT.GT.0 ) 01146 $ CALL CGERC( NX, KBT, -CONE, X( 1, I ), 1, BB( 2, I ), 01147 $ 1, X( 1, I+1 ), LDX ) 01148 END IF 01149 * 01150 * store a(i,i1) in RA1 for use in next loop over K 01151 * 01152 RA1 = AB( I-I1+1, I1 ) 01153 END IF 01154 * 01155 * Generate and apply vectors of rotations to chase all the 01156 * existing bulges KA positions up toward the top of the band 01157 * 01158 DO 840 K = 1, KB - 1 01159 IF( UPDATE ) THEN 01160 * 01161 * Determine the rotations which would annihilate the bulge 01162 * which has in theory just been created 01163 * 01164 IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN 01165 * 01166 * generate rotation to annihilate a(i,i+k-ka-1) 01167 * 01168 CALL CLARTG( AB( KA1-K, I+K-KA ), RA1, 01169 $ RWORK( I+K-KA ), WORK( I+K-KA ), RA ) 01170 * 01171 * create nonzero element a(i+k,i+k-ka-1) outside the 01172 * band and store it in WORK(m-kb+i+k) 01173 * 01174 T = -BB( K+1, I )*RA1 01175 WORK( M-KB+I+K ) = RWORK( I+K-KA )*T - 01176 $ CONJG( WORK( I+K-KA ) )* 01177 $ AB( KA1, I+K-KA ) 01178 AB( KA1, I+K-KA ) = WORK( I+K-KA )*T + 01179 $ RWORK( I+K-KA )*AB( KA1, I+K-KA ) 01180 RA1 = RA 01181 END IF 01182 END IF 01183 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1 01184 NR = ( J2+KA-1 ) / KA1 01185 J1 = J2 - ( NR-1 )*KA1 01186 IF( UPDATE ) THEN 01187 J2T = MIN( J2, I-2*KA+K-1 ) 01188 ELSE 01189 J2T = J2 01190 END IF 01191 NRT = ( J2T+KA-1 ) / KA1 01192 DO 800 J = J1, J2T, KA1 01193 * 01194 * create nonzero element a(j+ka,j-1) outside the band 01195 * and store it in WORK(j) 01196 * 01197 WORK( J ) = WORK( J )*AB( KA1, J-1 ) 01198 AB( KA1, J-1 ) = RWORK( J )*AB( KA1, J-1 ) 01199 800 CONTINUE 01200 * 01201 * generate rotations in 1st set to annihilate elements which 01202 * have been created outside the band 01203 * 01204 IF( NRT.GT.0 ) 01205 $ CALL CLARGV( NRT, AB( KA1, J1 ), INCA, WORK( J1 ), KA1, 01206 $ RWORK( J1 ), KA1 ) 01207 IF( NR.GT.0 ) THEN 01208 * 01209 * apply rotations in 1st set from the right 01210 * 01211 DO 810 L = 1, KA - 1 01212 CALL CLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ), 01213 $ INCA, RWORK( J1 ), WORK( J1 ), KA1 ) 01214 810 CONTINUE 01215 * 01216 * apply rotations in 1st set from both sides to diagonal 01217 * blocks 01218 * 01219 CALL CLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ), 01220 $ AB( 2, J1-1 ), INCA, RWORK( J1 ), 01221 $ WORK( J1 ), KA1 ) 01222 * 01223 CALL CLACGV( NR, WORK( J1 ), KA1 ) 01224 END IF 01225 * 01226 * start applying rotations in 1st set from the left 01227 * 01228 DO 820 L = KA - 1, KB - K + 1, -1 01229 NRT = ( J2+L-1 ) / KA1 01230 J1T = J2 - ( NRT-1 )*KA1 01231 IF( NRT.GT.0 ) 01232 $ CALL CLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA, 01233 $ AB( KA1-L, J1T-KA1+L ), INCA, 01234 $ RWORK( J1T ), WORK( J1T ), KA1 ) 01235 820 CONTINUE 01236 * 01237 IF( WANTX ) THEN 01238 * 01239 * post-multiply X by product of rotations in 1st set 01240 * 01241 DO 830 J = J1, J2, KA1 01242 CALL CROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1, 01243 $ RWORK( J ), CONJG( WORK( J ) ) ) 01244 830 CONTINUE 01245 END IF 01246 840 CONTINUE 01247 * 01248 IF( UPDATE ) THEN 01249 IF( I2.GT.0 .AND. KBT.GT.0 ) THEN 01250 * 01251 * create nonzero element a(i+kbt,i+kbt-ka-1) outside the 01252 * band and store it in WORK(m-kb+i+kbt) 01253 * 01254 WORK( M-KB+I+KBT ) = -BB( KBT+1, I )*RA1 01255 END IF 01256 END IF 01257 * 01258 DO 880 K = KB, 1, -1 01259 IF( UPDATE ) THEN 01260 J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1 01261 ELSE 01262 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1 01263 END IF 01264 * 01265 * finish applying rotations in 2nd set from the left 01266 * 01267 DO 850 L = KB - K, 1, -1 01268 NRT = ( J2+KA+L-1 ) / KA1 01269 J1T = J2 - ( NRT-1 )*KA1 01270 IF( NRT.GT.0 ) 01271 $ CALL CLARTV( NRT, AB( KA1-L+1, J1T+L-1 ), INCA, 01272 $ AB( KA1-L, J1T+L-1 ), INCA, 01273 $ RWORK( M-KB+J1T+KA ), 01274 $ WORK( M-KB+J1T+KA ), KA1 ) 01275 850 CONTINUE 01276 NR = ( J2+KA-1 ) / KA1 01277 J1 = J2 - ( NR-1 )*KA1 01278 DO 860 J = J1, J2, KA1 01279 WORK( M-KB+J ) = WORK( M-KB+J+KA ) 01280 RWORK( M-KB+J ) = RWORK( M-KB+J+KA ) 01281 860 CONTINUE 01282 DO 870 J = J1, J2, KA1 01283 * 01284 * create nonzero element a(j+ka,j-1) outside the band 01285 * and store it in WORK(m-kb+j) 01286 * 01287 WORK( M-KB+J ) = WORK( M-KB+J )*AB( KA1, J-1 ) 01288 AB( KA1, J-1 ) = RWORK( M-KB+J )*AB( KA1, J-1 ) 01289 870 CONTINUE 01290 IF( UPDATE ) THEN 01291 IF( I+K.GT.KA1 .AND. K.LE.KBT ) 01292 $ WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K ) 01293 END IF 01294 880 CONTINUE 01295 * 01296 DO 920 K = KB, 1, -1 01297 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1 01298 NR = ( J2+KA-1 ) / KA1 01299 J1 = J2 - ( NR-1 )*KA1 01300 IF( NR.GT.0 ) THEN 01301 * 01302 * generate rotations in 2nd set to annihilate elements 01303 * which have been created outside the band 01304 * 01305 CALL CLARGV( NR, AB( KA1, J1 ), INCA, WORK( M-KB+J1 ), 01306 $ KA1, RWORK( M-KB+J1 ), KA1 ) 01307 * 01308 * apply rotations in 2nd set from the right 01309 * 01310 DO 890 L = 1, KA - 1 01311 CALL CLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ), 01312 $ INCA, RWORK( M-KB+J1 ), WORK( M-KB+J1 ), 01313 $ KA1 ) 01314 890 CONTINUE 01315 * 01316 * apply rotations in 2nd set from both sides to diagonal 01317 * blocks 01318 * 01319 CALL CLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ), 01320 $ AB( 2, J1-1 ), INCA, RWORK( M-KB+J1 ), 01321 $ WORK( M-KB+J1 ), KA1 ) 01322 * 01323 CALL CLACGV( NR, WORK( M-KB+J1 ), KA1 ) 01324 END IF 01325 * 01326 * start applying rotations in 2nd set from the left 01327 * 01328 DO 900 L = KA - 1, KB - K + 1, -1 01329 NRT = ( J2+L-1 ) / KA1 01330 J1T = J2 - ( NRT-1 )*KA1 01331 IF( NRT.GT.0 ) 01332 $ CALL CLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA, 01333 $ AB( KA1-L, J1T-KA1+L ), INCA, 01334 $ RWORK( M-KB+J1T ), WORK( M-KB+J1T ), 01335 $ KA1 ) 01336 900 CONTINUE 01337 * 01338 IF( WANTX ) THEN 01339 * 01340 * post-multiply X by product of rotations in 2nd set 01341 * 01342 DO 910 J = J1, J2, KA1 01343 CALL CROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1, 01344 $ RWORK( M-KB+J ), CONJG( WORK( M-KB+J ) ) ) 01345 910 CONTINUE 01346 END IF 01347 920 CONTINUE 01348 * 01349 DO 940 K = 1, KB - 1 01350 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1 01351 * 01352 * finish applying rotations in 1st set from the left 01353 * 01354 DO 930 L = KB - K, 1, -1 01355 NRT = ( J2+L-1 ) / KA1 01356 J1T = J2 - ( NRT-1 )*KA1 01357 IF( NRT.GT.0 ) 01358 $ CALL CLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA, 01359 $ AB( KA1-L, J1T-KA1+L ), INCA, 01360 $ RWORK( J1T ), WORK( J1T ), KA1 ) 01361 930 CONTINUE 01362 940 CONTINUE 01363 * 01364 IF( KB.GT.1 ) THEN 01365 DO 950 J = 2, I2 - KA 01366 RWORK( J ) = RWORK( J+KA ) 01367 WORK( J ) = WORK( J+KA ) 01368 950 CONTINUE 01369 END IF 01370 * 01371 END IF 01372 * 01373 GO TO 490 01374 * 01375 * End of CHBGST 01376 * 01377 END
1.7.3 
