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