![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CGESDD 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CGESDD + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgesdd.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgesdd.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgesdd.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, 00022 * WORK, LWORK, RWORK, IWORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER JOBZ 00026 * INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N 00027 * .. 00028 * .. Array Arguments .. 00029 * INTEGER IWORK( * ) 00030 * REAL RWORK( * ), S( * ) 00031 * COMPLEX A( LDA, * ), U( LDU, * ), VT( LDVT, * ), 00032 * $ WORK( * ) 00033 * .. 00034 * 00035 * 00036 *> \par Purpose: 00037 * ============= 00038 *> 00039 *> \verbatim 00040 *> 00041 *> CGESDD computes the singular value decomposition (SVD) of a complex 00042 *> M-by-N matrix A, optionally computing the left and/or right singular 00043 *> vectors, by using divide-and-conquer method. The SVD is written 00044 *> 00045 *> A = U * SIGMA * conjugate-transpose(V) 00046 *> 00047 *> where SIGMA is an M-by-N matrix which is zero except for its 00048 *> min(m,n) diagonal elements, U is an M-by-M unitary matrix, and 00049 *> V is an N-by-N unitary matrix. The diagonal elements of SIGMA 00050 *> are the singular values of A; they are real and non-negative, and 00051 *> are returned in descending order. The first min(m,n) columns of 00052 *> U and V are the left and right singular vectors of A. 00053 *> 00054 *> Note that the routine returns VT = V**H, not V. 00055 *> 00056 *> The divide and conquer algorithm makes very mild assumptions about 00057 *> floating point arithmetic. It will work on machines with a guard 00058 *> digit in add/subtract, or on those binary machines without guard 00059 *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or 00060 *> Cray-2. It could conceivably fail on hexadecimal or decimal machines 00061 *> without guard digits, but we know of none. 00062 *> \endverbatim 00063 * 00064 * Arguments: 00065 * ========== 00066 * 00067 *> \param[in] JOBZ 00068 *> \verbatim 00069 *> JOBZ is CHARACTER*1 00070 *> Specifies options for computing all or part of the matrix U: 00071 *> = 'A': all M columns of U and all N rows of V**H are 00072 *> returned in the arrays U and VT; 00073 *> = 'S': the first min(M,N) columns of U and the first 00074 *> min(M,N) rows of V**H are returned in the arrays U 00075 *> and VT; 00076 *> = 'O': If M >= N, the first N columns of U are overwritten 00077 *> in the array A and all rows of V**H are returned in 00078 *> the array VT; 00079 *> otherwise, all columns of U are returned in the 00080 *> array U and the first M rows of V**H are overwritten 00081 *> in the array A; 00082 *> = 'N': no columns of U or rows of V**H are computed. 00083 *> \endverbatim 00084 *> 00085 *> \param[in] M 00086 *> \verbatim 00087 *> M is INTEGER 00088 *> The number of rows of the input matrix A. M >= 0. 00089 *> \endverbatim 00090 *> 00091 *> \param[in] N 00092 *> \verbatim 00093 *> N is INTEGER 00094 *> The number of columns of the input matrix A. N >= 0. 00095 *> \endverbatim 00096 *> 00097 *> \param[in,out] A 00098 *> \verbatim 00099 *> A is COMPLEX array, dimension (LDA,N) 00100 *> On entry, the M-by-N matrix A. 00101 *> On exit, 00102 *> if JOBZ = 'O', A is overwritten with the first N columns 00103 *> of U (the left singular vectors, stored 00104 *> columnwise) if M >= N; 00105 *> A is overwritten with the first M rows 00106 *> of V**H (the right singular vectors, stored 00107 *> rowwise) otherwise. 00108 *> if JOBZ .ne. 'O', the contents of A are destroyed. 00109 *> \endverbatim 00110 *> 00111 *> \param[in] LDA 00112 *> \verbatim 00113 *> LDA is INTEGER 00114 *> The leading dimension of the array A. LDA >= max(1,M). 00115 *> \endverbatim 00116 *> 00117 *> \param[out] S 00118 *> \verbatim 00119 *> S is REAL array, dimension (min(M,N)) 00120 *> The singular values of A, sorted so that S(i) >= S(i+1). 00121 *> \endverbatim 00122 *> 00123 *> \param[out] U 00124 *> \verbatim 00125 *> U is COMPLEX array, dimension (LDU,UCOL) 00126 *> UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N; 00127 *> UCOL = min(M,N) if JOBZ = 'S'. 00128 *> If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M 00129 *> unitary matrix U; 00130 *> if JOBZ = 'S', U contains the first min(M,N) columns of U 00131 *> (the left singular vectors, stored columnwise); 00132 *> if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced. 00133 *> \endverbatim 00134 *> 00135 *> \param[in] LDU 00136 *> \verbatim 00137 *> LDU is INTEGER 00138 *> The leading dimension of the array U. LDU >= 1; if 00139 *> JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M. 00140 *> \endverbatim 00141 *> 00142 *> \param[out] VT 00143 *> \verbatim 00144 *> VT is COMPLEX array, dimension (LDVT,N) 00145 *> If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the 00146 *> N-by-N unitary matrix V**H; 00147 *> if JOBZ = 'S', VT contains the first min(M,N) rows of 00148 *> V**H (the right singular vectors, stored rowwise); 00149 *> if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced. 00150 *> \endverbatim 00151 *> 00152 *> \param[in] LDVT 00153 *> \verbatim 00154 *> LDVT is INTEGER 00155 *> The leading dimension of the array VT. LDVT >= 1; if 00156 *> JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N; 00157 *> if JOBZ = 'S', LDVT >= min(M,N). 00158 *> \endverbatim 00159 *> 00160 *> \param[out] WORK 00161 *> \verbatim 00162 *> WORK is COMPLEX array, dimension (MAX(1,LWORK)) 00163 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00164 *> \endverbatim 00165 *> 00166 *> \param[in] LWORK 00167 *> \verbatim 00168 *> LWORK is INTEGER 00169 *> The dimension of the array WORK. LWORK >= 1. 00170 *> if JOBZ = 'N', LWORK >= 2*min(M,N)+max(M,N). 00171 *> if JOBZ = 'O', 00172 *> LWORK >= 2*min(M,N)*min(M,N)+2*min(M,N)+max(M,N). 00173 *> if JOBZ = 'S' or 'A', 00174 *> LWORK >= min(M,N)*min(M,N)+2*min(M,N)+max(M,N). 00175 *> For good performance, LWORK should generally be larger. 00176 *> 00177 *> If LWORK = -1, a workspace query is assumed. The optimal 00178 *> size for the WORK array is calculated and stored in WORK(1), 00179 *> and no other work except argument checking is performed. 00180 *> \endverbatim 00181 *> 00182 *> \param[out] RWORK 00183 *> \verbatim 00184 *> RWORK is REAL array, dimension (MAX(1,LRWORK)) 00185 *> If JOBZ = 'N', LRWORK >= 5*min(M,N). 00186 *> Otherwise, 00187 *> LRWORK >= min(M,N)*max(5*min(M,N)+7,2*max(M,N)+2*min(M,N)+1) 00188 *> \endverbatim 00189 *> 00190 *> \param[out] IWORK 00191 *> \verbatim 00192 *> IWORK is INTEGER array, dimension (8*min(M,N)) 00193 *> \endverbatim 00194 *> 00195 *> \param[out] INFO 00196 *> \verbatim 00197 *> INFO is INTEGER 00198 *> = 0: successful exit. 00199 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00200 *> > 0: The updating process of SBDSDC did not converge. 00201 *> \endverbatim 00202 * 00203 * Authors: 00204 * ======== 00205 * 00206 *> \author Univ. of Tennessee 00207 *> \author Univ. of California Berkeley 00208 *> \author Univ. of Colorado Denver 00209 *> \author NAG Ltd. 00210 * 00211 *> \date November 2011 00212 * 00213 *> \ingroup complexGEsing 00214 * 00215 *> \par Contributors: 00216 * ================== 00217 *> 00218 *> Ming Gu and Huan Ren, Computer Science Division, University of 00219 *> California at Berkeley, USA 00220 *> 00221 * ===================================================================== 00222 SUBROUTINE CGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, 00223 $ WORK, LWORK, RWORK, IWORK, INFO ) 00224 * 00225 * -- LAPACK driver routine (version 3.4.0) -- 00226 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00227 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00228 * November 2011 00229 * 00230 * .. Scalar Arguments .. 00231 CHARACTER JOBZ 00232 INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N 00233 * .. 00234 * .. Array Arguments .. 00235 INTEGER IWORK( * ) 00236 REAL RWORK( * ), S( * ) 00237 COMPLEX A( LDA, * ), U( LDU, * ), VT( LDVT, * ), 00238 $ WORK( * ) 00239 * .. 00240 * 00241 * ===================================================================== 00242 * 00243 * .. Parameters .. 00244 INTEGER LQUERV 00245 PARAMETER ( LQUERV = -1 ) 00246 COMPLEX CZERO, CONE 00247 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 00248 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 00249 REAL ZERO, ONE 00250 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00251 * .. 00252 * .. Local Scalars .. 00253 LOGICAL WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS 00254 INTEGER BLK, CHUNK, I, IE, IERR, IL, IR, IRU, IRVT, 00255 $ ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT, 00256 $ LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK, 00257 $ MNTHR1, MNTHR2, NRWORK, NWORK, WRKBL 00258 REAL ANRM, BIGNUM, EPS, SMLNUM 00259 * .. 00260 * .. Local Arrays .. 00261 INTEGER IDUM( 1 ) 00262 REAL DUM( 1 ) 00263 * .. 00264 * .. External Subroutines .. 00265 EXTERNAL CGEBRD, CGELQF, CGEMM, CGEQRF, CLACP2, CLACPY, 00266 $ CLACRM, CLARCM, CLASCL, CLASET, CUNGBR, CUNGLQ, 00267 $ CUNGQR, CUNMBR, SBDSDC, SLASCL, XERBLA 00268 * .. 00269 * .. External Functions .. 00270 LOGICAL LSAME 00271 INTEGER ILAENV 00272 REAL CLANGE, SLAMCH 00273 EXTERNAL CLANGE, SLAMCH, ILAENV, LSAME 00274 * .. 00275 * .. Intrinsic Functions .. 00276 INTRINSIC INT, MAX, MIN, SQRT 00277 * .. 00278 * .. Executable Statements .. 00279 * 00280 * Test the input arguments 00281 * 00282 INFO = 0 00283 MINMN = MIN( M, N ) 00284 MNTHR1 = INT( MINMN*17.0 / 9.0 ) 00285 MNTHR2 = INT( MINMN*5.0 / 3.0 ) 00286 WNTQA = LSAME( JOBZ, 'A' ) 00287 WNTQS = LSAME( JOBZ, 'S' ) 00288 WNTQAS = WNTQA .OR. WNTQS 00289 WNTQO = LSAME( JOBZ, 'O' ) 00290 WNTQN = LSAME( JOBZ, 'N' ) 00291 MINWRK = 1 00292 MAXWRK = 1 00293 * 00294 IF( .NOT.( WNTQA .OR. WNTQS .OR. WNTQO .OR. WNTQN ) ) THEN 00295 INFO = -1 00296 ELSE IF( M.LT.0 ) THEN 00297 INFO = -2 00298 ELSE IF( N.LT.0 ) THEN 00299 INFO = -3 00300 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00301 INFO = -5 00302 ELSE IF( LDU.LT.1 .OR. ( WNTQAS .AND. LDU.LT.M ) .OR. 00303 $ ( WNTQO .AND. M.LT.N .AND. LDU.LT.M ) ) THEN 00304 INFO = -8 00305 ELSE IF( LDVT.LT.1 .OR. ( WNTQA .AND. LDVT.LT.N ) .OR. 00306 $ ( WNTQS .AND. LDVT.LT.MINMN ) .OR. 00307 $ ( WNTQO .AND. M.GE.N .AND. LDVT.LT.N ) ) THEN 00308 INFO = -10 00309 END IF 00310 * 00311 * Compute workspace 00312 * (Note: Comments in the code beginning "Workspace:" describe the 00313 * minimal amount of workspace needed at that point in the code, 00314 * as well as the preferred amount for good performance. 00315 * CWorkspace refers to complex workspace, and RWorkspace to 00316 * real workspace. NB refers to the optimal block size for the 00317 * immediately following subroutine, as returned by ILAENV.) 00318 * 00319 IF( INFO.EQ.0 .AND. M.GT.0 .AND. N.GT.0 ) THEN 00320 IF( M.GE.N ) THEN 00321 * 00322 * There is no complex work space needed for bidiagonal SVD 00323 * The real work space needed for bidiagonal SVD is BDSPAC 00324 * for computing singular values and singular vectors; BDSPAN 00325 * for computing singular values only. 00326 * BDSPAC = 5*N*N + 7*N 00327 * BDSPAN = MAX(7*N+4, 3*N+2+SMLSIZ*(SMLSIZ+8)) 00328 * 00329 IF( M.GE.MNTHR1 ) THEN 00330 IF( WNTQN ) THEN 00331 * 00332 * Path 1 (M much larger than N, JOBZ='N') 00333 * 00334 MAXWRK = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, 00335 $ -1 ) 00336 MAXWRK = MAX( MAXWRK, 2*N+2*N* 00337 $ ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) ) 00338 MINWRK = 3*N 00339 ELSE IF( WNTQO ) THEN 00340 * 00341 * Path 2 (M much larger than N, JOBZ='O') 00342 * 00343 WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 ) 00344 WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'CUNGQR', ' ', M, 00345 $ N, N, -1 ) ) 00346 WRKBL = MAX( WRKBL, 2*N+2*N* 00347 $ ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) ) 00348 WRKBL = MAX( WRKBL, 2*N+N* 00349 $ ILAENV( 1, 'CUNMBR', 'QLN', N, N, N, -1 ) ) 00350 WRKBL = MAX( WRKBL, 2*N+N* 00351 $ ILAENV( 1, 'CUNMBR', 'PRC', N, N, N, -1 ) ) 00352 MAXWRK = M*N + N*N + WRKBL 00353 MINWRK = 2*N*N + 3*N 00354 ELSE IF( WNTQS ) THEN 00355 * 00356 * Path 3 (M much larger than N, JOBZ='S') 00357 * 00358 WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 ) 00359 WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'CUNGQR', ' ', M, 00360 $ N, N, -1 ) ) 00361 WRKBL = MAX( WRKBL, 2*N+2*N* 00362 $ ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) ) 00363 WRKBL = MAX( WRKBL, 2*N+N* 00364 $ ILAENV( 1, 'CUNMBR', 'QLN', N, N, N, -1 ) ) 00365 WRKBL = MAX( WRKBL, 2*N+N* 00366 $ ILAENV( 1, 'CUNMBR', 'PRC', N, N, N, -1 ) ) 00367 MAXWRK = N*N + WRKBL 00368 MINWRK = N*N + 3*N 00369 ELSE IF( WNTQA ) THEN 00370 * 00371 * Path 4 (M much larger than N, JOBZ='A') 00372 * 00373 WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 ) 00374 WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'CUNGQR', ' ', M, 00375 $ M, N, -1 ) ) 00376 WRKBL = MAX( WRKBL, 2*N+2*N* 00377 $ ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) ) 00378 WRKBL = MAX( WRKBL, 2*N+N* 00379 $ ILAENV( 1, 'CUNMBR', 'QLN', N, N, N, -1 ) ) 00380 WRKBL = MAX( WRKBL, 2*N+N* 00381 $ ILAENV( 1, 'CUNMBR', 'PRC', N, N, N, -1 ) ) 00382 MAXWRK = N*N + WRKBL 00383 MINWRK = N*N + 2*N + M 00384 END IF 00385 ELSE IF( M.GE.MNTHR2 ) THEN 00386 * 00387 * Path 5 (M much larger than N, but not as much as MNTHR1) 00388 * 00389 MAXWRK = 2*N + ( M+N )*ILAENV( 1, 'CGEBRD', ' ', M, N, 00390 $ -1, -1 ) 00391 MINWRK = 2*N + M 00392 IF( WNTQO ) THEN 00393 MAXWRK = MAX( MAXWRK, 2*N+N* 00394 $ ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) ) 00395 MAXWRK = MAX( MAXWRK, 2*N+N* 00396 $ ILAENV( 1, 'CUNGBR', 'Q', M, N, N, -1 ) ) 00397 MAXWRK = MAXWRK + M*N 00398 MINWRK = MINWRK + N*N 00399 ELSE IF( WNTQS ) THEN 00400 MAXWRK = MAX( MAXWRK, 2*N+N* 00401 $ ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) ) 00402 MAXWRK = MAX( MAXWRK, 2*N+N* 00403 $ ILAENV( 1, 'CUNGBR', 'Q', M, N, N, -1 ) ) 00404 ELSE IF( WNTQA ) THEN 00405 MAXWRK = MAX( MAXWRK, 2*N+N* 00406 $ ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) ) 00407 MAXWRK = MAX( MAXWRK, 2*N+M* 00408 $ ILAENV( 1, 'CUNGBR', 'Q', M, M, N, -1 ) ) 00409 END IF 00410 ELSE 00411 * 00412 * Path 6 (M at least N, but not much larger) 00413 * 00414 MAXWRK = 2*N + ( M+N )*ILAENV( 1, 'CGEBRD', ' ', M, N, 00415 $ -1, -1 ) 00416 MINWRK = 2*N + M 00417 IF( WNTQO ) THEN 00418 MAXWRK = MAX( MAXWRK, 2*N+N* 00419 $ ILAENV( 1, 'CUNMBR', 'PRC', N, N, N, -1 ) ) 00420 MAXWRK = MAX( MAXWRK, 2*N+N* 00421 $ ILAENV( 1, 'CUNMBR', 'QLN', M, N, N, -1 ) ) 00422 MAXWRK = MAXWRK + M*N 00423 MINWRK = MINWRK + N*N 00424 ELSE IF( WNTQS ) THEN 00425 MAXWRK = MAX( MAXWRK, 2*N+N* 00426 $ ILAENV( 1, 'CUNMBR', 'PRC', N, N, N, -1 ) ) 00427 MAXWRK = MAX( MAXWRK, 2*N+N* 00428 $ ILAENV( 1, 'CUNMBR', 'QLN', M, N, N, -1 ) ) 00429 ELSE IF( WNTQA ) THEN 00430 MAXWRK = MAX( MAXWRK, 2*N+N* 00431 $ ILAENV( 1, 'CUNGBR', 'PRC', N, N, N, -1 ) ) 00432 MAXWRK = MAX( MAXWRK, 2*N+M* 00433 $ ILAENV( 1, 'CUNGBR', 'QLN', M, M, N, -1 ) ) 00434 END IF 00435 END IF 00436 ELSE 00437 * 00438 * There is no complex work space needed for bidiagonal SVD 00439 * The real work space needed for bidiagonal SVD is BDSPAC 00440 * for computing singular values and singular vectors; BDSPAN 00441 * for computing singular values only. 00442 * BDSPAC = 5*M*M + 7*M 00443 * BDSPAN = MAX(7*M+4, 3*M+2+SMLSIZ*(SMLSIZ+8)) 00444 * 00445 IF( N.GE.MNTHR1 ) THEN 00446 IF( WNTQN ) THEN 00447 * 00448 * Path 1t (N much larger than M, JOBZ='N') 00449 * 00450 MAXWRK = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, 00451 $ -1 ) 00452 MAXWRK = MAX( MAXWRK, 2*M+2*M* 00453 $ ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) ) 00454 MINWRK = 3*M 00455 ELSE IF( WNTQO ) THEN 00456 * 00457 * Path 2t (N much larger than M, JOBZ='O') 00458 * 00459 WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 ) 00460 WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'CUNGLQ', ' ', M, 00461 $ N, M, -1 ) ) 00462 WRKBL = MAX( WRKBL, 2*M+2*M* 00463 $ ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) ) 00464 WRKBL = MAX( WRKBL, 2*M+M* 00465 $ ILAENV( 1, 'CUNMBR', 'PRC', M, M, M, -1 ) ) 00466 WRKBL = MAX( WRKBL, 2*M+M* 00467 $ ILAENV( 1, 'CUNMBR', 'QLN', M, M, M, -1 ) ) 00468 MAXWRK = M*N + M*M + WRKBL 00469 MINWRK = 2*M*M + 3*M 00470 ELSE IF( WNTQS ) THEN 00471 * 00472 * Path 3t (N much larger than M, JOBZ='S') 00473 * 00474 WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 ) 00475 WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'CUNGLQ', ' ', M, 00476 $ N, M, -1 ) ) 00477 WRKBL = MAX( WRKBL, 2*M+2*M* 00478 $ ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) ) 00479 WRKBL = MAX( WRKBL, 2*M+M* 00480 $ ILAENV( 1, 'CUNMBR', 'PRC', M, M, M, -1 ) ) 00481 WRKBL = MAX( WRKBL, 2*M+M* 00482 $ ILAENV( 1, 'CUNMBR', 'QLN', M, M, M, -1 ) ) 00483 MAXWRK = M*M + WRKBL 00484 MINWRK = M*M + 3*M 00485 ELSE IF( WNTQA ) THEN 00486 * 00487 * Path 4t (N much larger than M, JOBZ='A') 00488 * 00489 WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 ) 00490 WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'CUNGLQ', ' ', N, 00491 $ N, M, -1 ) ) 00492 WRKBL = MAX( WRKBL, 2*M+2*M* 00493 $ ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) ) 00494 WRKBL = MAX( WRKBL, 2*M+M* 00495 $ ILAENV( 1, 'CUNMBR', 'PRC', M, M, M, -1 ) ) 00496 WRKBL = MAX( WRKBL, 2*M+M* 00497 $ ILAENV( 1, 'CUNMBR', 'QLN', M, M, M, -1 ) ) 00498 MAXWRK = M*M + WRKBL 00499 MINWRK = M*M + 2*M + N 00500 END IF 00501 ELSE IF( N.GE.MNTHR2 ) THEN 00502 * 00503 * Path 5t (N much larger than M, but not as much as MNTHR1) 00504 * 00505 MAXWRK = 2*M + ( M+N )*ILAENV( 1, 'CGEBRD', ' ', M, N, 00506 $ -1, -1 ) 00507 MINWRK = 2*M + N 00508 IF( WNTQO ) THEN 00509 MAXWRK = MAX( MAXWRK, 2*M+M* 00510 $ ILAENV( 1, 'CUNGBR', 'P', M, N, M, -1 ) ) 00511 MAXWRK = MAX( MAXWRK, 2*M+M* 00512 $ ILAENV( 1, 'CUNGBR', 'Q', M, M, N, -1 ) ) 00513 MAXWRK = MAXWRK + M*N 00514 MINWRK = MINWRK + M*M 00515 ELSE IF( WNTQS ) THEN 00516 MAXWRK = MAX( MAXWRK, 2*M+M* 00517 $ ILAENV( 1, 'CUNGBR', 'P', M, N, M, -1 ) ) 00518 MAXWRK = MAX( MAXWRK, 2*M+M* 00519 $ ILAENV( 1, 'CUNGBR', 'Q', M, M, N, -1 ) ) 00520 ELSE IF( WNTQA ) THEN 00521 MAXWRK = MAX( MAXWRK, 2*M+N* 00522 $ ILAENV( 1, 'CUNGBR', 'P', N, N, M, -1 ) ) 00523 MAXWRK = MAX( MAXWRK, 2*M+M* 00524 $ ILAENV( 1, 'CUNGBR', 'Q', M, M, N, -1 ) ) 00525 END IF 00526 ELSE 00527 * 00528 * Path 6t (N greater than M, but not much larger) 00529 * 00530 MAXWRK = 2*M + ( M+N )*ILAENV( 1, 'CGEBRD', ' ', M, N, 00531 $ -1, -1 ) 00532 MINWRK = 2*M + N 00533 IF( WNTQO ) THEN 00534 MAXWRK = MAX( MAXWRK, 2*M+M* 00535 $ ILAENV( 1, 'CUNMBR', 'PRC', M, N, M, -1 ) ) 00536 MAXWRK = MAX( MAXWRK, 2*M+M* 00537 $ ILAENV( 1, 'CUNMBR', 'QLN', M, M, N, -1 ) ) 00538 MAXWRK = MAXWRK + M*N 00539 MINWRK = MINWRK + M*M 00540 ELSE IF( WNTQS ) THEN 00541 MAXWRK = MAX( MAXWRK, 2*M+M* 00542 $ ILAENV( 1, 'CUNGBR', 'PRC', M, N, M, -1 ) ) 00543 MAXWRK = MAX( MAXWRK, 2*M+M* 00544 $ ILAENV( 1, 'CUNGBR', 'QLN', M, M, N, -1 ) ) 00545 ELSE IF( WNTQA ) THEN 00546 MAXWRK = MAX( MAXWRK, 2*M+N* 00547 $ ILAENV( 1, 'CUNGBR', 'PRC', N, N, M, -1 ) ) 00548 MAXWRK = MAX( MAXWRK, 2*M+M* 00549 $ ILAENV( 1, 'CUNGBR', 'QLN', M, M, N, -1 ) ) 00550 END IF 00551 END IF 00552 END IF 00553 MAXWRK = MAX( MAXWRK, MINWRK ) 00554 END IF 00555 IF( INFO.EQ.0 ) THEN 00556 WORK( 1 ) = MAXWRK 00557 IF( LWORK.LT.MINWRK .AND. LWORK.NE.LQUERV ) 00558 $ INFO = -13 00559 END IF 00560 * 00561 * Quick returns 00562 * 00563 IF( INFO.NE.0 ) THEN 00564 CALL XERBLA( 'CGESDD', -INFO ) 00565 RETURN 00566 END IF 00567 IF( LWORK.EQ.LQUERV ) 00568 $ RETURN 00569 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 00570 RETURN 00571 END IF 00572 * 00573 * Get machine constants 00574 * 00575 EPS = SLAMCH( 'P' ) 00576 SMLNUM = SQRT( SLAMCH( 'S' ) ) / EPS 00577 BIGNUM = ONE / SMLNUM 00578 * 00579 * Scale A if max element outside range [SMLNUM,BIGNUM] 00580 * 00581 ANRM = CLANGE( 'M', M, N, A, LDA, DUM ) 00582 ISCL = 0 00583 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00584 ISCL = 1 00585 CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR ) 00586 ELSE IF( ANRM.GT.BIGNUM ) THEN 00587 ISCL = 1 00588 CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR ) 00589 END IF 00590 * 00591 IF( M.GE.N ) THEN 00592 * 00593 * A has at least as many rows as columns. If A has sufficiently 00594 * more rows than columns, first reduce using the QR 00595 * decomposition (if sufficient workspace available) 00596 * 00597 IF( M.GE.MNTHR1 ) THEN 00598 * 00599 IF( WNTQN ) THEN 00600 * 00601 * Path 1 (M much larger than N, JOBZ='N') 00602 * No singular vectors to be computed 00603 * 00604 ITAU = 1 00605 NWORK = ITAU + N 00606 * 00607 * Compute A=Q*R 00608 * (CWorkspace: need 2*N, prefer N+N*NB) 00609 * (RWorkspace: need 0) 00610 * 00611 CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 00612 $ LWORK-NWORK+1, IERR ) 00613 * 00614 * Zero out below R 00615 * 00616 CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ), 00617 $ LDA ) 00618 IE = 1 00619 ITAUQ = 1 00620 ITAUP = ITAUQ + N 00621 NWORK = ITAUP + N 00622 * 00623 * Bidiagonalize R in A 00624 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 00625 * (RWorkspace: need N) 00626 * 00627 CALL CGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 00628 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 00629 $ IERR ) 00630 NRWORK = IE + N 00631 * 00632 * Perform bidiagonal SVD, compute singular values only 00633 * (CWorkspace: 0) 00634 * (RWorkspace: need BDSPAN) 00635 * 00636 CALL SBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1, 00637 $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) 00638 * 00639 ELSE IF( WNTQO ) THEN 00640 * 00641 * Path 2 (M much larger than N, JOBZ='O') 00642 * N left singular vectors to be overwritten on A and 00643 * N right singular vectors to be computed in VT 00644 * 00645 IU = 1 00646 * 00647 * WORK(IU) is N by N 00648 * 00649 LDWRKU = N 00650 IR = IU + LDWRKU*N 00651 IF( LWORK.GE.M*N+N*N+3*N ) THEN 00652 * 00653 * WORK(IR) is M by N 00654 * 00655 LDWRKR = M 00656 ELSE 00657 LDWRKR = ( LWORK-N*N-3*N ) / N 00658 END IF 00659 ITAU = IR + LDWRKR*N 00660 NWORK = ITAU + N 00661 * 00662 * Compute A=Q*R 00663 * (CWorkspace: need N*N+2*N, prefer M*N+N+N*NB) 00664 * (RWorkspace: 0) 00665 * 00666 CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 00667 $ LWORK-NWORK+1, IERR ) 00668 * 00669 * Copy R to WORK( IR ), zeroing out below it 00670 * 00671 CALL CLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR ) 00672 CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, WORK( IR+1 ), 00673 $ LDWRKR ) 00674 * 00675 * Generate Q in A 00676 * (CWorkspace: need 2*N, prefer N+N*NB) 00677 * (RWorkspace: 0) 00678 * 00679 CALL CUNGQR( M, N, N, A, LDA, WORK( ITAU ), 00680 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 00681 IE = 1 00682 ITAUQ = ITAU 00683 ITAUP = ITAUQ + N 00684 NWORK = ITAUP + N 00685 * 00686 * Bidiagonalize R in WORK(IR) 00687 * (CWorkspace: need N*N+3*N, prefer M*N+2*N+2*N*NB) 00688 * (RWorkspace: need N) 00689 * 00690 CALL CGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ), 00691 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), 00692 $ LWORK-NWORK+1, IERR ) 00693 * 00694 * Perform bidiagonal SVD, computing left singular vectors 00695 * of R in WORK(IRU) and computing right singular vectors 00696 * of R in WORK(IRVT) 00697 * (CWorkspace: need 0) 00698 * (RWorkspace: need BDSPAC) 00699 * 00700 IRU = IE + N 00701 IRVT = IRU + N*N 00702 NRWORK = IRVT + N*N 00703 CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), 00704 $ N, RWORK( IRVT ), N, DUM, IDUM, 00705 $ RWORK( NRWORK ), IWORK, INFO ) 00706 * 00707 * Copy real matrix RWORK(IRU) to complex matrix WORK(IU) 00708 * Overwrite WORK(IU) by the left singular vectors of R 00709 * (CWorkspace: need 2*N*N+3*N, prefer M*N+N*N+2*N+N*NB) 00710 * (RWorkspace: 0) 00711 * 00712 CALL CLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ), 00713 $ LDWRKU ) 00714 CALL CUNMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR, 00715 $ WORK( ITAUQ ), WORK( IU ), LDWRKU, 00716 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 00717 * 00718 * Copy real matrix RWORK(IRVT) to complex matrix VT 00719 * Overwrite VT by the right singular vectors of R 00720 * (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB) 00721 * (RWorkspace: 0) 00722 * 00723 CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) 00724 CALL CUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR, 00725 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 00726 $ LWORK-NWORK+1, IERR ) 00727 * 00728 * Multiply Q in A by left singular vectors of R in 00729 * WORK(IU), storing result in WORK(IR) and copying to A 00730 * (CWorkspace: need 2*N*N, prefer N*N+M*N) 00731 * (RWorkspace: 0) 00732 * 00733 DO 10 I = 1, M, LDWRKR 00734 CHUNK = MIN( M-I+1, LDWRKR ) 00735 CALL CGEMM( 'N', 'N', CHUNK, N, N, CONE, A( I, 1 ), 00736 $ LDA, WORK( IU ), LDWRKU, CZERO, 00737 $ WORK( IR ), LDWRKR ) 00738 CALL CLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR, 00739 $ A( I, 1 ), LDA ) 00740 10 CONTINUE 00741 * 00742 ELSE IF( WNTQS ) THEN 00743 * 00744 * Path 3 (M much larger than N, JOBZ='S') 00745 * N left singular vectors to be computed in U and 00746 * N right singular vectors to be computed in VT 00747 * 00748 IR = 1 00749 * 00750 * WORK(IR) is N by N 00751 * 00752 LDWRKR = N 00753 ITAU = IR + LDWRKR*N 00754 NWORK = ITAU + N 00755 * 00756 * Compute A=Q*R 00757 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 00758 * (RWorkspace: 0) 00759 * 00760 CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 00761 $ LWORK-NWORK+1, IERR ) 00762 * 00763 * Copy R to WORK(IR), zeroing out below it 00764 * 00765 CALL CLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR ) 00766 CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, WORK( IR+1 ), 00767 $ LDWRKR ) 00768 * 00769 * Generate Q in A 00770 * (CWorkspace: need 2*N, prefer N+N*NB) 00771 * (RWorkspace: 0) 00772 * 00773 CALL CUNGQR( M, N, N, A, LDA, WORK( ITAU ), 00774 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 00775 IE = 1 00776 ITAUQ = ITAU 00777 ITAUP = ITAUQ + N 00778 NWORK = ITAUP + N 00779 * 00780 * Bidiagonalize R in WORK(IR) 00781 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) 00782 * (RWorkspace: need N) 00783 * 00784 CALL CGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ), 00785 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), 00786 $ LWORK-NWORK+1, IERR ) 00787 * 00788 * Perform bidiagonal SVD, computing left singular vectors 00789 * of bidiagonal matrix in RWORK(IRU) and computing right 00790 * singular vectors of bidiagonal matrix in RWORK(IRVT) 00791 * (CWorkspace: need 0) 00792 * (RWorkspace: need BDSPAC) 00793 * 00794 IRU = IE + N 00795 IRVT = IRU + N*N 00796 NRWORK = IRVT + N*N 00797 CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), 00798 $ N, RWORK( IRVT ), N, DUM, IDUM, 00799 $ RWORK( NRWORK ), IWORK, INFO ) 00800 * 00801 * Copy real matrix RWORK(IRU) to complex matrix U 00802 * Overwrite U by left singular vectors of R 00803 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) 00804 * (RWorkspace: 0) 00805 * 00806 CALL CLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU ) 00807 CALL CUNMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR, 00808 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 00809 $ LWORK-NWORK+1, IERR ) 00810 * 00811 * Copy real matrix RWORK(IRVT) to complex matrix VT 00812 * Overwrite VT by right singular vectors of R 00813 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) 00814 * (RWorkspace: 0) 00815 * 00816 CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) 00817 CALL CUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR, 00818 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 00819 $ LWORK-NWORK+1, IERR ) 00820 * 00821 * Multiply Q in A by left singular vectors of R in 00822 * WORK(IR), storing result in U 00823 * (CWorkspace: need N*N) 00824 * (RWorkspace: 0) 00825 * 00826 CALL CLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR ) 00827 CALL CGEMM( 'N', 'N', M, N, N, CONE, A, LDA, WORK( IR ), 00828 $ LDWRKR, CZERO, U, LDU ) 00829 * 00830 ELSE IF( WNTQA ) THEN 00831 * 00832 * Path 4 (M much larger than N, JOBZ='A') 00833 * M left singular vectors to be computed in U and 00834 * N right singular vectors to be computed in VT 00835 * 00836 IU = 1 00837 * 00838 * WORK(IU) is N by N 00839 * 00840 LDWRKU = N 00841 ITAU = IU + LDWRKU*N 00842 NWORK = ITAU + N 00843 * 00844 * Compute A=Q*R, copying result to U 00845 * (CWorkspace: need 2*N, prefer N+N*NB) 00846 * (RWorkspace: 0) 00847 * 00848 CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 00849 $ LWORK-NWORK+1, IERR ) 00850 CALL CLACPY( 'L', M, N, A, LDA, U, LDU ) 00851 * 00852 * Generate Q in U 00853 * (CWorkspace: need N+M, prefer N+M*NB) 00854 * (RWorkspace: 0) 00855 * 00856 CALL CUNGQR( M, M, N, U, LDU, WORK( ITAU ), 00857 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 00858 * 00859 * Produce R in A, zeroing out below it 00860 * 00861 CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ), 00862 $ LDA ) 00863 IE = 1 00864 ITAUQ = ITAU 00865 ITAUP = ITAUQ + N 00866 NWORK = ITAUP + N 00867 * 00868 * Bidiagonalize R in A 00869 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 00870 * (RWorkspace: need N) 00871 * 00872 CALL CGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 00873 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 00874 $ IERR ) 00875 IRU = IE + N 00876 IRVT = IRU + N*N 00877 NRWORK = IRVT + N*N 00878 * 00879 * Perform bidiagonal SVD, computing left singular vectors 00880 * of bidiagonal matrix in RWORK(IRU) and computing right 00881 * singular vectors of bidiagonal matrix in RWORK(IRVT) 00882 * (CWorkspace: need 0) 00883 * (RWorkspace: need BDSPAC) 00884 * 00885 CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), 00886 $ N, RWORK( IRVT ), N, DUM, IDUM, 00887 $ RWORK( NRWORK ), IWORK, INFO ) 00888 * 00889 * Copy real matrix RWORK(IRU) to complex matrix WORK(IU) 00890 * Overwrite WORK(IU) by left singular vectors of R 00891 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) 00892 * (RWorkspace: 0) 00893 * 00894 CALL CLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ), 00895 $ LDWRKU ) 00896 CALL CUNMBR( 'Q', 'L', 'N', N, N, N, A, LDA, 00897 $ WORK( ITAUQ ), WORK( IU ), LDWRKU, 00898 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 00899 * 00900 * Copy real matrix RWORK(IRVT) to complex matrix VT 00901 * Overwrite VT by right singular vectors of R 00902 * (CWorkspace: need 3*N, prefer 2*N+N*NB) 00903 * (RWorkspace: 0) 00904 * 00905 CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) 00906 CALL CUNMBR( 'P', 'R', 'C', N, N, N, A, LDA, 00907 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 00908 $ LWORK-NWORK+1, IERR ) 00909 * 00910 * Multiply Q in U by left singular vectors of R in 00911 * WORK(IU), storing result in A 00912 * (CWorkspace: need N*N) 00913 * (RWorkspace: 0) 00914 * 00915 CALL CGEMM( 'N', 'N', M, N, N, CONE, U, LDU, WORK( IU ), 00916 $ LDWRKU, CZERO, A, LDA ) 00917 * 00918 * Copy left singular vectors of A from A to U 00919 * 00920 CALL CLACPY( 'F', M, N, A, LDA, U, LDU ) 00921 * 00922 END IF 00923 * 00924 ELSE IF( M.GE.MNTHR2 ) THEN 00925 * 00926 * MNTHR2 <= M < MNTHR1 00927 * 00928 * Path 5 (M much larger than N, but not as much as MNTHR1) 00929 * Reduce to bidiagonal form without QR decomposition, use 00930 * CUNGBR and matrix multiplication to compute singular vectors 00931 * 00932 IE = 1 00933 NRWORK = IE + N 00934 ITAUQ = 1 00935 ITAUP = ITAUQ + N 00936 NWORK = ITAUP + N 00937 * 00938 * Bidiagonalize A 00939 * (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB) 00940 * (RWorkspace: need N) 00941 * 00942 CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 00943 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 00944 $ IERR ) 00945 IF( WNTQN ) THEN 00946 * 00947 * Compute singular values only 00948 * (Cworkspace: 0) 00949 * (Rworkspace: need BDSPAN) 00950 * 00951 CALL SBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1, 00952 $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) 00953 ELSE IF( WNTQO ) THEN 00954 IU = NWORK 00955 IRU = NRWORK 00956 IRVT = IRU + N*N 00957 NRWORK = IRVT + N*N 00958 * 00959 * Copy A to VT, generate P**H 00960 * (Cworkspace: need 2*N, prefer N+N*NB) 00961 * (Rworkspace: 0) 00962 * 00963 CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT ) 00964 CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 00965 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 00966 * 00967 * Generate Q in A 00968 * (CWorkspace: need 2*N, prefer N+N*NB) 00969 * (RWorkspace: 0) 00970 * 00971 CALL CUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ), 00972 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 00973 * 00974 IF( LWORK.GE.M*N+3*N ) THEN 00975 * 00976 * WORK( IU ) is M by N 00977 * 00978 LDWRKU = M 00979 ELSE 00980 * 00981 * WORK(IU) is LDWRKU by N 00982 * 00983 LDWRKU = ( LWORK-3*N ) / N 00984 END IF 00985 NWORK = IU + LDWRKU*N 00986 * 00987 * Perform bidiagonal SVD, computing left singular vectors 00988 * of bidiagonal matrix in RWORK(IRU) and computing right 00989 * singular vectors of bidiagonal matrix in RWORK(IRVT) 00990 * (CWorkspace: need 0) 00991 * (RWorkspace: need BDSPAC) 00992 * 00993 CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), 00994 $ N, RWORK( IRVT ), N, DUM, IDUM, 00995 $ RWORK( NRWORK ), IWORK, INFO ) 00996 * 00997 * Multiply real matrix RWORK(IRVT) by P**H in VT, 00998 * storing the result in WORK(IU), copying to VT 00999 * (Cworkspace: need 0) 01000 * (Rworkspace: need 3*N*N) 01001 * 01002 CALL CLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, 01003 $ WORK( IU ), LDWRKU, RWORK( NRWORK ) ) 01004 CALL CLACPY( 'F', N, N, WORK( IU ), LDWRKU, VT, LDVT ) 01005 * 01006 * Multiply Q in A by real matrix RWORK(IRU), storing the 01007 * result in WORK(IU), copying to A 01008 * (CWorkspace: need N*N, prefer M*N) 01009 * (Rworkspace: need 3*N*N, prefer N*N+2*M*N) 01010 * 01011 NRWORK = IRVT 01012 DO 20 I = 1, M, LDWRKU 01013 CHUNK = MIN( M-I+1, LDWRKU ) 01014 CALL CLACRM( CHUNK, N, A( I, 1 ), LDA, RWORK( IRU ), 01015 $ N, WORK( IU ), LDWRKU, RWORK( NRWORK ) ) 01016 CALL CLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU, 01017 $ A( I, 1 ), LDA ) 01018 20 CONTINUE 01019 * 01020 ELSE IF( WNTQS ) THEN 01021 * 01022 * Copy A to VT, generate P**H 01023 * (Cworkspace: need 2*N, prefer N+N*NB) 01024 * (Rworkspace: 0) 01025 * 01026 CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT ) 01027 CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 01028 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01029 * 01030 * Copy A to U, generate Q 01031 * (Cworkspace: need 2*N, prefer N+N*NB) 01032 * (Rworkspace: 0) 01033 * 01034 CALL CLACPY( 'L', M, N, A, LDA, U, LDU ) 01035 CALL CUNGBR( 'Q', M, N, N, U, LDU, WORK( ITAUQ ), 01036 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01037 * 01038 * Perform bidiagonal SVD, computing left singular vectors 01039 * of bidiagonal matrix in RWORK(IRU) and computing right 01040 * singular vectors of bidiagonal matrix in RWORK(IRVT) 01041 * (CWorkspace: need 0) 01042 * (RWorkspace: need BDSPAC) 01043 * 01044 IRU = NRWORK 01045 IRVT = IRU + N*N 01046 NRWORK = IRVT + N*N 01047 CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), 01048 $ N, RWORK( IRVT ), N, DUM, IDUM, 01049 $ RWORK( NRWORK ), IWORK, INFO ) 01050 * 01051 * Multiply real matrix RWORK(IRVT) by P**H in VT, 01052 * storing the result in A, copying to VT 01053 * (Cworkspace: need 0) 01054 * (Rworkspace: need 3*N*N) 01055 * 01056 CALL CLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA, 01057 $ RWORK( NRWORK ) ) 01058 CALL CLACPY( 'F', N, N, A, LDA, VT, LDVT ) 01059 * 01060 * Multiply Q in U by real matrix RWORK(IRU), storing the 01061 * result in A, copying to U 01062 * (CWorkspace: need 0) 01063 * (Rworkspace: need N*N+2*M*N) 01064 * 01065 NRWORK = IRVT 01066 CALL CLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA, 01067 $ RWORK( NRWORK ) ) 01068 CALL CLACPY( 'F', M, N, A, LDA, U, LDU ) 01069 ELSE 01070 * 01071 * Copy A to VT, generate P**H 01072 * (Cworkspace: need 2*N, prefer N+N*NB) 01073 * (Rworkspace: 0) 01074 * 01075 CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT ) 01076 CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 01077 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01078 * 01079 * Copy A to U, generate Q 01080 * (Cworkspace: need 2*N, prefer N+N*NB) 01081 * (Rworkspace: 0) 01082 * 01083 CALL CLACPY( 'L', M, N, A, LDA, U, LDU ) 01084 CALL CUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ), 01085 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01086 * 01087 * Perform bidiagonal SVD, computing left singular vectors 01088 * of bidiagonal matrix in RWORK(IRU) and computing right 01089 * singular vectors of bidiagonal matrix in RWORK(IRVT) 01090 * (CWorkspace: need 0) 01091 * (RWorkspace: need BDSPAC) 01092 * 01093 IRU = NRWORK 01094 IRVT = IRU + N*N 01095 NRWORK = IRVT + N*N 01096 CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), 01097 $ N, RWORK( IRVT ), N, DUM, IDUM, 01098 $ RWORK( NRWORK ), IWORK, INFO ) 01099 * 01100 * Multiply real matrix RWORK(IRVT) by P**H in VT, 01101 * storing the result in A, copying to VT 01102 * (Cworkspace: need 0) 01103 * (Rworkspace: need 3*N*N) 01104 * 01105 CALL CLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA, 01106 $ RWORK( NRWORK ) ) 01107 CALL CLACPY( 'F', N, N, A, LDA, VT, LDVT ) 01108 * 01109 * Multiply Q in U by real matrix RWORK(IRU), storing the 01110 * result in A, copying to U 01111 * (CWorkspace: 0) 01112 * (Rworkspace: need 3*N*N) 01113 * 01114 NRWORK = IRVT 01115 CALL CLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA, 01116 $ RWORK( NRWORK ) ) 01117 CALL CLACPY( 'F', M, N, A, LDA, U, LDU ) 01118 END IF 01119 * 01120 ELSE 01121 * 01122 * M .LT. MNTHR2 01123 * 01124 * Path 6 (M at least N, but not much larger) 01125 * Reduce to bidiagonal form without QR decomposition 01126 * Use CUNMBR to compute singular vectors 01127 * 01128 IE = 1 01129 NRWORK = IE + N 01130 ITAUQ = 1 01131 ITAUP = ITAUQ + N 01132 NWORK = ITAUP + N 01133 * 01134 * Bidiagonalize A 01135 * (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB) 01136 * (RWorkspace: need N) 01137 * 01138 CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 01139 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 01140 $ IERR ) 01141 IF( WNTQN ) THEN 01142 * 01143 * Compute singular values only 01144 * (Cworkspace: 0) 01145 * (Rworkspace: need BDSPAN) 01146 * 01147 CALL SBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1, 01148 $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) 01149 ELSE IF( WNTQO ) THEN 01150 IU = NWORK 01151 IRU = NRWORK 01152 IRVT = IRU + N*N 01153 NRWORK = IRVT + N*N 01154 IF( LWORK.GE.M*N+3*N ) THEN 01155 * 01156 * WORK( IU ) is M by N 01157 * 01158 LDWRKU = M 01159 ELSE 01160 * 01161 * WORK( IU ) is LDWRKU by N 01162 * 01163 LDWRKU = ( LWORK-3*N ) / N 01164 END IF 01165 NWORK = IU + LDWRKU*N 01166 * 01167 * Perform bidiagonal SVD, computing left singular vectors 01168 * of bidiagonal matrix in RWORK(IRU) and computing right 01169 * singular vectors of bidiagonal matrix in RWORK(IRVT) 01170 * (CWorkspace: need 0) 01171 * (RWorkspace: need BDSPAC) 01172 * 01173 CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), 01174 $ N, RWORK( IRVT ), N, DUM, IDUM, 01175 $ RWORK( NRWORK ), IWORK, INFO ) 01176 * 01177 * Copy real matrix RWORK(IRVT) to complex matrix VT 01178 * Overwrite VT by right singular vectors of A 01179 * (Cworkspace: need 2*N, prefer N+N*NB) 01180 * (Rworkspace: need 0) 01181 * 01182 CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) 01183 CALL CUNMBR( 'P', 'R', 'C', N, N, N, A, LDA, 01184 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 01185 $ LWORK-NWORK+1, IERR ) 01186 * 01187 IF( LWORK.GE.M*N+3*N ) THEN 01188 * 01189 * Copy real matrix RWORK(IRU) to complex matrix WORK(IU) 01190 * Overwrite WORK(IU) by left singular vectors of A, copying 01191 * to A 01192 * (Cworkspace: need M*N+2*N, prefer M*N+N+N*NB) 01193 * (Rworkspace: need 0) 01194 * 01195 CALL CLASET( 'F', M, N, CZERO, CZERO, WORK( IU ), 01196 $ LDWRKU ) 01197 CALL CLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ), 01198 $ LDWRKU ) 01199 CALL CUNMBR( 'Q', 'L', 'N', M, N, N, A, LDA, 01200 $ WORK( ITAUQ ), WORK( IU ), LDWRKU, 01201 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01202 CALL CLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA ) 01203 ELSE 01204 * 01205 * Generate Q in A 01206 * (Cworkspace: need 2*N, prefer N+N*NB) 01207 * (Rworkspace: need 0) 01208 * 01209 CALL CUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ), 01210 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01211 * 01212 * Multiply Q in A by real matrix RWORK(IRU), storing the 01213 * result in WORK(IU), copying to A 01214 * (CWorkspace: need N*N, prefer M*N) 01215 * (Rworkspace: need 3*N*N, prefer N*N+2*M*N) 01216 * 01217 NRWORK = IRVT 01218 DO 30 I = 1, M, LDWRKU 01219 CHUNK = MIN( M-I+1, LDWRKU ) 01220 CALL CLACRM( CHUNK, N, A( I, 1 ), LDA, 01221 $ RWORK( IRU ), N, WORK( IU ), LDWRKU, 01222 $ RWORK( NRWORK ) ) 01223 CALL CLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU, 01224 $ A( I, 1 ), LDA ) 01225 30 CONTINUE 01226 END IF 01227 * 01228 ELSE IF( WNTQS ) THEN 01229 * 01230 * Perform bidiagonal SVD, computing left singular vectors 01231 * of bidiagonal matrix in RWORK(IRU) and computing right 01232 * singular vectors of bidiagonal matrix in RWORK(IRVT) 01233 * (CWorkspace: need 0) 01234 * (RWorkspace: need BDSPAC) 01235 * 01236 IRU = NRWORK 01237 IRVT = IRU + N*N 01238 NRWORK = IRVT + N*N 01239 CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), 01240 $ N, RWORK( IRVT ), N, DUM, IDUM, 01241 $ RWORK( NRWORK ), IWORK, INFO ) 01242 * 01243 * Copy real matrix RWORK(IRU) to complex matrix U 01244 * Overwrite U by left singular vectors of A 01245 * (CWorkspace: need 3*N, prefer 2*N+N*NB) 01246 * (RWorkspace: 0) 01247 * 01248 CALL CLASET( 'F', M, N, CZERO, CZERO, U, LDU ) 01249 CALL CLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU ) 01250 CALL CUNMBR( 'Q', 'L', 'N', M, N, N, A, LDA, 01251 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 01252 $ LWORK-NWORK+1, IERR ) 01253 * 01254 * Copy real matrix RWORK(IRVT) to complex matrix VT 01255 * Overwrite VT by right singular vectors of A 01256 * (CWorkspace: need 3*N, prefer 2*N+N*NB) 01257 * (RWorkspace: 0) 01258 * 01259 CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) 01260 CALL CUNMBR( 'P', 'R', 'C', N, N, N, A, LDA, 01261 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 01262 $ LWORK-NWORK+1, IERR ) 01263 ELSE 01264 * 01265 * Perform bidiagonal SVD, computing left singular vectors 01266 * of bidiagonal matrix in RWORK(IRU) and computing right 01267 * singular vectors of bidiagonal matrix in RWORK(IRVT) 01268 * (CWorkspace: need 0) 01269 * (RWorkspace: need BDSPAC) 01270 * 01271 IRU = NRWORK 01272 IRVT = IRU + N*N 01273 NRWORK = IRVT + N*N 01274 CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), 01275 $ N, RWORK( IRVT ), N, DUM, IDUM, 01276 $ RWORK( NRWORK ), IWORK, INFO ) 01277 * 01278 * Set the right corner of U to identity matrix 01279 * 01280 CALL CLASET( 'F', M, M, CZERO, CZERO, U, LDU ) 01281 IF( M.GT.N ) THEN 01282 CALL CLASET( 'F', M-N, M-N, CZERO, CONE, 01283 $ U( N+1, N+1 ), LDU ) 01284 END IF 01285 * 01286 * Copy real matrix RWORK(IRU) to complex matrix U 01287 * Overwrite U by left singular vectors of A 01288 * (CWorkspace: need 2*N+M, prefer 2*N+M*NB) 01289 * (RWorkspace: 0) 01290 * 01291 CALL CLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU ) 01292 CALL CUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA, 01293 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 01294 $ LWORK-NWORK+1, IERR ) 01295 * 01296 * Copy real matrix RWORK(IRVT) to complex matrix VT 01297 * Overwrite VT by right singular vectors of A 01298 * (CWorkspace: need 3*N, prefer 2*N+N*NB) 01299 * (RWorkspace: 0) 01300 * 01301 CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) 01302 CALL CUNMBR( 'P', 'R', 'C', N, N, N, A, LDA, 01303 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 01304 $ LWORK-NWORK+1, IERR ) 01305 END IF 01306 * 01307 END IF 01308 * 01309 ELSE 01310 * 01311 * A has more columns than rows. If A has sufficiently more 01312 * columns than rows, first reduce using the LQ decomposition (if 01313 * sufficient workspace available) 01314 * 01315 IF( N.GE.MNTHR1 ) THEN 01316 * 01317 IF( WNTQN ) THEN 01318 * 01319 * Path 1t (N much larger than M, JOBZ='N') 01320 * No singular vectors to be computed 01321 * 01322 ITAU = 1 01323 NWORK = ITAU + M 01324 * 01325 * Compute A=L*Q 01326 * (CWorkspace: need 2*M, prefer M+M*NB) 01327 * (RWorkspace: 0) 01328 * 01329 CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 01330 $ LWORK-NWORK+1, IERR ) 01331 * 01332 * Zero out above L 01333 * 01334 CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ), 01335 $ LDA ) 01336 IE = 1 01337 ITAUQ = 1 01338 ITAUP = ITAUQ + M 01339 NWORK = ITAUP + M 01340 * 01341 * Bidiagonalize L in A 01342 * (CWorkspace: need 3*M, prefer 2*M+2*M*NB) 01343 * (RWorkspace: need M) 01344 * 01345 CALL CGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 01346 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 01347 $ IERR ) 01348 NRWORK = IE + M 01349 * 01350 * Perform bidiagonal SVD, compute singular values only 01351 * (CWorkspace: 0) 01352 * (RWorkspace: need BDSPAN) 01353 * 01354 CALL SBDSDC( 'U', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1, 01355 $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) 01356 * 01357 ELSE IF( WNTQO ) THEN 01358 * 01359 * Path 2t (N much larger than M, JOBZ='O') 01360 * M right singular vectors to be overwritten on A and 01361 * M left singular vectors to be computed in U 01362 * 01363 IVT = 1 01364 LDWKVT = M 01365 * 01366 * WORK(IVT) is M by M 01367 * 01368 IL = IVT + LDWKVT*M 01369 IF( LWORK.GE.M*N+M*M+3*M ) THEN 01370 * 01371 * WORK(IL) M by N 01372 * 01373 LDWRKL = M 01374 CHUNK = N 01375 ELSE 01376 * 01377 * WORK(IL) is M by CHUNK 01378 * 01379 LDWRKL = M 01380 CHUNK = ( LWORK-M*M-3*M ) / M 01381 END IF 01382 ITAU = IL + LDWRKL*CHUNK 01383 NWORK = ITAU + M 01384 * 01385 * Compute A=L*Q 01386 * (CWorkspace: need 2*M, prefer M+M*NB) 01387 * (RWorkspace: 0) 01388 * 01389 CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 01390 $ LWORK-NWORK+1, IERR ) 01391 * 01392 * Copy L to WORK(IL), zeroing about above it 01393 * 01394 CALL CLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL ) 01395 CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, 01396 $ WORK( IL+LDWRKL ), LDWRKL ) 01397 * 01398 * Generate Q in A 01399 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 01400 * (RWorkspace: 0) 01401 * 01402 CALL CUNGLQ( M, N, M, A, LDA, WORK( ITAU ), 01403 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01404 IE = 1 01405 ITAUQ = ITAU 01406 ITAUP = ITAUQ + M 01407 NWORK = ITAUP + M 01408 * 01409 * Bidiagonalize L in WORK(IL) 01410 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) 01411 * (RWorkspace: need M) 01412 * 01413 CALL CGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ), 01414 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), 01415 $ LWORK-NWORK+1, IERR ) 01416 * 01417 * Perform bidiagonal SVD, computing left singular vectors 01418 * of bidiagonal matrix in RWORK(IRU) and computing right 01419 * singular vectors of bidiagonal matrix in RWORK(IRVT) 01420 * (CWorkspace: need 0) 01421 * (RWorkspace: need BDSPAC) 01422 * 01423 IRU = IE + M 01424 IRVT = IRU + M*M 01425 NRWORK = IRVT + M*M 01426 CALL SBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ), 01427 $ M, RWORK( IRVT ), M, DUM, IDUM, 01428 $ RWORK( NRWORK ), IWORK, INFO ) 01429 * 01430 * Copy real matrix RWORK(IRU) to complex matrix WORK(IU) 01431 * Overwrite WORK(IU) by the left singular vectors of L 01432 * (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB) 01433 * (RWorkspace: 0) 01434 * 01435 CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) 01436 CALL CUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL, 01437 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 01438 $ LWORK-NWORK+1, IERR ) 01439 * 01440 * Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT) 01441 * Overwrite WORK(IVT) by the right singular vectors of L 01442 * (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB) 01443 * (RWorkspace: 0) 01444 * 01445 CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ), 01446 $ LDWKVT ) 01447 CALL CUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL, 01448 $ WORK( ITAUP ), WORK( IVT ), LDWKVT, 01449 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01450 * 01451 * Multiply right singular vectors of L in WORK(IL) by Q 01452 * in A, storing result in WORK(IL) and copying to A 01453 * (CWorkspace: need 2*M*M, prefer M*M+M*N)) 01454 * (RWorkspace: 0) 01455 * 01456 DO 40 I = 1, N, CHUNK 01457 BLK = MIN( N-I+1, CHUNK ) 01458 CALL CGEMM( 'N', 'N', M, BLK, M, CONE, WORK( IVT ), M, 01459 $ A( 1, I ), LDA, CZERO, WORK( IL ), 01460 $ LDWRKL ) 01461 CALL CLACPY( 'F', M, BLK, WORK( IL ), LDWRKL, 01462 $ A( 1, I ), LDA ) 01463 40 CONTINUE 01464 * 01465 ELSE IF( WNTQS ) THEN 01466 * 01467 * Path 3t (N much larger than M, JOBZ='S') 01468 * M right singular vectors to be computed in VT and 01469 * M left singular vectors to be computed in U 01470 * 01471 IL = 1 01472 * 01473 * WORK(IL) is M by M 01474 * 01475 LDWRKL = M 01476 ITAU = IL + LDWRKL*M 01477 NWORK = ITAU + M 01478 * 01479 * Compute A=L*Q 01480 * (CWorkspace: need 2*M, prefer M+M*NB) 01481 * (RWorkspace: 0) 01482 * 01483 CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 01484 $ LWORK-NWORK+1, IERR ) 01485 * 01486 * Copy L to WORK(IL), zeroing out above it 01487 * 01488 CALL CLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL ) 01489 CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, 01490 $ WORK( IL+LDWRKL ), LDWRKL ) 01491 * 01492 * Generate Q in A 01493 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 01494 * (RWorkspace: 0) 01495 * 01496 CALL CUNGLQ( M, N, M, A, LDA, WORK( ITAU ), 01497 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01498 IE = 1 01499 ITAUQ = ITAU 01500 ITAUP = ITAUQ + M 01501 NWORK = ITAUP + M 01502 * 01503 * Bidiagonalize L in WORK(IL) 01504 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) 01505 * (RWorkspace: need M) 01506 * 01507 CALL CGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ), 01508 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), 01509 $ LWORK-NWORK+1, IERR ) 01510 * 01511 * Perform bidiagonal SVD, computing left singular vectors 01512 * of bidiagonal matrix in RWORK(IRU) and computing right 01513 * singular vectors of bidiagonal matrix in RWORK(IRVT) 01514 * (CWorkspace: need 0) 01515 * (RWorkspace: need BDSPAC) 01516 * 01517 IRU = IE + M 01518 IRVT = IRU + M*M 01519 NRWORK = IRVT + M*M 01520 CALL SBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ), 01521 $ M, RWORK( IRVT ), M, DUM, IDUM, 01522 $ RWORK( NRWORK ), IWORK, INFO ) 01523 * 01524 * Copy real matrix RWORK(IRU) to complex matrix U 01525 * Overwrite U by left singular vectors of L 01526 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB) 01527 * (RWorkspace: 0) 01528 * 01529 CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) 01530 CALL CUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL, 01531 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 01532 $ LWORK-NWORK+1, IERR ) 01533 * 01534 * Copy real matrix RWORK(IRVT) to complex matrix VT 01535 * Overwrite VT by left singular vectors of L 01536 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB) 01537 * (RWorkspace: 0) 01538 * 01539 CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT ) 01540 CALL CUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL, 01541 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 01542 $ LWORK-NWORK+1, IERR ) 01543 * 01544 * Copy VT to WORK(IL), multiply right singular vectors of L 01545 * in WORK(IL) by Q in A, storing result in VT 01546 * (CWorkspace: need M*M) 01547 * (RWorkspace: 0) 01548 * 01549 CALL CLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL ) 01550 CALL CGEMM( 'N', 'N', M, N, M, CONE, WORK( IL ), LDWRKL, 01551 $ A, LDA, CZERO, VT, LDVT ) 01552 * 01553 ELSE IF( WNTQA ) THEN 01554 * 01555 * Path 9t (N much larger than M, JOBZ='A') 01556 * N right singular vectors to be computed in VT and 01557 * M left singular vectors to be computed in U 01558 * 01559 IVT = 1 01560 * 01561 * WORK(IVT) is M by M 01562 * 01563 LDWKVT = M 01564 ITAU = IVT + LDWKVT*M 01565 NWORK = ITAU + M 01566 * 01567 * Compute A=L*Q, copying result to VT 01568 * (CWorkspace: need 2*M, prefer M+M*NB) 01569 * (RWorkspace: 0) 01570 * 01571 CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 01572 $ LWORK-NWORK+1, IERR ) 01573 CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT ) 01574 * 01575 * Generate Q in VT 01576 * (CWorkspace: need M+N, prefer M+N*NB) 01577 * (RWorkspace: 0) 01578 * 01579 CALL CUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ), 01580 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01581 * 01582 * Produce L in A, zeroing out above it 01583 * 01584 CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ), 01585 $ LDA ) 01586 IE = 1 01587 ITAUQ = ITAU 01588 ITAUP = ITAUQ + M 01589 NWORK = ITAUP + M 01590 * 01591 * Bidiagonalize L in A 01592 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) 01593 * (RWorkspace: need M) 01594 * 01595 CALL CGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 01596 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 01597 $ IERR ) 01598 * 01599 * Perform bidiagonal SVD, computing left singular vectors 01600 * of bidiagonal matrix in RWORK(IRU) and computing right 01601 * singular vectors of bidiagonal matrix in RWORK(IRVT) 01602 * (CWorkspace: need 0) 01603 * (RWorkspace: need BDSPAC) 01604 * 01605 IRU = IE + M 01606 IRVT = IRU + M*M 01607 NRWORK = IRVT + M*M 01608 CALL SBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ), 01609 $ M, RWORK( IRVT ), M, DUM, IDUM, 01610 $ RWORK( NRWORK ), IWORK, INFO ) 01611 * 01612 * Copy real matrix RWORK(IRU) to complex matrix U 01613 * Overwrite U by left singular vectors of L 01614 * (CWorkspace: need 3*M, prefer 2*M+M*NB) 01615 * (RWorkspace: 0) 01616 * 01617 CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) 01618 CALL CUNMBR( 'Q', 'L', 'N', M, M, M, A, LDA, 01619 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 01620 $ LWORK-NWORK+1, IERR ) 01621 * 01622 * Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT) 01623 * Overwrite WORK(IVT) by right singular vectors of L 01624 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB) 01625 * (RWorkspace: 0) 01626 * 01627 CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ), 01628 $ LDWKVT ) 01629 CALL CUNMBR( 'P', 'R', 'C', M, M, M, A, LDA, 01630 $ WORK( ITAUP ), WORK( IVT ), LDWKVT, 01631 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01632 * 01633 * Multiply right singular vectors of L in WORK(IVT) by 01634 * Q in VT, storing result in A 01635 * (CWorkspace: need M*M) 01636 * (RWorkspace: 0) 01637 * 01638 CALL CGEMM( 'N', 'N', M, N, M, CONE, WORK( IVT ), 01639 $ LDWKVT, VT, LDVT, CZERO, A, LDA ) 01640 * 01641 * Copy right singular vectors of A from A to VT 01642 * 01643 CALL CLACPY( 'F', M, N, A, LDA, VT, LDVT ) 01644 * 01645 END IF 01646 * 01647 ELSE IF( N.GE.MNTHR2 ) THEN 01648 * 01649 * MNTHR2 <= N < MNTHR1 01650 * 01651 * Path 5t (N much larger than M, but not as much as MNTHR1) 01652 * Reduce to bidiagonal form without QR decomposition, use 01653 * CUNGBR and matrix multiplication to compute singular vectors 01654 * 01655 * 01656 IE = 1 01657 NRWORK = IE + M 01658 ITAUQ = 1 01659 ITAUP = ITAUQ + M 01660 NWORK = ITAUP + M 01661 * 01662 * Bidiagonalize A 01663 * (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB) 01664 * (RWorkspace: M) 01665 * 01666 CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 01667 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 01668 $ IERR ) 01669 * 01670 IF( WNTQN ) THEN 01671 * 01672 * Compute singular values only 01673 * (Cworkspace: 0) 01674 * (Rworkspace: need BDSPAN) 01675 * 01676 CALL SBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1, 01677 $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) 01678 ELSE IF( WNTQO ) THEN 01679 IRVT = NRWORK 01680 IRU = IRVT + M*M 01681 NRWORK = IRU + M*M 01682 IVT = NWORK 01683 * 01684 * Copy A to U, generate Q 01685 * (Cworkspace: need 2*M, prefer M+M*NB) 01686 * (Rworkspace: 0) 01687 * 01688 CALL CLACPY( 'L', M, M, A, LDA, U, LDU ) 01689 CALL CUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ), 01690 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01691 * 01692 * Generate P**H in A 01693 * (Cworkspace: need 2*M, prefer M+M*NB) 01694 * (Rworkspace: 0) 01695 * 01696 CALL CUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), 01697 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01698 * 01699 LDWKVT = M 01700 IF( LWORK.GE.M*N+3*M ) THEN 01701 * 01702 * WORK( IVT ) is M by N 01703 * 01704 NWORK = IVT + LDWKVT*N 01705 CHUNK = N 01706 ELSE 01707 * 01708 * WORK( IVT ) is M by CHUNK 01709 * 01710 CHUNK = ( LWORK-3*M ) / M 01711 NWORK = IVT + LDWKVT*CHUNK 01712 END IF 01713 * 01714 * Perform bidiagonal SVD, computing left singular vectors 01715 * of bidiagonal matrix in RWORK(IRU) and computing right 01716 * singular vectors of bidiagonal matrix in RWORK(IRVT) 01717 * (CWorkspace: need 0) 01718 * (RWorkspace: need BDSPAC) 01719 * 01720 CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ), 01721 $ M, RWORK( IRVT ), M, DUM, IDUM, 01722 $ RWORK( NRWORK ), IWORK, INFO ) 01723 * 01724 * Multiply Q in U by real matrix RWORK(IRVT) 01725 * storing the result in WORK(IVT), copying to U 01726 * (Cworkspace: need 0) 01727 * (Rworkspace: need 2*M*M) 01728 * 01729 CALL CLACRM( M, M, U, LDU, RWORK( IRU ), M, WORK( IVT ), 01730 $ LDWKVT, RWORK( NRWORK ) ) 01731 CALL CLACPY( 'F', M, M, WORK( IVT ), LDWKVT, U, LDU ) 01732 * 01733 * Multiply RWORK(IRVT) by P**H in A, storing the 01734 * result in WORK(IVT), copying to A 01735 * (CWorkspace: need M*M, prefer M*N) 01736 * (Rworkspace: need 2*M*M, prefer 2*M*N) 01737 * 01738 NRWORK = IRU 01739 DO 50 I = 1, N, CHUNK 01740 BLK = MIN( N-I+1, CHUNK ) 01741 CALL CLARCM( M, BLK, RWORK( IRVT ), M, A( 1, I ), LDA, 01742 $ WORK( IVT ), LDWKVT, RWORK( NRWORK ) ) 01743 CALL CLACPY( 'F', M, BLK, WORK( IVT ), LDWKVT, 01744 $ A( 1, I ), LDA ) 01745 50 CONTINUE 01746 ELSE IF( WNTQS ) THEN 01747 * 01748 * Copy A to U, generate Q 01749 * (Cworkspace: need 2*M, prefer M+M*NB) 01750 * (Rworkspace: 0) 01751 * 01752 CALL CLACPY( 'L', M, M, A, LDA, U, LDU ) 01753 CALL CUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ), 01754 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01755 * 01756 * Copy A to VT, generate P**H 01757 * (Cworkspace: need 2*M, prefer M+M*NB) 01758 * (Rworkspace: 0) 01759 * 01760 CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT ) 01761 CALL CUNGBR( 'P', M, N, M, VT, LDVT, WORK( ITAUP ), 01762 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01763 * 01764 * Perform bidiagonal SVD, computing left singular vectors 01765 * of bidiagonal matrix in RWORK(IRU) and computing right 01766 * singular vectors of bidiagonal matrix in RWORK(IRVT) 01767 * (CWorkspace: need 0) 01768 * (RWorkspace: need BDSPAC) 01769 * 01770 IRVT = NRWORK 01771 IRU = IRVT + M*M 01772 NRWORK = IRU + M*M 01773 CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ), 01774 $ M, RWORK( IRVT ), M, DUM, IDUM, 01775 $ RWORK( NRWORK ), IWORK, INFO ) 01776 * 01777 * Multiply Q in U by real matrix RWORK(IRU), storing the 01778 * result in A, copying to U 01779 * (CWorkspace: need 0) 01780 * (Rworkspace: need 3*M*M) 01781 * 01782 CALL CLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA, 01783 $ RWORK( NRWORK ) ) 01784 CALL CLACPY( 'F', M, M, A, LDA, U, LDU ) 01785 * 01786 * Multiply real matrix RWORK(IRVT) by P**H in VT, 01787 * storing the result in A, copying to VT 01788 * (Cworkspace: need 0) 01789 * (Rworkspace: need M*M+2*M*N) 01790 * 01791 NRWORK = IRU 01792 CALL CLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA, 01793 $ RWORK( NRWORK ) ) 01794 CALL CLACPY( 'F', M, N, A, LDA, VT, LDVT ) 01795 ELSE 01796 * 01797 * Copy A to U, generate Q 01798 * (Cworkspace: need 2*M, prefer M+M*NB) 01799 * (Rworkspace: 0) 01800 * 01801 CALL CLACPY( 'L', M, M, A, LDA, U, LDU ) 01802 CALL CUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ), 01803 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01804 * 01805 * Copy A to VT, generate P**H 01806 * (Cworkspace: need 2*M, prefer M+M*NB) 01807 * (Rworkspace: 0) 01808 * 01809 CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT ) 01810 CALL CUNGBR( 'P', N, N, M, VT, LDVT, WORK( ITAUP ), 01811 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01812 * 01813 * Perform bidiagonal SVD, computing left singular vectors 01814 * of bidiagonal matrix in RWORK(IRU) and computing right 01815 * singular vectors of bidiagonal matrix in RWORK(IRVT) 01816 * (CWorkspace: need 0) 01817 * (RWorkspace: need BDSPAC) 01818 * 01819 IRVT = NRWORK 01820 IRU = IRVT + M*M 01821 NRWORK = IRU + M*M 01822 CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ), 01823 $ M, RWORK( IRVT ), M, DUM, IDUM, 01824 $ RWORK( NRWORK ), IWORK, INFO ) 01825 * 01826 * Multiply Q in U by real matrix RWORK(IRU), storing the 01827 * result in A, copying to U 01828 * (CWorkspace: need 0) 01829 * (Rworkspace: need 3*M*M) 01830 * 01831 CALL CLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA, 01832 $ RWORK( NRWORK ) ) 01833 CALL CLACPY( 'F', M, M, A, LDA, U, LDU ) 01834 * 01835 * Multiply real matrix RWORK(IRVT) by P**H in VT, 01836 * storing the result in A, copying to VT 01837 * (Cworkspace: need 0) 01838 * (Rworkspace: need M*M+2*M*N) 01839 * 01840 CALL CLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA, 01841 $ RWORK( NRWORK ) ) 01842 CALL CLACPY( 'F', M, N, A, LDA, VT, LDVT ) 01843 END IF 01844 * 01845 ELSE 01846 * 01847 * N .LT. MNTHR2 01848 * 01849 * Path 6t (N greater than M, but not much larger) 01850 * Reduce to bidiagonal form without LQ decomposition 01851 * Use CUNMBR to compute singular vectors 01852 * 01853 IE = 1 01854 NRWORK = IE + M 01855 ITAUQ = 1 01856 ITAUP = ITAUQ + M 01857 NWORK = ITAUP + M 01858 * 01859 * Bidiagonalize A 01860 * (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB) 01861 * (RWorkspace: M) 01862 * 01863 CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 01864 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 01865 $ IERR ) 01866 IF( WNTQN ) THEN 01867 * 01868 * Compute singular values only 01869 * (Cworkspace: 0) 01870 * (Rworkspace: need BDSPAN) 01871 * 01872 CALL SBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1, 01873 $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) 01874 ELSE IF( WNTQO ) THEN 01875 LDWKVT = M 01876 IVT = NWORK 01877 IF( LWORK.GE.M*N+3*M ) THEN 01878 * 01879 * WORK( IVT ) is M by N 01880 * 01881 CALL CLASET( 'F', M, N, CZERO, CZERO, WORK( IVT ), 01882 $ LDWKVT ) 01883 NWORK = IVT + LDWKVT*N 01884 ELSE 01885 * 01886 * WORK( IVT ) is M by CHUNK 01887 * 01888 CHUNK = ( LWORK-3*M ) / M 01889 NWORK = IVT + LDWKVT*CHUNK 01890 END IF 01891 * 01892 * Perform bidiagonal SVD, computing left singular vectors 01893 * of bidiagonal matrix in RWORK(IRU) and computing right 01894 * singular vectors of bidiagonal matrix in RWORK(IRVT) 01895 * (CWorkspace: need 0) 01896 * (RWorkspace: need BDSPAC) 01897 * 01898 IRVT = NRWORK 01899 IRU = IRVT + M*M 01900 NRWORK = IRU + M*M 01901 CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ), 01902 $ M, RWORK( IRVT ), M, DUM, IDUM, 01903 $ RWORK( NRWORK ), IWORK, INFO ) 01904 * 01905 * Copy real matrix RWORK(IRU) to complex matrix U 01906 * Overwrite U by left singular vectors of A 01907 * (Cworkspace: need 2*M, prefer M+M*NB) 01908 * (Rworkspace: need 0) 01909 * 01910 CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) 01911 CALL CUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA, 01912 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 01913 $ LWORK-NWORK+1, IERR ) 01914 * 01915 IF( LWORK.GE.M*N+3*M ) THEN 01916 * 01917 * Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT) 01918 * Overwrite WORK(IVT) by right singular vectors of A, 01919 * copying to A 01920 * (Cworkspace: need M*N+2*M, prefer M*N+M+M*NB) 01921 * (Rworkspace: need 0) 01922 * 01923 CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ), 01924 $ LDWKVT ) 01925 CALL CUNMBR( 'P', 'R', 'C', M, N, M, A, LDA, 01926 $ WORK( ITAUP ), WORK( IVT ), LDWKVT, 01927 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01928 CALL CLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA ) 01929 ELSE 01930 * 01931 * Generate P**H in A 01932 * (Cworkspace: need 2*M, prefer M+M*NB) 01933 * (Rworkspace: need 0) 01934 * 01935 CALL CUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), 01936 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01937 * 01938 * Multiply Q in A by real matrix RWORK(IRU), storing the 01939 * result in WORK(IU), copying to A 01940 * (CWorkspace: need M*M, prefer M*N) 01941 * (Rworkspace: need 3*M*M, prefer M*M+2*M*N) 01942 * 01943 NRWORK = IRU 01944 DO 60 I = 1, N, CHUNK 01945 BLK = MIN( N-I+1, CHUNK ) 01946 CALL CLARCM( M, BLK, RWORK( IRVT ), M, A( 1, I ), 01947 $ LDA, WORK( IVT ), LDWKVT, 01948 $ RWORK( NRWORK ) ) 01949 CALL CLACPY( 'F', M, BLK, WORK( IVT ), LDWKVT, 01950 $ A( 1, I ), LDA ) 01951 60 CONTINUE 01952 END IF 01953 ELSE IF( WNTQS ) THEN 01954 * 01955 * Perform bidiagonal SVD, computing left singular vectors 01956 * of bidiagonal matrix in RWORK(IRU) and computing right 01957 * singular vectors of bidiagonal matrix in RWORK(IRVT) 01958 * (CWorkspace: need 0) 01959 * (RWorkspace: need BDSPAC) 01960 * 01961 IRVT = NRWORK 01962 IRU = IRVT + M*M 01963 NRWORK = IRU + M*M 01964 CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ), 01965 $ M, RWORK( IRVT ), M, DUM, IDUM, 01966 $ RWORK( NRWORK ), IWORK, INFO ) 01967 * 01968 * Copy real matrix RWORK(IRU) to complex matrix U 01969 * Overwrite U by left singular vectors of A 01970 * (CWorkspace: need 3*M, prefer 2*M+M*NB) 01971 * (RWorkspace: M*M) 01972 * 01973 CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) 01974 CALL CUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA, 01975 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 01976 $ LWORK-NWORK+1, IERR ) 01977 * 01978 * Copy real matrix RWORK(IRVT) to complex matrix VT 01979 * Overwrite VT by right singular vectors of A 01980 * (CWorkspace: need 3*M, prefer 2*M+M*NB) 01981 * (RWorkspace: M*M) 01982 * 01983 CALL CLASET( 'F', M, N, CZERO, CZERO, VT, LDVT ) 01984 CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT ) 01985 CALL CUNMBR( 'P', 'R', 'C', M, N, M, A, LDA, 01986 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 01987 $ LWORK-NWORK+1, IERR ) 01988 ELSE 01989 * 01990 * Perform bidiagonal SVD, computing left singular vectors 01991 * of bidiagonal matrix in RWORK(IRU) and computing right 01992 * singular vectors of bidiagonal matrix in RWORK(IRVT) 01993 * (CWorkspace: need 0) 01994 * (RWorkspace: need BDSPAC) 01995 * 01996 IRVT = NRWORK 01997 IRU = IRVT + M*M 01998 NRWORK = IRU + M*M 01999 * 02000 CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ), 02001 $ M, RWORK( IRVT ), M, DUM, IDUM, 02002 $ RWORK( NRWORK ), IWORK, INFO ) 02003 * 02004 * Copy real matrix RWORK(IRU) to complex matrix U 02005 * Overwrite U by left singular vectors of A 02006 * (CWorkspace: need 3*M, prefer 2*M+M*NB) 02007 * (RWorkspace: M*M) 02008 * 02009 CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) 02010 CALL CUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA, 02011 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 02012 $ LWORK-NWORK+1, IERR ) 02013 * 02014 * Set all of VT to identity matrix 02015 * 02016 CALL CLASET( 'F', N, N, CZERO, CONE, VT, LDVT ) 02017 * 02018 * Copy real matrix RWORK(IRVT) to complex matrix VT 02019 * Overwrite VT by right singular vectors of A 02020 * (CWorkspace: need 2*M+N, prefer 2*M+N*NB) 02021 * (RWorkspace: M*M) 02022 * 02023 CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT ) 02024 CALL CUNMBR( 'P', 'R', 'C', N, N, M, A, LDA, 02025 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 02026 $ LWORK-NWORK+1, IERR ) 02027 END IF 02028 * 02029 END IF 02030 * 02031 END IF 02032 * 02033 * Undo scaling if necessary 02034 * 02035 IF( ISCL.EQ.1 ) THEN 02036 IF( ANRM.GT.BIGNUM ) 02037 $ CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN, 02038 $ IERR ) 02039 IF( INFO.NE.0 .AND. ANRM.GT.BIGNUM ) 02040 $ CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN-1, 1, 02041 $ RWORK( IE ), MINMN, IERR ) 02042 IF( ANRM.LT.SMLNUM ) 02043 $ CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN, 02044 $ IERR ) 02045 IF( INFO.NE.0 .AND. ANRM.LT.SMLNUM ) 02046 $ CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN-1, 1, 02047 $ RWORK( IE ), MINMN, IERR ) 02048 END IF 02049 * 02050 * Return optimal workspace in WORK(1) 02051 * 02052 WORK( 1 ) = MAXWRK 02053 * 02054 RETURN 02055 * 02056 * End of CGESDD 02057 * 02058 END