![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZBBCSD 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZBBCSD + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zbbcsd.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zbbcsd.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zbbcsd.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, 00022 * THETA, PHI, U1, LDU1, U2, LDU2, V1T, LDV1T, 00023 * V2T, LDV2T, B11D, B11E, B12D, B12E, B21D, B21E, 00024 * B22D, B22E, RWORK, LRWORK, INFO ) 00025 * 00026 * .. Scalar Arguments .. 00027 * CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS 00028 * INTEGER INFO, LDU1, LDU2, LDV1T, LDV2T, LRWORK, M, P, Q 00029 * .. 00030 * .. Array Arguments .. 00031 * DOUBLE PRECISION B11D( * ), B11E( * ), B12D( * ), B12E( * ), 00032 * $ B21D( * ), B21E( * ), B22D( * ), B22E( * ), 00033 * $ PHI( * ), THETA( * ), RWORK( * ) 00034 * COMPLEX*16 U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ), 00035 * $ V2T( LDV2T, * ) 00036 * .. 00037 * 00038 * 00039 *> \par Purpose: 00040 * ============= 00041 *> 00042 *> \verbatim 00043 *> 00044 *> ZBBCSD computes the CS decomposition of a unitary matrix in 00045 *> bidiagonal-block form, 00046 *> 00047 *> 00048 *> [ B11 | B12 0 0 ] 00049 *> [ 0 | 0 -I 0 ] 00050 *> X = [----------------] 00051 *> [ B21 | B22 0 0 ] 00052 *> [ 0 | 0 0 I ] 00053 *> 00054 *> [ C | -S 0 0 ] 00055 *> [ U1 | ] [ 0 | 0 -I 0 ] [ V1 | ]**H 00056 *> = [---------] [---------------] [---------] . 00057 *> [ | U2 ] [ S | C 0 0 ] [ | V2 ] 00058 *> [ 0 | 0 0 I ] 00059 *> 00060 *> X is M-by-M, its top-left block is P-by-Q, and Q must be no larger 00061 *> than P, M-P, or M-Q. (If Q is not the smallest index, then X must be 00062 *> transposed and/or permuted. This can be done in constant time using 00063 *> the TRANS and SIGNS options. See ZUNCSD for details.) 00064 *> 00065 *> The bidiagonal matrices B11, B12, B21, and B22 are represented 00066 *> implicitly by angles THETA(1:Q) and PHI(1:Q-1). 00067 *> 00068 *> The unitary matrices U1, U2, V1T, and V2T are input/output. 00069 *> The input matrices are pre- or post-multiplied by the appropriate 00070 *> singular vector matrices. 00071 *> \endverbatim 00072 * 00073 * Arguments: 00074 * ========== 00075 * 00076 *> \param[in] JOBU1 00077 *> \verbatim 00078 *> JOBU1 is CHARACTER 00079 *> = 'Y': U1 is updated; 00080 *> otherwise: U1 is not updated. 00081 *> \endverbatim 00082 *> 00083 *> \param[in] JOBU2 00084 *> \verbatim 00085 *> JOBU2 is CHARACTER 00086 *> = 'Y': U2 is updated; 00087 *> otherwise: U2 is not updated. 00088 *> \endverbatim 00089 *> 00090 *> \param[in] JOBV1T 00091 *> \verbatim 00092 *> JOBV1T is CHARACTER 00093 *> = 'Y': V1T is updated; 00094 *> otherwise: V1T is not updated. 00095 *> \endverbatim 00096 *> 00097 *> \param[in] JOBV2T 00098 *> \verbatim 00099 *> JOBV2T is CHARACTER 00100 *> = 'Y': V2T is updated; 00101 *> otherwise: V2T is not updated. 00102 *> \endverbatim 00103 *> 00104 *> \param[in] TRANS 00105 *> \verbatim 00106 *> TRANS is CHARACTER 00107 *> = 'T': X, U1, U2, V1T, and V2T are stored in row-major 00108 *> order; 00109 *> otherwise: X, U1, U2, V1T, and V2T are stored in column- 00110 *> major order. 00111 *> \endverbatim 00112 *> 00113 *> \param[in] M 00114 *> \verbatim 00115 *> M is INTEGER 00116 *> The number of rows and columns in X, the unitary matrix in 00117 *> bidiagonal-block form. 00118 *> \endverbatim 00119 *> 00120 *> \param[in] P 00121 *> \verbatim 00122 *> P is INTEGER 00123 *> The number of rows in the top-left block of X. 0 <= P <= M. 00124 *> \endverbatim 00125 *> 00126 *> \param[in] Q 00127 *> \verbatim 00128 *> Q is INTEGER 00129 *> The number of columns in the top-left block of X. 00130 *> 0 <= Q <= MIN(P,M-P,M-Q). 00131 *> \endverbatim 00132 *> 00133 *> \param[in,out] THETA 00134 *> \verbatim 00135 *> THETA is DOUBLE PRECISION array, dimension (Q) 00136 *> On entry, the angles THETA(1),...,THETA(Q) that, along with 00137 *> PHI(1), ...,PHI(Q-1), define the matrix in bidiagonal-block 00138 *> form. On exit, the angles whose cosines and sines define the 00139 *> diagonal blocks in the CS decomposition. 00140 *> \endverbatim 00141 *> 00142 *> \param[in,out] PHI 00143 *> \verbatim 00144 *> PHI is DOUBLE PRECISION array, dimension (Q-1) 00145 *> The angles PHI(1),...,PHI(Q-1) that, along with THETA(1),..., 00146 *> THETA(Q), define the matrix in bidiagonal-block form. 00147 *> \endverbatim 00148 *> 00149 *> \param[in,out] U1 00150 *> \verbatim 00151 *> U1 is COMPLEX*16 array, dimension (LDU1,P) 00152 *> On entry, an LDU1-by-P matrix. On exit, U1 is postmultiplied 00153 *> by the left singular vector matrix common to [ B11 ; 0 ] and 00154 *> [ B12 0 0 ; 0 -I 0 0 ]. 00155 *> \endverbatim 00156 *> 00157 *> \param[in] LDU1 00158 *> \verbatim 00159 *> LDU1 is INTEGER 00160 *> The leading dimension of the array U1. 00161 *> \endverbatim 00162 *> 00163 *> \param[in,out] U2 00164 *> \verbatim 00165 *> U2 is COMPLEX*16 array, dimension (LDU2,M-P) 00166 *> On entry, an LDU2-by-(M-P) matrix. On exit, U2 is 00167 *> postmultiplied by the left singular vector matrix common to 00168 *> [ B21 ; 0 ] and [ B22 0 0 ; 0 0 I ]. 00169 *> \endverbatim 00170 *> 00171 *> \param[in] LDU2 00172 *> \verbatim 00173 *> LDU2 is INTEGER 00174 *> The leading dimension of the array U2. 00175 *> \endverbatim 00176 *> 00177 *> \param[in,out] V1T 00178 *> \verbatim 00179 *> V1T is COMPLEX*16 array, dimension (LDV1T,Q) 00180 *> On entry, a LDV1T-by-Q matrix. On exit, V1T is premultiplied 00181 *> by the conjugate transpose of the right singular vector 00182 *> matrix common to [ B11 ; 0 ] and [ B21 ; 0 ]. 00183 *> \endverbatim 00184 *> 00185 *> \param[in] LDV1T 00186 *> \verbatim 00187 *> LDV1T is INTEGER 00188 *> The leading dimension of the array V1T. 00189 *> \endverbatim 00190 *> 00191 *> \param[in,out] V2T 00192 *> \verbatim 00193 *> V2T is COMPLEX*16 array, dimenison (LDV2T,M-Q) 00194 *> On entry, a LDV2T-by-(M-Q) matrix. On exit, V2T is 00195 *> premultiplied by the conjugate transpose of the right 00196 *> singular vector matrix common to [ B12 0 0 ; 0 -I 0 ] and 00197 *> [ B22 0 0 ; 0 0 I ]. 00198 *> \endverbatim 00199 *> 00200 *> \param[in] LDV2T 00201 *> \verbatim 00202 *> LDV2T is INTEGER 00203 *> The leading dimension of the array V2T. 00204 *> \endverbatim 00205 *> 00206 *> \param[out] B11D 00207 *> \verbatim 00208 *> B11D is DOUBLE PRECISION array, dimension (Q) 00209 *> When ZBBCSD converges, B11D contains the cosines of THETA(1), 00210 *> ..., THETA(Q). If ZBBCSD fails to converge, then B11D 00211 *> contains the diagonal of the partially reduced top-left 00212 *> block. 00213 *> \endverbatim 00214 *> 00215 *> \param[out] B11E 00216 *> \verbatim 00217 *> B11E is DOUBLE PRECISION array, dimension (Q-1) 00218 *> When ZBBCSD converges, B11E contains zeros. If ZBBCSD fails 00219 *> to converge, then B11E contains the superdiagonal of the 00220 *> partially reduced top-left block. 00221 *> \endverbatim 00222 *> 00223 *> \param[out] B12D 00224 *> \verbatim 00225 *> B12D is DOUBLE PRECISION array, dimension (Q) 00226 *> When ZBBCSD converges, B12D contains the negative sines of 00227 *> THETA(1), ..., THETA(Q). If ZBBCSD fails to converge, then 00228 *> B12D contains the diagonal of the partially reduced top-right 00229 *> block. 00230 *> \endverbatim 00231 *> 00232 *> \param[out] B12E 00233 *> \verbatim 00234 *> B12E is DOUBLE PRECISION array, dimension (Q-1) 00235 *> When ZBBCSD converges, B12E contains zeros. If ZBBCSD fails 00236 *> to converge, then B12E contains the subdiagonal of the 00237 *> partially reduced top-right block. 00238 *> \endverbatim 00239 *> 00240 *> \param[out] B21D 00241 *> \verbatim 00242 *> B21D is DOUBLE PRECISION array, dimension (Q) 00243 *> When CBBCSD converges, B21D contains the negative sines of 00244 *> THETA(1), ..., THETA(Q). If CBBCSD fails to converge, then 00245 *> B21D contains the diagonal of the partially reduced bottom-left 00246 *> block. 00247 *> \endverbatim 00248 *> 00249 *> \param[out] B21E 00250 *> \verbatim 00251 *> B21E is DOUBLE PRECISION array, dimension (Q-1) 00252 *> When CBBCSD converges, B21E contains zeros. If CBBCSD fails 00253 *> to converge, then B21E contains the subdiagonal of the 00254 *> partially reduced bottom-left block. 00255 *> \endverbatim 00256 *> 00257 *> \param[out] B22D 00258 *> \verbatim 00259 *> B22D is DOUBLE PRECISION array, dimension (Q) 00260 *> When CBBCSD converges, B22D contains the negative sines of 00261 *> THETA(1), ..., THETA(Q). If CBBCSD fails to converge, then 00262 *> B22D contains the diagonal of the partially reduced bottom-right 00263 *> block. 00264 *> \endverbatim 00265 *> 00266 *> \param[out] B22E 00267 *> \verbatim 00268 *> B22E is DOUBLE PRECISION array, dimension (Q-1) 00269 *> When CBBCSD converges, B22E contains zeros. If CBBCSD fails 00270 *> to converge, then B22E contains the subdiagonal of the 00271 *> partially reduced bottom-right block. 00272 *> \endverbatim 00273 *> 00274 *> \param[out] RWORK 00275 *> \verbatim 00276 *> RWORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 00277 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00278 *> \endverbatim 00279 *> 00280 *> \param[in] LRWORK 00281 *> \verbatim 00282 *> LRWORK is INTEGER 00283 *> The dimension of the array RWORK. LRWORK >= MAX(1,8*Q). 00284 *> 00285 *> If LRWORK = -1, then a workspace query is assumed; the 00286 *> routine only calculates the optimal size of the RWORK array, 00287 *> returns this value as the first entry of the work array, and 00288 *> no error message related to LRWORK is issued by XERBLA. 00289 *> \endverbatim 00290 *> 00291 *> \param[out] INFO 00292 *> \verbatim 00293 *> INFO is INTEGER 00294 *> = 0: successful exit. 00295 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00296 *> > 0: if ZBBCSD did not converge, INFO specifies the number 00297 *> of nonzero entries in PHI, and B11D, B11E, etc., 00298 *> contain the partially reduced matrix. 00299 *> \endverbatim 00300 * 00301 *> \par Internal Parameters: 00302 * ========================= 00303 *> 00304 *> \verbatim 00305 *> TOLMUL DOUBLE PRECISION, default = MAX(10,MIN(100,EPS**(-1/8))) 00306 *> TOLMUL controls the convergence criterion of the QR loop. 00307 *> Angles THETA(i), PHI(i) are rounded to 0 or PI/2 when they 00308 *> are within TOLMUL*EPS of either bound. 00309 *> \endverbatim 00310 * 00311 *> \par References: 00312 * ================ 00313 *> 00314 *> [1] Brian D. Sutton. Computing the complete CS decomposition. Numer. 00315 *> Algorithms, 50(1):33-65, 2009. 00316 * 00317 * Authors: 00318 * ======== 00319 * 00320 *> \author Univ. of Tennessee 00321 *> \author Univ. of California Berkeley 00322 *> \author Univ. of Colorado Denver 00323 *> \author NAG Ltd. 00324 * 00325 *> \date November 2011 00326 * 00327 *> \ingroup complex16OTHERcomputational 00328 * 00329 * ===================================================================== 00330 SUBROUTINE ZBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, 00331 $ THETA, PHI, U1, LDU1, U2, LDU2, V1T, LDV1T, 00332 $ V2T, LDV2T, B11D, B11E, B12D, B12E, B21D, B21E, 00333 $ B22D, B22E, RWORK, LRWORK, INFO ) 00334 * 00335 * -- LAPACK computational routine (version 3.4.0) -- 00336 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00337 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00338 * November 2011 00339 * 00340 * .. Scalar Arguments .. 00341 CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS 00342 INTEGER INFO, LDU1, LDU2, LDV1T, LDV2T, LRWORK, M, P, Q 00343 * .. 00344 * .. Array Arguments .. 00345 DOUBLE PRECISION B11D( * ), B11E( * ), B12D( * ), B12E( * ), 00346 $ B21D( * ), B21E( * ), B22D( * ), B22E( * ), 00347 $ PHI( * ), THETA( * ), RWORK( * ) 00348 COMPLEX*16 U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ), 00349 $ V2T( LDV2T, * ) 00350 * .. 00351 * 00352 * =================================================================== 00353 * 00354 * .. Parameters .. 00355 INTEGER MAXITR 00356 PARAMETER ( MAXITR = 6 ) 00357 DOUBLE PRECISION HUNDRED, MEIGHTH, ONE, PIOVER2, TEN, ZERO 00358 PARAMETER ( HUNDRED = 100.0D0, MEIGHTH = -0.125D0, 00359 $ ONE = 1.0D0, PIOVER2 = 1.57079632679489662D0, 00360 $ TEN = 10.0D0, ZERO = 0.0D0 ) 00361 COMPLEX*16 NEGONECOMPLEX 00362 PARAMETER ( NEGONECOMPLEX = (-1.0D0,0.0D0) ) 00363 * .. 00364 * .. Local Scalars .. 00365 LOGICAL COLMAJOR, LQUERY, RESTART11, RESTART12, 00366 $ RESTART21, RESTART22, WANTU1, WANTU2, WANTV1T, 00367 $ WANTV2T 00368 INTEGER I, IMIN, IMAX, ITER, IU1CS, IU1SN, IU2CS, 00369 $ IU2SN, IV1TCS, IV1TSN, IV2TCS, IV2TSN, J, 00370 $ LRWORKMIN, LRWORKOPT, MAXIT, MINI 00371 DOUBLE PRECISION B11BULGE, B12BULGE, B21BULGE, B22BULGE, DUMMY, 00372 $ EPS, MU, NU, R, SIGMA11, SIGMA21, 00373 $ TEMP, THETAMAX, THETAMIN, THRESH, TOL, TOLMUL, 00374 $ UNFL, X1, X2, Y1, Y2 00375 * 00376 EXTERNAL DLARTGP, DLARTGS, DLAS2, XERBLA, ZLASR, ZSCAL, 00377 $ ZSWAP 00378 * .. 00379 * .. External Functions .. 00380 DOUBLE PRECISION DLAMCH 00381 LOGICAL LSAME 00382 EXTERNAL LSAME, DLAMCH 00383 * .. 00384 * .. Intrinsic Functions .. 00385 INTRINSIC ABS, ATAN2, COS, MAX, MIN, SIN, SQRT 00386 * .. 00387 * .. Executable Statements .. 00388 * 00389 * Test input arguments 00390 * 00391 INFO = 0 00392 LQUERY = LRWORK .EQ. -1 00393 WANTU1 = LSAME( JOBU1, 'Y' ) 00394 WANTU2 = LSAME( JOBU2, 'Y' ) 00395 WANTV1T = LSAME( JOBV1T, 'Y' ) 00396 WANTV2T = LSAME( JOBV2T, 'Y' ) 00397 COLMAJOR = .NOT. LSAME( TRANS, 'T' ) 00398 * 00399 IF( M .LT. 0 ) THEN 00400 INFO = -6 00401 ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN 00402 INFO = -7 00403 ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN 00404 INFO = -8 00405 ELSE IF( Q .GT. P .OR. Q .GT. M-P .OR. Q .GT. M-Q ) THEN 00406 INFO = -8 00407 ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN 00408 INFO = -12 00409 ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN 00410 INFO = -14 00411 ELSE IF( WANTV1T .AND. LDV1T .LT. Q ) THEN 00412 INFO = -16 00413 ELSE IF( WANTV2T .AND. LDV2T .LT. M-Q ) THEN 00414 INFO = -18 00415 END IF 00416 * 00417 * Quick return if Q = 0 00418 * 00419 IF( INFO .EQ. 0 .AND. Q .EQ. 0 ) THEN 00420 LRWORKMIN = 1 00421 RWORK(1) = LRWORKMIN 00422 RETURN 00423 END IF 00424 * 00425 * Compute workspace 00426 * 00427 IF( INFO .EQ. 0 ) THEN 00428 IU1CS = 1 00429 IU1SN = IU1CS + Q 00430 IU2CS = IU1SN + Q 00431 IU2SN = IU2CS + Q 00432 IV1TCS = IU2SN + Q 00433 IV1TSN = IV1TCS + Q 00434 IV2TCS = IV1TSN + Q 00435 IV2TSN = IV2TCS + Q 00436 LRWORKOPT = IV2TSN + Q - 1 00437 LRWORKMIN = LRWORKOPT 00438 RWORK(1) = LRWORKOPT 00439 IF( LRWORK .LT. LRWORKMIN .AND. .NOT. LQUERY ) THEN 00440 INFO = -28 00441 END IF 00442 END IF 00443 * 00444 IF( INFO .NE. 0 ) THEN 00445 CALL XERBLA( 'ZBBCSD', -INFO ) 00446 RETURN 00447 ELSE IF( LQUERY ) THEN 00448 RETURN 00449 END IF 00450 * 00451 * Get machine constants 00452 * 00453 EPS = DLAMCH( 'Epsilon' ) 00454 UNFL = DLAMCH( 'Safe minimum' ) 00455 TOLMUL = MAX( TEN, MIN( HUNDRED, EPS**MEIGHTH ) ) 00456 TOL = TOLMUL*EPS 00457 THRESH = MAX( TOL, MAXITR*Q*Q*UNFL ) 00458 * 00459 * Test for negligible sines or cosines 00460 * 00461 DO I = 1, Q 00462 IF( THETA(I) .LT. THRESH ) THEN 00463 THETA(I) = ZERO 00464 ELSE IF( THETA(I) .GT. PIOVER2-THRESH ) THEN 00465 THETA(I) = PIOVER2 00466 END IF 00467 END DO 00468 DO I = 1, Q-1 00469 IF( PHI(I) .LT. THRESH ) THEN 00470 PHI(I) = ZERO 00471 ELSE IF( PHI(I) .GT. PIOVER2-THRESH ) THEN 00472 PHI(I) = PIOVER2 00473 END IF 00474 END DO 00475 * 00476 * Initial deflation 00477 * 00478 IMAX = Q 00479 DO WHILE( ( IMAX .GT. 1 ) .AND. ( PHI(IMAX-1) .EQ. ZERO ) ) 00480 IMAX = IMAX - 1 00481 END DO 00482 IMIN = IMAX - 1 00483 IF ( IMIN .GT. 1 ) THEN 00484 DO WHILE( PHI(IMIN-1) .NE. ZERO ) 00485 IMIN = IMIN - 1 00486 IF ( IMIN .LE. 1 ) EXIT 00487 END DO 00488 END IF 00489 * 00490 * Initialize iteration counter 00491 * 00492 MAXIT = MAXITR*Q*Q 00493 ITER = 0 00494 * 00495 * Begin main iteration loop 00496 * 00497 DO WHILE( IMAX .GT. 1 ) 00498 * 00499 * Compute the matrix entries 00500 * 00501 B11D(IMIN) = COS( THETA(IMIN) ) 00502 B21D(IMIN) = -SIN( THETA(IMIN) ) 00503 DO I = IMIN, IMAX - 1 00504 B11E(I) = -SIN( THETA(I) ) * SIN( PHI(I) ) 00505 B11D(I+1) = COS( THETA(I+1) ) * COS( PHI(I) ) 00506 B12D(I) = SIN( THETA(I) ) * COS( PHI(I) ) 00507 B12E(I) = COS( THETA(I+1) ) * SIN( PHI(I) ) 00508 B21E(I) = -COS( THETA(I) ) * SIN( PHI(I) ) 00509 B21D(I+1) = -SIN( THETA(I+1) ) * COS( PHI(I) ) 00510 B22D(I) = COS( THETA(I) ) * COS( PHI(I) ) 00511 B22E(I) = -SIN( THETA(I+1) ) * SIN( PHI(I) ) 00512 END DO 00513 B12D(IMAX) = SIN( THETA(IMAX) ) 00514 B22D(IMAX) = COS( THETA(IMAX) ) 00515 * 00516 * Abort if not converging; otherwise, increment ITER 00517 * 00518 IF( ITER .GT. MAXIT ) THEN 00519 INFO = 0 00520 DO I = 1, Q 00521 IF( PHI(I) .NE. ZERO ) 00522 $ INFO = INFO + 1 00523 END DO 00524 RETURN 00525 END IF 00526 * 00527 ITER = ITER + IMAX - IMIN 00528 * 00529 * Compute shifts 00530 * 00531 THETAMAX = THETA(IMIN) 00532 THETAMIN = THETA(IMIN) 00533 DO I = IMIN+1, IMAX 00534 IF( THETA(I) > THETAMAX ) 00535 $ THETAMAX = THETA(I) 00536 IF( THETA(I) < THETAMIN ) 00537 $ THETAMIN = THETA(I) 00538 END DO 00539 * 00540 IF( THETAMAX .GT. PIOVER2 - THRESH ) THEN 00541 * 00542 * Zero on diagonals of B11 and B22; induce deflation with a 00543 * zero shift 00544 * 00545 MU = ZERO 00546 NU = ONE 00547 * 00548 ELSE IF( THETAMIN .LT. THRESH ) THEN 00549 * 00550 * Zero on diagonals of B12 and B22; induce deflation with a 00551 * zero shift 00552 * 00553 MU = ONE 00554 NU = ZERO 00555 * 00556 ELSE 00557 * 00558 * Compute shifts for B11 and B21 and use the lesser 00559 * 00560 CALL DLAS2( B11D(IMAX-1), B11E(IMAX-1), B11D(IMAX), SIGMA11, 00561 $ DUMMY ) 00562 CALL DLAS2( B21D(IMAX-1), B21E(IMAX-1), B21D(IMAX), SIGMA21, 00563 $ DUMMY ) 00564 * 00565 IF( SIGMA11 .LE. SIGMA21 ) THEN 00566 MU = SIGMA11 00567 NU = SQRT( ONE - MU**2 ) 00568 IF( MU .LT. THRESH ) THEN 00569 MU = ZERO 00570 NU = ONE 00571 END IF 00572 ELSE 00573 NU = SIGMA21 00574 MU = SQRT( 1.0 - NU**2 ) 00575 IF( NU .LT. THRESH ) THEN 00576 MU = ONE 00577 NU = ZERO 00578 END IF 00579 END IF 00580 END IF 00581 * 00582 * Rotate to produce bulges in B11 and B21 00583 * 00584 IF( MU .LE. NU ) THEN 00585 CALL DLARTGS( B11D(IMIN), B11E(IMIN), MU, 00586 $ RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1) ) 00587 ELSE 00588 CALL DLARTGS( B21D(IMIN), B21E(IMIN), NU, 00589 $ RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1) ) 00590 END IF 00591 * 00592 TEMP = RWORK(IV1TCS+IMIN-1)*B11D(IMIN) + 00593 $ RWORK(IV1TSN+IMIN-1)*B11E(IMIN) 00594 B11E(IMIN) = RWORK(IV1TCS+IMIN-1)*B11E(IMIN) - 00595 $ RWORK(IV1TSN+IMIN-1)*B11D(IMIN) 00596 B11D(IMIN) = TEMP 00597 B11BULGE = RWORK(IV1TSN+IMIN-1)*B11D(IMIN+1) 00598 B11D(IMIN+1) = RWORK(IV1TCS+IMIN-1)*B11D(IMIN+1) 00599 TEMP = RWORK(IV1TCS+IMIN-1)*B21D(IMIN) + 00600 $ RWORK(IV1TSN+IMIN-1)*B21E(IMIN) 00601 B21E(IMIN) = RWORK(IV1TCS+IMIN-1)*B21E(IMIN) - 00602 $ RWORK(IV1TSN+IMIN-1)*B21D(IMIN) 00603 B21D(IMIN) = TEMP 00604 B21BULGE = RWORK(IV1TSN+IMIN-1)*B21D(IMIN+1) 00605 B21D(IMIN+1) = RWORK(IV1TCS+IMIN-1)*B21D(IMIN+1) 00606 * 00607 * Compute THETA(IMIN) 00608 * 00609 THETA( IMIN ) = ATAN2( SQRT( B21D(IMIN)**2+B21BULGE**2 ), 00610 $ SQRT( B11D(IMIN)**2+B11BULGE**2 ) ) 00611 * 00612 * Chase the bulges in B11(IMIN+1,IMIN) and B21(IMIN+1,IMIN) 00613 * 00614 IF( B11D(IMIN)**2+B11BULGE**2 .GT. THRESH**2 ) THEN 00615 CALL DLARTGP( B11BULGE, B11D(IMIN), RWORK(IU1SN+IMIN-1), 00616 $ RWORK(IU1CS+IMIN-1), R ) 00617 ELSE IF( MU .LE. NU ) THEN 00618 CALL DLARTGS( B11E( IMIN ), B11D( IMIN + 1 ), MU, 00619 $ RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1) ) 00620 ELSE 00621 CALL DLARTGS( B12D( IMIN ), B12E( IMIN ), NU, 00622 $ RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1) ) 00623 END IF 00624 IF( B21D(IMIN)**2+B21BULGE**2 .GT. THRESH**2 ) THEN 00625 CALL DLARTGP( B21BULGE, B21D(IMIN), RWORK(IU2SN+IMIN-1), 00626 $ RWORK(IU2CS+IMIN-1), R ) 00627 ELSE IF( NU .LT. MU ) THEN 00628 CALL DLARTGS( B21E( IMIN ), B21D( IMIN + 1 ), NU, 00629 $ RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1) ) 00630 ELSE 00631 CALL DLARTGS( B22D(IMIN), B22E(IMIN), MU, 00632 $ RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1) ) 00633 END IF 00634 RWORK(IU2CS+IMIN-1) = -RWORK(IU2CS+IMIN-1) 00635 RWORK(IU2SN+IMIN-1) = -RWORK(IU2SN+IMIN-1) 00636 * 00637 TEMP = RWORK(IU1CS+IMIN-1)*B11E(IMIN) + 00638 $ RWORK(IU1SN+IMIN-1)*B11D(IMIN+1) 00639 B11D(IMIN+1) = RWORK(IU1CS+IMIN-1)*B11D(IMIN+1) - 00640 $ RWORK(IU1SN+IMIN-1)*B11E(IMIN) 00641 B11E(IMIN) = TEMP 00642 IF( IMAX .GT. IMIN+1 ) THEN 00643 B11BULGE = RWORK(IU1SN+IMIN-1)*B11E(IMIN+1) 00644 B11E(IMIN+1) = RWORK(IU1CS+IMIN-1)*B11E(IMIN+1) 00645 END IF 00646 TEMP = RWORK(IU1CS+IMIN-1)*B12D(IMIN) + 00647 $ RWORK(IU1SN+IMIN-1)*B12E(IMIN) 00648 B12E(IMIN) = RWORK(IU1CS+IMIN-1)*B12E(IMIN) - 00649 $ RWORK(IU1SN+IMIN-1)*B12D(IMIN) 00650 B12D(IMIN) = TEMP 00651 B12BULGE = RWORK(IU1SN+IMIN-1)*B12D(IMIN+1) 00652 B12D(IMIN+1) = RWORK(IU1CS+IMIN-1)*B12D(IMIN+1) 00653 TEMP = RWORK(IU2CS+IMIN-1)*B21E(IMIN) + 00654 $ RWORK(IU2SN+IMIN-1)*B21D(IMIN+1) 00655 B21D(IMIN+1) = RWORK(IU2CS+IMIN-1)*B21D(IMIN+1) - 00656 $ RWORK(IU2SN+IMIN-1)*B21E(IMIN) 00657 B21E(IMIN) = TEMP 00658 IF( IMAX .GT. IMIN+1 ) THEN 00659 B21BULGE = RWORK(IU2SN+IMIN-1)*B21E(IMIN+1) 00660 B21E(IMIN+1) = RWORK(IU2CS+IMIN-1)*B21E(IMIN+1) 00661 END IF 00662 TEMP = RWORK(IU2CS+IMIN-1)*B22D(IMIN) + 00663 $ RWORK(IU2SN+IMIN-1)*B22E(IMIN) 00664 B22E(IMIN) = RWORK(IU2CS+IMIN-1)*B22E(IMIN) - 00665 $ RWORK(IU2SN+IMIN-1)*B22D(IMIN) 00666 B22D(IMIN) = TEMP 00667 B22BULGE = RWORK(IU2SN+IMIN-1)*B22D(IMIN+1) 00668 B22D(IMIN+1) = RWORK(IU2CS+IMIN-1)*B22D(IMIN+1) 00669 * 00670 * Inner loop: chase bulges from B11(IMIN,IMIN+2), 00671 * B12(IMIN,IMIN+1), B21(IMIN,IMIN+2), and B22(IMIN,IMIN+1) to 00672 * bottom-right 00673 * 00674 DO I = IMIN+1, IMAX-1 00675 * 00676 * Compute PHI(I-1) 00677 * 00678 X1 = SIN(THETA(I-1))*B11E(I-1) + COS(THETA(I-1))*B21E(I-1) 00679 X2 = SIN(THETA(I-1))*B11BULGE + COS(THETA(I-1))*B21BULGE 00680 Y1 = SIN(THETA(I-1))*B12D(I-1) + COS(THETA(I-1))*B22D(I-1) 00681 Y2 = SIN(THETA(I-1))*B12BULGE + COS(THETA(I-1))*B22BULGE 00682 * 00683 PHI(I-1) = ATAN2( SQRT(X1**2+X2**2), SQRT(Y1**2+Y2**2) ) 00684 * 00685 * Determine if there are bulges to chase or if a new direct 00686 * summand has been reached 00687 * 00688 RESTART11 = B11E(I-1)**2 + B11BULGE**2 .LE. THRESH**2 00689 RESTART21 = B21E(I-1)**2 + B21BULGE**2 .LE. THRESH**2 00690 RESTART12 = B12D(I-1)**2 + B12BULGE**2 .LE. THRESH**2 00691 RESTART22 = B22D(I-1)**2 + B22BULGE**2 .LE. THRESH**2 00692 * 00693 * If possible, chase bulges from B11(I-1,I+1), B12(I-1,I), 00694 * B21(I-1,I+1), and B22(I-1,I). If necessary, restart bulge- 00695 * chasing by applying the original shift again. 00696 * 00697 IF( .NOT. RESTART11 .AND. .NOT. RESTART21 ) THEN 00698 CALL DLARTGP( X2, X1, RWORK(IV1TSN+I-1), 00699 $ RWORK(IV1TCS+I-1), R ) 00700 ELSE IF( .NOT. RESTART11 .AND. RESTART21 ) THEN 00701 CALL DLARTGP( B11BULGE, B11E(I-1), RWORK(IV1TSN+I-1), 00702 $ RWORK(IV1TCS+I-1), R ) 00703 ELSE IF( RESTART11 .AND. .NOT. RESTART21 ) THEN 00704 CALL DLARTGP( B21BULGE, B21E(I-1), RWORK(IV1TSN+I-1), 00705 $ RWORK(IV1TCS+I-1), R ) 00706 ELSE IF( MU .LE. NU ) THEN 00707 CALL DLARTGS( B11D(I), B11E(I), MU, RWORK(IV1TCS+I-1), 00708 $ RWORK(IV1TSN+I-1) ) 00709 ELSE 00710 CALL DLARTGS( B21D(I), B21E(I), NU, RWORK(IV1TCS+I-1), 00711 $ RWORK(IV1TSN+I-1) ) 00712 END IF 00713 RWORK(IV1TCS+I-1) = -RWORK(IV1TCS+I-1) 00714 RWORK(IV1TSN+I-1) = -RWORK(IV1TSN+I-1) 00715 IF( .NOT. RESTART12 .AND. .NOT. RESTART22 ) THEN 00716 CALL DLARTGP( Y2, Y1, RWORK(IV2TSN+I-1-1), 00717 $ RWORK(IV2TCS+I-1-1), R ) 00718 ELSE IF( .NOT. RESTART12 .AND. RESTART22 ) THEN 00719 CALL DLARTGP( B12BULGE, B12D(I-1), RWORK(IV2TSN+I-1-1), 00720 $ RWORK(IV2TCS+I-1-1), R ) 00721 ELSE IF( RESTART12 .AND. .NOT. RESTART22 ) THEN 00722 CALL DLARTGP( B22BULGE, B22D(I-1), RWORK(IV2TSN+I-1-1), 00723 $ RWORK(IV2TCS+I-1-1), R ) 00724 ELSE IF( NU .LT. MU ) THEN 00725 CALL DLARTGS( B12E(I-1), B12D(I), NU, 00726 $ RWORK(IV2TCS+I-1-1), RWORK(IV2TSN+I-1-1) ) 00727 ELSE 00728 CALL DLARTGS( B22E(I-1), B22D(I), MU, 00729 $ RWORK(IV2TCS+I-1-1), RWORK(IV2TSN+I-1-1) ) 00730 END IF 00731 * 00732 TEMP = RWORK(IV1TCS+I-1)*B11D(I) + RWORK(IV1TSN+I-1)*B11E(I) 00733 B11E(I) = RWORK(IV1TCS+I-1)*B11E(I) - 00734 $ RWORK(IV1TSN+I-1)*B11D(I) 00735 B11D(I) = TEMP 00736 B11BULGE = RWORK(IV1TSN+I-1)*B11D(I+1) 00737 B11D(I+1) = RWORK(IV1TCS+I-1)*B11D(I+1) 00738 TEMP = RWORK(IV1TCS+I-1)*B21D(I) + RWORK(IV1TSN+I-1)*B21E(I) 00739 B21E(I) = RWORK(IV1TCS+I-1)*B21E(I) - 00740 $ RWORK(IV1TSN+I-1)*B21D(I) 00741 B21D(I) = TEMP 00742 B21BULGE = RWORK(IV1TSN+I-1)*B21D(I+1) 00743 B21D(I+1) = RWORK(IV1TCS+I-1)*B21D(I+1) 00744 TEMP = RWORK(IV2TCS+I-1-1)*B12E(I-1) + 00745 $ RWORK(IV2TSN+I-1-1)*B12D(I) 00746 B12D(I) = RWORK(IV2TCS+I-1-1)*B12D(I) - 00747 $ RWORK(IV2TSN+I-1-1)*B12E(I-1) 00748 B12E(I-1) = TEMP 00749 B12BULGE = RWORK(IV2TSN+I-1-1)*B12E(I) 00750 B12E(I) = RWORK(IV2TCS+I-1-1)*B12E(I) 00751 TEMP = RWORK(IV2TCS+I-1-1)*B22E(I-1) + 00752 $ RWORK(IV2TSN+I-1-1)*B22D(I) 00753 B22D(I) = RWORK(IV2TCS+I-1-1)*B22D(I) - 00754 $ RWORK(IV2TSN+I-1-1)*B22E(I-1) 00755 B22E(I-1) = TEMP 00756 B22BULGE = RWORK(IV2TSN+I-1-1)*B22E(I) 00757 B22E(I) = RWORK(IV2TCS+I-1-1)*B22E(I) 00758 * 00759 * Compute THETA(I) 00760 * 00761 X1 = COS(PHI(I-1))*B11D(I) + SIN(PHI(I-1))*B12E(I-1) 00762 X2 = COS(PHI(I-1))*B11BULGE + SIN(PHI(I-1))*B12BULGE 00763 Y1 = COS(PHI(I-1))*B21D(I) + SIN(PHI(I-1))*B22E(I-1) 00764 Y2 = COS(PHI(I-1))*B21BULGE + SIN(PHI(I-1))*B22BULGE 00765 * 00766 THETA(I) = ATAN2( SQRT(Y1**2+Y2**2), SQRT(X1**2+X2**2) ) 00767 * 00768 * Determine if there are bulges to chase or if a new direct 00769 * summand has been reached 00770 * 00771 RESTART11 = B11D(I)**2 + B11BULGE**2 .LE. THRESH**2 00772 RESTART12 = B12E(I-1)**2 + B12BULGE**2 .LE. THRESH**2 00773 RESTART21 = B21D(I)**2 + B21BULGE**2 .LE. THRESH**2 00774 RESTART22 = B22E(I-1)**2 + B22BULGE**2 .LE. THRESH**2 00775 * 00776 * If possible, chase bulges from B11(I+1,I), B12(I+1,I-1), 00777 * B21(I+1,I), and B22(I+1,I-1). If necessary, restart bulge- 00778 * chasing by applying the original shift again. 00779 * 00780 IF( .NOT. RESTART11 .AND. .NOT. RESTART12 ) THEN 00781 CALL DLARTGP( X2, X1, RWORK(IU1SN+I-1), RWORK(IU1CS+I-1), 00782 $ R ) 00783 ELSE IF( .NOT. RESTART11 .AND. RESTART12 ) THEN 00784 CALL DLARTGP( B11BULGE, B11D(I), RWORK(IU1SN+I-1), 00785 $ RWORK(IU1CS+I-1), R ) 00786 ELSE IF( RESTART11 .AND. .NOT. RESTART12 ) THEN 00787 CALL DLARTGP( B12BULGE, B12E(I-1), RWORK(IU1SN+I-1), 00788 $ RWORK(IU1CS+I-1), R ) 00789 ELSE IF( MU .LE. NU ) THEN 00790 CALL DLARTGS( B11E(I), B11D(I+1), MU, RWORK(IU1CS+I-1), 00791 $ RWORK(IU1SN+I-1) ) 00792 ELSE 00793 CALL DLARTGS( B12D(I), B12E(I), NU, RWORK(IU1CS+I-1), 00794 $ RWORK(IU1SN+I-1) ) 00795 END IF 00796 IF( .NOT. RESTART21 .AND. .NOT. RESTART22 ) THEN 00797 CALL DLARTGP( Y2, Y1, RWORK(IU2SN+I-1), RWORK(IU2CS+I-1), 00798 $ R ) 00799 ELSE IF( .NOT. RESTART21 .AND. RESTART22 ) THEN 00800 CALL DLARTGP( B21BULGE, B21D(I), RWORK(IU2SN+I-1), 00801 $ RWORK(IU2CS+I-1), R ) 00802 ELSE IF( RESTART21 .AND. .NOT. RESTART22 ) THEN 00803 CALL DLARTGP( B22BULGE, B22E(I-1), RWORK(IU2SN+I-1), 00804 $ RWORK(IU2CS+I-1), R ) 00805 ELSE IF( NU .LT. MU ) THEN 00806 CALL DLARTGS( B21E(I), B21E(I+1), NU, RWORK(IU2CS+I-1), 00807 $ RWORK(IU2SN+I-1) ) 00808 ELSE 00809 CALL DLARTGS( B22D(I), B22E(I), MU, RWORK(IU2CS+I-1), 00810 $ RWORK(IU2SN+I-1) ) 00811 END IF 00812 RWORK(IU2CS+I-1) = -RWORK(IU2CS+I-1) 00813 RWORK(IU2SN+I-1) = -RWORK(IU2SN+I-1) 00814 * 00815 TEMP = RWORK(IU1CS+I-1)*B11E(I) + RWORK(IU1SN+I-1)*B11D(I+1) 00816 B11D(I+1) = RWORK(IU1CS+I-1)*B11D(I+1) - 00817 $ RWORK(IU1SN+I-1)*B11E(I) 00818 B11E(I) = TEMP 00819 IF( I .LT. IMAX - 1 ) THEN 00820 B11BULGE = RWORK(IU1SN+I-1)*B11E(I+1) 00821 B11E(I+1) = RWORK(IU1CS+I-1)*B11E(I+1) 00822 END IF 00823 TEMP = RWORK(IU2CS+I-1)*B21E(I) + RWORK(IU2SN+I-1)*B21D(I+1) 00824 B21D(I+1) = RWORK(IU2CS+I-1)*B21D(I+1) - 00825 $ RWORK(IU2SN+I-1)*B21E(I) 00826 B21E(I) = TEMP 00827 IF( I .LT. IMAX - 1 ) THEN 00828 B21BULGE = RWORK(IU2SN+I-1)*B21E(I+1) 00829 B21E(I+1) = RWORK(IU2CS+I-1)*B21E(I+1) 00830 END IF 00831 TEMP = RWORK(IU1CS+I-1)*B12D(I) + RWORK(IU1SN+I-1)*B12E(I) 00832 B12E(I) = RWORK(IU1CS+I-1)*B12E(I) - 00833 $ RWORK(IU1SN+I-1)*B12D(I) 00834 B12D(I) = TEMP 00835 B12BULGE = RWORK(IU1SN+I-1)*B12D(I+1) 00836 B12D(I+1) = RWORK(IU1CS+I-1)*B12D(I+1) 00837 TEMP = RWORK(IU2CS+I-1)*B22D(I) + RWORK(IU2SN+I-1)*B22E(I) 00838 B22E(I) = RWORK(IU2CS+I-1)*B22E(I) - 00839 $ RWORK(IU2SN+I-1)*B22D(I) 00840 B22D(I) = TEMP 00841 B22BULGE = RWORK(IU2SN+I-1)*B22D(I+1) 00842 B22D(I+1) = RWORK(IU2CS+I-1)*B22D(I+1) 00843 * 00844 END DO 00845 * 00846 * Compute PHI(IMAX-1) 00847 * 00848 X1 = SIN(THETA(IMAX-1))*B11E(IMAX-1) + 00849 $ COS(THETA(IMAX-1))*B21E(IMAX-1) 00850 Y1 = SIN(THETA(IMAX-1))*B12D(IMAX-1) + 00851 $ COS(THETA(IMAX-1))*B22D(IMAX-1) 00852 Y2 = SIN(THETA(IMAX-1))*B12BULGE + COS(THETA(IMAX-1))*B22BULGE 00853 * 00854 PHI(IMAX-1) = ATAN2( ABS(X1), SQRT(Y1**2+Y2**2) ) 00855 * 00856 * Chase bulges from B12(IMAX-1,IMAX) and B22(IMAX-1,IMAX) 00857 * 00858 RESTART12 = B12D(IMAX-1)**2 + B12BULGE**2 .LE. THRESH**2 00859 RESTART22 = B22D(IMAX-1)**2 + B22BULGE**2 .LE. THRESH**2 00860 * 00861 IF( .NOT. RESTART12 .AND. .NOT. RESTART22 ) THEN 00862 CALL DLARTGP( Y2, Y1, RWORK(IV2TSN+IMAX-1-1), 00863 $ RWORK(IV2TCS+IMAX-1-1), R ) 00864 ELSE IF( .NOT. RESTART12 .AND. RESTART22 ) THEN 00865 CALL DLARTGP( B12BULGE, B12D(IMAX-1), 00866 $ RWORK(IV2TSN+IMAX-1-1), 00867 $ RWORK(IV2TCS+IMAX-1-1), R ) 00868 ELSE IF( RESTART12 .AND. .NOT. RESTART22 ) THEN 00869 CALL DLARTGP( B22BULGE, B22D(IMAX-1), 00870 $ RWORK(IV2TSN+IMAX-1-1), 00871 $ RWORK(IV2TCS+IMAX-1-1), R ) 00872 ELSE IF( NU .LT. MU ) THEN 00873 CALL DLARTGS( B12E(IMAX-1), B12D(IMAX), NU, 00874 $ RWORK(IV2TCS+IMAX-1-1), 00875 $ RWORK(IV2TSN+IMAX-1-1) ) 00876 ELSE 00877 CALL DLARTGS( B22E(IMAX-1), B22D(IMAX), MU, 00878 $ RWORK(IV2TCS+IMAX-1-1), 00879 $ RWORK(IV2TSN+IMAX-1-1) ) 00880 END IF 00881 * 00882 TEMP = RWORK(IV2TCS+IMAX-1-1)*B12E(IMAX-1) + 00883 $ RWORK(IV2TSN+IMAX-1-1)*B12D(IMAX) 00884 B12D(IMAX) = RWORK(IV2TCS+IMAX-1-1)*B12D(IMAX) - 00885 $ RWORK(IV2TSN+IMAX-1-1)*B12E(IMAX-1) 00886 B12E(IMAX-1) = TEMP 00887 TEMP = RWORK(IV2TCS+IMAX-1-1)*B22E(IMAX-1) + 00888 $ RWORK(IV2TSN+IMAX-1-1)*B22D(IMAX) 00889 B22D(IMAX) = RWORK(IV2TCS+IMAX-1-1)*B22D(IMAX) - 00890 $ RWORK(IV2TSN+IMAX-1-1)*B22E(IMAX-1) 00891 B22E(IMAX-1) = TEMP 00892 * 00893 * Update singular vectors 00894 * 00895 IF( WANTU1 ) THEN 00896 IF( COLMAJOR ) THEN 00897 CALL ZLASR( 'R', 'V', 'F', P, IMAX-IMIN+1, 00898 $ RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1), 00899 $ U1(1,IMIN), LDU1 ) 00900 ELSE 00901 CALL ZLASR( 'L', 'V', 'F', IMAX-IMIN+1, P, 00902 $ RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1), 00903 $ U1(IMIN,1), LDU1 ) 00904 END IF 00905 END IF 00906 IF( WANTU2 ) THEN 00907 IF( COLMAJOR ) THEN 00908 CALL ZLASR( 'R', 'V', 'F', M-P, IMAX-IMIN+1, 00909 $ RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1), 00910 $ U2(1,IMIN), LDU2 ) 00911 ELSE 00912 CALL ZLASR( 'L', 'V', 'F', IMAX-IMIN+1, M-P, 00913 $ RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1), 00914 $ U2(IMIN,1), LDU2 ) 00915 END IF 00916 END IF 00917 IF( WANTV1T ) THEN 00918 IF( COLMAJOR ) THEN 00919 CALL ZLASR( 'L', 'V', 'F', IMAX-IMIN+1, Q, 00920 $ RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1), 00921 $ V1T(IMIN,1), LDV1T ) 00922 ELSE 00923 CALL ZLASR( 'R', 'V', 'F', Q, IMAX-IMIN+1, 00924 $ RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1), 00925 $ V1T(1,IMIN), LDV1T ) 00926 END IF 00927 END IF 00928 IF( WANTV2T ) THEN 00929 IF( COLMAJOR ) THEN 00930 CALL ZLASR( 'L', 'V', 'F', IMAX-IMIN+1, M-Q, 00931 $ RWORK(IV2TCS+IMIN-1), RWORK(IV2TSN+IMIN-1), 00932 $ V2T(IMIN,1), LDV2T ) 00933 ELSE 00934 CALL ZLASR( 'R', 'V', 'F', M-Q, IMAX-IMIN+1, 00935 $ RWORK(IV2TCS+IMIN-1), RWORK(IV2TSN+IMIN-1), 00936 $ V2T(1,IMIN), LDV2T ) 00937 END IF 00938 END IF 00939 * 00940 * Fix signs on B11(IMAX-1,IMAX) and B21(IMAX-1,IMAX) 00941 * 00942 IF( B11E(IMAX-1)+B21E(IMAX-1) .GT. 0 ) THEN 00943 B11D(IMAX) = -B11D(IMAX) 00944 B21D(IMAX) = -B21D(IMAX) 00945 IF( WANTV1T ) THEN 00946 IF( COLMAJOR ) THEN 00947 CALL ZSCAL( Q, NEGONECOMPLEX, V1T(IMAX,1), LDV1T ) 00948 ELSE 00949 CALL ZSCAL( Q, NEGONECOMPLEX, V1T(1,IMAX), 1 ) 00950 END IF 00951 END IF 00952 END IF 00953 * 00954 * Compute THETA(IMAX) 00955 * 00956 X1 = COS(PHI(IMAX-1))*B11D(IMAX) + 00957 $ SIN(PHI(IMAX-1))*B12E(IMAX-1) 00958 Y1 = COS(PHI(IMAX-1))*B21D(IMAX) + 00959 $ SIN(PHI(IMAX-1))*B22E(IMAX-1) 00960 * 00961 THETA(IMAX) = ATAN2( ABS(Y1), ABS(X1) ) 00962 * 00963 * Fix signs on B11(IMAX,IMAX), B12(IMAX,IMAX-1), B21(IMAX,IMAX), 00964 * and B22(IMAX,IMAX-1) 00965 * 00966 IF( B11D(IMAX)+B12E(IMAX-1) .LT. 0 ) THEN 00967 B12D(IMAX) = -B12D(IMAX) 00968 IF( WANTU1 ) THEN 00969 IF( COLMAJOR ) THEN 00970 CALL ZSCAL( P, NEGONECOMPLEX, U1(1,IMAX), 1 ) 00971 ELSE 00972 CALL ZSCAL( P, NEGONECOMPLEX, U1(IMAX,1), LDU1 ) 00973 END IF 00974 END IF 00975 END IF 00976 IF( B21D(IMAX)+B22E(IMAX-1) .GT. 0 ) THEN 00977 B22D(IMAX) = -B22D(IMAX) 00978 IF( WANTU2 ) THEN 00979 IF( COLMAJOR ) THEN 00980 CALL ZSCAL( M-P, NEGONECOMPLEX, U2(1,IMAX), 1 ) 00981 ELSE 00982 CALL ZSCAL( M-P, NEGONECOMPLEX, U2(IMAX,1), LDU2 ) 00983 END IF 00984 END IF 00985 END IF 00986 * 00987 * Fix signs on B12(IMAX,IMAX) and B22(IMAX,IMAX) 00988 * 00989 IF( B12D(IMAX)+B22D(IMAX) .LT. 0 ) THEN 00990 IF( WANTV2T ) THEN 00991 IF( COLMAJOR ) THEN 00992 CALL ZSCAL( M-Q, NEGONECOMPLEX, V2T(IMAX,1), LDV2T ) 00993 ELSE 00994 CALL ZSCAL( M-Q, NEGONECOMPLEX, V2T(1,IMAX), 1 ) 00995 END IF 00996 END IF 00997 END IF 00998 * 00999 * Test for negligible sines or cosines 01000 * 01001 DO I = IMIN, IMAX 01002 IF( THETA(I) .LT. THRESH ) THEN 01003 THETA(I) = ZERO 01004 ELSE IF( THETA(I) .GT. PIOVER2-THRESH ) THEN 01005 THETA(I) = PIOVER2 01006 END IF 01007 END DO 01008 DO I = IMIN, IMAX-1 01009 IF( PHI(I) .LT. THRESH ) THEN 01010 PHI(I) = ZERO 01011 ELSE IF( PHI(I) .GT. PIOVER2-THRESH ) THEN 01012 PHI(I) = PIOVER2 01013 END IF 01014 END DO 01015 * 01016 * Deflate 01017 * 01018 IF (IMAX .GT. 1) THEN 01019 DO WHILE( PHI(IMAX-1) .EQ. ZERO ) 01020 IMAX = IMAX - 1 01021 IF (IMAX .LE. 1) EXIT 01022 END DO 01023 END IF 01024 IF( IMIN .GT. IMAX - 1 ) 01025 $ IMIN = IMAX - 1 01026 IF (IMIN .GT. 1) THEN 01027 DO WHILE (PHI(IMIN-1) .NE. ZERO) 01028 IMIN = IMIN - 1 01029 IF (IMIN .LE. 1) EXIT 01030 END DO 01031 END IF 01032 * 01033 * Repeat main iteration loop 01034 * 01035 END DO 01036 * 01037 * Postprocessing: order THETA from least to greatest 01038 * 01039 DO I = 1, Q 01040 * 01041 MINI = I 01042 THETAMIN = THETA(I) 01043 DO J = I+1, Q 01044 IF( THETA(J) .LT. THETAMIN ) THEN 01045 MINI = J 01046 THETAMIN = THETA(J) 01047 END IF 01048 END DO 01049 * 01050 IF( MINI .NE. I ) THEN 01051 THETA(MINI) = THETA(I) 01052 THETA(I) = THETAMIN 01053 IF( COLMAJOR ) THEN 01054 IF( WANTU1 ) 01055 $ CALL ZSWAP( P, U1(1,I), 1, U1(1,MINI), 1 ) 01056 IF( WANTU2 ) 01057 $ CALL ZSWAP( M-P, U2(1,I), 1, U2(1,MINI), 1 ) 01058 IF( WANTV1T ) 01059 $ CALL ZSWAP( Q, V1T(I,1), LDV1T, V1T(MINI,1), LDV1T ) 01060 IF( WANTV2T ) 01061 $ CALL ZSWAP( M-Q, V2T(I,1), LDV2T, V2T(MINI,1), 01062 $ LDV2T ) 01063 ELSE 01064 IF( WANTU1 ) 01065 $ CALL ZSWAP( P, U1(I,1), LDU1, U1(MINI,1), LDU1 ) 01066 IF( WANTU2 ) 01067 $ CALL ZSWAP( M-P, U2(I,1), LDU2, U2(MINI,1), LDU2 ) 01068 IF( WANTV1T ) 01069 $ CALL ZSWAP( Q, V1T(1,I), 1, V1T(1,MINI), 1 ) 01070 IF( WANTV2T ) 01071 $ CALL ZSWAP( M-Q, V2T(1,I), 1, V2T(1,MINI), 1 ) 01072 END IF 01073 END IF 01074 * 01075 END DO 01076 * 01077 RETURN 01078 * 01079 * End of ZBBCSD 01080 * 01081 END 01082