 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE ZBBCSD( 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 DOUBLE PRECISION B11D( * ), B11E( * ), B12D( * ), B12E( * ), 00021 $ B21D( * ), B21E( * ), B22D( * ), B22E( * ), 00022 $ PHI( * ), THETA( * ), RWORK( * ) 00023 COMPLEX*16 U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ), 00024 $ V2T( LDV2T, * ) 00025 * .. 00026 * 00027 * Purpose 00028 * ======= 00029 * 00030 * ZBBCSD 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 ZUNCSD 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) 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) COMPLEX*16 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*16 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*16 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*16 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) DOUBLE PRECISION array, dimension (Q) 00138 * When ZBBCSD converges, B11D contains the cosines of THETA(1), 00139 * ..., THETA(Q). If ZBBCSD 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 ZBBCSD converges, B11E contains zeros. If ZBBCSD 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 ZBBCSD converges, B12D contains the negative sines of 00150 * THETA(1), ..., THETA(Q). If ZBBCSD 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 ZBBCSD converges, B12E contains zeros. If ZBBCSD fails 00156 * to converge, then B12E contains the subdiagonal of the 00157 * partially reduced top-right block. 00158 * 00159 * RWORK (workspace) DOUBLE PRECISION 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 ZBBCSD 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 COMPLEX*16 NEGONECOMPLEX 00201 PARAMETER ( NEGONECOMPLEX = (-1.0D0,0.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 $ LRWORKMIN, LRWORKOPT, 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 DLARTGP, DLARTGS, DLAS2, XERBLA, ZLASR, ZSCAL, 00216 $ ZSWAP 00217 * .. 00218 * .. External Functions .. 00219 DOUBLE PRECISION DLAMCH 00220 LOGICAL LSAME 00221 EXTERNAL LSAME, DLAMCH 00222 * .. 00223 * .. Intrinsic Functions .. 00224 INTRINSIC ABS, ATAN2, COS, MAX, MIN, SIN, SQRT 00225 * .. 00226 * .. Executable Statements .. 00227 * 00228 * Test input arguments 00229 * 00230 INFO = 0 00231 LQUERY = LRWORK .EQ. -1 00232 WANTU1 = LSAME( JOBU1, 'Y' ) 00233 WANTU2 = LSAME( JOBU2, 'Y' ) 00234 WANTV1T = LSAME( JOBV1T, 'Y' ) 00235 WANTV2T = LSAME( JOBV2T, 'Y' ) 00236 COLMAJOR = .NOT. LSAME( TRANS, 'T' ) 00237 * 00238 IF( M .LT. 0 ) THEN 00239 INFO = -6 00240 ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN 00241 INFO = -7 00242 ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN 00243 INFO = -8 00244 ELSE IF( Q .GT. P .OR. Q .GT. M-P .OR. Q .GT. M-Q ) THEN 00245 INFO = -8 00246 ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN 00247 INFO = -12 00248 ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN 00249 INFO = -14 00250 ELSE IF( WANTV1T .AND. LDV1T .LT. Q ) THEN 00251 INFO = -16 00252 ELSE IF( WANTV2T .AND. LDV2T .LT. M-Q ) THEN 00253 INFO = -18 00254 END IF 00255 * 00256 * Quick return if Q = 0 00257 * 00258 IF( INFO .EQ. 0 .AND. Q .EQ. 0 ) THEN 00259 LRWORKMIN = 1 00260 RWORK(1) = LRWORKMIN 00261 RETURN 00262 END IF 00263 * 00264 * Compute workspace 00265 * 00266 IF( INFO .EQ. 0 ) THEN 00267 IU1CS = 1 00268 IU1SN = IU1CS + Q 00269 IU2CS = IU1SN + Q 00270 IU2SN = IU2CS + Q 00271 IV1TCS = IU2SN + Q 00272 IV1TSN = IV1TCS + Q 00273 IV2TCS = IV1TSN + Q 00274 IV2TSN = IV2TCS + Q 00275 LRWORKOPT = IV2TSN + Q - 1 00276 LRWORKMIN = LRWORKOPT 00277 RWORK(1) = LRWORKOPT 00278 IF( LRWORK .LT. LRWORKMIN .AND. .NOT. LQUERY ) THEN 00279 INFO = -28 00280 END IF 00281 END IF 00282 * 00283 IF( INFO .NE. 0 ) THEN 00284 CALL XERBLA( 'ZBBCSD', -INFO ) 00285 RETURN 00286 ELSE IF( LQUERY ) THEN 00287 RETURN 00288 END IF 00289 * 00290 * Get machine constants 00291 * 00292 EPS = DLAMCH( 'Epsilon' ) 00293 UNFL = DLAMCH( 'Safe minimum' ) 00294 TOLMUL = MAX( TEN, MIN( HUNDRED, EPS**MEIGHTH ) ) 00295 TOL = TOLMUL*EPS 00296 THRESH = MAX( TOL, MAXITR*Q*Q*UNFL ) 00297 * 00298 * Test for negligible sines or cosines 00299 * 00300 DO I = 1, Q 00301 IF( THETA(I) .LT. THRESH ) THEN 00302 THETA(I) = ZERO 00303 ELSE IF( THETA(I) .GT. PIOVER2-THRESH ) THEN 00304 THETA(I) = PIOVER2 00305 END IF 00306 END DO 00307 DO I = 1, Q-1 00308 IF( PHI(I) .LT. THRESH ) THEN 00309 PHI(I) = ZERO 00310 ELSE IF( PHI(I) .GT. PIOVER2-THRESH ) THEN 00311 PHI(I) = PIOVER2 00312 END IF 00313 END DO 00314 * 00315 * Initial deflation 00316 * 00317 IMAX = Q 00318 DO WHILE( ( IMAX .GT. 1 ) .AND. ( PHI(IMAX-1) .EQ. ZERO ) ) 00319 IMAX = IMAX - 1 00320 END DO 00321 IMIN = IMAX - 1 00322 IF ( IMIN .GT. 1 ) THEN 00323 DO WHILE( PHI(IMIN-1) .NE. ZERO ) 00324 IMIN = IMIN - 1 00325 IF ( IMIN .LE. 1 ) EXIT 00326 END DO 00327 END IF 00328 * 00329 * Initialize iteration counter 00330 * 00331 MAXIT = MAXITR*Q*Q 00332 ITER = 0 00333 * 00334 * Begin main iteration loop 00335 * 00336 DO WHILE( IMAX .GT. 1 ) 00337 * 00338 * Compute the matrix entries 00339 * 00340 B11D(IMIN) = COS( THETA(IMIN) ) 00341 B21D(IMIN) = -SIN( THETA(IMIN) ) 00342 DO I = IMIN, IMAX - 1 00343 B11E(I) = -SIN( THETA(I) ) * SIN( PHI(I) ) 00344 B11D(I+1) = COS( THETA(I+1) ) * COS( PHI(I) ) 00345 B12D(I) = SIN( THETA(I) ) * COS( PHI(I) ) 00346 B12E(I) = COS( THETA(I+1) ) * SIN( PHI(I) ) 00347 B21E(I) = -COS( THETA(I) ) * SIN( PHI(I) ) 00348 B21D(I+1) = -SIN( THETA(I+1) ) * COS( PHI(I) ) 00349 B22D(I) = COS( THETA(I) ) * COS( PHI(I) ) 00350 B22E(I) = -SIN( THETA(I+1) ) * SIN( PHI(I) ) 00351 END DO 00352 B12D(IMAX) = SIN( THETA(IMAX) ) 00353 B22D(IMAX) = COS( THETA(IMAX) ) 00354 * 00355 * Abort if not converging; otherwise, increment ITER 00356 * 00357 IF( ITER .GT. MAXIT ) THEN 00358 INFO = 0 00359 DO I = 1, Q 00360 IF( PHI(I) .NE. ZERO ) 00361 $ INFO = INFO + 1 00362 END DO 00363 RETURN 00364 END IF 00365 * 00366 ITER = ITER + IMAX - IMIN 00367 * 00368 * Compute shifts 00369 * 00370 THETAMAX = THETA(IMIN) 00371 THETAMIN = THETA(IMIN) 00372 DO I = IMIN+1, IMAX 00373 IF( THETA(I) > THETAMAX ) 00374 $ THETAMAX = THETA(I) 00375 IF( THETA(I) < THETAMIN ) 00376 $ THETAMIN = THETA(I) 00377 END DO 00378 * 00379 IF( THETAMAX .GT. PIOVER2 - THRESH ) THEN 00380 * 00381 * Zero on diagonals of B11 and B22; induce deflation with a 00382 * zero shift 00383 * 00384 MU = ZERO 00385 NU = ONE 00386 * 00387 ELSE IF( THETAMIN .LT. THRESH ) THEN 00388 * 00389 * Zero on diagonals of B12 and B22; induce deflation with a 00390 * zero shift 00391 * 00392 MU = ONE 00393 NU = ZERO 00394 * 00395 ELSE 00396 * 00397 * Compute shifts for B11 and B21 and use the lesser 00398 * 00399 CALL DLAS2( B11D(IMAX-1), B11E(IMAX-1), B11D(IMAX), SIGMA11, 00400 $ DUMMY ) 00401 CALL DLAS2( B21D(IMAX-1), B21E(IMAX-1), B21D(IMAX), SIGMA21, 00402 $ DUMMY ) 00403 * 00404 IF( SIGMA11 .LE. SIGMA21 ) THEN 00405 MU = SIGMA11 00406 NU = SQRT( ONE - MU**2 ) 00407 IF( MU .LT. THRESH ) THEN 00408 MU = ZERO 00409 NU = ONE 00410 END IF 00411 ELSE 00412 NU = SIGMA21 00413 MU = SQRT( 1.0 - NU**2 ) 00414 IF( NU .LT. THRESH ) THEN 00415 MU = ONE 00416 NU = ZERO 00417 END IF 00418 END IF 00419 END IF 00420 * 00421 * Rotate to produce bulges in B11 and B21 00422 * 00423 IF( MU .LE. NU ) THEN 00424 CALL DLARTGS( B11D(IMIN), B11E(IMIN), MU, 00425 $ RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1) ) 00426 ELSE 00427 CALL DLARTGS( B21D(IMIN), B21E(IMIN), NU, 00428 $ RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1) ) 00429 END IF 00430 * 00431 TEMP = RWORK(IV1TCS+IMIN-1)*B11D(IMIN) + 00432 $ RWORK(IV1TSN+IMIN-1)*B11E(IMIN) 00433 B11E(IMIN) = RWORK(IV1TCS+IMIN-1)*B11E(IMIN) - 00434 $ RWORK(IV1TSN+IMIN-1)*B11D(IMIN) 00435 B11D(IMIN) = TEMP 00436 B11BULGE = RWORK(IV1TSN+IMIN-1)*B11D(IMIN+1) 00437 B11D(IMIN+1) = RWORK(IV1TCS+IMIN-1)*B11D(IMIN+1) 00438 TEMP = RWORK(IV1TCS+IMIN-1)*B21D(IMIN) + 00439 $ RWORK(IV1TSN+IMIN-1)*B21E(IMIN) 00440 B21E(IMIN) = RWORK(IV1TCS+IMIN-1)*B21E(IMIN) - 00441 $ RWORK(IV1TSN+IMIN-1)*B21D(IMIN) 00442 B21D(IMIN) = TEMP 00443 B21BULGE = RWORK(IV1TSN+IMIN-1)*B21D(IMIN+1) 00444 B21D(IMIN+1) = RWORK(IV1TCS+IMIN-1)*B21D(IMIN+1) 00445 * 00446 * Compute THETA(IMIN) 00447 * 00448 THETA( IMIN ) = ATAN2( SQRT( B21D(IMIN)**2+B21BULGE**2 ), 00449 $ SQRT( B11D(IMIN)**2+B11BULGE**2 ) ) 00450 * 00451 * Chase the bulges in B11(IMIN+1,IMIN) and B21(IMIN+1,IMIN) 00452 * 00453 IF( B11D(IMIN)**2+B11BULGE**2 .GT. THRESH**2 ) THEN 00454 CALL DLARTGP( B11BULGE, B11D(IMIN), RWORK(IU1SN+IMIN-1), 00455 $ RWORK(IU1CS+IMIN-1), R ) 00456 ELSE IF( MU .LE. NU ) THEN 00457 CALL DLARTGS( B11E( IMIN ), B11D( IMIN + 1 ), MU, 00458 $ RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1) ) 00459 ELSE 00460 CALL DLARTGS( B12D( IMIN ), B12E( IMIN ), NU, 00461 $ RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1) ) 00462 END IF 00463 IF( B21D(IMIN)**2+B21BULGE**2 .GT. THRESH**2 ) THEN 00464 CALL DLARTGP( B21BULGE, B21D(IMIN), RWORK(IU2SN+IMIN-1), 00465 $ RWORK(IU2CS+IMIN-1), R ) 00466 ELSE IF( NU .LT. MU ) THEN 00467 CALL DLARTGS( B21E( IMIN ), B21D( IMIN + 1 ), NU, 00468 $ RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1) ) 00469 ELSE 00470 CALL DLARTGS( B22D(IMIN), B22E(IMIN), MU, 00471 $ RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1) ) 00472 END IF 00473 RWORK(IU2CS+IMIN-1) = -RWORK(IU2CS+IMIN-1) 00474 RWORK(IU2SN+IMIN-1) = -RWORK(IU2SN+IMIN-1) 00475 * 00476 TEMP = RWORK(IU1CS+IMIN-1)*B11E(IMIN) + 00477 $ RWORK(IU1SN+IMIN-1)*B11D(IMIN+1) 00478 B11D(IMIN+1) = RWORK(IU1CS+IMIN-1)*B11D(IMIN+1) - 00479 $ RWORK(IU1SN+IMIN-1)*B11E(IMIN) 00480 B11E(IMIN) = TEMP 00481 IF( IMAX .GT. IMIN+1 ) THEN 00482 B11BULGE = RWORK(IU1SN+IMIN-1)*B11E(IMIN+1) 00483 B11E(IMIN+1) = RWORK(IU1CS+IMIN-1)*B11E(IMIN+1) 00484 END IF 00485 TEMP = RWORK(IU1CS+IMIN-1)*B12D(IMIN) + 00486 $ RWORK(IU1SN+IMIN-1)*B12E(IMIN) 00487 B12E(IMIN) = RWORK(IU1CS+IMIN-1)*B12E(IMIN) - 00488 $ RWORK(IU1SN+IMIN-1)*B12D(IMIN) 00489 B12D(IMIN) = TEMP 00490 B12BULGE = RWORK(IU1SN+IMIN-1)*B12D(IMIN+1) 00491 B12D(IMIN+1) = RWORK(IU1CS+IMIN-1)*B12D(IMIN+1) 00492 TEMP = RWORK(IU2CS+IMIN-1)*B21E(IMIN) + 00493 $ RWORK(IU2SN+IMIN-1)*B21D(IMIN+1) 00494 B21D(IMIN+1) = RWORK(IU2CS+IMIN-1)*B21D(IMIN+1) - 00495 $ RWORK(IU2SN+IMIN-1)*B21E(IMIN) 00496 B21E(IMIN) = TEMP 00497 IF( IMAX .GT. IMIN+1 ) THEN 00498 B21BULGE = RWORK(IU2SN+IMIN-1)*B21E(IMIN+1) 00499 B21E(IMIN+1) = RWORK(IU2CS+IMIN-1)*B21E(IMIN+1) 00500 END IF 00501 TEMP = RWORK(IU2CS+IMIN-1)*B22D(IMIN) + 00502 $ RWORK(IU2SN+IMIN-1)*B22E(IMIN) 00503 B22E(IMIN) = RWORK(IU2CS+IMIN-1)*B22E(IMIN) - 00504 $ RWORK(IU2SN+IMIN-1)*B22D(IMIN) 00505 B22D(IMIN) = TEMP 00506 B22BULGE = RWORK(IU2SN+IMIN-1)*B22D(IMIN+1) 00507 B22D(IMIN+1) = RWORK(IU2CS+IMIN-1)*B22D(IMIN+1) 00508 * 00509 * Inner loop: chase bulges from B11(IMIN,IMIN+2), 00510 * B12(IMIN,IMIN+1), B21(IMIN,IMIN+2), and B22(IMIN,IMIN+1) to 00511 * bottom-right 00512 * 00513 DO I = IMIN+1, IMAX-1 00514 * 00515 * Compute PHI(I-1) 00516 * 00517 X1 = SIN(THETA(I-1))*B11E(I-1) + COS(THETA(I-1))*B21E(I-1) 00518 X2 = SIN(THETA(I-1))*B11BULGE + COS(THETA(I-1))*B21BULGE 00519 Y1 = SIN(THETA(I-1))*B12D(I-1) + COS(THETA(I-1))*B22D(I-1) 00520 Y2 = SIN(THETA(I-1))*B12BULGE + COS(THETA(I-1))*B22BULGE 00521 * 00522 PHI(I-1) = ATAN2( SQRT(X1**2+X2**2), SQRT(Y1**2+Y2**2) ) 00523 * 00524 * Determine if there are bulges to chase or if a new direct 00525 * summand has been reached 00526 * 00527 RESTART11 = B11E(I-1)**2 + B11BULGE**2 .LE. THRESH**2 00528 RESTART21 = B21E(I-1)**2 + B21BULGE**2 .LE. THRESH**2 00529 RESTART12 = B12D(I-1)**2 + B12BULGE**2 .LE. THRESH**2 00530 RESTART22 = B22D(I-1)**2 + B22BULGE**2 .LE. THRESH**2 00531 * 00532 * If possible, chase bulges from B11(I-1,I+1), B12(I-1,I), 00533 * B21(I-1,I+1), and B22(I-1,I). If necessary, restart bulge- 00534 * chasing by applying the original shift again. 00535 * 00536 IF( .NOT. RESTART11 .AND. .NOT. RESTART21 ) THEN 00537 CALL DLARTGP( X2, X1, RWORK(IV1TSN+I-1), 00538 $ RWORK(IV1TCS+I-1), R ) 00539 ELSE IF( .NOT. RESTART11 .AND. RESTART21 ) THEN 00540 CALL DLARTGP( B11BULGE, B11E(I-1), RWORK(IV1TSN+I-1), 00541 $ RWORK(IV1TCS+I-1), R ) 00542 ELSE IF( RESTART11 .AND. .NOT. RESTART21 ) THEN 00543 CALL DLARTGP( B21BULGE, B21E(I-1), RWORK(IV1TSN+I-1), 00544 $ RWORK(IV1TCS+I-1), R ) 00545 ELSE IF( MU .LE. NU ) THEN 00546 CALL DLARTGS( B11D(I), B11E(I), MU, RWORK(IV1TCS+I-1), 00547 $ RWORK(IV1TSN+I-1) ) 00548 ELSE 00549 CALL DLARTGS( B21D(I), B21E(I), NU, RWORK(IV1TCS+I-1), 00550 $ RWORK(IV1TSN+I-1) ) 00551 END IF 00552 RWORK(IV1TCS+I-1) = -RWORK(IV1TCS+I-1) 00553 RWORK(IV1TSN+I-1) = -RWORK(IV1TSN+I-1) 00554 IF( .NOT. RESTART12 .AND. .NOT. RESTART22 ) THEN 00555 CALL DLARTGP( Y2, Y1, RWORK(IV2TSN+I-1-1), 00556 $ RWORK(IV2TCS+I-1-1), R ) 00557 ELSE IF( .NOT. RESTART12 .AND. RESTART22 ) THEN 00558 CALL DLARTGP( B12BULGE, B12D(I-1), RWORK(IV2TSN+I-1-1), 00559 $ RWORK(IV2TCS+I-1-1), R ) 00560 ELSE IF( RESTART12 .AND. .NOT. RESTART22 ) THEN 00561 CALL DLARTGP( B22BULGE, B22D(I-1), RWORK(IV2TSN+I-1-1), 00562 $ RWORK(IV2TCS+I-1-1), R ) 00563 ELSE IF( NU .LT. MU ) THEN 00564 CALL DLARTGS( B12E(I-1), B12D(I), NU, 00565 $ RWORK(IV2TCS+I-1-1), RWORK(IV2TSN+I-1-1) ) 00566 ELSE 00567 CALL DLARTGS( B22E(I-1), B22D(I), MU, 00568 $ RWORK(IV2TCS+I-1-1), RWORK(IV2TSN+I-1-1) ) 00569 END IF 00570 * 00571 TEMP = RWORK(IV1TCS+I-1)*B11D(I) + RWORK(IV1TSN+I-1)*B11E(I) 00572 B11E(I) = RWORK(IV1TCS+I-1)*B11E(I) - 00573 $ RWORK(IV1TSN+I-1)*B11D(I) 00574 B11D(I) = TEMP 00575 B11BULGE = RWORK(IV1TSN+I-1)*B11D(I+1) 00576 B11D(I+1) = RWORK(IV1TCS+I-1)*B11D(I+1) 00577 TEMP = RWORK(IV1TCS+I-1)*B21D(I) + RWORK(IV1TSN+I-1)*B21E(I) 00578 B21E(I) = RWORK(IV1TCS+I-1)*B21E(I) - 00579 $ RWORK(IV1TSN+I-1)*B21D(I) 00580 B21D(I) = TEMP 00581 B21BULGE = RWORK(IV1TSN+I-1)*B21D(I+1) 00582 B21D(I+1) = RWORK(IV1TCS+I-1)*B21D(I+1) 00583 TEMP = RWORK(IV2TCS+I-1-1)*B12E(I-1) + 00584 $ RWORK(IV2TSN+I-1-1)*B12D(I) 00585 B12D(I) = RWORK(IV2TCS+I-1-1)*B12D(I) - 00586 $ RWORK(IV2TSN+I-1-1)*B12E(I-1) 00587 B12E(I-1) = TEMP 00588 B12BULGE = RWORK(IV2TSN+I-1-1)*B12E(I) 00589 B12E(I) = RWORK(IV2TCS+I-1-1)*B12E(I) 00590 TEMP = RWORK(IV2TCS+I-1-1)*B22E(I-1) + 00591 $ RWORK(IV2TSN+I-1-1)*B22D(I) 00592 B22D(I) = RWORK(IV2TCS+I-1-1)*B22D(I) - 00593 $ RWORK(IV2TSN+I-1-1)*B22E(I-1) 00594 B22E(I-1) = TEMP 00595 B22BULGE = RWORK(IV2TSN+I-1-1)*B22E(I) 00596 B22E(I) = RWORK(IV2TCS+I-1-1)*B22E(I) 00597 * 00598 * Compute THETA(I) 00599 * 00600 X1 = COS(PHI(I-1))*B11D(I) + SIN(PHI(I-1))*B12E(I-1) 00601 X2 = COS(PHI(I-1))*B11BULGE + SIN(PHI(I-1))*B12BULGE 00602 Y1 = COS(PHI(I-1))*B21D(I) + SIN(PHI(I-1))*B22E(I-1) 00603 Y2 = COS(PHI(I-1))*B21BULGE + SIN(PHI(I-1))*B22BULGE 00604 * 00605 THETA(I) = ATAN2( SQRT(Y1**2+Y2**2), SQRT(X1**2+X2**2) ) 00606 * 00607 * Determine if there are bulges to chase or if a new direct 00608 * summand has been reached 00609 * 00610 RESTART11 = B11D(I)**2 + B11BULGE**2 .LE. THRESH**2 00611 RESTART12 = B12E(I-1)**2 + B12BULGE**2 .LE. THRESH**2 00612 RESTART21 = B21D(I)**2 + B21BULGE**2 .LE. THRESH**2 00613 RESTART22 = B22E(I-1)**2 + B22BULGE**2 .LE. THRESH**2 00614 * 00615 * If possible, chase bulges from B11(I+1,I), B12(I+1,I-1), 00616 * B21(I+1,I), and B22(I+1,I-1). If necessary, restart bulge- 00617 * chasing by applying the original shift again. 00618 * 00619 IF( .NOT. RESTART11 .AND. .NOT. RESTART12 ) THEN 00620 CALL DLARTGP( X2, X1, RWORK(IU1SN+I-1), RWORK(IU1CS+I-1), 00621 $ R ) 00622 ELSE IF( .NOT. RESTART11 .AND. RESTART12 ) THEN 00623 CALL DLARTGP( B11BULGE, B11D(I), RWORK(IU1SN+I-1), 00624 $ RWORK(IU1CS+I-1), R ) 00625 ELSE IF( RESTART11 .AND. .NOT. RESTART12 ) THEN 00626 CALL DLARTGP( B12BULGE, B12E(I-1), RWORK(IU1SN+I-1), 00627 $ RWORK(IU1CS+I-1), R ) 00628 ELSE IF( MU .LE. NU ) THEN 00629 CALL DLARTGS( B11E(I), B11D(I+1), MU, RWORK(IU1CS+I-1), 00630 $ RWORK(IU1SN+I-1) ) 00631 ELSE 00632 CALL DLARTGS( B12D(I), B12E(I), NU, RWORK(IU1CS+I-1), 00633 $ RWORK(IU1SN+I-1) ) 00634 END IF 00635 IF( .NOT. RESTART21 .AND. .NOT. RESTART22 ) THEN 00636 CALL DLARTGP( Y2, Y1, RWORK(IU2SN+I-1), RWORK(IU2CS+I-1), 00637 $ R ) 00638 ELSE IF( .NOT. RESTART21 .AND. RESTART22 ) THEN 00639 CALL DLARTGP( B21BULGE, B21D(I), RWORK(IU2SN+I-1), 00640 $ RWORK(IU2CS+I-1), R ) 00641 ELSE IF( RESTART21 .AND. .NOT. RESTART22 ) THEN 00642 CALL DLARTGP( B22BULGE, B22E(I-1), RWORK(IU2SN+I-1), 00643 $ RWORK(IU2CS+I-1), R ) 00644 ELSE IF( NU .LT. MU ) THEN 00645 CALL DLARTGS( B21E(I), B21E(I+1), NU, RWORK(IU2CS+I-1), 00646 $ RWORK(IU2SN+I-1) ) 00647 ELSE 00648 CALL DLARTGS( B22D(I), B22E(I), MU, RWORK(IU2CS+I-1), 00649 $ RWORK(IU2SN+I-1) ) 00650 END IF 00651 RWORK(IU2CS+I-1) = -RWORK(IU2CS+I-1) 00652 RWORK(IU2SN+I-1) = -RWORK(IU2SN+I-1) 00653 * 00654 TEMP = RWORK(IU1CS+I-1)*B11E(I) + RWORK(IU1SN+I-1)*B11D(I+1) 00655 B11D(I+1) = RWORK(IU1CS+I-1)*B11D(I+1) - 00656 $ RWORK(IU1SN+I-1)*B11E(I) 00657 B11E(I) = TEMP 00658 IF( I .LT. IMAX - 1 ) THEN 00659 B11BULGE = RWORK(IU1SN+I-1)*B11E(I+1) 00660 B11E(I+1) = RWORK(IU1CS+I-1)*B11E(I+1) 00661 END IF 00662 TEMP = RWORK(IU2CS+I-1)*B21E(I) + RWORK(IU2SN+I-1)*B21D(I+1) 00663 B21D(I+1) = RWORK(IU2CS+I-1)*B21D(I+1) - 00664 $ RWORK(IU2SN+I-1)*B21E(I) 00665 B21E(I) = TEMP 00666 IF( I .LT. IMAX - 1 ) THEN 00667 B21BULGE = RWORK(IU2SN+I-1)*B21E(I+1) 00668 B21E(I+1) = RWORK(IU2CS+I-1)*B21E(I+1) 00669 END IF 00670 TEMP = RWORK(IU1CS+I-1)*B12D(I) + RWORK(IU1SN+I-1)*B12E(I) 00671 B12E(I) = RWORK(IU1CS+I-1)*B12E(I) - 00672 $ RWORK(IU1SN+I-1)*B12D(I) 00673 B12D(I) = TEMP 00674 B12BULGE = RWORK(IU1SN+I-1)*B12D(I+1) 00675 B12D(I+1) = RWORK(IU1CS+I-1)*B12D(I+1) 00676 TEMP = RWORK(IU2CS+I-1)*B22D(I) + RWORK(IU2SN+I-1)*B22E(I) 00677 B22E(I) = RWORK(IU2CS+I-1)*B22E(I) - 00678 $ RWORK(IU2SN+I-1)*B22D(I) 00679 B22D(I) = TEMP 00680 B22BULGE = RWORK(IU2SN+I-1)*B22D(I+1) 00681 B22D(I+1) = RWORK(IU2CS+I-1)*B22D(I+1) 00682 * 00683 END DO 00684 * 00685 * Compute PHI(IMAX-1) 00686 * 00687 X1 = SIN(THETA(IMAX-1))*B11E(IMAX-1) + 00688 $ COS(THETA(IMAX-1))*B21E(IMAX-1) 00689 Y1 = SIN(THETA(IMAX-1))*B12D(IMAX-1) + 00690 $ COS(THETA(IMAX-1))*B22D(IMAX-1) 00691 Y2 = SIN(THETA(IMAX-1))*B12BULGE + COS(THETA(IMAX-1))*B22BULGE 00692 * 00693 PHI(IMAX-1) = ATAN2( ABS(X1), SQRT(Y1**2+Y2**2) ) 00694 * 00695 * Chase bulges from B12(IMAX-1,IMAX) and B22(IMAX-1,IMAX) 00696 * 00697 RESTART12 = B12D(IMAX-1)**2 + B12BULGE**2 .LE. THRESH**2 00698 RESTART22 = B22D(IMAX-1)**2 + B22BULGE**2 .LE. THRESH**2 00699 * 00700 IF( .NOT. RESTART12 .AND. .NOT. RESTART22 ) THEN 00701 CALL DLARTGP( Y2, Y1, RWORK(IV2TSN+IMAX-1-1), 00702 $ RWORK(IV2TCS+IMAX-1-1), R ) 00703 ELSE IF( .NOT. RESTART12 .AND. RESTART22 ) THEN 00704 CALL DLARTGP( B12BULGE, B12D(IMAX-1), 00705 $ RWORK(IV2TSN+IMAX-1-1), 00706 $ RWORK(IV2TCS+IMAX-1-1), R ) 00707 ELSE IF( RESTART12 .AND. .NOT. RESTART22 ) THEN 00708 CALL DLARTGP( B22BULGE, B22D(IMAX-1), 00709 $ RWORK(IV2TSN+IMAX-1-1), 00710 $ RWORK(IV2TCS+IMAX-1-1), R ) 00711 ELSE IF( NU .LT. MU ) THEN 00712 CALL DLARTGS( B12E(IMAX-1), B12D(IMAX), NU, 00713 $ RWORK(IV2TCS+IMAX-1-1), 00714 $ RWORK(IV2TSN+IMAX-1-1) ) 00715 ELSE 00716 CALL DLARTGS( B22E(IMAX-1), B22D(IMAX), MU, 00717 $ RWORK(IV2TCS+IMAX-1-1), 00718 $ RWORK(IV2TSN+IMAX-1-1) ) 00719 END IF 00720 * 00721 TEMP = RWORK(IV2TCS+IMAX-1-1)*B12E(IMAX-1) + 00722 $ RWORK(IV2TSN+IMAX-1-1)*B12D(IMAX) 00723 B12D(IMAX) = RWORK(IV2TCS+IMAX-1-1)*B12D(IMAX) - 00724 $ RWORK(IV2TSN+IMAX-1-1)*B12E(IMAX-1) 00725 B12E(IMAX-1) = TEMP 00726 TEMP = RWORK(IV2TCS+IMAX-1-1)*B22E(IMAX-1) + 00727 $ RWORK(IV2TSN+IMAX-1-1)*B22D(IMAX) 00728 B22D(IMAX) = RWORK(IV2TCS+IMAX-1-1)*B22D(IMAX) - 00729 $ RWORK(IV2TSN+IMAX-1-1)*B22E(IMAX-1) 00730 B22E(IMAX-1) = TEMP 00731 * 00732 * Update singular vectors 00733 * 00734 IF( WANTU1 ) THEN 00735 IF( COLMAJOR ) THEN 00736 CALL ZLASR( 'R', 'V', 'F', P, IMAX-IMIN+1, 00737 $ RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1), 00738 $ U1(1,IMIN), LDU1 ) 00739 ELSE 00740 CALL ZLASR( 'L', 'V', 'F', IMAX-IMIN+1, P, 00741 $ RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1), 00742 $ U1(IMIN,1), LDU1 ) 00743 END IF 00744 END IF 00745 IF( WANTU2 ) THEN 00746 IF( COLMAJOR ) THEN 00747 CALL ZLASR( 'R', 'V', 'F', M-P, IMAX-IMIN+1, 00748 $ RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1), 00749 $ U2(1,IMIN), LDU2 ) 00750 ELSE 00751 CALL ZLASR( 'L', 'V', 'F', IMAX-IMIN+1, M-P, 00752 $ RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1), 00753 $ U2(IMIN,1), LDU2 ) 00754 END IF 00755 END IF 00756 IF( WANTV1T ) THEN 00757 IF( COLMAJOR ) THEN 00758 CALL ZLASR( 'L', 'V', 'F', IMAX-IMIN+1, Q, 00759 $ RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1), 00760 $ V1T(IMIN,1), LDV1T ) 00761 ELSE 00762 CALL ZLASR( 'R', 'V', 'F', Q, IMAX-IMIN+1, 00763 $ RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1), 00764 $ V1T(1,IMIN), LDV1T ) 00765 END IF 00766 END IF 00767 IF( WANTV2T ) THEN 00768 IF( COLMAJOR ) THEN 00769 CALL ZLASR( 'L', 'V', 'F', IMAX-IMIN+1, M-Q, 00770 $ RWORK(IV2TCS+IMIN-1), RWORK(IV2TSN+IMIN-1), 00771 $ V2T(IMIN,1), LDV2T ) 00772 ELSE 00773 CALL ZLASR( 'R', 'V', 'F', M-Q, IMAX-IMIN+1, 00774 $ RWORK(IV2TCS+IMIN-1), RWORK(IV2TSN+IMIN-1), 00775 $ V2T(1,IMIN), LDV2T ) 00776 END IF 00777 END IF 00778 * 00779 * Fix signs on B11(IMAX-1,IMAX) and B21(IMAX-1,IMAX) 00780 * 00781 IF( B11E(IMAX-1)+B21E(IMAX-1) .GT. 0 ) THEN 00782 B11D(IMAX) = -B11D(IMAX) 00783 B21D(IMAX) = -B21D(IMAX) 00784 IF( WANTV1T ) THEN 00785 IF( COLMAJOR ) THEN 00786 CALL ZSCAL( Q, NEGONECOMPLEX, V1T(IMAX,1), LDV1T ) 00787 ELSE 00788 CALL ZSCAL( Q, NEGONECOMPLEX, V1T(1,IMAX), 1 ) 00789 END IF 00790 END IF 00791 END IF 00792 * 00793 * Compute THETA(IMAX) 00794 * 00795 X1 = COS(PHI(IMAX-1))*B11D(IMAX) + 00796 $ SIN(PHI(IMAX-1))*B12E(IMAX-1) 00797 Y1 = COS(PHI(IMAX-1))*B21D(IMAX) + 00798 $ SIN(PHI(IMAX-1))*B22E(IMAX-1) 00799 * 00800 THETA(IMAX) = ATAN2( ABS(Y1), ABS(X1) ) 00801 * 00802 * Fix signs on B11(IMAX,IMAX), B12(IMAX,IMAX-1), B21(IMAX,IMAX), 00803 * and B22(IMAX,IMAX-1) 00804 * 00805 IF( B11D(IMAX)+B12E(IMAX-1) .LT. 0 ) THEN 00806 B12D(IMAX) = -B12D(IMAX) 00807 IF( WANTU1 ) THEN 00808 IF( COLMAJOR ) THEN 00809 CALL ZSCAL( P, NEGONECOMPLEX, U1(1,IMAX), 1 ) 00810 ELSE 00811 CALL ZSCAL( P, NEGONECOMPLEX, U1(IMAX,1), LDU1 ) 00812 END IF 00813 END IF 00814 END IF 00815 IF( B21D(IMAX)+B22E(IMAX-1) .GT. 0 ) THEN 00816 B22D(IMAX) = -B22D(IMAX) 00817 IF( WANTU2 ) THEN 00818 IF( COLMAJOR ) THEN 00819 CALL ZSCAL( M-P, NEGONECOMPLEX, U2(1,IMAX), 1 ) 00820 ELSE 00821 CALL ZSCAL( M-P, NEGONECOMPLEX, U2(IMAX,1), LDU2 ) 00822 END IF 00823 END IF 00824 END IF 00825 * 00826 * Fix signs on B12(IMAX,IMAX) and B22(IMAX,IMAX) 00827 * 00828 IF( B12D(IMAX)+B22D(IMAX) .LT. 0 ) THEN 00829 IF( WANTV2T ) THEN 00830 IF( COLMAJOR ) THEN 00831 CALL ZSCAL( M-Q, NEGONECOMPLEX, V2T(IMAX,1), LDV2T ) 00832 ELSE 00833 CALL ZSCAL( M-Q, NEGONECOMPLEX, V2T(1,IMAX), 1 ) 00834 END IF 00835 END IF 00836 END IF 00837 * 00838 * Test for negligible sines or cosines 00839 * 00840 DO I = IMIN, IMAX 00841 IF( THETA(I) .LT. THRESH ) THEN 00842 THETA(I) = ZERO 00843 ELSE IF( THETA(I) .GT. PIOVER2-THRESH ) THEN 00844 THETA(I) = PIOVER2 00845 END IF 00846 END DO 00847 DO I = IMIN, IMAX-1 00848 IF( PHI(I) .LT. THRESH ) THEN 00849 PHI(I) = ZERO 00850 ELSE IF( PHI(I) .GT. PIOVER2-THRESH ) THEN 00851 PHI(I) = PIOVER2 00852 END IF 00853 END DO 00854 * 00855 * Deflate 00856 * 00857 IF (IMAX .GT. 1) THEN 00858 DO WHILE( PHI(IMAX-1) .EQ. ZERO ) 00859 IMAX = IMAX - 1 00860 IF (IMAX .LE. 1) EXIT 00861 END DO 00862 END IF 00863 IF( IMIN .GT. IMAX - 1 ) 00864 $ IMIN = IMAX - 1 00865 IF (IMIN .GT. 1) THEN 00866 DO WHILE (PHI(IMIN-1) .NE. ZERO) 00867 IMIN = IMIN - 1 00868 IF (IMIN .LE. 1) EXIT 00869 END DO 00870 END IF 00871 * 00872 * Repeat main iteration loop 00873 * 00874 END DO 00875 * 00876 * Postprocessing: order THETA from least to greatest 00877 * 00878 DO I = 1, Q 00879 * 00880 MINI = I 00881 THETAMIN = THETA(I) 00882 DO J = I+1, Q 00883 IF( THETA(J) .LT. THETAMIN ) THEN 00884 MINI = J 00885 THETAMIN = THETA(J) 00886 END IF 00887 END DO 00888 * 00889 IF( MINI .NE. I ) THEN 00890 THETA(MINI) = THETA(I) 00891 THETA(I) = THETAMIN 00892 IF( COLMAJOR ) THEN 00893 IF( WANTU1 ) 00894 $ CALL ZSWAP( P, U1(1,I), 1, U1(1,MINI), 1 ) 00895 IF( WANTU2 ) 00896 $ CALL ZSWAP( M-P, U2(1,I), 1, U2(1,MINI), 1 ) 00897 IF( WANTV1T ) 00898 $ CALL ZSWAP( Q, V1T(I,1), LDV1T, V1T(MINI,1), LDV1T ) 00899 IF( WANTV2T ) 00900 $ CALL ZSWAP( M-Q, V2T(I,1), LDV2T, V2T(MINI,1), 00901 $ LDV2T ) 00902 ELSE 00903 IF( WANTU1 ) 00904 $ CALL ZSWAP( P, U1(I,1), LDU1, U1(MINI,1), LDU1 ) 00905 IF( WANTU2 ) 00906 $ CALL ZSWAP( M-P, U2(I,1), LDU2, U2(MINI,1), LDU2 ) 00907 IF( WANTV1T ) 00908 $ CALL ZSWAP( Q, V1T(1,I), 1, V1T(1,MINI), 1 ) 00909 IF( WANTV2T ) 00910 $ CALL ZSWAP( M-Q, V2T(1,I), 1, V2T(1,MINI), 1 ) 00911 END IF 00912 END IF 00913 * 00914 END DO 00915 * 00916 RETURN 00917 * 00918 * End of ZBBCSD 00919 * 00920 END 00921
1.7.3 
