 |
LAPACK
3.4.0
LAPACK: Linear Algebra PACKage
|
|
LAPACK
3.4.0
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DORCSD 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DORCSD + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dorcsd.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dorcsd.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dorcsd.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * RECURSIVE SUBROUTINE DORCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, 00022 * SIGNS, M, P, Q, X11, LDX11, X12, 00023 * LDX12, X21, LDX21, X22, LDX22, THETA, 00024 * U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, 00025 * LDV2T, WORK, LWORK, IWORK, INFO ) 00026 * 00027 * .. Scalar Arguments .. 00028 * CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, SIGNS, TRANS 00029 * INTEGER INFO, LDU1, LDU2, LDV1T, LDV2T, LDX11, LDX12, 00030 * $ LDX21, LDX22, LWORK, M, P, Q 00031 * .. 00032 * .. Array Arguments .. 00033 * INTEGER IWORK( * ) 00034 * DOUBLE PRECISION THETA( * ) 00035 * DOUBLE PRECISION U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ), 00036 * $ V2T( LDV2T, * ), WORK( * ), X11( LDX11, * ), 00037 * $ X12( LDX12, * ), X21( LDX21, * ), X22( LDX22, 00038 * $ * ) 00039 * .. 00040 * 00041 * 00042 *> \par Purpose: 00043 * ============= 00044 *> 00045 *> \verbatim 00046 *> 00047 *> DORCSD computes the CS decomposition of an M-by-M partitioned 00048 *> orthogonal matrix X: 00049 *> 00050 *> [ I 0 0 | 0 0 0 ] 00051 *> [ 0 C 0 | 0 -S 0 ] 00052 *> [ X11 | X12 ] [ U1 | ] [ 0 0 0 | 0 0 -I ] [ V1 | ]**T 00053 *> X = [-----------] = [---------] [---------------------] [---------] . 00054 *> [ X21 | X22 ] [ | U2 ] [ 0 0 0 | I 0 0 ] [ | V2 ] 00055 *> [ 0 S 0 | 0 C 0 ] 00056 *> [ 0 0 I | 0 0 0 ] 00057 *> 00058 *> X11 is P-by-Q. The orthogonal matrices U1, U2, V1, and V2 are P-by-P, 00059 *> (M-P)-by-(M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. C and S are 00060 *> R-by-R nonnegative diagonal matrices satisfying C^2 + S^2 = I, in 00061 *> which R = MIN(P,M-P,Q,M-Q). 00062 *> \endverbatim 00063 * 00064 * Arguments: 00065 * ========== 00066 * 00067 *> \param[in] JOBU1 00068 *> \verbatim 00069 *> JOBU1 is CHARACTER 00070 *> = 'Y': U1 is computed; 00071 *> otherwise: U1 is not computed. 00072 *> \endverbatim 00073 *> 00074 *> \param[in] JOBU2 00075 *> \verbatim 00076 *> JOBU2 is CHARACTER 00077 *> = 'Y': U2 is computed; 00078 *> otherwise: U2 is not computed. 00079 *> \endverbatim 00080 *> 00081 *> \param[in] JOBV1T 00082 *> \verbatim 00083 *> JOBV1T is CHARACTER 00084 *> = 'Y': V1T is computed; 00085 *> otherwise: V1T is not computed. 00086 *> \endverbatim 00087 *> 00088 *> \param[in] JOBV2T 00089 *> \verbatim 00090 *> JOBV2T is CHARACTER 00091 *> = 'Y': V2T is computed; 00092 *> otherwise: V2T is not computed. 00093 *> \endverbatim 00094 *> 00095 *> \param[in] TRANS 00096 *> \verbatim 00097 *> TRANS is CHARACTER 00098 *> = 'T': X, U1, U2, V1T, and V2T are stored in row-major 00099 *> order; 00100 *> otherwise: X, U1, U2, V1T, and V2T are stored in column- 00101 *> major order. 00102 *> \endverbatim 00103 *> 00104 *> \param[in] SIGNS 00105 *> \verbatim 00106 *> SIGNS is CHARACTER 00107 *> = 'O': The lower-left block is made nonpositive (the 00108 *> "other" convention); 00109 *> otherwise: The upper-right block is made nonpositive (the 00110 *> "default" convention). 00111 *> \endverbatim 00112 *> 00113 *> \param[in] M 00114 *> \verbatim 00115 *> M is INTEGER 00116 *> The number of rows and columns in X. 00117 *> \endverbatim 00118 *> 00119 *> \param[in] P 00120 *> \verbatim 00121 *> P is INTEGER 00122 *> The number of rows in X11 and X12. 0 <= P <= M. 00123 *> \endverbatim 00124 *> 00125 *> \param[in] Q 00126 *> \verbatim 00127 *> Q is INTEGER 00128 *> The number of columns in X11 and X21. 0 <= Q <= M. 00129 *> \endverbatim 00130 *> 00131 *> \param[in,out] X11 00132 *> \verbatim 00133 *> X11 is DOUBLE PRECISION array, dimension (LDX11,Q) 00134 *> On entry, part of the orthogonal matrix whose CSD is desired. 00135 *> \endverbatim 00136 *> 00137 *> \param[in] LDX11 00138 *> \verbatim 00139 *> LDX11 is INTEGER 00140 *> The leading dimension of X11. LDX11 >= MAX(1,P). 00141 *> \endverbatim 00142 *> 00143 *> \param[in,out] X12 00144 *> \verbatim 00145 *> X12 is DOUBLE PRECISION array, dimension (LDX12,M-Q) 00146 *> On entry, part of the orthogonal matrix whose CSD is desired. 00147 *> \endverbatim 00148 *> 00149 *> \param[in] LDX12 00150 *> \verbatim 00151 *> LDX12 is INTEGER 00152 *> The leading dimension of X12. LDX12 >= MAX(1,P). 00153 *> \endverbatim 00154 *> 00155 *> \param[in,out] X21 00156 *> \verbatim 00157 *> X21 is DOUBLE PRECISION array, dimension (LDX21,Q) 00158 *> On entry, part of the orthogonal matrix whose CSD is desired. 00159 *> \endverbatim 00160 *> 00161 *> \param[in] LDX21 00162 *> \verbatim 00163 *> LDX21 is INTEGER 00164 *> The leading dimension of X11. LDX21 >= MAX(1,M-P). 00165 *> \endverbatim 00166 *> 00167 *> \param[in,out] X22 00168 *> \verbatim 00169 *> X22 is DOUBLE PRECISION array, dimension (LDX22,M-Q) 00170 *> On entry, part of the orthogonal matrix whose CSD is desired. 00171 *> \endverbatim 00172 *> 00173 *> \param[in] LDX22 00174 *> \verbatim 00175 *> LDX22 is INTEGER 00176 *> The leading dimension of X11. LDX22 >= MAX(1,M-P). 00177 *> \endverbatim 00178 *> 00179 *> \param[out] THETA 00180 *> \verbatim 00181 *> THETA is DOUBLE PRECISION array, dimension (R), in which R = 00182 *> MIN(P,M-P,Q,M-Q). 00183 *> C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and 00184 *> S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ). 00185 *> \endverbatim 00186 *> 00187 *> \param[out] U1 00188 *> \verbatim 00189 *> U1 is DOUBLE PRECISION array, dimension (P) 00190 *> If JOBU1 = 'Y', U1 contains the P-by-P orthogonal matrix U1. 00191 *> \endverbatim 00192 *> 00193 *> \param[in] LDU1 00194 *> \verbatim 00195 *> LDU1 is INTEGER 00196 *> The leading dimension of U1. If JOBU1 = 'Y', LDU1 >= 00197 *> MAX(1,P). 00198 *> \endverbatim 00199 *> 00200 *> \param[out] U2 00201 *> \verbatim 00202 *> U2 is DOUBLE PRECISION array, dimension (M-P) 00203 *> If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) orthogonal 00204 *> matrix U2. 00205 *> \endverbatim 00206 *> 00207 *> \param[in] LDU2 00208 *> \verbatim 00209 *> LDU2 is INTEGER 00210 *> The leading dimension of U2. If JOBU2 = 'Y', LDU2 >= 00211 *> MAX(1,M-P). 00212 *> \endverbatim 00213 *> 00214 *> \param[out] V1T 00215 *> \verbatim 00216 *> V1T is DOUBLE PRECISION array, dimension (Q) 00217 *> If JOBV1T = 'Y', V1T contains the Q-by-Q matrix orthogonal 00218 *> matrix V1**T. 00219 *> \endverbatim 00220 *> 00221 *> \param[in] LDV1T 00222 *> \verbatim 00223 *> LDV1T is INTEGER 00224 *> The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >= 00225 *> MAX(1,Q). 00226 *> \endverbatim 00227 *> 00228 *> \param[out] V2T 00229 *> \verbatim 00230 *> V2T is DOUBLE PRECISION array, dimension (M-Q) 00231 *> If JOBV2T = 'Y', V2T contains the (M-Q)-by-(M-Q) orthogonal 00232 *> matrix V2**T. 00233 *> \endverbatim 00234 *> 00235 *> \param[in] LDV2T 00236 *> \verbatim 00237 *> LDV2T is INTEGER 00238 *> The leading dimension of V2T. If JOBV2T = 'Y', LDV2T >= 00239 *> MAX(1,M-Q). 00240 *> \endverbatim 00241 *> 00242 *> \param[out] WORK 00243 *> \verbatim 00244 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 00245 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00246 *> If INFO > 0 on exit, WORK(2:R) contains the values PHI(1), 00247 *> ..., PHI(R-1) that, together with THETA(1), ..., THETA(R), 00248 *> define the matrix in intermediate bidiagonal-block form 00249 *> remaining after nonconvergence. INFO specifies the number 00250 *> of nonzero PHI's. 00251 *> \endverbatim 00252 *> 00253 *> \param[in] LWORK 00254 *> \verbatim 00255 *> LWORK is INTEGER 00256 *> The dimension of the array WORK. 00257 *> 00258 *> If LWORK = -1, then a workspace query is assumed; the routine 00259 *> only calculates the optimal size of the WORK array, returns 00260 *> this value as the first entry of the work array, and no error 00261 *> message related to LWORK is issued by XERBLA. 00262 *> \endverbatim 00263 *> 00264 *> \param[out] IWORK 00265 *> \verbatim 00266 *> IWORK is INTEGER array, dimension (M-MIN(P, M-P, Q, M-Q)) 00267 *> \endverbatim 00268 *> 00269 *> \param[out] INFO 00270 *> \verbatim 00271 *> INFO is INTEGER 00272 *> = 0: successful exit. 00273 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00274 *> > 0: DBBCSD did not converge. See the description of WORK 00275 *> above for details. 00276 *> \endverbatim 00277 * 00278 *> \par References: 00279 * ================ 00280 *> 00281 *> [1] Brian D. Sutton. Computing the complete CS decomposition. Numer. 00282 *> Algorithms, 50(1):33-65, 2009. 00283 * 00284 * Authors: 00285 * ======== 00286 * 00287 *> \author Univ. of Tennessee 00288 *> \author Univ. of California Berkeley 00289 *> \author Univ. of Colorado Denver 00290 *> \author NAG Ltd. 00291 * 00292 *> \date November 2011 00293 * 00294 *> \ingroup doubleOTHERcomputational 00295 * 00296 * ===================================================================== 00297 RECURSIVE SUBROUTINE DORCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, 00298 $ SIGNS, M, P, Q, X11, LDX11, X12, 00299 $ LDX12, X21, LDX21, X22, LDX22, THETA, 00300 $ U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, 00301 $ LDV2T, WORK, LWORK, IWORK, INFO ) 00302 * 00303 * -- LAPACK computational routine (version 3.4.0) -- 00304 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00305 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00306 * November 2011 00307 * 00308 * .. Scalar Arguments .. 00309 CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, SIGNS, TRANS 00310 INTEGER INFO, LDU1, LDU2, LDV1T, LDV2T, LDX11, LDX12, 00311 $ LDX21, LDX22, LWORK, M, P, Q 00312 * .. 00313 * .. Array Arguments .. 00314 INTEGER IWORK( * ) 00315 DOUBLE PRECISION THETA( * ) 00316 DOUBLE PRECISION U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ), 00317 $ V2T( LDV2T, * ), WORK( * ), X11( LDX11, * ), 00318 $ X12( LDX12, * ), X21( LDX21, * ), X22( LDX22, 00319 $ * ) 00320 * .. 00321 * 00322 * =================================================================== 00323 * 00324 * .. Parameters .. 00325 DOUBLE PRECISION REALONE 00326 PARAMETER ( REALONE = 1.0D0 ) 00327 DOUBLE PRECISION NEGONE, ONE, PIOVER2, ZERO 00328 PARAMETER ( NEGONE = -1.0D0, ONE = 1.0D0, 00329 $ PIOVER2 = 1.57079632679489662D0, 00330 $ ZERO = 0.0D0 ) 00331 * .. 00332 * .. Local Scalars .. 00333 CHARACTER TRANST, SIGNST 00334 INTEGER CHILDINFO, I, IB11D, IB11E, IB12D, IB12E, 00335 $ IB21D, IB21E, IB22D, IB22E, IBBCSD, IORBDB, 00336 $ IORGLQ, IORGQR, IPHI, ITAUP1, ITAUP2, ITAUQ1, 00337 $ ITAUQ2, J, LBBCSDWORK, LBBCSDWORKMIN, 00338 $ LBBCSDWORKOPT, LORBDBWORK, LORBDBWORKMIN, 00339 $ LORBDBWORKOPT, LORGLQWORK, LORGLQWORKMIN, 00340 $ LORGLQWORKOPT, LORGQRWORK, LORGQRWORKMIN, 00341 $ LORGQRWORKOPT, LWORKMIN, LWORKOPT 00342 LOGICAL COLMAJOR, DEFAULTSIGNS, LQUERY, WANTU1, WANTU2, 00343 $ WANTV1T, WANTV2T 00344 * .. 00345 * .. External Subroutines .. 00346 EXTERNAL DBBCSD, DLACPY, DLAPMR, DLAPMT, DLASCL, DLASET, 00347 $ DORBDB, DORGLQ, DORGQR, XERBLA 00348 * .. 00349 * .. External Functions .. 00350 LOGICAL LSAME 00351 EXTERNAL LSAME 00352 * .. 00353 * .. Intrinsic Functions 00354 INTRINSIC INT, MAX, MIN 00355 * .. 00356 * .. Executable Statements .. 00357 * 00358 * Test input arguments 00359 * 00360 INFO = 0 00361 WANTU1 = LSAME( JOBU1, 'Y' ) 00362 WANTU2 = LSAME( JOBU2, 'Y' ) 00363 WANTV1T = LSAME( JOBV1T, 'Y' ) 00364 WANTV2T = LSAME( JOBV2T, 'Y' ) 00365 COLMAJOR = .NOT. LSAME( TRANS, 'T' ) 00366 DEFAULTSIGNS = .NOT. LSAME( SIGNS, 'O' ) 00367 LQUERY = LWORK .EQ. -1 00368 IF( M .LT. 0 ) THEN 00369 INFO = -7 00370 ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN 00371 INFO = -8 00372 ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN 00373 INFO = -9 00374 ELSE IF( ( COLMAJOR .AND. LDX11 .LT. MAX(1,P) ) .OR. 00375 $ ( .NOT.COLMAJOR .AND. LDX11 .LT. MAX(1,Q) ) ) THEN 00376 INFO = -11 00377 ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN 00378 INFO = -20 00379 ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN 00380 INFO = -22 00381 ELSE IF( WANTV1T .AND. LDV1T .LT. Q ) THEN 00382 INFO = -24 00383 ELSE IF( WANTV2T .AND. LDV2T .LT. M-Q ) THEN 00384 INFO = -26 00385 END IF 00386 * 00387 * Work with transpose if convenient 00388 * 00389 IF( INFO .EQ. 0 .AND. MIN( P, M-P ) .LT. MIN( Q, M-Q ) ) THEN 00390 IF( COLMAJOR ) THEN 00391 TRANST = 'T' 00392 ELSE 00393 TRANST = 'N' 00394 END IF 00395 IF( DEFAULTSIGNS ) THEN 00396 SIGNST = 'O' 00397 ELSE 00398 SIGNST = 'D' 00399 END IF 00400 CALL DORCSD( JOBV1T, JOBV2T, JOBU1, JOBU2, TRANST, SIGNST, M, 00401 $ Q, P, X11, LDX11, X21, LDX21, X12, LDX12, X22, 00402 $ LDX22, THETA, V1T, LDV1T, V2T, LDV2T, U1, LDU1, 00403 $ U2, LDU2, WORK, LWORK, IWORK, INFO ) 00404 RETURN 00405 END IF 00406 * 00407 * Work with permutation [ 0 I; I 0 ] * X * [ 0 I; I 0 ] if 00408 * convenient 00409 * 00410 IF( INFO .EQ. 0 .AND. M-Q .LT. Q ) THEN 00411 IF( DEFAULTSIGNS ) THEN 00412 SIGNST = 'O' 00413 ELSE 00414 SIGNST = 'D' 00415 END IF 00416 CALL DORCSD( JOBU2, JOBU1, JOBV2T, JOBV1T, TRANS, SIGNST, M, 00417 $ M-P, M-Q, X22, LDX22, X21, LDX21, X12, LDX12, X11, 00418 $ LDX11, THETA, U2, LDU2, U1, LDU1, V2T, LDV2T, V1T, 00419 $ LDV1T, WORK, LWORK, IWORK, INFO ) 00420 RETURN 00421 END IF 00422 * 00423 * Compute workspace 00424 * 00425 IF( INFO .EQ. 0 ) THEN 00426 * 00427 IPHI = 2 00428 ITAUP1 = IPHI + MAX( 1, Q - 1 ) 00429 ITAUP2 = ITAUP1 + MAX( 1, P ) 00430 ITAUQ1 = ITAUP2 + MAX( 1, M - P ) 00431 ITAUQ2 = ITAUQ1 + MAX( 1, Q ) 00432 IORGQR = ITAUQ2 + MAX( 1, M - Q ) 00433 CALL DORGQR( M-Q, M-Q, M-Q, 0, MAX(1,M-Q), 0, WORK, -1, 00434 $ CHILDINFO ) 00435 LORGQRWORKOPT = INT( WORK(1) ) 00436 LORGQRWORKMIN = MAX( 1, M - Q ) 00437 IORGLQ = ITAUQ2 + MAX( 1, M - Q ) 00438 CALL DORGLQ( M-Q, M-Q, M-Q, 0, MAX(1,M-Q), 0, WORK, -1, 00439 $ CHILDINFO ) 00440 LORGLQWORKOPT = INT( WORK(1) ) 00441 LORGLQWORKMIN = MAX( 1, M - Q ) 00442 IORBDB = ITAUQ2 + MAX( 1, M - Q ) 00443 CALL DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, 00444 $ X21, LDX21, X22, LDX22, 0, 0, 0, 0, 0, 0, WORK, 00445 $ -1, CHILDINFO ) 00446 LORBDBWORKOPT = INT( WORK(1) ) 00447 LORBDBWORKMIN = LORBDBWORKOPT 00448 IB11D = ITAUQ2 + MAX( 1, M - Q ) 00449 IB11E = IB11D + MAX( 1, Q ) 00450 IB12D = IB11E + MAX( 1, Q - 1 ) 00451 IB12E = IB12D + MAX( 1, Q ) 00452 IB21D = IB12E + MAX( 1, Q - 1 ) 00453 IB21E = IB21D + MAX( 1, Q ) 00454 IB22D = IB21E + MAX( 1, Q - 1 ) 00455 IB22E = IB22D + MAX( 1, Q ) 00456 IBBCSD = IB22E + MAX( 1, Q - 1 ) 00457 CALL DBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, 0, 00458 $ 0, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, 0, 00459 $ 0, 0, 0, 0, 0, 0, 0, WORK, -1, CHILDINFO ) 00460 LBBCSDWORKOPT = INT( WORK(1) ) 00461 LBBCSDWORKMIN = LBBCSDWORKOPT 00462 LWORKOPT = MAX( IORGQR + LORGQRWORKOPT, IORGLQ + LORGLQWORKOPT, 00463 $ IORBDB + LORBDBWORKOPT, IBBCSD + LBBCSDWORKOPT ) - 1 00464 LWORKMIN = MAX( IORGQR + LORGQRWORKMIN, IORGLQ + LORGLQWORKMIN, 00465 $ IORBDB + LORBDBWORKOPT, IBBCSD + LBBCSDWORKMIN ) - 1 00466 WORK(1) = MAX(LWORKOPT,LWORKMIN) 00467 * 00468 IF( LWORK .LT. LWORKMIN .AND. .NOT. LQUERY ) THEN 00469 INFO = -22 00470 ELSE 00471 LORGQRWORK = LWORK - IORGQR + 1 00472 LORGLQWORK = LWORK - IORGLQ + 1 00473 LORBDBWORK = LWORK - IORBDB + 1 00474 LBBCSDWORK = LWORK - IBBCSD + 1 00475 END IF 00476 END IF 00477 * 00478 * Abort if any illegal arguments 00479 * 00480 IF( INFO .NE. 0 ) THEN 00481 CALL XERBLA( 'DORCSD', -INFO ) 00482 RETURN 00483 ELSE IF( LQUERY ) THEN 00484 RETURN 00485 END IF 00486 * 00487 * Transform to bidiagonal block form 00488 * 00489 CALL DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21, 00490 $ LDX21, X22, LDX22, THETA, WORK(IPHI), WORK(ITAUP1), 00491 $ WORK(ITAUP2), WORK(ITAUQ1), WORK(ITAUQ2), 00492 $ WORK(IORBDB), LORBDBWORK, CHILDINFO ) 00493 * 00494 * Accumulate Householder reflectors 00495 * 00496 IF( COLMAJOR ) THEN 00497 IF( WANTU1 .AND. P .GT. 0 ) THEN 00498 CALL DLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 ) 00499 CALL DORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR), 00500 $ LORGQRWORK, INFO) 00501 END IF 00502 IF( WANTU2 .AND. M-P .GT. 0 ) THEN 00503 CALL DLACPY( 'L', M-P, Q, X21, LDX21, U2, LDU2 ) 00504 CALL DORGQR( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2), 00505 $ WORK(IORGQR), LORGQRWORK, INFO ) 00506 END IF 00507 IF( WANTV1T .AND. Q .GT. 0 ) THEN 00508 CALL DLACPY( 'U', Q-1, Q-1, X11(1,2), LDX11, V1T(2,2), 00509 $ LDV1T ) 00510 V1T(1, 1) = ONE 00511 DO J = 2, Q 00512 V1T(1,J) = ZERO 00513 V1T(J,1) = ZERO 00514 END DO 00515 CALL DORGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1), 00516 $ WORK(IORGLQ), LORGLQWORK, INFO ) 00517 END IF 00518 IF( WANTV2T .AND. M-Q .GT. 0 ) THEN 00519 CALL DLACPY( 'U', P, M-Q, X12, LDX12, V2T, LDV2T ) 00520 CALL DLACPY( 'U', M-P-Q, M-P-Q, X22(Q+1,P+1), LDX22, 00521 $ V2T(P+1,P+1), LDV2T ) 00522 CALL DORGLQ( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2), 00523 $ WORK(IORGLQ), LORGLQWORK, INFO ) 00524 END IF 00525 ELSE 00526 IF( WANTU1 .AND. P .GT. 0 ) THEN 00527 CALL DLACPY( 'U', Q, P, X11, LDX11, U1, LDU1 ) 00528 CALL DORGLQ( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGLQ), 00529 $ LORGLQWORK, INFO) 00530 END IF 00531 IF( WANTU2 .AND. M-P .GT. 0 ) THEN 00532 CALL DLACPY( 'U', Q, M-P, X21, LDX21, U2, LDU2 ) 00533 CALL DORGLQ( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2), 00534 $ WORK(IORGLQ), LORGLQWORK, INFO ) 00535 END IF 00536 IF( WANTV1T .AND. Q .GT. 0 ) THEN 00537 CALL DLACPY( 'L', Q-1, Q-1, X11(2,1), LDX11, V1T(2,2), 00538 $ LDV1T ) 00539 V1T(1, 1) = ONE 00540 DO J = 2, Q 00541 V1T(1,J) = ZERO 00542 V1T(J,1) = ZERO 00543 END DO 00544 CALL DORGQR( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1), 00545 $ WORK(IORGQR), LORGQRWORK, INFO ) 00546 END IF 00547 IF( WANTV2T .AND. M-Q .GT. 0 ) THEN 00548 CALL DLACPY( 'L', M-Q, P, X12, LDX12, V2T, LDV2T ) 00549 CALL DLACPY( 'L', M-P-Q, M-P-Q, X22(P+1,Q+1), LDX22, 00550 $ V2T(P+1,P+1), LDV2T ) 00551 CALL DORGQR( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2), 00552 $ WORK(IORGQR), LORGQRWORK, INFO ) 00553 END IF 00554 END IF 00555 * 00556 * Compute the CSD of the matrix in bidiagonal-block form 00557 * 00558 CALL DBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, THETA, 00559 $ WORK(IPHI), U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, 00560 $ LDV2T, WORK(IB11D), WORK(IB11E), WORK(IB12D), 00561 $ WORK(IB12E), WORK(IB21D), WORK(IB21E), WORK(IB22D), 00562 $ WORK(IB22E), WORK(IBBCSD), LBBCSDWORK, INFO ) 00563 * 00564 * Permute rows and columns to place identity submatrices in top- 00565 * left corner of (1,1)-block and/or bottom-right corner of (1,2)- 00566 * block and/or bottom-right corner of (2,1)-block and/or top-left 00567 * corner of (2,2)-block 00568 * 00569 IF( Q .GT. 0 .AND. WANTU2 ) THEN 00570 DO I = 1, Q 00571 IWORK(I) = M - P - Q + I 00572 END DO 00573 DO I = Q + 1, M - P 00574 IWORK(I) = I - Q 00575 END DO 00576 IF( COLMAJOR ) THEN 00577 CALL DLAPMT( .FALSE., M-P, M-P, U2, LDU2, IWORK ) 00578 ELSE 00579 CALL DLAPMR( .FALSE., M-P, M-P, U2, LDU2, IWORK ) 00580 END IF 00581 END IF 00582 IF( M .GT. 0 .AND. WANTV2T ) THEN 00583 DO I = 1, P 00584 IWORK(I) = M - P - Q + I 00585 END DO 00586 DO I = P + 1, M - Q 00587 IWORK(I) = I - P 00588 END DO 00589 IF( .NOT. COLMAJOR ) THEN 00590 CALL DLAPMT( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK ) 00591 ELSE 00592 CALL DLAPMR( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK ) 00593 END IF 00594 END IF 00595 * 00596 RETURN 00597 * 00598 * End DORCSD 00599 * 00600 END 00601
1.7.6.1 
