![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SORCSD 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SORCSD + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sorcsd.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sorcsd.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sorcsd.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * RECURSIVE SUBROUTINE SORCSD( 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 * REAL THETA( * ) 00035 * REAL 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 *> SORCSD 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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: SBBCSD 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 realOTHERcomputational 00295 * 00296 * ===================================================================== 00297 RECURSIVE SUBROUTINE SORCSD( 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 REAL THETA( * ) 00316 REAL 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 REAL ONE, ZERO 00326 PARAMETER ( ONE = 1.0E+0, 00327 $ ZERO = 0.0E+0 ) 00328 * .. 00329 * .. Local Arrays .. 00330 REAL DUMMY(1) 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 SBBCSD, SLACPY, SLAPMR, SLAPMT, SLASCL, SLASET, 00347 $ SORBDB, SORGLQ, SORGQR, 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 SORCSD( 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 SORCSD( 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 SORGQR( M-Q, M-Q, M-Q, DUMMY, MAX(1,M-Q), DUMMY, WORK, -1, 00434 $ CHILDINFO ) 00435 LORGQRWORKOPT = INT( WORK(1) ) 00436 LORGQRWORKMIN = MAX( 1, M - Q ) 00437 IORGLQ = ITAUQ2 + MAX( 1, M - Q ) 00438 CALL SORGLQ( M-Q, M-Q, M-Q, DUMMY, MAX(1,M-Q), DUMMY, WORK, -1, 00439 $ CHILDINFO ) 00440 LORGLQWORKOPT = INT( WORK(1) ) 00441 LORGLQWORKMIN = MAX( 1, M - Q ) 00442 IORBDB = ITAUQ2 + MAX( 1, M - Q ) 00443 CALL SORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, 00444 $ X21, LDX21, X22, LDX22, DUMMY, DUMMY, DUMMY, DUMMY, DUMMY, 00445 $ DUMMY,WORK,-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 SBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, 00458 $ DUMMY, DUMMY, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, 00459 $ LDV2T, DUMMY, DUMMY, DUMMY, DUMMY, DUMMY, DUMMY, 00460 $ DUMMY, DUMMY, WORK, -1, CHILDINFO ) 00461 LBBCSDWORKOPT = INT( WORK(1) ) 00462 LBBCSDWORKMIN = LBBCSDWORKOPT 00463 LWORKOPT = MAX( IORGQR + LORGQRWORKOPT, IORGLQ + LORGLQWORKOPT, 00464 $ IORBDB + LORBDBWORKOPT, IBBCSD + LBBCSDWORKOPT ) - 1 00465 LWORKMIN = MAX( IORGQR + LORGQRWORKMIN, IORGLQ + LORGLQWORKMIN, 00466 $ IORBDB + LORBDBWORKOPT, IBBCSD + LBBCSDWORKMIN ) - 1 00467 WORK(1) = MAX(LWORKOPT,LWORKMIN) 00468 * 00469 IF( LWORK .LT. LWORKMIN .AND. .NOT. LQUERY ) THEN 00470 INFO = -22 00471 ELSE 00472 LORGQRWORK = LWORK - IORGQR + 1 00473 LORGLQWORK = LWORK - IORGLQ + 1 00474 LORBDBWORK = LWORK - IORBDB + 1 00475 LBBCSDWORK = LWORK - IBBCSD + 1 00476 END IF 00477 END IF 00478 * 00479 * Abort if any illegal arguments 00480 * 00481 IF( INFO .NE. 0 ) THEN 00482 CALL XERBLA( 'SORCSD', -INFO ) 00483 RETURN 00484 ELSE IF( LQUERY ) THEN 00485 RETURN 00486 END IF 00487 * 00488 * Transform to bidiagonal block form 00489 * 00490 CALL SORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21, 00491 $ LDX21, X22, LDX22, THETA, WORK(IPHI), WORK(ITAUP1), 00492 $ WORK(ITAUP2), WORK(ITAUQ1), WORK(ITAUQ2), 00493 $ WORK(IORBDB), LORBDBWORK, CHILDINFO ) 00494 * 00495 * Accumulate Householder reflectors 00496 * 00497 IF( COLMAJOR ) THEN 00498 IF( WANTU1 .AND. P .GT. 0 ) THEN 00499 CALL SLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 ) 00500 CALL SORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR), 00501 $ LORGQRWORK, INFO) 00502 END IF 00503 IF( WANTU2 .AND. M-P .GT. 0 ) THEN 00504 CALL SLACPY( 'L', M-P, Q, X21, LDX21, U2, LDU2 ) 00505 CALL SORGQR( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2), 00506 $ WORK(IORGQR), LORGQRWORK, INFO ) 00507 END IF 00508 IF( WANTV1T .AND. Q .GT. 0 ) THEN 00509 CALL SLACPY( 'U', Q-1, Q-1, X11(1,2), LDX11, V1T(2,2), 00510 $ LDV1T ) 00511 V1T(1, 1) = ONE 00512 DO J = 2, Q 00513 V1T(1,J) = ZERO 00514 V1T(J,1) = ZERO 00515 END DO 00516 CALL SORGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1), 00517 $ WORK(IORGLQ), LORGLQWORK, INFO ) 00518 END IF 00519 IF( WANTV2T .AND. M-Q .GT. 0 ) THEN 00520 CALL SLACPY( 'U', P, M-Q, X12, LDX12, V2T, LDV2T ) 00521 CALL SLACPY( 'U', M-P-Q, M-P-Q, X22(Q+1,P+1), LDX22, 00522 $ V2T(P+1,P+1), LDV2T ) 00523 CALL SORGLQ( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2), 00524 $ WORK(IORGLQ), LORGLQWORK, INFO ) 00525 END IF 00526 ELSE 00527 IF( WANTU1 .AND. P .GT. 0 ) THEN 00528 CALL SLACPY( 'U', Q, P, X11, LDX11, U1, LDU1 ) 00529 CALL SORGLQ( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGLQ), 00530 $ LORGLQWORK, INFO) 00531 END IF 00532 IF( WANTU2 .AND. M-P .GT. 0 ) THEN 00533 CALL SLACPY( 'U', Q, M-P, X21, LDX21, U2, LDU2 ) 00534 CALL SORGLQ( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2), 00535 $ WORK(IORGLQ), LORGLQWORK, INFO ) 00536 END IF 00537 IF( WANTV1T .AND. Q .GT. 0 ) THEN 00538 CALL SLACPY( 'L', Q-1, Q-1, X11(2,1), LDX11, V1T(2,2), 00539 $ LDV1T ) 00540 V1T(1, 1) = ONE 00541 DO J = 2, Q 00542 V1T(1,J) = ZERO 00543 V1T(J,1) = ZERO 00544 END DO 00545 CALL SORGQR( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1), 00546 $ WORK(IORGQR), LORGQRWORK, INFO ) 00547 END IF 00548 IF( WANTV2T .AND. M-Q .GT. 0 ) THEN 00549 CALL SLACPY( 'L', M-Q, P, X12, LDX12, V2T, LDV2T ) 00550 CALL SLACPY( 'L', M-P-Q, M-P-Q, X22(P+1,Q+1), LDX22, 00551 $ V2T(P+1,P+1), LDV2T ) 00552 CALL SORGQR( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2), 00553 $ WORK(IORGQR), LORGQRWORK, INFO ) 00554 END IF 00555 END IF 00556 * 00557 * Compute the CSD of the matrix in bidiagonal-block form 00558 * 00559 CALL SBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, THETA, 00560 $ WORK(IPHI), U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, 00561 $ LDV2T, WORK(IB11D), WORK(IB11E), WORK(IB12D), 00562 $ WORK(IB12E), WORK(IB21D), WORK(IB21E), WORK(IB22D), 00563 $ WORK(IB22E), WORK(IBBCSD), LBBCSDWORK, INFO ) 00564 * 00565 * Permute rows and columns to place identity submatrices in top- 00566 * left corner of (1,1)-block and/or bottom-right corner of (1,2)- 00567 * block and/or bottom-right corner of (2,1)-block and/or top-left 00568 * corner of (2,2)-block 00569 * 00570 IF( Q .GT. 0 .AND. WANTU2 ) THEN 00571 DO I = 1, Q 00572 IWORK(I) = M - P - Q + I 00573 END DO 00574 DO I = Q + 1, M - P 00575 IWORK(I) = I - Q 00576 END DO 00577 IF( COLMAJOR ) THEN 00578 CALL SLAPMT( .FALSE., M-P, M-P, U2, LDU2, IWORK ) 00579 ELSE 00580 CALL SLAPMR( .FALSE., M-P, M-P, U2, LDU2, IWORK ) 00581 END IF 00582 END IF 00583 IF( M .GT. 0 .AND. WANTV2T ) THEN 00584 DO I = 1, P 00585 IWORK(I) = M - P - Q + I 00586 END DO 00587 DO I = P + 1, M - Q 00588 IWORK(I) = I - P 00589 END DO 00590 IF( .NOT. COLMAJOR ) THEN 00591 CALL SLAPMT( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK ) 00592 ELSE 00593 CALL SLAPMR( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK ) 00594 END IF 00595 END IF 00596 * 00597 RETURN 00598 * 00599 * End SORCSD 00600 * 00601 END 00602