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