![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZHBGST 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZHBGST + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhbgst.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhbgst.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhbgst.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZHBGST( 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 * DOUBLE PRECISION RWORK( * ) 00030 * COMPLEX*16 AB( LDAB, * ), BB( LDBB, * ), WORK( * ), 00031 * $ X( LDX, * ) 00032 * .. 00033 * 00034 * 00035 *> \par Purpose: 00036 * ============= 00037 *> 00038 *> \verbatim 00039 *> 00040 *> ZHBGST 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 ZPBSTF, 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*16 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*16 array, dimension (LDBB,N) 00110 *> The banded factor S from the split Cholesky factorization of 00111 *> B, as returned by ZPBSTF, 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*16 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*16 array, dimension (N) 00138 *> \endverbatim 00139 *> 00140 *> \param[out] RWORK 00141 *> \verbatim 00142 *> RWORK is DOUBLE PRECISION 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 complex16OTHERcomputational 00163 * 00164 * ===================================================================== 00165 SUBROUTINE ZHBGST( 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 DOUBLE PRECISION RWORK( * ) 00179 COMPLEX*16 AB( LDAB, * ), BB( LDBB, * ), WORK( * ), 00180 $ X( LDX, * ) 00181 * .. 00182 * 00183 * ===================================================================== 00184 * 00185 * .. Parameters .. 00186 COMPLEX*16 CZERO, CONE 00187 DOUBLE PRECISION ONE 00188 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 00189 $ CONE = ( 1.0D+0, 0.0D+0 ), ONE = 1.0D+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 DOUBLE PRECISION BII 00196 COMPLEX*16 RA, RA1, T 00197 * .. 00198 * .. External Functions .. 00199 LOGICAL LSAME 00200 EXTERNAL LSAME 00201 * .. 00202 * .. External Subroutines .. 00203 EXTERNAL XERBLA, ZDSCAL, ZGERC, ZGERU, ZLACGV, ZLAR2V, 00204 $ ZLARGV, ZLARTG, ZLARTV, ZLASET, ZROT 00205 * .. 00206 * .. Intrinsic Functions .. 00207 INTRINSIC DBLE, DCONJG, MAX, MIN 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( 'ZHBGST', -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 ZLASET( '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 ZPBSTF. 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 = DBLE( BB( KB1, I ) ) 00349 AB( KA1, I ) = ( DBLE( 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 $ DCONJG( AB( K-I+KA1, I ) ) - 00361 $ DCONJG( BB( K-I+KB1, I ) )* 00362 $ AB( J-I+KA1, I ) + 00363 $ DBLE( AB( KA1, I ) )* 00364 $ BB( J-I+KB1, I )* 00365 $ DCONJG( 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 $ DCONJG( 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 ZDSCAL( N-M, ONE / BII, X( M+1, I ), 1 ) 00385 IF( KBT.GT.0 ) 00386 $ CALL ZGERC( 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 ZLARTG( 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 $ DCONJG( 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 ZLARGV( 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 ZLARTV( 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 ZLAR2V( 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 ZLACGV( 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 ZLARTV( 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 ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1, 00485 $ RWORK( J-M ), DCONJG( 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 ZLARTV( 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 ZLARGV( 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 ZLARTV( 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 ZLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ), 00560 $ AB( KA, J2+1 ), INCA, RWORK( J2 ), 00561 $ WORK( J2 ), KA1 ) 00562 * 00563 CALL ZLACGV( 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 ZLARTV( 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 ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1, 00582 $ RWORK( J ), DCONJG( 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 ZLARTV( 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 = DBLE( BB( 1, I ) ) 00617 AB( 1, I ) = ( DBLE( 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 )*DCONJG( AB( I-K+1, 00628 $ K ) ) - DCONJG( BB( I-K+1, K ) )* 00629 $ AB( I-J+1, J ) + DBLE( AB( 1, I ) )* 00630 $ BB( I-J+1, J )*DCONJG( 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 $ DCONJG( 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 ZDSCAL( N-M, ONE / BII, X( M+1, I ), 1 ) 00651 IF( KBT.GT.0 ) 00652 $ CALL ZGERU( 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 ZLARTG( 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 $ DCONJG( WORK( I-K+KA-M ) )* 00685 $ AB( KA1, I-K ) 00686 AB( KA1, I-K ) = WORK( I-K+KA-M )*T + 00687 $ RWORK( I-K+KA-M )*AB( KA1, I-K ) 00688 RA1 = RA 00689 END IF 00690 END IF 00691 J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1 00692 NR = ( N-J2+KA ) / KA1 00693 J1 = J2 + ( NR-1 )*KA1 00694 IF( UPDATE ) THEN 00695 J2T = MAX( J2, I+2*KA-K+1 ) 00696 ELSE 00697 J2T = J2 00698 END IF 00699 NRT = ( N-J2T+KA ) / KA1 00700 DO 320 J = J2T, J1, KA1 00701 * 00702 * create nonzero element a(j+1,j-ka) outside the band 00703 * and store it in WORK(j-m) 00704 * 00705 WORK( J-M ) = WORK( J-M )*AB( KA1, J-KA+1 ) 00706 AB( KA1, J-KA+1 ) = RWORK( J-M )*AB( KA1, J-KA+1 ) 00707 320 CONTINUE 00708 * 00709 * generate rotations in 1st set to annihilate elements which 00710 * have been created outside the band 00711 * 00712 IF( NRT.GT.0 ) 00713 $ CALL ZLARGV( NRT, AB( KA1, J2T-KA ), INCA, WORK( J2T-M ), 00714 $ KA1, RWORK( J2T-M ), KA1 ) 00715 IF( NR.GT.0 ) THEN 00716 * 00717 * apply rotations in 1st set from the left 00718 * 00719 DO 330 L = 1, KA - 1 00720 CALL ZLARTV( NR, AB( L+1, J2-L ), INCA, 00721 $ AB( L+2, J2-L ), INCA, RWORK( J2-M ), 00722 $ WORK( J2-M ), KA1 ) 00723 330 CONTINUE 00724 * 00725 * apply rotations in 1st set from both sides to diagonal 00726 * blocks 00727 * 00728 CALL ZLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ), 00729 $ INCA, RWORK( J2-M ), WORK( J2-M ), KA1 ) 00730 * 00731 CALL ZLACGV( NR, WORK( J2-M ), KA1 ) 00732 END IF 00733 * 00734 * start applying rotations in 1st set from the right 00735 * 00736 DO 340 L = KA - 1, KB - K + 1, -1 00737 NRT = ( N-J2+L ) / KA1 00738 IF( NRT.GT.0 ) 00739 $ CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA, 00740 $ AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ), 00741 $ WORK( J2-M ), KA1 ) 00742 340 CONTINUE 00743 * 00744 IF( WANTX ) THEN 00745 * 00746 * post-multiply X by product of rotations in 1st set 00747 * 00748 DO 350 J = J2, J1, KA1 00749 CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1, 00750 $ RWORK( J-M ), WORK( J-M ) ) 00751 350 CONTINUE 00752 END IF 00753 360 CONTINUE 00754 * 00755 IF( UPDATE ) THEN 00756 IF( I2.LE.N .AND. KBT.GT.0 ) THEN 00757 * 00758 * create nonzero element a(i-kbt+ka+1,i-kbt) outside the 00759 * band and store it in WORK(i-kbt) 00760 * 00761 WORK( I-KBT ) = -BB( KBT+1, I-KBT )*RA1 00762 END IF 00763 END IF 00764 * 00765 DO 400 K = KB, 1, -1 00766 IF( UPDATE ) THEN 00767 J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1 00768 ELSE 00769 J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1 00770 END IF 00771 * 00772 * finish applying rotations in 2nd set from the right 00773 * 00774 DO 370 L = KB - K, 1, -1 00775 NRT = ( N-J2+KA+L ) / KA1 00776 IF( NRT.GT.0 ) 00777 $ CALL ZLARTV( NRT, AB( KA1-L+1, J2-KA ), INCA, 00778 $ AB( KA1-L, J2-KA+1 ), INCA, 00779 $ RWORK( J2-KA ), WORK( J2-KA ), KA1 ) 00780 370 CONTINUE 00781 NR = ( N-J2+KA ) / KA1 00782 J1 = J2 + ( NR-1 )*KA1 00783 DO 380 J = J1, J2, -KA1 00784 WORK( J ) = WORK( J-KA ) 00785 RWORK( J ) = RWORK( J-KA ) 00786 380 CONTINUE 00787 DO 390 J = J2, J1, KA1 00788 * 00789 * create nonzero element a(j+1,j-ka) outside the band 00790 * and store it in WORK(j) 00791 * 00792 WORK( J ) = WORK( J )*AB( KA1, J-KA+1 ) 00793 AB( KA1, J-KA+1 ) = RWORK( J )*AB( KA1, J-KA+1 ) 00794 390 CONTINUE 00795 IF( UPDATE ) THEN 00796 IF( I-K.LT.N-KA .AND. K.LE.KBT ) 00797 $ WORK( I-K+KA ) = WORK( I-K ) 00798 END IF 00799 400 CONTINUE 00800 * 00801 DO 440 K = KB, 1, -1 00802 J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1 00803 NR = ( N-J2+KA ) / KA1 00804 J1 = J2 + ( NR-1 )*KA1 00805 IF( NR.GT.0 ) THEN 00806 * 00807 * generate rotations in 2nd set to annihilate elements 00808 * which have been created outside the band 00809 * 00810 CALL ZLARGV( NR, AB( KA1, J2-KA ), INCA, WORK( J2 ), KA1, 00811 $ RWORK( J2 ), KA1 ) 00812 * 00813 * apply rotations in 2nd set from the left 00814 * 00815 DO 410 L = 1, KA - 1 00816 CALL ZLARTV( NR, AB( L+1, J2-L ), INCA, 00817 $ AB( L+2, J2-L ), INCA, RWORK( J2 ), 00818 $ WORK( J2 ), KA1 ) 00819 410 CONTINUE 00820 * 00821 * apply rotations in 2nd set from both sides to diagonal 00822 * blocks 00823 * 00824 CALL ZLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ), 00825 $ INCA, RWORK( J2 ), WORK( J2 ), KA1 ) 00826 * 00827 CALL ZLACGV( NR, WORK( J2 ), KA1 ) 00828 END IF 00829 * 00830 * start applying rotations in 2nd set from the right 00831 * 00832 DO 420 L = KA - 1, KB - K + 1, -1 00833 NRT = ( N-J2+L ) / KA1 00834 IF( NRT.GT.0 ) 00835 $ CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA, 00836 $ AB( KA1-L, J2+1 ), INCA, RWORK( J2 ), 00837 $ WORK( J2 ), KA1 ) 00838 420 CONTINUE 00839 * 00840 IF( WANTX ) THEN 00841 * 00842 * post-multiply X by product of rotations in 2nd set 00843 * 00844 DO 430 J = J2, J1, KA1 00845 CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1, 00846 $ RWORK( J ), WORK( J ) ) 00847 430 CONTINUE 00848 END IF 00849 440 CONTINUE 00850 * 00851 DO 460 K = 1, KB - 1 00852 J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1 00853 * 00854 * finish applying rotations in 1st set from the right 00855 * 00856 DO 450 L = KB - K, 1, -1 00857 NRT = ( N-J2+L ) / KA1 00858 IF( NRT.GT.0 ) 00859 $ CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA, 00860 $ AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ), 00861 $ WORK( J2-M ), KA1 ) 00862 450 CONTINUE 00863 460 CONTINUE 00864 * 00865 IF( KB.GT.1 ) THEN 00866 DO 470 J = N - 1, J2 + KA, -1 00867 RWORK( J-M ) = RWORK( J-KA-M ) 00868 WORK( J-M ) = WORK( J-KA-M ) 00869 470 CONTINUE 00870 END IF 00871 * 00872 END IF 00873 * 00874 GO TO 10 00875 * 00876 480 CONTINUE 00877 * 00878 * **************************** Phase 2 ***************************** 00879 * 00880 * The logical structure of this phase is: 00881 * 00882 * UPDATE = .TRUE. 00883 * DO I = 1, M 00884 * use S(i) to update A and create a new bulge 00885 * apply rotations to push all bulges KA positions upward 00886 * END DO 00887 * UPDATE = .FALSE. 00888 * DO I = M - KA - 1, 2, -1 00889 * apply rotations to push all bulges KA positions upward 00890 * END DO 00891 * 00892 * To avoid duplicating code, the two loops are merged. 00893 * 00894 UPDATE = .TRUE. 00895 I = 0 00896 490 CONTINUE 00897 IF( UPDATE ) THEN 00898 I = I + 1 00899 KBT = MIN( KB, M-I ) 00900 I0 = I + 1 00901 I1 = MAX( 1, I-KA ) 00902 I2 = I + KBT - KA1 00903 IF( I.GT.M ) THEN 00904 UPDATE = .FALSE. 00905 I = I - 1 00906 I0 = M + 1 00907 IF( KA.EQ.0 ) 00908 $ RETURN 00909 GO TO 490 00910 END IF 00911 ELSE 00912 I = I - KA 00913 IF( I.LT.2 ) 00914 $ RETURN 00915 END IF 00916 * 00917 IF( I.LT.M-KBT ) THEN 00918 NX = M 00919 ELSE 00920 NX = N 00921 END IF 00922 * 00923 IF( UPPER ) THEN 00924 * 00925 * Transform A, working with the upper triangle 00926 * 00927 IF( UPDATE ) THEN 00928 * 00929 * Form inv(S(i))**H * A * inv(S(i)) 00930 * 00931 BII = DBLE( BB( KB1, I ) ) 00932 AB( KA1, I ) = ( DBLE( AB( KA1, I ) ) / BII ) / BII 00933 DO 500 J = I1, I - 1 00934 AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII 00935 500 CONTINUE 00936 DO 510 J = I + 1, MIN( N, I+KA ) 00937 AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII 00938 510 CONTINUE 00939 DO 540 K = I + 1, I + KBT 00940 DO 520 J = K, I + KBT 00941 AB( K-J+KA1, J ) = AB( K-J+KA1, J ) - 00942 $ BB( I-J+KB1, J )* 00943 $ DCONJG( AB( I-K+KA1, K ) ) - 00944 $ DCONJG( BB( I-K+KB1, K ) )* 00945 $ AB( I-J+KA1, J ) + 00946 $ DBLE( AB( KA1, I ) )* 00947 $ BB( I-J+KB1, J )* 00948 $ DCONJG( BB( I-K+KB1, K ) ) 00949 520 CONTINUE 00950 DO 530 J = I + KBT + 1, MIN( N, I+KA ) 00951 AB( K-J+KA1, J ) = AB( K-J+KA1, J ) - 00952 $ DCONJG( BB( I-K+KB1, K ) )* 00953 $ AB( I-J+KA1, J ) 00954 530 CONTINUE 00955 540 CONTINUE 00956 DO 560 J = I1, I 00957 DO 550 K = I + 1, MIN( J+KA, I+KBT ) 00958 AB( J-K+KA1, K ) = AB( J-K+KA1, K ) - 00959 $ BB( I-K+KB1, K )*AB( J-I+KA1, I ) 00960 550 CONTINUE 00961 560 CONTINUE 00962 * 00963 IF( WANTX ) THEN 00964 * 00965 * post-multiply X by inv(S(i)) 00966 * 00967 CALL ZDSCAL( NX, ONE / BII, X( 1, I ), 1 ) 00968 IF( KBT.GT.0 ) 00969 $ CALL ZGERU( NX, KBT, -CONE, X( 1, I ), 1, 00970 $ BB( KB, I+1 ), LDBB-1, X( 1, I+1 ), LDX ) 00971 END IF 00972 * 00973 * store a(i1,i) in RA1 for use in next loop over K 00974 * 00975 RA1 = AB( I1-I+KA1, I ) 00976 END IF 00977 * 00978 * Generate and apply vectors of rotations to chase all the 00979 * existing bulges KA positions up toward the top of the band 00980 * 00981 DO 610 K = 1, KB - 1 00982 IF( UPDATE ) THEN 00983 * 00984 * Determine the rotations which would annihilate the bulge 00985 * which has in theory just been created 00986 * 00987 IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN 00988 * 00989 * generate rotation to annihilate a(i+k-ka-1,i) 00990 * 00991 CALL ZLARTG( AB( K+1, I ), RA1, RWORK( I+K-KA ), 00992 $ WORK( I+K-KA ), RA ) 00993 * 00994 * create nonzero element a(i+k-ka-1,i+k) outside the 00995 * band and store it in WORK(m-kb+i+k) 00996 * 00997 T = -BB( KB1-K, I+K )*RA1 00998 WORK( M-KB+I+K ) = RWORK( I+K-KA )*T - 00999 $ DCONJG( WORK( I+K-KA ) )* 01000 $ AB( 1, I+K ) 01001 AB( 1, I+K ) = WORK( I+K-KA )*T + 01002 $ RWORK( I+K-KA )*AB( 1, I+K ) 01003 RA1 = RA 01004 END IF 01005 END IF 01006 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1 01007 NR = ( J2+KA-1 ) / KA1 01008 J1 = J2 - ( NR-1 )*KA1 01009 IF( UPDATE ) THEN 01010 J2T = MIN( J2, I-2*KA+K-1 ) 01011 ELSE 01012 J2T = J2 01013 END IF 01014 NRT = ( J2T+KA-1 ) / KA1 01015 DO 570 J = J1, J2T, KA1 01016 * 01017 * create nonzero element a(j-1,j+ka) outside the band 01018 * and store it in WORK(j) 01019 * 01020 WORK( J ) = WORK( J )*AB( 1, J+KA-1 ) 01021 AB( 1, J+KA-1 ) = RWORK( J )*AB( 1, J+KA-1 ) 01022 570 CONTINUE 01023 * 01024 * generate rotations in 1st set to annihilate elements which 01025 * have been created outside the band 01026 * 01027 IF( NRT.GT.0 ) 01028 $ CALL ZLARGV( NRT, AB( 1, J1+KA ), INCA, WORK( J1 ), KA1, 01029 $ RWORK( J1 ), KA1 ) 01030 IF( NR.GT.0 ) THEN 01031 * 01032 * apply rotations in 1st set from the left 01033 * 01034 DO 580 L = 1, KA - 1 01035 CALL ZLARTV( NR, AB( KA1-L, J1+L ), INCA, 01036 $ AB( KA-L, J1+L ), INCA, RWORK( J1 ), 01037 $ WORK( J1 ), KA1 ) 01038 580 CONTINUE 01039 * 01040 * apply rotations in 1st set from both sides to diagonal 01041 * blocks 01042 * 01043 CALL ZLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ), 01044 $ AB( KA, J1 ), INCA, RWORK( J1 ), WORK( J1 ), 01045 $ KA1 ) 01046 * 01047 CALL ZLACGV( NR, WORK( J1 ), KA1 ) 01048 END IF 01049 * 01050 * start applying rotations in 1st set from the right 01051 * 01052 DO 590 L = KA - 1, KB - K + 1, -1 01053 NRT = ( J2+L-1 ) / KA1 01054 J1T = J2 - ( NRT-1 )*KA1 01055 IF( NRT.GT.0 ) 01056 $ CALL ZLARTV( NRT, AB( L, J1T ), INCA, 01057 $ AB( L+1, J1T-1 ), INCA, RWORK( J1T ), 01058 $ WORK( J1T ), KA1 ) 01059 590 CONTINUE 01060 * 01061 IF( WANTX ) THEN 01062 * 01063 * post-multiply X by product of rotations in 1st set 01064 * 01065 DO 600 J = J1, J2, KA1 01066 CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1, 01067 $ RWORK( J ), WORK( J ) ) 01068 600 CONTINUE 01069 END IF 01070 610 CONTINUE 01071 * 01072 IF( UPDATE ) THEN 01073 IF( I2.GT.0 .AND. KBT.GT.0 ) THEN 01074 * 01075 * create nonzero element a(i+kbt-ka-1,i+kbt) outside the 01076 * band and store it in WORK(m-kb+i+kbt) 01077 * 01078 WORK( M-KB+I+KBT ) = -BB( KB1-KBT, I+KBT )*RA1 01079 END IF 01080 END IF 01081 * 01082 DO 650 K = KB, 1, -1 01083 IF( UPDATE ) THEN 01084 J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1 01085 ELSE 01086 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1 01087 END IF 01088 * 01089 * finish applying rotations in 2nd set from the right 01090 * 01091 DO 620 L = KB - K, 1, -1 01092 NRT = ( J2+KA+L-1 ) / KA1 01093 J1T = J2 - ( NRT-1 )*KA1 01094 IF( NRT.GT.0 ) 01095 $ CALL ZLARTV( NRT, AB( L, J1T+KA ), INCA, 01096 $ AB( L+1, J1T+KA-1 ), INCA, 01097 $ RWORK( M-KB+J1T+KA ), 01098 $ WORK( M-KB+J1T+KA ), KA1 ) 01099 620 CONTINUE 01100 NR = ( J2+KA-1 ) / KA1 01101 J1 = J2 - ( NR-1 )*KA1 01102 DO 630 J = J1, J2, KA1 01103 WORK( M-KB+J ) = WORK( M-KB+J+KA ) 01104 RWORK( M-KB+J ) = RWORK( M-KB+J+KA ) 01105 630 CONTINUE 01106 DO 640 J = J1, J2, KA1 01107 * 01108 * create nonzero element a(j-1,j+ka) outside the band 01109 * and store it in WORK(m-kb+j) 01110 * 01111 WORK( M-KB+J ) = WORK( M-KB+J )*AB( 1, J+KA-1 ) 01112 AB( 1, J+KA-1 ) = RWORK( M-KB+J )*AB( 1, J+KA-1 ) 01113 640 CONTINUE 01114 IF( UPDATE ) THEN 01115 IF( I+K.GT.KA1 .AND. K.LE.KBT ) 01116 $ WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K ) 01117 END IF 01118 650 CONTINUE 01119 * 01120 DO 690 K = KB, 1, -1 01121 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1 01122 NR = ( J2+KA-1 ) / KA1 01123 J1 = J2 - ( NR-1 )*KA1 01124 IF( NR.GT.0 ) THEN 01125 * 01126 * generate rotations in 2nd set to annihilate elements 01127 * which have been created outside the band 01128 * 01129 CALL ZLARGV( NR, AB( 1, J1+KA ), INCA, WORK( M-KB+J1 ), 01130 $ KA1, RWORK( M-KB+J1 ), KA1 ) 01131 * 01132 * apply rotations in 2nd set from the left 01133 * 01134 DO 660 L = 1, KA - 1 01135 CALL ZLARTV( NR, AB( KA1-L, J1+L ), INCA, 01136 $ AB( KA-L, J1+L ), INCA, RWORK( M-KB+J1 ), 01137 $ WORK( M-KB+J1 ), KA1 ) 01138 660 CONTINUE 01139 * 01140 * apply rotations in 2nd set from both sides to diagonal 01141 * blocks 01142 * 01143 CALL ZLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ), 01144 $ AB( KA, J1 ), INCA, RWORK( M-KB+J1 ), 01145 $ WORK( M-KB+J1 ), KA1 ) 01146 * 01147 CALL ZLACGV( NR, WORK( M-KB+J1 ), KA1 ) 01148 END IF 01149 * 01150 * start applying rotations in 2nd set from the right 01151 * 01152 DO 670 L = KA - 1, KB - K + 1, -1 01153 NRT = ( J2+L-1 ) / KA1 01154 J1T = J2 - ( NRT-1 )*KA1 01155 IF( NRT.GT.0 ) 01156 $ CALL ZLARTV( NRT, AB( L, J1T ), INCA, 01157 $ AB( L+1, J1T-1 ), INCA, 01158 $ RWORK( M-KB+J1T ), WORK( M-KB+J1T ), 01159 $ KA1 ) 01160 670 CONTINUE 01161 * 01162 IF( WANTX ) THEN 01163 * 01164 * post-multiply X by product of rotations in 2nd set 01165 * 01166 DO 680 J = J1, J2, KA1 01167 CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1, 01168 $ RWORK( M-KB+J ), WORK( M-KB+J ) ) 01169 680 CONTINUE 01170 END IF 01171 690 CONTINUE 01172 * 01173 DO 710 K = 1, KB - 1 01174 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1 01175 * 01176 * finish applying rotations in 1st set from the right 01177 * 01178 DO 700 L = KB - K, 1, -1 01179 NRT = ( J2+L-1 ) / KA1 01180 J1T = J2 - ( NRT-1 )*KA1 01181 IF( NRT.GT.0 ) 01182 $ CALL ZLARTV( NRT, AB( L, J1T ), INCA, 01183 $ AB( L+1, J1T-1 ), INCA, RWORK( J1T ), 01184 $ WORK( J1T ), KA1 ) 01185 700 CONTINUE 01186 710 CONTINUE 01187 * 01188 IF( KB.GT.1 ) THEN 01189 DO 720 J = 2, I2 - KA 01190 RWORK( J ) = RWORK( J+KA ) 01191 WORK( J ) = WORK( J+KA ) 01192 720 CONTINUE 01193 END IF 01194 * 01195 ELSE 01196 * 01197 * Transform A, working with the lower triangle 01198 * 01199 IF( UPDATE ) THEN 01200 * 01201 * Form inv(S(i))**H * A * inv(S(i)) 01202 * 01203 BII = DBLE( BB( 1, I ) ) 01204 AB( 1, I ) = ( DBLE( AB( 1, I ) ) / BII ) / BII 01205 DO 730 J = I1, I - 1 01206 AB( I-J+1, J ) = AB( I-J+1, J ) / BII 01207 730 CONTINUE 01208 DO 740 J = I + 1, MIN( N, I+KA ) 01209 AB( J-I+1, I ) = AB( J-I+1, I ) / BII 01210 740 CONTINUE 01211 DO 770 K = I + 1, I + KBT 01212 DO 750 J = K, I + KBT 01213 AB( J-K+1, K ) = AB( J-K+1, K ) - 01214 $ BB( J-I+1, I )*DCONJG( AB( K-I+1, 01215 $ I ) ) - DCONJG( BB( K-I+1, I ) )* 01216 $ AB( J-I+1, I ) + DBLE( AB( 1, I ) )* 01217 $ BB( J-I+1, I )*DCONJG( BB( K-I+1, 01218 $ I ) ) 01219 750 CONTINUE 01220 DO 760 J = I + KBT + 1, MIN( N, I+KA ) 01221 AB( J-K+1, K ) = AB( J-K+1, K ) - 01222 $ DCONJG( BB( K-I+1, I ) )* 01223 $ AB( J-I+1, I ) 01224 760 CONTINUE 01225 770 CONTINUE 01226 DO 790 J = I1, I 01227 DO 780 K = I + 1, MIN( J+KA, I+KBT ) 01228 AB( K-J+1, J ) = AB( K-J+1, J ) - 01229 $ BB( K-I+1, I )*AB( I-J+1, J ) 01230 780 CONTINUE 01231 790 CONTINUE 01232 * 01233 IF( WANTX ) THEN 01234 * 01235 * post-multiply X by inv(S(i)) 01236 * 01237 CALL ZDSCAL( NX, ONE / BII, X( 1, I ), 1 ) 01238 IF( KBT.GT.0 ) 01239 $ CALL ZGERC( NX, KBT, -CONE, X( 1, I ), 1, BB( 2, I ), 01240 $ 1, X( 1, I+1 ), LDX ) 01241 END IF 01242 * 01243 * store a(i,i1) in RA1 for use in next loop over K 01244 * 01245 RA1 = AB( I-I1+1, I1 ) 01246 END IF 01247 * 01248 * Generate and apply vectors of rotations to chase all the 01249 * existing bulges KA positions up toward the top of the band 01250 * 01251 DO 840 K = 1, KB - 1 01252 IF( UPDATE ) THEN 01253 * 01254 * Determine the rotations which would annihilate the bulge 01255 * which has in theory just been created 01256 * 01257 IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN 01258 * 01259 * generate rotation to annihilate a(i,i+k-ka-1) 01260 * 01261 CALL ZLARTG( AB( KA1-K, I+K-KA ), RA1, 01262 $ RWORK( I+K-KA ), WORK( I+K-KA ), RA ) 01263 * 01264 * create nonzero element a(i+k,i+k-ka-1) outside the 01265 * band and store it in WORK(m-kb+i+k) 01266 * 01267 T = -BB( K+1, I )*RA1 01268 WORK( M-KB+I+K ) = RWORK( I+K-KA )*T - 01269 $ DCONJG( WORK( I+K-KA ) )* 01270 $ AB( KA1, I+K-KA ) 01271 AB( KA1, I+K-KA ) = WORK( I+K-KA )*T + 01272 $ RWORK( I+K-KA )*AB( KA1, I+K-KA ) 01273 RA1 = RA 01274 END IF 01275 END IF 01276 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1 01277 NR = ( J2+KA-1 ) / KA1 01278 J1 = J2 - ( NR-1 )*KA1 01279 IF( UPDATE ) THEN 01280 J2T = MIN( J2, I-2*KA+K-1 ) 01281 ELSE 01282 J2T = J2 01283 END IF 01284 NRT = ( J2T+KA-1 ) / KA1 01285 DO 800 J = J1, J2T, KA1 01286 * 01287 * create nonzero element a(j+ka,j-1) outside the band 01288 * and store it in WORK(j) 01289 * 01290 WORK( J ) = WORK( J )*AB( KA1, J-1 ) 01291 AB( KA1, J-1 ) = RWORK( J )*AB( KA1, J-1 ) 01292 800 CONTINUE 01293 * 01294 * generate rotations in 1st set to annihilate elements which 01295 * have been created outside the band 01296 * 01297 IF( NRT.GT.0 ) 01298 $ CALL ZLARGV( NRT, AB( KA1, J1 ), INCA, WORK( J1 ), KA1, 01299 $ RWORK( J1 ), KA1 ) 01300 IF( NR.GT.0 ) THEN 01301 * 01302 * apply rotations in 1st set from the right 01303 * 01304 DO 810 L = 1, KA - 1 01305 CALL ZLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ), 01306 $ INCA, RWORK( J1 ), WORK( J1 ), KA1 ) 01307 810 CONTINUE 01308 * 01309 * apply rotations in 1st set from both sides to diagonal 01310 * blocks 01311 * 01312 CALL ZLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ), 01313 $ AB( 2, J1-1 ), INCA, RWORK( J1 ), 01314 $ WORK( J1 ), KA1 ) 01315 * 01316 CALL ZLACGV( NR, WORK( J1 ), KA1 ) 01317 END IF 01318 * 01319 * start applying rotations in 1st set from the left 01320 * 01321 DO 820 L = KA - 1, KB - K + 1, -1 01322 NRT = ( J2+L-1 ) / KA1 01323 J1T = J2 - ( NRT-1 )*KA1 01324 IF( NRT.GT.0 ) 01325 $ CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA, 01326 $ AB( KA1-L, J1T-KA1+L ), INCA, 01327 $ RWORK( J1T ), WORK( J1T ), KA1 ) 01328 820 CONTINUE 01329 * 01330 IF( WANTX ) THEN 01331 * 01332 * post-multiply X by product of rotations in 1st set 01333 * 01334 DO 830 J = J1, J2, KA1 01335 CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1, 01336 $ RWORK( J ), DCONJG( WORK( J ) ) ) 01337 830 CONTINUE 01338 END IF 01339 840 CONTINUE 01340 * 01341 IF( UPDATE ) THEN 01342 IF( I2.GT.0 .AND. KBT.GT.0 ) THEN 01343 * 01344 * create nonzero element a(i+kbt,i+kbt-ka-1) outside the 01345 * band and store it in WORK(m-kb+i+kbt) 01346 * 01347 WORK( M-KB+I+KBT ) = -BB( KBT+1, I )*RA1 01348 END IF 01349 END IF 01350 * 01351 DO 880 K = KB, 1, -1 01352 IF( UPDATE ) THEN 01353 J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1 01354 ELSE 01355 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1 01356 END IF 01357 * 01358 * finish applying rotations in 2nd set from the left 01359 * 01360 DO 850 L = KB - K, 1, -1 01361 NRT = ( J2+KA+L-1 ) / KA1 01362 J1T = J2 - ( NRT-1 )*KA1 01363 IF( NRT.GT.0 ) 01364 $ CALL ZLARTV( NRT, AB( KA1-L+1, J1T+L-1 ), INCA, 01365 $ AB( KA1-L, J1T+L-1 ), INCA, 01366 $ RWORK( M-KB+J1T+KA ), 01367 $ WORK( M-KB+J1T+KA ), KA1 ) 01368 850 CONTINUE 01369 NR = ( J2+KA-1 ) / KA1 01370 J1 = J2 - ( NR-1 )*KA1 01371 DO 860 J = J1, J2, KA1 01372 WORK( M-KB+J ) = WORK( M-KB+J+KA ) 01373 RWORK( M-KB+J ) = RWORK( M-KB+J+KA ) 01374 860 CONTINUE 01375 DO 870 J = J1, J2, KA1 01376 * 01377 * create nonzero element a(j+ka,j-1) outside the band 01378 * and store it in WORK(m-kb+j) 01379 * 01380 WORK( M-KB+J ) = WORK( M-KB+J )*AB( KA1, J-1 ) 01381 AB( KA1, J-1 ) = RWORK( M-KB+J )*AB( KA1, J-1 ) 01382 870 CONTINUE 01383 IF( UPDATE ) THEN 01384 IF( I+K.GT.KA1 .AND. K.LE.KBT ) 01385 $ WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K ) 01386 END IF 01387 880 CONTINUE 01388 * 01389 DO 920 K = KB, 1, -1 01390 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1 01391 NR = ( J2+KA-1 ) / KA1 01392 J1 = J2 - ( NR-1 )*KA1 01393 IF( NR.GT.0 ) THEN 01394 * 01395 * generate rotations in 2nd set to annihilate elements 01396 * which have been created outside the band 01397 * 01398 CALL ZLARGV( NR, AB( KA1, J1 ), INCA, WORK( M-KB+J1 ), 01399 $ KA1, RWORK( M-KB+J1 ), KA1 ) 01400 * 01401 * apply rotations in 2nd set from the right 01402 * 01403 DO 890 L = 1, KA - 1 01404 CALL ZLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ), 01405 $ INCA, RWORK( M-KB+J1 ), WORK( M-KB+J1 ), 01406 $ KA1 ) 01407 890 CONTINUE 01408 * 01409 * apply rotations in 2nd set from both sides to diagonal 01410 * blocks 01411 * 01412 CALL ZLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ), 01413 $ AB( 2, J1-1 ), INCA, RWORK( M-KB+J1 ), 01414 $ WORK( M-KB+J1 ), KA1 ) 01415 * 01416 CALL ZLACGV( NR, WORK( M-KB+J1 ), KA1 ) 01417 END IF 01418 * 01419 * start applying rotations in 2nd set from the left 01420 * 01421 DO 900 L = KA - 1, KB - K + 1, -1 01422 NRT = ( J2+L-1 ) / KA1 01423 J1T = J2 - ( NRT-1 )*KA1 01424 IF( NRT.GT.0 ) 01425 $ CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA, 01426 $ AB( KA1-L, J1T-KA1+L ), INCA, 01427 $ RWORK( M-KB+J1T ), WORK( M-KB+J1T ), 01428 $ KA1 ) 01429 900 CONTINUE 01430 * 01431 IF( WANTX ) THEN 01432 * 01433 * post-multiply X by product of rotations in 2nd set 01434 * 01435 DO 910 J = J1, J2, KA1 01436 CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1, 01437 $ RWORK( M-KB+J ), DCONJG( WORK( M-KB+J ) ) ) 01438 910 CONTINUE 01439 END IF 01440 920 CONTINUE 01441 * 01442 DO 940 K = 1, KB - 1 01443 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1 01444 * 01445 * finish applying rotations in 1st set from the left 01446 * 01447 DO 930 L = KB - K, 1, -1 01448 NRT = ( J2+L-1 ) / KA1 01449 J1T = J2 - ( NRT-1 )*KA1 01450 IF( NRT.GT.0 ) 01451 $ CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA, 01452 $ AB( KA1-L, J1T-KA1+L ), INCA, 01453 $ RWORK( J1T ), WORK( J1T ), KA1 ) 01454 930 CONTINUE 01455 940 CONTINUE 01456 * 01457 IF( KB.GT.1 ) THEN 01458 DO 950 J = 2, I2 - KA 01459 RWORK( J ) = RWORK( J+KA ) 01460 WORK( J ) = WORK( J+KA ) 01461 950 CONTINUE 01462 END IF 01463 * 01464 END IF 01465 * 01466 GO TO 490 01467 * 01468 * End of ZHBGST 01469 * 01470 END