![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CBBCSD 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CBBCSD + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cbbcsd.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cbbcsd.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cbbcsd.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CBBCSD( 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 * REAL B11D( * ), B11E( * ), B12D( * ), B12E( * ), 00032 * $ B21D( * ), B21E( * ), B22D( * ), B22E( * ), 00033 * $ PHI( * ), THETA( * ), RWORK( * ) 00034 * COMPLEX U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ), 00035 * $ V2T( LDV2T, * ) 00036 * .. 00037 * 00038 * 00039 *> \par Purpose: 00040 * ============= 00041 *> 00042 *> \verbatim 00043 *> 00044 *> CBBCSD 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 CUNCSD 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 REAL 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 REAL 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 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 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 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 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 REAL array, dimension (Q) 00209 *> When CBBCSD converges, B11D contains the cosines of THETA(1), 00210 *> ..., THETA(Q). If CBBCSD 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 REAL array, dimension (Q-1) 00218 *> When CBBCSD converges, B11E contains zeros. If CBBCSD 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 REAL array, dimension (Q) 00226 *> When CBBCSD converges, B12D contains the negative sines of 00227 *> THETA(1), ..., THETA(Q). If CBBCSD 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 REAL array, dimension (Q-1) 00235 *> When CBBCSD converges, B12E contains zeros. If CBBCSD 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 CBBCSD 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 REAL, 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 complexOTHERcomputational 00328 * 00329 * ===================================================================== 00330 SUBROUTINE CBBCSD( 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 REAL B11D( * ), B11E( * ), B12D( * ), B12E( * ), 00346 $ B21D( * ), B21E( * ), B22D( * ), B22E( * ), 00347 $ PHI( * ), THETA( * ), RWORK( * ) 00348 COMPLEX U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ), 00349 $ V2T( LDV2T, * ) 00350 * .. 00351 * 00352 * =================================================================== 00353 * 00354 * .. Parameters .. 00355 INTEGER MAXITR 00356 PARAMETER ( MAXITR = 6 ) 00357 REAL HUNDRED, MEIGHTH, ONE, PIOVER2, TEN, ZERO 00358 PARAMETER ( HUNDRED = 100.0E0, MEIGHTH = -0.125E0, 00359 $ ONE = 1.0E0, PIOVER2 = 1.57079632679489662E0, 00360 $ TEN = 10.0E0, ZERO = 0.0E0 ) 00361 COMPLEX NEGONECOMPLEX 00362 PARAMETER ( NEGONECOMPLEX = (-1.0E0,0.0E0) ) 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 REAL 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 Subroutines .. 00377 EXTERNAL CLASR, CSCAL, CSWAP, SLARTGP, SLARTGS, SLAS2, 00378 $ XERBLA 00379 * .. 00380 * .. External Functions .. 00381 REAL SLAMCH 00382 LOGICAL LSAME 00383 EXTERNAL LSAME, SLAMCH 00384 * .. 00385 * .. Intrinsic Functions .. 00386 INTRINSIC ABS, ATAN2, COS, MAX, MIN, SIN, SQRT 00387 * .. 00388 * .. Executable Statements .. 00389 * 00390 * Test input arguments 00391 * 00392 INFO = 0 00393 LQUERY = LRWORK .EQ. -1 00394 WANTU1 = LSAME( JOBU1, 'Y' ) 00395 WANTU2 = LSAME( JOBU2, 'Y' ) 00396 WANTV1T = LSAME( JOBV1T, 'Y' ) 00397 WANTV2T = LSAME( JOBV2T, 'Y' ) 00398 COLMAJOR = .NOT. LSAME( TRANS, 'T' ) 00399 * 00400 IF( M .LT. 0 ) THEN 00401 INFO = -6 00402 ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN 00403 INFO = -7 00404 ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN 00405 INFO = -8 00406 ELSE IF( Q .GT. P .OR. Q .GT. M-P .OR. Q .GT. M-Q ) THEN 00407 INFO = -8 00408 ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN 00409 INFO = -12 00410 ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN 00411 INFO = -14 00412 ELSE IF( WANTV1T .AND. LDV1T .LT. Q ) THEN 00413 INFO = -16 00414 ELSE IF( WANTV2T .AND. LDV2T .LT. M-Q ) THEN 00415 INFO = -18 00416 END IF 00417 * 00418 * Quick return if Q = 0 00419 * 00420 IF( INFO .EQ. 0 .AND. Q .EQ. 0 ) THEN 00421 LRWORKMIN = 1 00422 RWORK(1) = LRWORKMIN 00423 RETURN 00424 END IF 00425 * 00426 * Compute workspace 00427 * 00428 IF( INFO .EQ. 0 ) THEN 00429 IU1CS = 1 00430 IU1SN = IU1CS + Q 00431 IU2CS = IU1SN + Q 00432 IU2SN = IU2CS + Q 00433 IV1TCS = IU2SN + Q 00434 IV1TSN = IV1TCS + Q 00435 IV2TCS = IV1TSN + Q 00436 IV2TSN = IV2TCS + Q 00437 LRWORKOPT = IV2TSN + Q - 1 00438 LRWORKMIN = LRWORKOPT 00439 RWORK(1) = LRWORKOPT 00440 IF( LRWORK .LT. LRWORKMIN .AND. .NOT. LQUERY ) THEN 00441 INFO = -28 00442 END IF 00443 END IF 00444 * 00445 IF( INFO .NE. 0 ) THEN 00446 CALL XERBLA( 'CBBCSD', -INFO ) 00447 RETURN 00448 ELSE IF( LQUERY ) THEN 00449 RETURN 00450 END IF 00451 * 00452 * Get machine constants 00453 * 00454 EPS = SLAMCH( 'Epsilon' ) 00455 UNFL = SLAMCH( 'Safe minimum' ) 00456 TOLMUL = MAX( TEN, MIN( HUNDRED, EPS**MEIGHTH ) ) 00457 TOL = TOLMUL*EPS 00458 THRESH = MAX( TOL, MAXITR*Q*Q*UNFL ) 00459 * 00460 * Test for negligible sines or cosines 00461 * 00462 DO I = 1, Q 00463 IF( THETA(I) .LT. THRESH ) THEN 00464 THETA(I) = ZERO 00465 ELSE IF( THETA(I) .GT. PIOVER2-THRESH ) THEN 00466 THETA(I) = PIOVER2 00467 END IF 00468 END DO 00469 DO I = 1, Q-1 00470 IF( PHI(I) .LT. THRESH ) THEN 00471 PHI(I) = ZERO 00472 ELSE IF( PHI(I) .GT. PIOVER2-THRESH ) THEN 00473 PHI(I) = PIOVER2 00474 END IF 00475 END DO 00476 * 00477 * Initial deflation 00478 * 00479 IMAX = Q 00480 DO WHILE( ( IMAX .GT. 1 ) .AND. ( PHI(IMAX-1) .EQ. ZERO ) ) 00481 IMAX = IMAX - 1 00482 END DO 00483 IMIN = IMAX - 1 00484 IF ( IMIN .GT. 1 ) THEN 00485 DO WHILE( PHI(IMIN-1) .NE. ZERO ) 00486 IMIN = IMIN - 1 00487 IF ( IMIN .LE. 1 ) EXIT 00488 END DO 00489 END IF 00490 * 00491 * Initialize iteration counter 00492 * 00493 MAXIT = MAXITR*Q*Q 00494 ITER = 0 00495 * 00496 * Begin main iteration loop 00497 * 00498 DO WHILE( IMAX .GT. 1 ) 00499 * 00500 * Compute the matrix entries 00501 * 00502 B11D(IMIN) = COS( THETA(IMIN) ) 00503 B21D(IMIN) = -SIN( THETA(IMIN) ) 00504 DO I = IMIN, IMAX - 1 00505 B11E(I) = -SIN( THETA(I) ) * SIN( PHI(I) ) 00506 B11D(I+1) = COS( THETA(I+1) ) * COS( PHI(I) ) 00507 B12D(I) = SIN( THETA(I) ) * COS( PHI(I) ) 00508 B12E(I) = COS( THETA(I+1) ) * SIN( PHI(I) ) 00509 B21E(I) = -COS( THETA(I) ) * SIN( PHI(I) ) 00510 B21D(I+1) = -SIN( THETA(I+1) ) * COS( PHI(I) ) 00511 B22D(I) = COS( THETA(I) ) * COS( PHI(I) ) 00512 B22E(I) = -SIN( THETA(I+1) ) * SIN( PHI(I) ) 00513 END DO 00514 B12D(IMAX) = SIN( THETA(IMAX) ) 00515 B22D(IMAX) = COS( THETA(IMAX) ) 00516 * 00517 * Abort if not converging; otherwise, increment ITER 00518 * 00519 IF( ITER .GT. MAXIT ) THEN 00520 INFO = 0 00521 DO I = 1, Q 00522 IF( PHI(I) .NE. ZERO ) 00523 $ INFO = INFO + 1 00524 END DO 00525 RETURN 00526 END IF 00527 * 00528 ITER = ITER + IMAX - IMIN 00529 * 00530 * Compute shifts 00531 * 00532 THETAMAX = THETA(IMIN) 00533 THETAMIN = THETA(IMIN) 00534 DO I = IMIN+1, IMAX 00535 IF( THETA(I) > THETAMAX ) 00536 $ THETAMAX = THETA(I) 00537 IF( THETA(I) < THETAMIN ) 00538 $ THETAMIN = THETA(I) 00539 END DO 00540 * 00541 IF( THETAMAX .GT. PIOVER2 - THRESH ) THEN 00542 * 00543 * Zero on diagonals of B11 and B22; induce deflation with a 00544 * zero shift 00545 * 00546 MU = ZERO 00547 NU = ONE 00548 * 00549 ELSE IF( THETAMIN .LT. THRESH ) THEN 00550 * 00551 * Zero on diagonals of B12 and B22; induce deflation with a 00552 * zero shift 00553 * 00554 MU = ONE 00555 NU = ZERO 00556 * 00557 ELSE 00558 * 00559 * Compute shifts for B11 and B21 and use the lesser 00560 * 00561 CALL SLAS2( B11D(IMAX-1), B11E(IMAX-1), B11D(IMAX), SIGMA11, 00562 $ DUMMY ) 00563 CALL SLAS2( B21D(IMAX-1), B21E(IMAX-1), B21D(IMAX), SIGMA21, 00564 $ DUMMY ) 00565 * 00566 IF( SIGMA11 .LE. SIGMA21 ) THEN 00567 MU = SIGMA11 00568 NU = SQRT( ONE - MU**2 ) 00569 IF( MU .LT. THRESH ) THEN 00570 MU = ZERO 00571 NU = ONE 00572 END IF 00573 ELSE 00574 NU = SIGMA21 00575 MU = SQRT( 1.0 - NU**2 ) 00576 IF( NU .LT. THRESH ) THEN 00577 MU = ONE 00578 NU = ZERO 00579 END IF 00580 END IF 00581 END IF 00582 * 00583 * Rotate to produce bulges in B11 and B21 00584 * 00585 IF( MU .LE. NU ) THEN 00586 CALL SLARTGS( B11D(IMIN), B11E(IMIN), MU, 00587 $ RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1) ) 00588 ELSE 00589 CALL SLARTGS( B21D(IMIN), B21E(IMIN), NU, 00590 $ RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1) ) 00591 END IF 00592 * 00593 TEMP = RWORK(IV1TCS+IMIN-1)*B11D(IMIN) + 00594 $ RWORK(IV1TSN+IMIN-1)*B11E(IMIN) 00595 B11E(IMIN) = RWORK(IV1TCS+IMIN-1)*B11E(IMIN) - 00596 $ RWORK(IV1TSN+IMIN-1)*B11D(IMIN) 00597 B11D(IMIN) = TEMP 00598 B11BULGE = RWORK(IV1TSN+IMIN-1)*B11D(IMIN+1) 00599 B11D(IMIN+1) = RWORK(IV1TCS+IMIN-1)*B11D(IMIN+1) 00600 TEMP = RWORK(IV1TCS+IMIN-1)*B21D(IMIN) + 00601 $ RWORK(IV1TSN+IMIN-1)*B21E(IMIN) 00602 B21E(IMIN) = RWORK(IV1TCS+IMIN-1)*B21E(IMIN) - 00603 $ RWORK(IV1TSN+IMIN-1)*B21D(IMIN) 00604 B21D(IMIN) = TEMP 00605 B21BULGE = RWORK(IV1TSN+IMIN-1)*B21D(IMIN+1) 00606 B21D(IMIN+1) = RWORK(IV1TCS+IMIN-1)*B21D(IMIN+1) 00607 * 00608 * Compute THETA(IMIN) 00609 * 00610 THETA( IMIN ) = ATAN2( SQRT( B21D(IMIN)**2+B21BULGE**2 ), 00611 $ SQRT( B11D(IMIN)**2+B11BULGE**2 ) ) 00612 * 00613 * Chase the bulges in B11(IMIN+1,IMIN) and B21(IMIN+1,IMIN) 00614 * 00615 IF( B11D(IMIN)**2+B11BULGE**2 .GT. THRESH**2 ) THEN 00616 CALL SLARTGP( B11BULGE, B11D(IMIN), RWORK(IU1SN+IMIN-1), 00617 $ RWORK(IU1CS+IMIN-1), R ) 00618 ELSE IF( MU .LE. NU ) THEN 00619 CALL SLARTGS( B11E( IMIN ), B11D( IMIN + 1 ), MU, 00620 $ RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1) ) 00621 ELSE 00622 CALL SLARTGS( B12D( IMIN ), B12E( IMIN ), NU, 00623 $ RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1) ) 00624 END IF 00625 IF( B21D(IMIN)**2+B21BULGE**2 .GT. THRESH**2 ) THEN 00626 CALL SLARTGP( B21BULGE, B21D(IMIN), RWORK(IU2SN+IMIN-1), 00627 $ RWORK(IU2CS+IMIN-1), R ) 00628 ELSE IF( NU .LT. MU ) THEN 00629 CALL SLARTGS( B21E( IMIN ), B21D( IMIN + 1 ), NU, 00630 $ RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1) ) 00631 ELSE 00632 CALL SLARTGS( B22D(IMIN), B22E(IMIN), MU, 00633 $ RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1) ) 00634 END IF 00635 RWORK(IU2CS+IMIN-1) = -RWORK(IU2CS+IMIN-1) 00636 RWORK(IU2SN+IMIN-1) = -RWORK(IU2SN+IMIN-1) 00637 * 00638 TEMP = RWORK(IU1CS+IMIN-1)*B11E(IMIN) + 00639 $ RWORK(IU1SN+IMIN-1)*B11D(IMIN+1) 00640 B11D(IMIN+1) = RWORK(IU1CS+IMIN-1)*B11D(IMIN+1) - 00641 $ RWORK(IU1SN+IMIN-1)*B11E(IMIN) 00642 B11E(IMIN) = TEMP 00643 IF( IMAX .GT. IMIN+1 ) THEN 00644 B11BULGE = RWORK(IU1SN+IMIN-1)*B11E(IMIN+1) 00645 B11E(IMIN+1) = RWORK(IU1CS+IMIN-1)*B11E(IMIN+1) 00646 END IF 00647 TEMP = RWORK(IU1CS+IMIN-1)*B12D(IMIN) + 00648 $ RWORK(IU1SN+IMIN-1)*B12E(IMIN) 00649 B12E(IMIN) = RWORK(IU1CS+IMIN-1)*B12E(IMIN) - 00650 $ RWORK(IU1SN+IMIN-1)*B12D(IMIN) 00651 B12D(IMIN) = TEMP 00652 B12BULGE = RWORK(IU1SN+IMIN-1)*B12D(IMIN+1) 00653 B12D(IMIN+1) = RWORK(IU1CS+IMIN-1)*B12D(IMIN+1) 00654 TEMP = RWORK(IU2CS+IMIN-1)*B21E(IMIN) + 00655 $ RWORK(IU2SN+IMIN-1)*B21D(IMIN+1) 00656 B21D(IMIN+1) = RWORK(IU2CS+IMIN-1)*B21D(IMIN+1) - 00657 $ RWORK(IU2SN+IMIN-1)*B21E(IMIN) 00658 B21E(IMIN) = TEMP 00659 IF( IMAX .GT. IMIN+1 ) THEN 00660 B21BULGE = RWORK(IU2SN+IMIN-1)*B21E(IMIN+1) 00661 B21E(IMIN+1) = RWORK(IU2CS+IMIN-1)*B21E(IMIN+1) 00662 END IF 00663 TEMP = RWORK(IU2CS+IMIN-1)*B22D(IMIN) + 00664 $ RWORK(IU2SN+IMIN-1)*B22E(IMIN) 00665 B22E(IMIN) = RWORK(IU2CS+IMIN-1)*B22E(IMIN) - 00666 $ RWORK(IU2SN+IMIN-1)*B22D(IMIN) 00667 B22D(IMIN) = TEMP 00668 B22BULGE = RWORK(IU2SN+IMIN-1)*B22D(IMIN+1) 00669 B22D(IMIN+1) = RWORK(IU2CS+IMIN-1)*B22D(IMIN+1) 00670 * 00671 * Inner loop: chase bulges from B11(IMIN,IMIN+2), 00672 * B12(IMIN,IMIN+1), B21(IMIN,IMIN+2), and B22(IMIN,IMIN+1) to 00673 * bottom-right 00674 * 00675 DO I = IMIN+1, IMAX-1 00676 * 00677 * Compute PHI(I-1) 00678 * 00679 X1 = SIN(THETA(I-1))*B11E(I-1) + COS(THETA(I-1))*B21E(I-1) 00680 X2 = SIN(THETA(I-1))*B11BULGE + COS(THETA(I-1))*B21BULGE 00681 Y1 = SIN(THETA(I-1))*B12D(I-1) + COS(THETA(I-1))*B22D(I-1) 00682 Y2 = SIN(THETA(I-1))*B12BULGE + COS(THETA(I-1))*B22BULGE 00683 * 00684 PHI(I-1) = ATAN2( SQRT(X1**2+X2**2), SQRT(Y1**2+Y2**2) ) 00685 * 00686 * Determine if there are bulges to chase or if a new direct 00687 * summand has been reached 00688 * 00689 RESTART11 = B11E(I-1)**2 + B11BULGE**2 .LE. THRESH**2 00690 RESTART21 = B21E(I-1)**2 + B21BULGE**2 .LE. THRESH**2 00691 RESTART12 = B12D(I-1)**2 + B12BULGE**2 .LE. THRESH**2 00692 RESTART22 = B22D(I-1)**2 + B22BULGE**2 .LE. THRESH**2 00693 * 00694 * If possible, chase bulges from B11(I-1,I+1), B12(I-1,I), 00695 * B21(I-1,I+1), and B22(I-1,I). If necessary, restart bulge- 00696 * chasing by applying the original shift again. 00697 * 00698 IF( .NOT. RESTART11 .AND. .NOT. RESTART21 ) THEN 00699 CALL SLARTGP( X2, X1, RWORK(IV1TSN+I-1), 00700 $ RWORK(IV1TCS+I-1), R ) 00701 ELSE IF( .NOT. RESTART11 .AND. RESTART21 ) THEN 00702 CALL SLARTGP( B11BULGE, B11E(I-1), RWORK(IV1TSN+I-1), 00703 $ RWORK(IV1TCS+I-1), R ) 00704 ELSE IF( RESTART11 .AND. .NOT. RESTART21 ) THEN 00705 CALL SLARTGP( B21BULGE, B21E(I-1), RWORK(IV1TSN+I-1), 00706 $ RWORK(IV1TCS+I-1), R ) 00707 ELSE IF( MU .LE. NU ) THEN 00708 CALL SLARTGS( B11D(I), B11E(I), MU, RWORK(IV1TCS+I-1), 00709 $ RWORK(IV1TSN+I-1) ) 00710 ELSE 00711 CALL SLARTGS( B21D(I), B21E(I), NU, RWORK(IV1TCS+I-1), 00712 $ RWORK(IV1TSN+I-1) ) 00713 END IF 00714 RWORK(IV1TCS+I-1) = -RWORK(IV1TCS+I-1) 00715 RWORK(IV1TSN+I-1) = -RWORK(IV1TSN+I-1) 00716 IF( .NOT. RESTART12 .AND. .NOT. RESTART22 ) THEN 00717 CALL SLARTGP( Y2, Y1, RWORK(IV2TSN+I-1-1), 00718 $ RWORK(IV2TCS+I-1-1), R ) 00719 ELSE IF( .NOT. RESTART12 .AND. RESTART22 ) THEN 00720 CALL SLARTGP( B12BULGE, B12D(I-1), RWORK(IV2TSN+I-1-1), 00721 $ RWORK(IV2TCS+I-1-1), R ) 00722 ELSE IF( RESTART12 .AND. .NOT. RESTART22 ) THEN 00723 CALL SLARTGP( B22BULGE, B22D(I-1), RWORK(IV2TSN+I-1-1), 00724 $ RWORK(IV2TCS+I-1-1), R ) 00725 ELSE IF( NU .LT. MU ) THEN 00726 CALL SLARTGS( B12E(I-1), B12D(I), NU, 00727 $ RWORK(IV2TCS+I-1-1), RWORK(IV2TSN+I-1-1) ) 00728 ELSE 00729 CALL SLARTGS( B22E(I-1), B22D(I), MU, 00730 $ RWORK(IV2TCS+I-1-1), RWORK(IV2TSN+I-1-1) ) 00731 END IF 00732 * 00733 TEMP = RWORK(IV1TCS+I-1)*B11D(I) + RWORK(IV1TSN+I-1)*B11E(I) 00734 B11E(I) = RWORK(IV1TCS+I-1)*B11E(I) - 00735 $ RWORK(IV1TSN+I-1)*B11D(I) 00736 B11D(I) = TEMP 00737 B11BULGE = RWORK(IV1TSN+I-1)*B11D(I+1) 00738 B11D(I+1) = RWORK(IV1TCS+I-1)*B11D(I+1) 00739 TEMP = RWORK(IV1TCS+I-1)*B21D(I) + RWORK(IV1TSN+I-1)*B21E(I) 00740 B21E(I) = RWORK(IV1TCS+I-1)*B21E(I) - 00741 $ RWORK(IV1TSN+I-1)*B21D(I) 00742 B21D(I) = TEMP 00743 B21BULGE = RWORK(IV1TSN+I-1)*B21D(I+1) 00744 B21D(I+1) = RWORK(IV1TCS+I-1)*B21D(I+1) 00745 TEMP = RWORK(IV2TCS+I-1-1)*B12E(I-1) + 00746 $ RWORK(IV2TSN+I-1-1)*B12D(I) 00747 B12D(I) = RWORK(IV2TCS+I-1-1)*B12D(I) - 00748 $ RWORK(IV2TSN+I-1-1)*B12E(I-1) 00749 B12E(I-1) = TEMP 00750 B12BULGE = RWORK(IV2TSN+I-1-1)*B12E(I) 00751 B12E(I) = RWORK(IV2TCS+I-1-1)*B12E(I) 00752 TEMP = RWORK(IV2TCS+I-1-1)*B22E(I-1) + 00753 $ RWORK(IV2TSN+I-1-1)*B22D(I) 00754 B22D(I) = RWORK(IV2TCS+I-1-1)*B22D(I) - 00755 $ RWORK(IV2TSN+I-1-1)*B22E(I-1) 00756 B22E(I-1) = TEMP 00757 B22BULGE = RWORK(IV2TSN+I-1-1)*B22E(I) 00758 B22E(I) = RWORK(IV2TCS+I-1-1)*B22E(I) 00759 * 00760 * Compute THETA(I) 00761 * 00762 X1 = COS(PHI(I-1))*B11D(I) + SIN(PHI(I-1))*B12E(I-1) 00763 X2 = COS(PHI(I-1))*B11BULGE + SIN(PHI(I-1))*B12BULGE 00764 Y1 = COS(PHI(I-1))*B21D(I) + SIN(PHI(I-1))*B22E(I-1) 00765 Y2 = COS(PHI(I-1))*B21BULGE + SIN(PHI(I-1))*B22BULGE 00766 * 00767 THETA(I) = ATAN2( SQRT(Y1**2+Y2**2), SQRT(X1**2+X2**2) ) 00768 * 00769 * Determine if there are bulges to chase or if a new direct 00770 * summand has been reached 00771 * 00772 RESTART11 = B11D(I)**2 + B11BULGE**2 .LE. THRESH**2 00773 RESTART12 = B12E(I-1)**2 + B12BULGE**2 .LE. THRESH**2 00774 RESTART21 = B21D(I)**2 + B21BULGE**2 .LE. THRESH**2 00775 RESTART22 = B22E(I-1)**2 + B22BULGE**2 .LE. THRESH**2 00776 * 00777 * If possible, chase bulges from B11(I+1,I), B12(I+1,I-1), 00778 * B21(I+1,I), and B22(I+1,I-1). If necessary, restart bulge- 00779 * chasing by applying the original shift again. 00780 * 00781 IF( .NOT. RESTART11 .AND. .NOT. RESTART12 ) THEN 00782 CALL SLARTGP( X2, X1, RWORK(IU1SN+I-1), RWORK(IU1CS+I-1), 00783 $ R ) 00784 ELSE IF( .NOT. RESTART11 .AND. RESTART12 ) THEN 00785 CALL SLARTGP( B11BULGE, B11D(I), RWORK(IU1SN+I-1), 00786 $ RWORK(IU1CS+I-1), R ) 00787 ELSE IF( RESTART11 .AND. .NOT. RESTART12 ) THEN 00788 CALL SLARTGP( B12BULGE, B12E(I-1), RWORK(IU1SN+I-1), 00789 $ RWORK(IU1CS+I-1), R ) 00790 ELSE IF( MU .LE. NU ) THEN 00791 CALL SLARTGS( B11E(I), B11D(I+1), MU, RWORK(IU1CS+I-1), 00792 $ RWORK(IU1SN+I-1) ) 00793 ELSE 00794 CALL SLARTGS( B12D(I), B12E(I), NU, RWORK(IU1CS+I-1), 00795 $ RWORK(IU1SN+I-1) ) 00796 END IF 00797 IF( .NOT. RESTART21 .AND. .NOT. RESTART22 ) THEN 00798 CALL SLARTGP( Y2, Y1, RWORK(IU2SN+I-1), RWORK(IU2CS+I-1), 00799 $ R ) 00800 ELSE IF( .NOT. RESTART21 .AND. RESTART22 ) THEN 00801 CALL SLARTGP( B21BULGE, B21D(I), RWORK(IU2SN+I-1), 00802 $ RWORK(IU2CS+I-1), R ) 00803 ELSE IF( RESTART21 .AND. .NOT. RESTART22 ) THEN 00804 CALL SLARTGP( B22BULGE, B22E(I-1), RWORK(IU2SN+I-1), 00805 $ RWORK(IU2CS+I-1), R ) 00806 ELSE IF( NU .LT. MU ) THEN 00807 CALL SLARTGS( B21E(I), B21E(I+1), NU, RWORK(IU2CS+I-1), 00808 $ RWORK(IU2SN+I-1) ) 00809 ELSE 00810 CALL SLARTGS( B22D(I), B22E(I), MU, RWORK(IU2CS+I-1), 00811 $ RWORK(IU2SN+I-1) ) 00812 END IF 00813 RWORK(IU2CS+I-1) = -RWORK(IU2CS+I-1) 00814 RWORK(IU2SN+I-1) = -RWORK(IU2SN+I-1) 00815 * 00816 TEMP = RWORK(IU1CS+I-1)*B11E(I) + RWORK(IU1SN+I-1)*B11D(I+1) 00817 B11D(I+1) = RWORK(IU1CS+I-1)*B11D(I+1) - 00818 $ RWORK(IU1SN+I-1)*B11E(I) 00819 B11E(I) = TEMP 00820 IF( I .LT. IMAX - 1 ) THEN 00821 B11BULGE = RWORK(IU1SN+I-1)*B11E(I+1) 00822 B11E(I+1) = RWORK(IU1CS+I-1)*B11E(I+1) 00823 END IF 00824 TEMP = RWORK(IU2CS+I-1)*B21E(I) + RWORK(IU2SN+I-1)*B21D(I+1) 00825 B21D(I+1) = RWORK(IU2CS+I-1)*B21D(I+1) - 00826 $ RWORK(IU2SN+I-1)*B21E(I) 00827 B21E(I) = TEMP 00828 IF( I .LT. IMAX - 1 ) THEN 00829 B21BULGE = RWORK(IU2SN+I-1)*B21E(I+1) 00830 B21E(I+1) = RWORK(IU2CS+I-1)*B21E(I+1) 00831 END IF 00832 TEMP = RWORK(IU1CS+I-1)*B12D(I) + RWORK(IU1SN+I-1)*B12E(I) 00833 B12E(I) = RWORK(IU1CS+I-1)*B12E(I) - 00834 $ RWORK(IU1SN+I-1)*B12D(I) 00835 B12D(I) = TEMP 00836 B12BULGE = RWORK(IU1SN+I-1)*B12D(I+1) 00837 B12D(I+1) = RWORK(IU1CS+I-1)*B12D(I+1) 00838 TEMP = RWORK(IU2CS+I-1)*B22D(I) + RWORK(IU2SN+I-1)*B22E(I) 00839 B22E(I) = RWORK(IU2CS+I-1)*B22E(I) - 00840 $ RWORK(IU2SN+I-1)*B22D(I) 00841 B22D(I) = TEMP 00842 B22BULGE = RWORK(IU2SN+I-1)*B22D(I+1) 00843 B22D(I+1) = RWORK(IU2CS+I-1)*B22D(I+1) 00844 * 00845 END DO 00846 * 00847 * Compute PHI(IMAX-1) 00848 * 00849 X1 = SIN(THETA(IMAX-1))*B11E(IMAX-1) + 00850 $ COS(THETA(IMAX-1))*B21E(IMAX-1) 00851 Y1 = SIN(THETA(IMAX-1))*B12D(IMAX-1) + 00852 $ COS(THETA(IMAX-1))*B22D(IMAX-1) 00853 Y2 = SIN(THETA(IMAX-1))*B12BULGE + COS(THETA(IMAX-1))*B22BULGE 00854 * 00855 PHI(IMAX-1) = ATAN2( ABS(X1), SQRT(Y1**2+Y2**2) ) 00856 * 00857 * Chase bulges from B12(IMAX-1,IMAX) and B22(IMAX-1,IMAX) 00858 * 00859 RESTART12 = B12D(IMAX-1)**2 + B12BULGE**2 .LE. THRESH**2 00860 RESTART22 = B22D(IMAX-1)**2 + B22BULGE**2 .LE. THRESH**2 00861 * 00862 IF( .NOT. RESTART12 .AND. .NOT. RESTART22 ) THEN 00863 CALL SLARTGP( Y2, Y1, RWORK(IV2TSN+IMAX-1-1), 00864 $ RWORK(IV2TCS+IMAX-1-1), R ) 00865 ELSE IF( .NOT. RESTART12 .AND. RESTART22 ) THEN 00866 CALL SLARTGP( B12BULGE, B12D(IMAX-1), 00867 $ RWORK(IV2TSN+IMAX-1-1), 00868 $ RWORK(IV2TCS+IMAX-1-1), R ) 00869 ELSE IF( RESTART12 .AND. .NOT. RESTART22 ) THEN 00870 CALL SLARTGP( B22BULGE, B22D(IMAX-1), 00871 $ RWORK(IV2TSN+IMAX-1-1), 00872 $ RWORK(IV2TCS+IMAX-1-1), R ) 00873 ELSE IF( NU .LT. MU ) THEN 00874 CALL SLARTGS( B12E(IMAX-1), B12D(IMAX), NU, 00875 $ RWORK(IV2TCS+IMAX-1-1), 00876 $ RWORK(IV2TSN+IMAX-1-1) ) 00877 ELSE 00878 CALL SLARTGS( B22E(IMAX-1), B22D(IMAX), MU, 00879 $ RWORK(IV2TCS+IMAX-1-1), 00880 $ RWORK(IV2TSN+IMAX-1-1) ) 00881 END IF 00882 * 00883 TEMP = RWORK(IV2TCS+IMAX-1-1)*B12E(IMAX-1) + 00884 $ RWORK(IV2TSN+IMAX-1-1)*B12D(IMAX) 00885 B12D(IMAX) = RWORK(IV2TCS+IMAX-1-1)*B12D(IMAX) - 00886 $ RWORK(IV2TSN+IMAX-1-1)*B12E(IMAX-1) 00887 B12E(IMAX-1) = TEMP 00888 TEMP = RWORK(IV2TCS+IMAX-1-1)*B22E(IMAX-1) + 00889 $ RWORK(IV2TSN+IMAX-1-1)*B22D(IMAX) 00890 B22D(IMAX) = RWORK(IV2TCS+IMAX-1-1)*B22D(IMAX) - 00891 $ RWORK(IV2TSN+IMAX-1-1)*B22E(IMAX-1) 00892 B22E(IMAX-1) = TEMP 00893 * 00894 * Update singular vectors 00895 * 00896 IF( WANTU1 ) THEN 00897 IF( COLMAJOR ) THEN 00898 CALL CLASR( 'R', 'V', 'F', P, IMAX-IMIN+1, 00899 $ RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1), 00900 $ U1(1,IMIN), LDU1 ) 00901 ELSE 00902 CALL CLASR( 'L', 'V', 'F', IMAX-IMIN+1, P, 00903 $ RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1), 00904 $ U1(IMIN,1), LDU1 ) 00905 END IF 00906 END IF 00907 IF( WANTU2 ) THEN 00908 IF( COLMAJOR ) THEN 00909 CALL CLASR( 'R', 'V', 'F', M-P, IMAX-IMIN+1, 00910 $ RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1), 00911 $ U2(1,IMIN), LDU2 ) 00912 ELSE 00913 CALL CLASR( 'L', 'V', 'F', IMAX-IMIN+1, M-P, 00914 $ RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1), 00915 $ U2(IMIN,1), LDU2 ) 00916 END IF 00917 END IF 00918 IF( WANTV1T ) THEN 00919 IF( COLMAJOR ) THEN 00920 CALL CLASR( 'L', 'V', 'F', IMAX-IMIN+1, Q, 00921 $ RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1), 00922 $ V1T(IMIN,1), LDV1T ) 00923 ELSE 00924 CALL CLASR( 'R', 'V', 'F', Q, IMAX-IMIN+1, 00925 $ RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1), 00926 $ V1T(1,IMIN), LDV1T ) 00927 END IF 00928 END IF 00929 IF( WANTV2T ) THEN 00930 IF( COLMAJOR ) THEN 00931 CALL CLASR( 'L', 'V', 'F', IMAX-IMIN+1, M-Q, 00932 $ RWORK(IV2TCS+IMIN-1), RWORK(IV2TSN+IMIN-1), 00933 $ V2T(IMIN,1), LDV2T ) 00934 ELSE 00935 CALL CLASR( 'R', 'V', 'F', M-Q, IMAX-IMIN+1, 00936 $ RWORK(IV2TCS+IMIN-1), RWORK(IV2TSN+IMIN-1), 00937 $ V2T(1,IMIN), LDV2T ) 00938 END IF 00939 END IF 00940 * 00941 * Fix signs on B11(IMAX-1,IMAX) and B21(IMAX-1,IMAX) 00942 * 00943 IF( B11E(IMAX-1)+B21E(IMAX-1) .GT. 0 ) THEN 00944 B11D(IMAX) = -B11D(IMAX) 00945 B21D(IMAX) = -B21D(IMAX) 00946 IF( WANTV1T ) THEN 00947 IF( COLMAJOR ) THEN 00948 CALL CSCAL( Q, NEGONECOMPLEX, V1T(IMAX,1), LDV1T ) 00949 ELSE 00950 CALL CSCAL( Q, NEGONECOMPLEX, V1T(1,IMAX), 1 ) 00951 END IF 00952 END IF 00953 END IF 00954 * 00955 * Compute THETA(IMAX) 00956 * 00957 X1 = COS(PHI(IMAX-1))*B11D(IMAX) + 00958 $ SIN(PHI(IMAX-1))*B12E(IMAX-1) 00959 Y1 = COS(PHI(IMAX-1))*B21D(IMAX) + 00960 $ SIN(PHI(IMAX-1))*B22E(IMAX-1) 00961 * 00962 THETA(IMAX) = ATAN2( ABS(Y1), ABS(X1) ) 00963 * 00964 * Fix signs on B11(IMAX,IMAX), B12(IMAX,IMAX-1), B21(IMAX,IMAX), 00965 * and B22(IMAX,IMAX-1) 00966 * 00967 IF( B11D(IMAX)+B12E(IMAX-1) .LT. 0 ) THEN 00968 B12D(IMAX) = -B12D(IMAX) 00969 IF( WANTU1 ) THEN 00970 IF( COLMAJOR ) THEN 00971 CALL CSCAL( P, NEGONECOMPLEX, U1(1,IMAX), 1 ) 00972 ELSE 00973 CALL CSCAL( P, NEGONECOMPLEX, U1(IMAX,1), LDU1 ) 00974 END IF 00975 END IF 00976 END IF 00977 IF( B21D(IMAX)+B22E(IMAX-1) .GT. 0 ) THEN 00978 B22D(IMAX) = -B22D(IMAX) 00979 IF( WANTU2 ) THEN 00980 IF( COLMAJOR ) THEN 00981 CALL CSCAL( M-P, NEGONECOMPLEX, U2(1,IMAX), 1 ) 00982 ELSE 00983 CALL CSCAL( M-P, NEGONECOMPLEX, U2(IMAX,1), LDU2 ) 00984 END IF 00985 END IF 00986 END IF 00987 * 00988 * Fix signs on B12(IMAX,IMAX) and B22(IMAX,IMAX) 00989 * 00990 IF( B12D(IMAX)+B22D(IMAX) .LT. 0 ) THEN 00991 IF( WANTV2T ) THEN 00992 IF( COLMAJOR ) THEN 00993 CALL CSCAL( M-Q, NEGONECOMPLEX, V2T(IMAX,1), LDV2T ) 00994 ELSE 00995 CALL CSCAL( M-Q, NEGONECOMPLEX, V2T(1,IMAX), 1 ) 00996 END IF 00997 END IF 00998 END IF 00999 * 01000 * Test for negligible sines or cosines 01001 * 01002 DO I = IMIN, IMAX 01003 IF( THETA(I) .LT. THRESH ) THEN 01004 THETA(I) = ZERO 01005 ELSE IF( THETA(I) .GT. PIOVER2-THRESH ) THEN 01006 THETA(I) = PIOVER2 01007 END IF 01008 END DO 01009 DO I = IMIN, IMAX-1 01010 IF( PHI(I) .LT. THRESH ) THEN 01011 PHI(I) = ZERO 01012 ELSE IF( PHI(I) .GT. PIOVER2-THRESH ) THEN 01013 PHI(I) = PIOVER2 01014 END IF 01015 END DO 01016 * 01017 * Deflate 01018 * 01019 IF (IMAX .GT. 1) THEN 01020 DO WHILE( PHI(IMAX-1) .EQ. ZERO ) 01021 IMAX = IMAX - 1 01022 IF (IMAX .LE. 1) EXIT 01023 END DO 01024 END IF 01025 IF( IMIN .GT. IMAX - 1 ) 01026 $ IMIN = IMAX - 1 01027 IF (IMIN .GT. 1) THEN 01028 DO WHILE (PHI(IMIN-1) .NE. ZERO) 01029 IMIN = IMIN - 1 01030 IF (IMIN .LE. 1) EXIT 01031 END DO 01032 END IF 01033 * 01034 * Repeat main iteration loop 01035 * 01036 END DO 01037 * 01038 * Postprocessing: order THETA from least to greatest 01039 * 01040 DO I = 1, Q 01041 * 01042 MINI = I 01043 THETAMIN = THETA(I) 01044 DO J = I+1, Q 01045 IF( THETA(J) .LT. THETAMIN ) THEN 01046 MINI = J 01047 THETAMIN = THETA(J) 01048 END IF 01049 END DO 01050 * 01051 IF( MINI .NE. I ) THEN 01052 THETA(MINI) = THETA(I) 01053 THETA(I) = THETAMIN 01054 IF( COLMAJOR ) THEN 01055 IF( WANTU1 ) 01056 $ CALL CSWAP( P, U1(1,I), 1, U1(1,MINI), 1 ) 01057 IF( WANTU2 ) 01058 $ CALL CSWAP( M-P, U2(1,I), 1, U2(1,MINI), 1 ) 01059 IF( WANTV1T ) 01060 $ CALL CSWAP( Q, V1T(I,1), LDV1T, V1T(MINI,1), LDV1T ) 01061 IF( WANTV2T ) 01062 $ CALL CSWAP( M-Q, V2T(I,1), LDV2T, V2T(MINI,1), 01063 $ LDV2T ) 01064 ELSE 01065 IF( WANTU1 ) 01066 $ CALL CSWAP( P, U1(I,1), LDU1, U1(MINI,1), LDU1 ) 01067 IF( WANTU2 ) 01068 $ CALL CSWAP( M-P, U2(I,1), LDU2, U2(MINI,1), LDU2 ) 01069 IF( WANTV1T ) 01070 $ CALL CSWAP( Q, V1T(1,I), 1, V1T(1,MINI), 1 ) 01071 IF( WANTV2T ) 01072 $ CALL CSWAP( M-Q, V2T(1,I), 1, V2T(1,MINI), 1 ) 01073 END IF 01074 END IF 01075 * 01076 END DO 01077 * 01078 RETURN 01079 * 01080 * End of CBBCSD 01081 * 01082 END 01083