 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DBBCSD( 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, WORK, LWORK, 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, LWORK, M, P, Q 00018 * .. 00019 * .. Array Arguments .. 00020 DOUBLE PRECISION B11D( * ), B11E( * ), B12D( * ), B12E( * ), 00021 $ B21D( * ), B21E( * ), B22D( * ), B22E( * ), 00022 $ PHI( * ), THETA( * ), WORK( * ) 00023 DOUBLE PRECISION U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ), 00024 $ V2T( LDV2T, * ) 00025 * .. 00026 * 00027 * Purpose 00028 * ======= 00029 * 00030 * DBBCSD computes the CS decomposition of an orthogonal 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 | ]**T 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 DORCSD 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 orthogonal 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 orthogonal 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (LDV1T,Q) 00121 * On entry, a LDV1T-by-Q matrix. On exit, V1T is premultiplied 00122 * by the 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) DOUBLE PRECISION array, dimenison (LDV2T,M-Q) 00129 * On entry, a LDV2T-by-(M-Q) matrix. On exit, V2T is 00130 * premultiplied by the 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) DOUBLE PRECISION array, dimension (Q) 00138 * When DBBCSD converges, B11D contains the cosines of THETA(1), 00139 * ..., THETA(Q). If DBBCSD fails to converge, then B11D 00140 * contains the diagonal of the partially reduced top-left 00141 * block. 00142 * 00143 * B11E (output) DOUBLE PRECISION array, dimension (Q-1) 00144 * When DBBCSD converges, B11E contains zeros. If DBBCSD fails 00145 * to converge, then B11E contains the superdiagonal of the 00146 * partially reduced top-left block. 00147 * 00148 * B12D (output) DOUBLE PRECISION array, dimension (Q) 00149 * When DBBCSD converges, B12D contains the negative sines of 00150 * THETA(1), ..., THETA(Q). If DBBCSD fails to converge, then 00151 * B12D contains the diagonal of the partially reduced top-right 00152 * block. 00153 * 00154 * B12E (output) DOUBLE PRECISION array, dimension (Q-1) 00155 * When DBBCSD converges, B12E contains zeros. If DBBCSD fails 00156 * to converge, then B12E contains the subdiagonal of the 00157 * partially reduced top-right block. 00158 * 00159 * WORK (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 00160 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00161 * 00162 * LWORK (input) INTEGER 00163 * The dimension of the array WORK. LWORK >= MAX(1,8*Q). 00164 * 00165 * If LWORK = -1, then a workspace query is assumed; the 00166 * routine only calculates the optimal size of the WORK array, 00167 * returns this value as the first entry of the work array, and 00168 * no error message related to LWORK 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 DBBCSD 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 DOUBLE PRECISION, 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 DOUBLE PRECISION HUNDRED, MEIGHTH, ONE, PIOVER2, TEN, ZERO 00197 PARAMETER ( HUNDRED = 100.0D0, MEIGHTH = -0.125D0, 00198 $ ONE = 1.0D0, PIOVER2 = 1.57079632679489662D0, 00199 $ TEN = 10.0D0, ZERO = 0.0D0 ) 00200 DOUBLE PRECISION NEGONECOMPLEX 00201 PARAMETER ( NEGONECOMPLEX = -1.0D0 ) 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 $ LWORKMIN, LWORKOPT, MAXIT, MINI 00210 DOUBLE PRECISION 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 DLASR, DSCAL, DSWAP, DLARTGP, DLARTGS, DLAS2, 00217 $ XERBLA 00218 * .. 00219 * .. External Functions .. 00220 DOUBLE PRECISION DLAMCH 00221 LOGICAL LSAME 00222 EXTERNAL LSAME, DLAMCH 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 = LWORK .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 LWORKMIN = 1 00261 WORK(1) = LWORKMIN 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 LWORKOPT = IV2TSN + Q - 1 00277 LWORKMIN = LWORKOPT 00278 WORK(1) = LWORKOPT 00279 IF( LWORK .LT. LWORKMIN .AND. .NOT. LQUERY ) THEN 00280 INFO = -28 00281 END IF 00282 END IF 00283 * 00284 IF( INFO .NE. 0 ) THEN 00285 CALL XERBLA( 'DBBCSD', -INFO ) 00286 RETURN 00287 ELSE IF( LQUERY ) THEN 00288 RETURN 00289 END IF 00290 * 00291 * Get machine constants 00292 * 00293 EPS = DLAMCH( 'Epsilon' ) 00294 UNFL = DLAMCH( '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 DLAS2( B11D(IMAX-1), B11E(IMAX-1), B11D(IMAX), SIGMA11, 00401 $ DUMMY ) 00402 CALL DLAS2( 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 DLARTGS( B11D(IMIN), B11E(IMIN), MU, 00426 $ WORK(IV1TCS+IMIN-1), WORK(IV1TSN+IMIN-1) ) 00427 ELSE 00428 CALL DLARTGS( B21D(IMIN), B21E(IMIN), NU, 00429 $ WORK(IV1TCS+IMIN-1), WORK(IV1TSN+IMIN-1) ) 00430 END IF 00431 * 00432 TEMP = WORK(IV1TCS+IMIN-1)*B11D(IMIN) + 00433 $ WORK(IV1TSN+IMIN-1)*B11E(IMIN) 00434 B11E(IMIN) = WORK(IV1TCS+IMIN-1)*B11E(IMIN) - 00435 $ WORK(IV1TSN+IMIN-1)*B11D(IMIN) 00436 B11D(IMIN) = TEMP 00437 B11BULGE = WORK(IV1TSN+IMIN-1)*B11D(IMIN+1) 00438 B11D(IMIN+1) = WORK(IV1TCS+IMIN-1)*B11D(IMIN+1) 00439 TEMP = WORK(IV1TCS+IMIN-1)*B21D(IMIN) + 00440 $ WORK(IV1TSN+IMIN-1)*B21E(IMIN) 00441 B21E(IMIN) = WORK(IV1TCS+IMIN-1)*B21E(IMIN) - 00442 $ WORK(IV1TSN+IMIN-1)*B21D(IMIN) 00443 B21D(IMIN) = TEMP 00444 B21BULGE = WORK(IV1TSN+IMIN-1)*B21D(IMIN+1) 00445 B21D(IMIN+1) = WORK(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 DLARTGP( B11BULGE, B11D(IMIN), WORK(IU1SN+IMIN-1), 00456 $ WORK(IU1CS+IMIN-1), R ) 00457 ELSE IF( MU .LE. NU ) THEN 00458 CALL DLARTGS( B11E( IMIN ), B11D( IMIN + 1 ), MU, 00459 $ WORK(IU1CS+IMIN-1), WORK(IU1SN+IMIN-1) ) 00460 ELSE 00461 CALL DLARTGS( B12D( IMIN ), B12E( IMIN ), NU, 00462 $ WORK(IU1CS+IMIN-1), WORK(IU1SN+IMIN-1) ) 00463 END IF 00464 IF( B21D(IMIN)**2+B21BULGE**2 .GT. THRESH**2 ) THEN 00465 CALL DLARTGP( B21BULGE, B21D(IMIN), WORK(IU2SN+IMIN-1), 00466 $ WORK(IU2CS+IMIN-1), R ) 00467 ELSE IF( NU .LT. MU ) THEN 00468 CALL DLARTGS( B21E( IMIN ), B21D( IMIN + 1 ), NU, 00469 $ WORK(IU2CS+IMIN-1), WORK(IU2SN+IMIN-1) ) 00470 ELSE 00471 CALL DLARTGS( B22D(IMIN), B22E(IMIN), MU, 00472 $ WORK(IU2CS+IMIN-1), WORK(IU2SN+IMIN-1) ) 00473 END IF 00474 WORK(IU2CS+IMIN-1) = -WORK(IU2CS+IMIN-1) 00475 WORK(IU2SN+IMIN-1) = -WORK(IU2SN+IMIN-1) 00476 * 00477 TEMP = WORK(IU1CS+IMIN-1)*B11E(IMIN) + 00478 $ WORK(IU1SN+IMIN-1)*B11D(IMIN+1) 00479 B11D(IMIN+1) = WORK(IU1CS+IMIN-1)*B11D(IMIN+1) - 00480 $ WORK(IU1SN+IMIN-1)*B11E(IMIN) 00481 B11E(IMIN) = TEMP 00482 IF( IMAX .GT. IMIN+1 ) THEN 00483 B11BULGE = WORK(IU1SN+IMIN-1)*B11E(IMIN+1) 00484 B11E(IMIN+1) = WORK(IU1CS+IMIN-1)*B11E(IMIN+1) 00485 END IF 00486 TEMP = WORK(IU1CS+IMIN-1)*B12D(IMIN) + 00487 $ WORK(IU1SN+IMIN-1)*B12E(IMIN) 00488 B12E(IMIN) = WORK(IU1CS+IMIN-1)*B12E(IMIN) - 00489 $ WORK(IU1SN+IMIN-1)*B12D(IMIN) 00490 B12D(IMIN) = TEMP 00491 B12BULGE = WORK(IU1SN+IMIN-1)*B12D(IMIN+1) 00492 B12D(IMIN+1) = WORK(IU1CS+IMIN-1)*B12D(IMIN+1) 00493 TEMP = WORK(IU2CS+IMIN-1)*B21E(IMIN) + 00494 $ WORK(IU2SN+IMIN-1)*B21D(IMIN+1) 00495 B21D(IMIN+1) = WORK(IU2CS+IMIN-1)*B21D(IMIN+1) - 00496 $ WORK(IU2SN+IMIN-1)*B21E(IMIN) 00497 B21E(IMIN) = TEMP 00498 IF( IMAX .GT. IMIN+1 ) THEN 00499 B21BULGE = WORK(IU2SN+IMIN-1)*B21E(IMIN+1) 00500 B21E(IMIN+1) = WORK(IU2CS+IMIN-1)*B21E(IMIN+1) 00501 END IF 00502 TEMP = WORK(IU2CS+IMIN-1)*B22D(IMIN) + 00503 $ WORK(IU2SN+IMIN-1)*B22E(IMIN) 00504 B22E(IMIN) = WORK(IU2CS+IMIN-1)*B22E(IMIN) - 00505 $ WORK(IU2SN+IMIN-1)*B22D(IMIN) 00506 B22D(IMIN) = TEMP 00507 B22BULGE = WORK(IU2SN+IMIN-1)*B22D(IMIN+1) 00508 B22D(IMIN+1) = WORK(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 DLARTGP( X2, X1, WORK(IV1TSN+I-1), WORK(IV1TCS+I-1), 00539 $ R ) 00540 ELSE IF( .NOT. RESTART11 .AND. RESTART21 ) THEN 00541 CALL DLARTGP( B11BULGE, B11E(I-1), WORK(IV1TSN+I-1), 00542 $ WORK(IV1TCS+I-1), R ) 00543 ELSE IF( RESTART11 .AND. .NOT. RESTART21 ) THEN 00544 CALL DLARTGP( B21BULGE, B21E(I-1), WORK(IV1TSN+I-1), 00545 $ WORK(IV1TCS+I-1), R ) 00546 ELSE IF( MU .LE. NU ) THEN 00547 CALL DLARTGS( B11D(I), B11E(I), MU, WORK(IV1TCS+I-1), 00548 $ WORK(IV1TSN+I-1) ) 00549 ELSE 00550 CALL DLARTGS( B21D(I), B21E(I), NU, WORK(IV1TCS+I-1), 00551 $ WORK(IV1TSN+I-1) ) 00552 END IF 00553 WORK(IV1TCS+I-1) = -WORK(IV1TCS+I-1) 00554 WORK(IV1TSN+I-1) = -WORK(IV1TSN+I-1) 00555 IF( .NOT. RESTART12 .AND. .NOT. RESTART22 ) THEN 00556 CALL DLARTGP( Y2, Y1, WORK(IV2TSN+I-1-1), 00557 $ WORK(IV2TCS+I-1-1), R ) 00558 ELSE IF( .NOT. RESTART12 .AND. RESTART22 ) THEN 00559 CALL DLARTGP( B12BULGE, B12D(I-1), WORK(IV2TSN+I-1-1), 00560 $ WORK(IV2TCS+I-1-1), R ) 00561 ELSE IF( RESTART12 .AND. .NOT. RESTART22 ) THEN 00562 CALL DLARTGP( B22BULGE, B22D(I-1), WORK(IV2TSN+I-1-1), 00563 $ WORK(IV2TCS+I-1-1), R ) 00564 ELSE IF( NU .LT. MU ) THEN 00565 CALL DLARTGS( B12E(I-1), B12D(I), NU, WORK(IV2TCS+I-1-1), 00566 $ WORK(IV2TSN+I-1-1) ) 00567 ELSE 00568 CALL DLARTGS( B22E(I-1), B22D(I), MU, WORK(IV2TCS+I-1-1), 00569 $ WORK(IV2TSN+I-1-1) ) 00570 END IF 00571 * 00572 TEMP = WORK(IV1TCS+I-1)*B11D(I) + WORK(IV1TSN+I-1)*B11E(I) 00573 B11E(I) = WORK(IV1TCS+I-1)*B11E(I) - 00574 $ WORK(IV1TSN+I-1)*B11D(I) 00575 B11D(I) = TEMP 00576 B11BULGE = WORK(IV1TSN+I-1)*B11D(I+1) 00577 B11D(I+1) = WORK(IV1TCS+I-1)*B11D(I+1) 00578 TEMP = WORK(IV1TCS+I-1)*B21D(I) + WORK(IV1TSN+I-1)*B21E(I) 00579 B21E(I) = WORK(IV1TCS+I-1)*B21E(I) - 00580 $ WORK(IV1TSN+I-1)*B21D(I) 00581 B21D(I) = TEMP 00582 B21BULGE = WORK(IV1TSN+I-1)*B21D(I+1) 00583 B21D(I+1) = WORK(IV1TCS+I-1)*B21D(I+1) 00584 TEMP = WORK(IV2TCS+I-1-1)*B12E(I-1) + 00585 $ WORK(IV2TSN+I-1-1)*B12D(I) 00586 B12D(I) = WORK(IV2TCS+I-1-1)*B12D(I) - 00587 $ WORK(IV2TSN+I-1-1)*B12E(I-1) 00588 B12E(I-1) = TEMP 00589 B12BULGE = WORK(IV2TSN+I-1-1)*B12E(I) 00590 B12E(I) = WORK(IV2TCS+I-1-1)*B12E(I) 00591 TEMP = WORK(IV2TCS+I-1-1)*B22E(I-1) + 00592 $ WORK(IV2TSN+I-1-1)*B22D(I) 00593 B22D(I) = WORK(IV2TCS+I-1-1)*B22D(I) - 00594 $ WORK(IV2TSN+I-1-1)*B22E(I-1) 00595 B22E(I-1) = TEMP 00596 B22BULGE = WORK(IV2TSN+I-1-1)*B22E(I) 00597 B22E(I) = WORK(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 DLARTGP( X2, X1, WORK(IU1SN+I-1), WORK(IU1CS+I-1), 00622 $ R ) 00623 ELSE IF( .NOT. RESTART11 .AND. RESTART12 ) THEN 00624 CALL DLARTGP( B11BULGE, B11D(I), WORK(IU1SN+I-1), 00625 $ WORK(IU1CS+I-1), R ) 00626 ELSE IF( RESTART11 .AND. .NOT. RESTART12 ) THEN 00627 CALL DLARTGP( B12BULGE, B12E(I-1), WORK(IU1SN+I-1), 00628 $ WORK(IU1CS+I-1), R ) 00629 ELSE IF( MU .LE. NU ) THEN 00630 CALL DLARTGS( B11E(I), B11D(I+1), MU, WORK(IU1CS+I-1), 00631 $ WORK(IU1SN+I-1) ) 00632 ELSE 00633 CALL DLARTGS( B12D(I), B12E(I), NU, WORK(IU1CS+I-1), 00634 $ WORK(IU1SN+I-1) ) 00635 END IF 00636 IF( .NOT. RESTART21 .AND. .NOT. RESTART22 ) THEN 00637 CALL DLARTGP( Y2, Y1, WORK(IU2SN+I-1), WORK(IU2CS+I-1), 00638 $ R ) 00639 ELSE IF( .NOT. RESTART21 .AND. RESTART22 ) THEN 00640 CALL DLARTGP( B21BULGE, B21D(I), WORK(IU2SN+I-1), 00641 $ WORK(IU2CS+I-1), R ) 00642 ELSE IF( RESTART21 .AND. .NOT. RESTART22 ) THEN 00643 CALL DLARTGP( B22BULGE, B22E(I-1), WORK(IU2SN+I-1), 00644 $ WORK(IU2CS+I-1), R ) 00645 ELSE IF( NU .LT. MU ) THEN 00646 CALL DLARTGS( B21E(I), B21E(I+1), NU, WORK(IU2CS+I-1), 00647 $ WORK(IU2SN+I-1) ) 00648 ELSE 00649 CALL DLARTGS( B22D(I), B22E(I), MU, WORK(IU2CS+I-1), 00650 $ WORK(IU2SN+I-1) ) 00651 END IF 00652 WORK(IU2CS+I-1) = -WORK(IU2CS+I-1) 00653 WORK(IU2SN+I-1) = -WORK(IU2SN+I-1) 00654 * 00655 TEMP = WORK(IU1CS+I-1)*B11E(I) + WORK(IU1SN+I-1)*B11D(I+1) 00656 B11D(I+1) = WORK(IU1CS+I-1)*B11D(I+1) - 00657 $ WORK(IU1SN+I-1)*B11E(I) 00658 B11E(I) = TEMP 00659 IF( I .LT. IMAX - 1 ) THEN 00660 B11BULGE = WORK(IU1SN+I-1)*B11E(I+1) 00661 B11E(I+1) = WORK(IU1CS+I-1)*B11E(I+1) 00662 END IF 00663 TEMP = WORK(IU2CS+I-1)*B21E(I) + WORK(IU2SN+I-1)*B21D(I+1) 00664 B21D(I+1) = WORK(IU2CS+I-1)*B21D(I+1) - 00665 $ WORK(IU2SN+I-1)*B21E(I) 00666 B21E(I) = TEMP 00667 IF( I .LT. IMAX - 1 ) THEN 00668 B21BULGE = WORK(IU2SN+I-1)*B21E(I+1) 00669 B21E(I+1) = WORK(IU2CS+I-1)*B21E(I+1) 00670 END IF 00671 TEMP = WORK(IU1CS+I-1)*B12D(I) + WORK(IU1SN+I-1)*B12E(I) 00672 B12E(I) = WORK(IU1CS+I-1)*B12E(I) - WORK(IU1SN+I-1)*B12D(I) 00673 B12D(I) = TEMP 00674 B12BULGE = WORK(IU1SN+I-1)*B12D(I+1) 00675 B12D(I+1) = WORK(IU1CS+I-1)*B12D(I+1) 00676 TEMP = WORK(IU2CS+I-1)*B22D(I) + WORK(IU2SN+I-1)*B22E(I) 00677 B22E(I) = WORK(IU2CS+I-1)*B22E(I) - WORK(IU2SN+I-1)*B22D(I) 00678 B22D(I) = TEMP 00679 B22BULGE = WORK(IU2SN+I-1)*B22D(I+1) 00680 B22D(I+1) = WORK(IU2CS+I-1)*B22D(I+1) 00681 * 00682 END DO 00683 * 00684 * Compute PHI(IMAX-1) 00685 * 00686 X1 = SIN(THETA(IMAX-1))*B11E(IMAX-1) + 00687 $ COS(THETA(IMAX-1))*B21E(IMAX-1) 00688 Y1 = SIN(THETA(IMAX-1))*B12D(IMAX-1) + 00689 $ COS(THETA(IMAX-1))*B22D(IMAX-1) 00690 Y2 = SIN(THETA(IMAX-1))*B12BULGE + COS(THETA(IMAX-1))*B22BULGE 00691 * 00692 PHI(IMAX-1) = ATAN2( ABS(X1), SQRT(Y1**2+Y2**2) ) 00693 * 00694 * Chase bulges from B12(IMAX-1,IMAX) and B22(IMAX-1,IMAX) 00695 * 00696 RESTART12 = B12D(IMAX-1)**2 + B12BULGE**2 .LE. THRESH**2 00697 RESTART22 = B22D(IMAX-1)**2 + B22BULGE**2 .LE. THRESH**2 00698 * 00699 IF( .NOT. RESTART12 .AND. .NOT. RESTART22 ) THEN 00700 CALL DLARTGP( Y2, Y1, WORK(IV2TSN+IMAX-1-1), 00701 $ WORK(IV2TCS+IMAX-1-1), R ) 00702 ELSE IF( .NOT. RESTART12 .AND. RESTART22 ) THEN 00703 CALL DLARTGP( B12BULGE, B12D(IMAX-1), WORK(IV2TSN+IMAX-1-1), 00704 $ WORK(IV2TCS+IMAX-1-1), R ) 00705 ELSE IF( RESTART12 .AND. .NOT. RESTART22 ) THEN 00706 CALL DLARTGP( B22BULGE, B22D(IMAX-1), WORK(IV2TSN+IMAX-1-1), 00707 $ WORK(IV2TCS+IMAX-1-1), R ) 00708 ELSE IF( NU .LT. MU ) THEN 00709 CALL DLARTGS( B12E(IMAX-1), B12D(IMAX), NU, 00710 $ WORK(IV2TCS+IMAX-1-1), WORK(IV2TSN+IMAX-1-1) ) 00711 ELSE 00712 CALL DLARTGS( B22E(IMAX-1), B22D(IMAX), MU, 00713 $ WORK(IV2TCS+IMAX-1-1), WORK(IV2TSN+IMAX-1-1) ) 00714 END IF 00715 * 00716 TEMP = WORK(IV2TCS+IMAX-1-1)*B12E(IMAX-1) + 00717 $ WORK(IV2TSN+IMAX-1-1)*B12D(IMAX) 00718 B12D(IMAX) = WORK(IV2TCS+IMAX-1-1)*B12D(IMAX) - 00719 $ WORK(IV2TSN+IMAX-1-1)*B12E(IMAX-1) 00720 B12E(IMAX-1) = TEMP 00721 TEMP = WORK(IV2TCS+IMAX-1-1)*B22E(IMAX-1) + 00722 $ WORK(IV2TSN+IMAX-1-1)*B22D(IMAX) 00723 B22D(IMAX) = WORK(IV2TCS+IMAX-1-1)*B22D(IMAX) - 00724 $ WORK(IV2TSN+IMAX-1-1)*B22E(IMAX-1) 00725 B22E(IMAX-1) = TEMP 00726 * 00727 * Update singular vectors 00728 * 00729 IF( WANTU1 ) THEN 00730 IF( COLMAJOR ) THEN 00731 CALL DLASR( 'R', 'V', 'F', P, IMAX-IMIN+1, 00732 $ WORK(IU1CS+IMIN-1), WORK(IU1SN+IMIN-1), 00733 $ U1(1,IMIN), LDU1 ) 00734 ELSE 00735 CALL DLASR( 'L', 'V', 'F', IMAX-IMIN+1, P, 00736 $ WORK(IU1CS+IMIN-1), WORK(IU1SN+IMIN-1), 00737 $ U1(IMIN,1), LDU1 ) 00738 END IF 00739 END IF 00740 IF( WANTU2 ) THEN 00741 IF( COLMAJOR ) THEN 00742 CALL DLASR( 'R', 'V', 'F', M-P, IMAX-IMIN+1, 00743 $ WORK(IU2CS+IMIN-1), WORK(IU2SN+IMIN-1), 00744 $ U2(1,IMIN), LDU2 ) 00745 ELSE 00746 CALL DLASR( 'L', 'V', 'F', IMAX-IMIN+1, M-P, 00747 $ WORK(IU2CS+IMIN-1), WORK(IU2SN+IMIN-1), 00748 $ U2(IMIN,1), LDU2 ) 00749 END IF 00750 END IF 00751 IF( WANTV1T ) THEN 00752 IF( COLMAJOR ) THEN 00753 CALL DLASR( 'L', 'V', 'F', IMAX-IMIN+1, Q, 00754 $ WORK(IV1TCS+IMIN-1), WORK(IV1TSN+IMIN-1), 00755 $ V1T(IMIN,1), LDV1T ) 00756 ELSE 00757 CALL DLASR( 'R', 'V', 'F', Q, IMAX-IMIN+1, 00758 $ WORK(IV1TCS+IMIN-1), WORK(IV1TSN+IMIN-1), 00759 $ V1T(1,IMIN), LDV1T ) 00760 END IF 00761 END IF 00762 IF( WANTV2T ) THEN 00763 IF( COLMAJOR ) THEN 00764 CALL DLASR( 'L', 'V', 'F', IMAX-IMIN+1, M-Q, 00765 $ WORK(IV2TCS+IMIN-1), WORK(IV2TSN+IMIN-1), 00766 $ V2T(IMIN,1), LDV2T ) 00767 ELSE 00768 CALL DLASR( 'R', 'V', 'F', M-Q, IMAX-IMIN+1, 00769 $ WORK(IV2TCS+IMIN-1), WORK(IV2TSN+IMIN-1), 00770 $ V2T(1,IMIN), LDV2T ) 00771 END IF 00772 END IF 00773 * 00774 * Fix signs on B11(IMAX-1,IMAX) and B21(IMAX-1,IMAX) 00775 * 00776 IF( B11E(IMAX-1)+B21E(IMAX-1) .GT. 0 ) THEN 00777 B11D(IMAX) = -B11D(IMAX) 00778 B21D(IMAX) = -B21D(IMAX) 00779 IF( WANTV1T ) THEN 00780 IF( COLMAJOR ) THEN 00781 CALL DSCAL( Q, NEGONECOMPLEX, V1T(IMAX,1), LDV1T ) 00782 ELSE 00783 CALL DSCAL( Q, NEGONECOMPLEX, V1T(1,IMAX), 1 ) 00784 END IF 00785 END IF 00786 END IF 00787 * 00788 * Compute THETA(IMAX) 00789 * 00790 X1 = COS(PHI(IMAX-1))*B11D(IMAX) + 00791 $ SIN(PHI(IMAX-1))*B12E(IMAX-1) 00792 Y1 = COS(PHI(IMAX-1))*B21D(IMAX) + 00793 $ SIN(PHI(IMAX-1))*B22E(IMAX-1) 00794 * 00795 THETA(IMAX) = ATAN2( ABS(Y1), ABS(X1) ) 00796 * 00797 * Fix signs on B11(IMAX,IMAX), B12(IMAX,IMAX-1), B21(IMAX,IMAX), 00798 * and B22(IMAX,IMAX-1) 00799 * 00800 IF( B11D(IMAX)+B12E(IMAX-1) .LT. 0 ) THEN 00801 B12D(IMAX) = -B12D(IMAX) 00802 IF( WANTU1 ) THEN 00803 IF( COLMAJOR ) THEN 00804 CALL DSCAL( P, NEGONECOMPLEX, U1(1,IMAX), 1 ) 00805 ELSE 00806 CALL DSCAL( P, NEGONECOMPLEX, U1(IMAX,1), LDU1 ) 00807 END IF 00808 END IF 00809 END IF 00810 IF( B21D(IMAX)+B22E(IMAX-1) .GT. 0 ) THEN 00811 B22D(IMAX) = -B22D(IMAX) 00812 IF( WANTU2 ) THEN 00813 IF( COLMAJOR ) THEN 00814 CALL DSCAL( M-P, NEGONECOMPLEX, U2(1,IMAX), 1 ) 00815 ELSE 00816 CALL DSCAL( M-P, NEGONECOMPLEX, U2(IMAX,1), LDU2 ) 00817 END IF 00818 END IF 00819 END IF 00820 * 00821 * Fix signs on B12(IMAX,IMAX) and B22(IMAX,IMAX) 00822 * 00823 IF( B12D(IMAX)+B22D(IMAX) .LT. 0 ) THEN 00824 IF( WANTV2T ) THEN 00825 IF( COLMAJOR ) THEN 00826 CALL DSCAL( M-Q, NEGONECOMPLEX, V2T(IMAX,1), LDV2T ) 00827 ELSE 00828 CALL DSCAL( M-Q, NEGONECOMPLEX, V2T(1,IMAX), 1 ) 00829 END IF 00830 END IF 00831 END IF 00832 * 00833 * Test for negligible sines or cosines 00834 * 00835 DO I = IMIN, IMAX 00836 IF( THETA(I) .LT. THRESH ) THEN 00837 THETA(I) = ZERO 00838 ELSE IF( THETA(I) .GT. PIOVER2-THRESH ) THEN 00839 THETA(I) = PIOVER2 00840 END IF 00841 END DO 00842 DO I = IMIN, IMAX-1 00843 IF( PHI(I) .LT. THRESH ) THEN 00844 PHI(I) = ZERO 00845 ELSE IF( PHI(I) .GT. PIOVER2-THRESH ) THEN 00846 PHI(I) = PIOVER2 00847 END IF 00848 END DO 00849 * 00850 * Deflate 00851 * 00852 IF (IMAX .GT. 1) THEN 00853 DO WHILE( PHI(IMAX-1) .EQ. ZERO ) 00854 IMAX = IMAX - 1 00855 IF (IMAX .LE. 1) EXIT 00856 END DO 00857 END IF 00858 IF( IMIN .GT. IMAX - 1 ) 00859 $ IMIN = IMAX - 1 00860 IF (IMIN .GT. 1) THEN 00861 DO WHILE (PHI(IMIN-1) .NE. ZERO) 00862 IMIN = IMIN - 1 00863 IF (IMIN .LE. 1) EXIT 00864 END DO 00865 END IF 00866 * 00867 * Repeat main iteration loop 00868 * 00869 END DO 00870 * 00871 * Postprocessing: order THETA from least to greatest 00872 * 00873 DO I = 1, Q 00874 * 00875 MINI = I 00876 THETAMIN = THETA(I) 00877 DO J = I+1, Q 00878 IF( THETA(J) .LT. THETAMIN ) THEN 00879 MINI = J 00880 THETAMIN = THETA(J) 00881 END IF 00882 END DO 00883 * 00884 IF( MINI .NE. I ) THEN 00885 THETA(MINI) = THETA(I) 00886 THETA(I) = THETAMIN 00887 IF( COLMAJOR ) THEN 00888 IF( WANTU1 ) 00889 $ CALL DSWAP( P, U1(1,I), 1, U1(1,MINI), 1 ) 00890 IF( WANTU2 ) 00891 $ CALL DSWAP( M-P, U2(1,I), 1, U2(1,MINI), 1 ) 00892 IF( WANTV1T ) 00893 $ CALL DSWAP( Q, V1T(I,1), LDV1T, V1T(MINI,1), LDV1T ) 00894 IF( WANTV2T ) 00895 $ CALL DSWAP( M-Q, V2T(I,1), LDV2T, V2T(MINI,1), 00896 $ LDV2T ) 00897 ELSE 00898 IF( WANTU1 ) 00899 $ CALL DSWAP( P, U1(I,1), LDU1, U1(MINI,1), LDU1 ) 00900 IF( WANTU2 ) 00901 $ CALL DSWAP( M-P, U2(I,1), LDU2, U2(MINI,1), LDU2 ) 00902 IF( WANTV1T ) 00903 $ CALL DSWAP( Q, V1T(1,I), 1, V1T(1,MINI), 1 ) 00904 IF( WANTV2T ) 00905 $ CALL DSWAP( M-Q, V2T(1,I), 1, V2T(1,MINI), 1 ) 00906 END IF 00907 END IF 00908 * 00909 END DO 00910 * 00911 RETURN 00912 * 00913 * End of DBBCSD 00914 * 00915 END 00916
1.7.3 
