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