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