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