![]() |
LAPACK
3.4.1
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 ONE, ZERO 00326 PARAMETER ( ONE = 1.0D0, 00327 $ ZERO = 0.0D0 ) 00328 * .. 00329 * .. Local Scalars .. 00330 CHARACTER TRANST, SIGNST 00331 INTEGER CHILDINFO, I, IB11D, IB11E, IB12D, IB12E, 00332 $ IB21D, IB21E, IB22D, IB22E, IBBCSD, IORBDB, 00333 $ IORGLQ, IORGQR, IPHI, ITAUP1, ITAUP2, ITAUQ1, 00334 $ ITAUQ2, J, LBBCSDWORK, LBBCSDWORKMIN, 00335 $ LBBCSDWORKOPT, LORBDBWORK, LORBDBWORKMIN, 00336 $ LORBDBWORKOPT, LORGLQWORK, LORGLQWORKMIN, 00337 $ LORGLQWORKOPT, LORGQRWORK, LORGQRWORKMIN, 00338 $ LORGQRWORKOPT, LWORKMIN, LWORKOPT 00339 LOGICAL COLMAJOR, DEFAULTSIGNS, LQUERY, WANTU1, WANTU2, 00340 $ WANTV1T, WANTV2T 00341 * .. 00342 * .. External Subroutines .. 00343 EXTERNAL DBBCSD, DLACPY, DLAPMR, DLAPMT, DLASCL, DLASET, 00344 $ DORBDB, DORGLQ, DORGQR, XERBLA 00345 * .. 00346 * .. External Functions .. 00347 LOGICAL LSAME 00348 EXTERNAL LSAME 00349 * .. 00350 * .. Intrinsic Functions 00351 INTRINSIC INT, MAX, MIN 00352 * .. 00353 * .. Executable Statements .. 00354 * 00355 * Test input arguments 00356 * 00357 INFO = 0 00358 WANTU1 = LSAME( JOBU1, 'Y' ) 00359 WANTU2 = LSAME( JOBU2, 'Y' ) 00360 WANTV1T = LSAME( JOBV1T, 'Y' ) 00361 WANTV2T = LSAME( JOBV2T, 'Y' ) 00362 COLMAJOR = .NOT. LSAME( TRANS, 'T' ) 00363 DEFAULTSIGNS = .NOT. LSAME( SIGNS, 'O' ) 00364 LQUERY = LWORK .EQ. -1 00365 IF( M .LT. 0 ) THEN 00366 INFO = -7 00367 ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN 00368 INFO = -8 00369 ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN 00370 INFO = -9 00371 ELSE IF( ( COLMAJOR .AND. LDX11 .LT. MAX(1,P) ) .OR. 00372 $ ( .NOT.COLMAJOR .AND. LDX11 .LT. MAX(1,Q) ) ) THEN 00373 INFO = -11 00374 ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN 00375 INFO = -20 00376 ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN 00377 INFO = -22 00378 ELSE IF( WANTV1T .AND. LDV1T .LT. Q ) THEN 00379 INFO = -24 00380 ELSE IF( WANTV2T .AND. LDV2T .LT. M-Q ) THEN 00381 INFO = -26 00382 END IF 00383 * 00384 * Work with transpose if convenient 00385 * 00386 IF( INFO .EQ. 0 .AND. MIN( P, M-P ) .LT. MIN( Q, M-Q ) ) THEN 00387 IF( COLMAJOR ) THEN 00388 TRANST = 'T' 00389 ELSE 00390 TRANST = 'N' 00391 END IF 00392 IF( DEFAULTSIGNS ) THEN 00393 SIGNST = 'O' 00394 ELSE 00395 SIGNST = 'D' 00396 END IF 00397 CALL DORCSD( JOBV1T, JOBV2T, JOBU1, JOBU2, TRANST, SIGNST, M, 00398 $ Q, P, X11, LDX11, X21, LDX21, X12, LDX12, X22, 00399 $ LDX22, THETA, V1T, LDV1T, V2T, LDV2T, U1, LDU1, 00400 $ U2, LDU2, WORK, LWORK, IWORK, INFO ) 00401 RETURN 00402 END IF 00403 * 00404 * Work with permutation [ 0 I; I 0 ] * X * [ 0 I; I 0 ] if 00405 * convenient 00406 * 00407 IF( INFO .EQ. 0 .AND. M-Q .LT. Q ) THEN 00408 IF( DEFAULTSIGNS ) THEN 00409 SIGNST = 'O' 00410 ELSE 00411 SIGNST = 'D' 00412 END IF 00413 CALL DORCSD( JOBU2, JOBU1, JOBV2T, JOBV1T, TRANS, SIGNST, M, 00414 $ M-P, M-Q, X22, LDX22, X21, LDX21, X12, LDX12, X11, 00415 $ LDX11, THETA, U2, LDU2, U1, LDU1, V2T, LDV2T, V1T, 00416 $ LDV1T, WORK, LWORK, IWORK, INFO ) 00417 RETURN 00418 END IF 00419 * 00420 * Compute workspace 00421 * 00422 IF( INFO .EQ. 0 ) THEN 00423 * 00424 IPHI = 2 00425 ITAUP1 = IPHI + MAX( 1, Q - 1 ) 00426 ITAUP2 = ITAUP1 + MAX( 1, P ) 00427 ITAUQ1 = ITAUP2 + MAX( 1, M - P ) 00428 ITAUQ2 = ITAUQ1 + MAX( 1, Q ) 00429 IORGQR = ITAUQ2 + MAX( 1, M - Q ) 00430 CALL DORGQR( M-Q, M-Q, M-Q, 0, MAX(1,M-Q), 0, WORK, -1, 00431 $ CHILDINFO ) 00432 LORGQRWORKOPT = INT( WORK(1) ) 00433 LORGQRWORKMIN = MAX( 1, M - Q ) 00434 IORGLQ = ITAUQ2 + MAX( 1, M - Q ) 00435 CALL DORGLQ( M-Q, M-Q, M-Q, 0, MAX(1,M-Q), 0, WORK, -1, 00436 $ CHILDINFO ) 00437 LORGLQWORKOPT = INT( WORK(1) ) 00438 LORGLQWORKMIN = MAX( 1, M - Q ) 00439 IORBDB = ITAUQ2 + MAX( 1, M - Q ) 00440 CALL DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, 00441 $ X21, LDX21, X22, LDX22, 0, 0, 0, 0, 0, 0, WORK, 00442 $ -1, CHILDINFO ) 00443 LORBDBWORKOPT = INT( WORK(1) ) 00444 LORBDBWORKMIN = LORBDBWORKOPT 00445 IB11D = ITAUQ2 + MAX( 1, M - Q ) 00446 IB11E = IB11D + MAX( 1, Q ) 00447 IB12D = IB11E + MAX( 1, Q - 1 ) 00448 IB12E = IB12D + MAX( 1, Q ) 00449 IB21D = IB12E + MAX( 1, Q - 1 ) 00450 IB21E = IB21D + MAX( 1, Q ) 00451 IB22D = IB21E + MAX( 1, Q - 1 ) 00452 IB22E = IB22D + MAX( 1, Q ) 00453 IBBCSD = IB22E + MAX( 1, Q - 1 ) 00454 CALL DBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, 0, 00455 $ 0, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, 0, 00456 $ 0, 0, 0, 0, 0, 0, 0, WORK, -1, CHILDINFO ) 00457 LBBCSDWORKOPT = INT( WORK(1) ) 00458 LBBCSDWORKMIN = LBBCSDWORKOPT 00459 LWORKOPT = MAX( IORGQR + LORGQRWORKOPT, IORGLQ + LORGLQWORKOPT, 00460 $ IORBDB + LORBDBWORKOPT, IBBCSD + LBBCSDWORKOPT ) - 1 00461 LWORKMIN = MAX( IORGQR + LORGQRWORKMIN, IORGLQ + LORGLQWORKMIN, 00462 $ IORBDB + LORBDBWORKOPT, IBBCSD + LBBCSDWORKMIN ) - 1 00463 WORK(1) = MAX(LWORKOPT,LWORKMIN) 00464 * 00465 IF( LWORK .LT. LWORKMIN .AND. .NOT. LQUERY ) THEN 00466 INFO = -22 00467 ELSE 00468 LORGQRWORK = LWORK - IORGQR + 1 00469 LORGLQWORK = LWORK - IORGLQ + 1 00470 LORBDBWORK = LWORK - IORBDB + 1 00471 LBBCSDWORK = LWORK - IBBCSD + 1 00472 END IF 00473 END IF 00474 * 00475 * Abort if any illegal arguments 00476 * 00477 IF( INFO .NE. 0 ) THEN 00478 CALL XERBLA( 'DORCSD', -INFO ) 00479 RETURN 00480 ELSE IF( LQUERY ) THEN 00481 RETURN 00482 END IF 00483 * 00484 * Transform to bidiagonal block form 00485 * 00486 CALL DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21, 00487 $ LDX21, X22, LDX22, THETA, WORK(IPHI), WORK(ITAUP1), 00488 $ WORK(ITAUP2), WORK(ITAUQ1), WORK(ITAUQ2), 00489 $ WORK(IORBDB), LORBDBWORK, CHILDINFO ) 00490 * 00491 * Accumulate Householder reflectors 00492 * 00493 IF( COLMAJOR ) THEN 00494 IF( WANTU1 .AND. P .GT. 0 ) THEN 00495 CALL DLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 ) 00496 CALL DORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR), 00497 $ LORGQRWORK, INFO) 00498 END IF 00499 IF( WANTU2 .AND. M-P .GT. 0 ) THEN 00500 CALL DLACPY( 'L', M-P, Q, X21, LDX21, U2, LDU2 ) 00501 CALL DORGQR( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2), 00502 $ WORK(IORGQR), LORGQRWORK, INFO ) 00503 END IF 00504 IF( WANTV1T .AND. Q .GT. 0 ) THEN 00505 CALL DLACPY( 'U', Q-1, Q-1, X11(1,2), LDX11, V1T(2,2), 00506 $ LDV1T ) 00507 V1T(1, 1) = ONE 00508 DO J = 2, Q 00509 V1T(1,J) = ZERO 00510 V1T(J,1) = ZERO 00511 END DO 00512 CALL DORGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1), 00513 $ WORK(IORGLQ), LORGLQWORK, INFO ) 00514 END IF 00515 IF( WANTV2T .AND. M-Q .GT. 0 ) THEN 00516 CALL DLACPY( 'U', P, M-Q, X12, LDX12, V2T, LDV2T ) 00517 CALL DLACPY( 'U', M-P-Q, M-P-Q, X22(Q+1,P+1), LDX22, 00518 $ V2T(P+1,P+1), LDV2T ) 00519 CALL DORGLQ( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2), 00520 $ WORK(IORGLQ), LORGLQWORK, INFO ) 00521 END IF 00522 ELSE 00523 IF( WANTU1 .AND. P .GT. 0 ) THEN 00524 CALL DLACPY( 'U', Q, P, X11, LDX11, U1, LDU1 ) 00525 CALL DORGLQ( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGLQ), 00526 $ LORGLQWORK, INFO) 00527 END IF 00528 IF( WANTU2 .AND. M-P .GT. 0 ) THEN 00529 CALL DLACPY( 'U', Q, M-P, X21, LDX21, U2, LDU2 ) 00530 CALL DORGLQ( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2), 00531 $ WORK(IORGLQ), LORGLQWORK, INFO ) 00532 END IF 00533 IF( WANTV1T .AND. Q .GT. 0 ) THEN 00534 CALL DLACPY( 'L', Q-1, Q-1, X11(2,1), LDX11, V1T(2,2), 00535 $ LDV1T ) 00536 V1T(1, 1) = ONE 00537 DO J = 2, Q 00538 V1T(1,J) = ZERO 00539 V1T(J,1) = ZERO 00540 END DO 00541 CALL DORGQR( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1), 00542 $ WORK(IORGQR), LORGQRWORK, INFO ) 00543 END IF 00544 IF( WANTV2T .AND. M-Q .GT. 0 ) THEN 00545 CALL DLACPY( 'L', M-Q, P, X12, LDX12, V2T, LDV2T ) 00546 CALL DLACPY( 'L', M-P-Q, M-P-Q, X22(P+1,Q+1), LDX22, 00547 $ V2T(P+1,P+1), LDV2T ) 00548 CALL DORGQR( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2), 00549 $ WORK(IORGQR), LORGQRWORK, INFO ) 00550 END IF 00551 END IF 00552 * 00553 * Compute the CSD of the matrix in bidiagonal-block form 00554 * 00555 CALL DBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, THETA, 00556 $ WORK(IPHI), U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, 00557 $ LDV2T, WORK(IB11D), WORK(IB11E), WORK(IB12D), 00558 $ WORK(IB12E), WORK(IB21D), WORK(IB21E), WORK(IB22D), 00559 $ WORK(IB22E), WORK(IBBCSD), LBBCSDWORK, INFO ) 00560 * 00561 * Permute rows and columns to place identity submatrices in top- 00562 * left corner of (1,1)-block and/or bottom-right corner of (1,2)- 00563 * block and/or bottom-right corner of (2,1)-block and/or top-left 00564 * corner of (2,2)-block 00565 * 00566 IF( Q .GT. 0 .AND. WANTU2 ) THEN 00567 DO I = 1, Q 00568 IWORK(I) = M - P - Q + I 00569 END DO 00570 DO I = Q + 1, M - P 00571 IWORK(I) = I - Q 00572 END DO 00573 IF( COLMAJOR ) THEN 00574 CALL DLAPMT( .FALSE., M-P, M-P, U2, LDU2, IWORK ) 00575 ELSE 00576 CALL DLAPMR( .FALSE., M-P, M-P, U2, LDU2, IWORK ) 00577 END IF 00578 END IF 00579 IF( M .GT. 0 .AND. WANTV2T ) THEN 00580 DO I = 1, P 00581 IWORK(I) = M - P - Q + I 00582 END DO 00583 DO I = P + 1, M - Q 00584 IWORK(I) = I - P 00585 END DO 00586 IF( .NOT. COLMAJOR ) THEN 00587 CALL DLAPMT( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK ) 00588 ELSE 00589 CALL DLAPMR( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK ) 00590 END IF 00591 END IF 00592 * 00593 RETURN 00594 * 00595 * End DORCSD 00596 * 00597 END 00598