![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZUNCSD 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZUNCSD + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zuncsd.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zuncsd.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zuncsd.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * RECURSIVE SUBROUTINE ZUNCSD( 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 * DOUBLE PRECISION THETA( * ) 00036 * DOUBLE PRECISION RWORK( * ) 00037 * COMPLEX*16 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 *> ZUNCSD 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*16 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*16 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*16 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*16 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 DOUBLE PRECISION 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*16 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*16 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*16 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*16 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*16 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 DOUBLE PRECISION 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: ZBBCSD 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 complex16OTHERcomputational 00314 * 00315 * ===================================================================== 00316 RECURSIVE SUBROUTINE ZUNCSD( 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 DOUBLE PRECISION THETA( * ) 00336 DOUBLE PRECISION RWORK( * ) 00337 COMPLEX*16 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 COMPLEX*16 ONE, ZERO 00347 PARAMETER ( ONE = (1.0D0,0.0D0), 00348 $ ZERO = (0.0D0,0.0D0) ) 00349 * .. 00350 * .. Local Scalars .. 00351 CHARACTER TRANST, SIGNST 00352 INTEGER CHILDINFO, I, IB11D, IB11E, IB12D, IB12E, 00353 $ IB21D, IB21E, IB22D, IB22E, IBBCSD, IORBDB, 00354 $ IORGLQ, IORGQR, IPHI, ITAUP1, ITAUP2, ITAUQ1, 00355 $ ITAUQ2, J, LBBCSDWORK, LBBCSDWORKMIN, 00356 $ LBBCSDWORKOPT, LORBDBWORK, LORBDBWORKMIN, 00357 $ LORBDBWORKOPT, LORGLQWORK, LORGLQWORKMIN, 00358 $ LORGLQWORKOPT, LORGQRWORK, LORGQRWORKMIN, 00359 $ LORGQRWORKOPT, LWORKMIN, LWORKOPT 00360 LOGICAL COLMAJOR, DEFAULTSIGNS, LQUERY, WANTU1, WANTU2, 00361 $ WANTV1T, WANTV2T 00362 INTEGER LRWORKMIN, LRWORKOPT 00363 LOGICAL LRQUERY 00364 * .. 00365 * .. External Subroutines .. 00366 EXTERNAL XERBLA, ZBBCSD, ZLACPY, ZLAPMR, ZLAPMT, ZLASCL, 00367 $ ZLASET, ZUNBDB, ZUNGLQ, ZUNGQR 00368 * .. 00369 * .. External Functions .. 00370 LOGICAL LSAME 00371 EXTERNAL LSAME 00372 * .. 00373 * .. Intrinsic Functions 00374 INTRINSIC COS, INT, MAX, MIN, SIN 00375 * .. 00376 * .. Executable Statements .. 00377 * 00378 * Test input arguments 00379 * 00380 INFO = 0 00381 WANTU1 = LSAME( JOBU1, 'Y' ) 00382 WANTU2 = LSAME( JOBU2, 'Y' ) 00383 WANTV1T = LSAME( JOBV1T, 'Y' ) 00384 WANTV2T = LSAME( JOBV2T, 'Y' ) 00385 COLMAJOR = .NOT. LSAME( TRANS, 'T' ) 00386 DEFAULTSIGNS = .NOT. LSAME( SIGNS, 'O' ) 00387 LQUERY = LWORK .EQ. -1 00388 LRQUERY = LRWORK .EQ. -1 00389 IF( M .LT. 0 ) THEN 00390 INFO = -7 00391 ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN 00392 INFO = -8 00393 ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN 00394 INFO = -9 00395 ELSE IF( ( COLMAJOR .AND. LDX11 .LT. MAX(1,P) ) .OR. 00396 $ ( .NOT.COLMAJOR .AND. LDX11 .LT. MAX(1,Q) ) ) THEN 00397 INFO = -11 00398 ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN 00399 INFO = -20 00400 ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN 00401 INFO = -22 00402 ELSE IF( WANTV1T .AND. LDV1T .LT. Q ) THEN 00403 INFO = -24 00404 ELSE IF( WANTV2T .AND. LDV2T .LT. M-Q ) THEN 00405 INFO = -26 00406 END IF 00407 * 00408 * Work with transpose if convenient 00409 * 00410 IF( INFO .EQ. 0 .AND. MIN( P, M-P ) .LT. MIN( Q, M-Q ) ) THEN 00411 IF( COLMAJOR ) THEN 00412 TRANST = 'T' 00413 ELSE 00414 TRANST = 'N' 00415 END IF 00416 IF( DEFAULTSIGNS ) THEN 00417 SIGNST = 'O' 00418 ELSE 00419 SIGNST = 'D' 00420 END IF 00421 CALL ZUNCSD( JOBV1T, JOBV2T, JOBU1, JOBU2, TRANST, SIGNST, M, 00422 $ Q, P, X11, LDX11, X21, LDX21, X12, LDX12, X22, 00423 $ LDX22, THETA, V1T, LDV1T, V2T, LDV2T, U1, LDU1, 00424 $ U2, LDU2, WORK, LWORK, RWORK, LRWORK, IWORK, 00425 $ INFO ) 00426 RETURN 00427 END IF 00428 * 00429 * Work with permutation [ 0 I; I 0 ] * X * [ 0 I; I 0 ] if 00430 * convenient 00431 * 00432 IF( INFO .EQ. 0 .AND. M-Q .LT. Q ) THEN 00433 IF( DEFAULTSIGNS ) THEN 00434 SIGNST = 'O' 00435 ELSE 00436 SIGNST = 'D' 00437 END IF 00438 CALL ZUNCSD( JOBU2, JOBU1, JOBV2T, JOBV1T, TRANS, SIGNST, M, 00439 $ M-P, M-Q, X22, LDX22, X21, LDX21, X12, LDX12, X11, 00440 $ LDX11, THETA, U2, LDU2, U1, LDU1, V2T, LDV2T, V1T, 00441 $ LDV1T, WORK, LWORK, RWORK, LRWORK, IWORK, INFO ) 00442 RETURN 00443 END IF 00444 * 00445 * Compute workspace 00446 * 00447 IF( INFO .EQ. 0 ) THEN 00448 * 00449 * Real workspace 00450 * 00451 IPHI = 2 00452 IB11D = IPHI + MAX( 1, Q - 1 ) 00453 IB11E = IB11D + MAX( 1, Q ) 00454 IB12D = IB11E + MAX( 1, Q - 1 ) 00455 IB12E = IB12D + MAX( 1, Q ) 00456 IB21D = IB12E + MAX( 1, Q - 1 ) 00457 IB21E = IB21D + MAX( 1, Q ) 00458 IB22D = IB21E + MAX( 1, Q - 1 ) 00459 IB22E = IB22D + MAX( 1, Q ) 00460 IBBCSD = IB22E + MAX( 1, Q - 1 ) 00461 CALL ZBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, 0, 00462 $ 0, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, 0, 00463 $ 0, 0, 0, 0, 0, 0, 0, RWORK, -1, CHILDINFO ) 00464 LBBCSDWORKOPT = INT( RWORK(1) ) 00465 LBBCSDWORKMIN = LBBCSDWORKOPT 00466 LRWORKOPT = IBBCSD + LBBCSDWORKOPT - 1 00467 LRWORKMIN = IBBCSD + LBBCSDWORKMIN - 1 00468 RWORK(1) = LRWORKOPT 00469 * 00470 * Complex workspace 00471 * 00472 ITAUP1 = 2 00473 ITAUP2 = ITAUP1 + MAX( 1, P ) 00474 ITAUQ1 = ITAUP2 + MAX( 1, M - P ) 00475 ITAUQ2 = ITAUQ1 + MAX( 1, Q ) 00476 IORGQR = ITAUQ2 + MAX( 1, M - Q ) 00477 CALL ZUNGQR( M-Q, M-Q, M-Q, 0, MAX(1,M-Q), 0, WORK, -1, 00478 $ CHILDINFO ) 00479 LORGQRWORKOPT = INT( WORK(1) ) 00480 LORGQRWORKMIN = MAX( 1, M - Q ) 00481 IORGLQ = ITAUQ2 + MAX( 1, M - Q ) 00482 CALL ZUNGLQ( M-Q, M-Q, M-Q, 0, MAX(1,M-Q), 0, WORK, -1, 00483 $ CHILDINFO ) 00484 LORGLQWORKOPT = INT( WORK(1) ) 00485 LORGLQWORKMIN = MAX( 1, M - Q ) 00486 IORBDB = ITAUQ2 + MAX( 1, M - Q ) 00487 CALL ZUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, 00488 $ X21, LDX21, X22, LDX22, 0, 0, 0, 0, 0, 0, WORK, 00489 $ -1, CHILDINFO ) 00490 LORBDBWORKOPT = INT( WORK(1) ) 00491 LORBDBWORKMIN = LORBDBWORKOPT 00492 LWORKOPT = MAX( IORGQR + LORGQRWORKOPT, IORGLQ + LORGLQWORKOPT, 00493 $ IORBDB + LORBDBWORKOPT ) - 1 00494 LWORKMIN = MAX( IORGQR + LORGQRWORKMIN, IORGLQ + LORGLQWORKMIN, 00495 $ IORBDB + LORBDBWORKMIN ) - 1 00496 WORK(1) = MAX(LWORKOPT,LWORKMIN) 00497 * 00498 IF( LWORK .LT. LWORKMIN 00499 $ .AND. .NOT. ( LQUERY .OR. LRQUERY ) ) THEN 00500 INFO = -22 00501 ELSE IF( LRWORK .LT. LRWORKMIN 00502 $ .AND. .NOT. ( LQUERY .OR. LRQUERY ) ) THEN 00503 INFO = -24 00504 ELSE 00505 LORGQRWORK = LWORK - IORGQR + 1 00506 LORGLQWORK = LWORK - IORGLQ + 1 00507 LORBDBWORK = LWORK - IORBDB + 1 00508 LBBCSDWORK = LRWORK - IBBCSD + 1 00509 END IF 00510 END IF 00511 * 00512 * Abort if any illegal arguments 00513 * 00514 IF( INFO .NE. 0 ) THEN 00515 CALL XERBLA( 'ZUNCSD', -INFO ) 00516 RETURN 00517 ELSE IF( LQUERY .OR. LRQUERY ) THEN 00518 RETURN 00519 END IF 00520 * 00521 * Transform to bidiagonal block form 00522 * 00523 CALL ZUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21, 00524 $ LDX21, X22, LDX22, THETA, RWORK(IPHI), WORK(ITAUP1), 00525 $ WORK(ITAUP2), WORK(ITAUQ1), WORK(ITAUQ2), 00526 $ WORK(IORBDB), LORBDBWORK, CHILDINFO ) 00527 * 00528 * Accumulate Householder reflectors 00529 * 00530 IF( COLMAJOR ) THEN 00531 IF( WANTU1 .AND. P .GT. 0 ) THEN 00532 CALL ZLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 ) 00533 CALL ZUNGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR), 00534 $ LORGQRWORK, INFO) 00535 END IF 00536 IF( WANTU2 .AND. M-P .GT. 0 ) THEN 00537 CALL ZLACPY( 'L', M-P, Q, X21, LDX21, U2, LDU2 ) 00538 CALL ZUNGQR( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2), 00539 $ WORK(IORGQR), LORGQRWORK, INFO ) 00540 END IF 00541 IF( WANTV1T .AND. Q .GT. 0 ) THEN 00542 CALL ZLACPY( 'U', Q-1, Q-1, X11(1,2), LDX11, V1T(2,2), 00543 $ LDV1T ) 00544 V1T(1, 1) = ONE 00545 DO J = 2, Q 00546 V1T(1,J) = ZERO 00547 V1T(J,1) = ZERO 00548 END DO 00549 CALL ZUNGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1), 00550 $ WORK(IORGLQ), LORGLQWORK, INFO ) 00551 END IF 00552 IF( WANTV2T .AND. M-Q .GT. 0 ) THEN 00553 CALL ZLACPY( 'U', P, M-Q, X12, LDX12, V2T, LDV2T ) 00554 CALL ZLACPY( 'U', M-P-Q, M-P-Q, X22(Q+1,P+1), LDX22, 00555 $ V2T(P+1,P+1), LDV2T ) 00556 CALL ZUNGLQ( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2), 00557 $ WORK(IORGLQ), LORGLQWORK, INFO ) 00558 END IF 00559 ELSE 00560 IF( WANTU1 .AND. P .GT. 0 ) THEN 00561 CALL ZLACPY( 'U', Q, P, X11, LDX11, U1, LDU1 ) 00562 CALL ZUNGLQ( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGLQ), 00563 $ LORGLQWORK, INFO) 00564 END IF 00565 IF( WANTU2 .AND. M-P .GT. 0 ) THEN 00566 CALL ZLACPY( 'U', Q, M-P, X21, LDX21, U2, LDU2 ) 00567 CALL ZUNGLQ( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2), 00568 $ WORK(IORGLQ), LORGLQWORK, INFO ) 00569 END IF 00570 IF( WANTV1T .AND. Q .GT. 0 ) THEN 00571 CALL ZLACPY( 'L', Q-1, Q-1, X11(2,1), LDX11, V1T(2,2), 00572 $ LDV1T ) 00573 V1T(1, 1) = ONE 00574 DO J = 2, Q 00575 V1T(1,J) = ZERO 00576 V1T(J,1) = ZERO 00577 END DO 00578 CALL ZUNGQR( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1), 00579 $ WORK(IORGQR), LORGQRWORK, INFO ) 00580 END IF 00581 IF( WANTV2T .AND. M-Q .GT. 0 ) THEN 00582 CALL ZLACPY( 'L', M-Q, P, X12, LDX12, V2T, LDV2T ) 00583 CALL ZLACPY( 'L', M-P-Q, M-P-Q, X22(P+1,Q+1), LDX22, 00584 $ V2T(P+1,P+1), LDV2T ) 00585 CALL ZUNGQR( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2), 00586 $ WORK(IORGQR), LORGQRWORK, INFO ) 00587 END IF 00588 END IF 00589 * 00590 * Compute the CSD of the matrix in bidiagonal-block form 00591 * 00592 CALL ZBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, THETA, 00593 $ RWORK(IPHI), U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, 00594 $ LDV2T, RWORK(IB11D), RWORK(IB11E), RWORK(IB12D), 00595 $ RWORK(IB12E), RWORK(IB21D), RWORK(IB21E), 00596 $ RWORK(IB22D), RWORK(IB22E), RWORK(IBBCSD), 00597 $ LBBCSDWORK, INFO ) 00598 * 00599 * Permute rows and columns to place identity submatrices in top- 00600 * left corner of (1,1)-block and/or bottom-right corner of (1,2)- 00601 * block and/or bottom-right corner of (2,1)-block and/or top-left 00602 * corner of (2,2)-block 00603 * 00604 IF( Q .GT. 0 .AND. WANTU2 ) THEN 00605 DO I = 1, Q 00606 IWORK(I) = M - P - Q + I 00607 END DO 00608 DO I = Q + 1, M - P 00609 IWORK(I) = I - Q 00610 END DO 00611 IF( COLMAJOR ) THEN 00612 CALL ZLAPMT( .FALSE., M-P, M-P, U2, LDU2, IWORK ) 00613 ELSE 00614 CALL ZLAPMR( .FALSE., M-P, M-P, U2, LDU2, IWORK ) 00615 END IF 00616 END IF 00617 IF( M .GT. 0 .AND. WANTV2T ) THEN 00618 DO I = 1, P 00619 IWORK(I) = M - P - Q + I 00620 END DO 00621 DO I = P + 1, M - Q 00622 IWORK(I) = I - P 00623 END DO 00624 IF( .NOT. COLMAJOR ) THEN 00625 CALL ZLAPMT( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK ) 00626 ELSE 00627 CALL ZLAPMR( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK ) 00628 END IF 00629 END IF 00630 * 00631 RETURN 00632 * 00633 * End ZUNCSD 00634 * 00635 END 00636