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