![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> CGESVD computes the singular value decomposition (SVD) for GE matrices</b> 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CGESVD + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgesvd.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgesvd.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgesvd.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, 00022 * WORK, LWORK, RWORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER JOBU, JOBVT 00026 * INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N 00027 * .. 00028 * .. Array Arguments .. 00029 * REAL RWORK( * ), S( * ) 00030 * COMPLEX A( LDA, * ), U( LDU, * ), VT( LDVT, * ), 00031 * $ WORK( * ) 00032 * .. 00033 * 00034 * 00035 *> \par Purpose: 00036 * ============= 00037 *> 00038 *> \verbatim 00039 *> 00040 *> CGESVD computes the singular value decomposition (SVD) of a complex 00041 *> M-by-N matrix A, optionally computing the left and/or right singular 00042 *> vectors. The SVD is written 00043 *> 00044 *> A = U * SIGMA * conjugate-transpose(V) 00045 *> 00046 *> where SIGMA is an M-by-N matrix which is zero except for its 00047 *> min(m,n) diagonal elements, U is an M-by-M unitary matrix, and 00048 *> V is an N-by-N unitary matrix. The diagonal elements of SIGMA 00049 *> are the singular values of A; they are real and non-negative, and 00050 *> are returned in descending order. The first min(m,n) columns of 00051 *> U and V are the left and right singular vectors of A. 00052 *> 00053 *> Note that the routine returns V**H, not V. 00054 *> \endverbatim 00055 * 00056 * Arguments: 00057 * ========== 00058 * 00059 *> \param[in] JOBU 00060 *> \verbatim 00061 *> JOBU is CHARACTER*1 00062 *> Specifies options for computing all or part of the matrix U: 00063 *> = 'A': all M columns of U are returned in array U: 00064 *> = 'S': the first min(m,n) columns of U (the left singular 00065 *> vectors) are returned in the array U; 00066 *> = 'O': the first min(m,n) columns of U (the left singular 00067 *> vectors) are overwritten on the array A; 00068 *> = 'N': no columns of U (no left singular vectors) are 00069 *> computed. 00070 *> \endverbatim 00071 *> 00072 *> \param[in] JOBVT 00073 *> \verbatim 00074 *> JOBVT is CHARACTER*1 00075 *> Specifies options for computing all or part of the matrix 00076 *> V**H: 00077 *> = 'A': all N rows of V**H are returned in the array VT; 00078 *> = 'S': the first min(m,n) rows of V**H (the right singular 00079 *> vectors) are returned in the array VT; 00080 *> = 'O': the first min(m,n) rows of V**H (the right singular 00081 *> vectors) are overwritten on the array A; 00082 *> = 'N': no rows of V**H (no right singular vectors) are 00083 *> computed. 00084 *> 00085 *> JOBVT and JOBU cannot both be 'O'. 00086 *> \endverbatim 00087 *> 00088 *> \param[in] M 00089 *> \verbatim 00090 *> M is INTEGER 00091 *> The number of rows of the input matrix A. M >= 0. 00092 *> \endverbatim 00093 *> 00094 *> \param[in] N 00095 *> \verbatim 00096 *> N is INTEGER 00097 *> The number of columns of the input matrix A. N >= 0. 00098 *> \endverbatim 00099 *> 00100 *> \param[in,out] A 00101 *> \verbatim 00102 *> A is COMPLEX array, dimension (LDA,N) 00103 *> On entry, the M-by-N matrix A. 00104 *> On exit, 00105 *> if JOBU = 'O', A is overwritten with the first min(m,n) 00106 *> columns of U (the left singular vectors, 00107 *> stored columnwise); 00108 *> if JOBVT = 'O', A is overwritten with the first min(m,n) 00109 *> rows of V**H (the right singular vectors, 00110 *> stored rowwise); 00111 *> if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A 00112 *> are destroyed. 00113 *> \endverbatim 00114 *> 00115 *> \param[in] LDA 00116 *> \verbatim 00117 *> LDA is INTEGER 00118 *> The leading dimension of the array A. LDA >= max(1,M). 00119 *> \endverbatim 00120 *> 00121 *> \param[out] S 00122 *> \verbatim 00123 *> S is REAL array, dimension (min(M,N)) 00124 *> The singular values of A, sorted so that S(i) >= S(i+1). 00125 *> \endverbatim 00126 *> 00127 *> \param[out] U 00128 *> \verbatim 00129 *> U is COMPLEX array, dimension (LDU,UCOL) 00130 *> (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'. 00131 *> If JOBU = 'A', U contains the M-by-M unitary matrix U; 00132 *> if JOBU = 'S', U contains the first min(m,n) columns of U 00133 *> (the left singular vectors, stored columnwise); 00134 *> if JOBU = 'N' or 'O', U is not referenced. 00135 *> \endverbatim 00136 *> 00137 *> \param[in] LDU 00138 *> \verbatim 00139 *> LDU is INTEGER 00140 *> The leading dimension of the array U. LDU >= 1; if 00141 *> JOBU = 'S' or 'A', LDU >= M. 00142 *> \endverbatim 00143 *> 00144 *> \param[out] VT 00145 *> \verbatim 00146 *> VT is COMPLEX array, dimension (LDVT,N) 00147 *> If JOBVT = 'A', VT contains the N-by-N unitary matrix 00148 *> V**H; 00149 *> if JOBVT = 'S', VT contains the first min(m,n) rows of 00150 *> V**H (the right singular vectors, stored rowwise); 00151 *> if JOBVT = 'N' or 'O', VT is not referenced. 00152 *> \endverbatim 00153 *> 00154 *> \param[in] LDVT 00155 *> \verbatim 00156 *> LDVT is INTEGER 00157 *> The leading dimension of the array VT. LDVT >= 1; if 00158 *> JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N). 00159 *> \endverbatim 00160 *> 00161 *> \param[out] WORK 00162 *> \verbatim 00163 *> WORK is COMPLEX array, dimension (MAX(1,LWORK)) 00164 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00165 *> \endverbatim 00166 *> 00167 *> \param[in] LWORK 00168 *> \verbatim 00169 *> LWORK is INTEGER 00170 *> The dimension of the array WORK. 00171 *> LWORK >= MAX(1,2*MIN(M,N)+MAX(M,N)). 00172 *> For good performance, LWORK should generally be larger. 00173 *> 00174 *> If LWORK = -1, then a workspace query is assumed; the routine 00175 *> only calculates the optimal size of the WORK array, returns 00176 *> this value as the first entry of the WORK array, and no error 00177 *> message related to LWORK is issued by XERBLA. 00178 *> \endverbatim 00179 *> 00180 *> \param[out] RWORK 00181 *> \verbatim 00182 *> RWORK is REAL array, dimension (5*min(M,N)) 00183 *> On exit, if INFO > 0, RWORK(1:MIN(M,N)-1) contains the 00184 *> unconverged superdiagonal elements of an upper bidiagonal 00185 *> matrix B whose diagonal is in S (not necessarily sorted). 00186 *> B satisfies A = U * B * VT, so it has the same singular 00187 *> values as A, and singular vectors related by U and VT. 00188 *> \endverbatim 00189 *> 00190 *> \param[out] INFO 00191 *> \verbatim 00192 *> INFO is INTEGER 00193 *> = 0: successful exit. 00194 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00195 *> > 0: if CBDSQR did not converge, INFO specifies how many 00196 *> superdiagonals of an intermediate bidiagonal form B 00197 *> did not converge to zero. See the description of RWORK 00198 *> above for details. 00199 *> \endverbatim 00200 * 00201 * Authors: 00202 * ======== 00203 * 00204 *> \author Univ. of Tennessee 00205 *> \author Univ. of California Berkeley 00206 *> \author Univ. of Colorado Denver 00207 *> \author NAG Ltd. 00208 * 00209 *> \date April 2012 00210 * 00211 *> \ingroup complexGEsing 00212 * 00213 * ===================================================================== 00214 SUBROUTINE CGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, 00215 $ WORK, LWORK, RWORK, INFO ) 00216 * 00217 * -- LAPACK driver routine (version 3.4.1) -- 00218 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00219 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00220 * April 2012 00221 * 00222 * .. Scalar Arguments .. 00223 CHARACTER JOBU, JOBVT 00224 INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N 00225 * .. 00226 * .. Array Arguments .. 00227 REAL RWORK( * ), S( * ) 00228 COMPLEX A( LDA, * ), U( LDU, * ), VT( LDVT, * ), 00229 $ WORK( * ) 00230 * .. 00231 * 00232 * ===================================================================== 00233 * 00234 * .. Parameters .. 00235 COMPLEX CZERO, CONE 00236 PARAMETER ( CZERO = ( 0.0E0, 0.0E0 ), 00237 $ CONE = ( 1.0E0, 0.0E0 ) ) 00238 REAL ZERO, ONE 00239 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 00240 * .. 00241 * .. Local Scalars .. 00242 LOGICAL LQUERY, WNTUA, WNTUAS, WNTUN, WNTUO, WNTUS, 00243 $ WNTVA, WNTVAS, WNTVN, WNTVO, WNTVS 00244 INTEGER BLK, CHUNK, I, IE, IERR, IR, IRWORK, ISCL, 00245 $ ITAU, ITAUP, ITAUQ, IU, IWORK, LDWRKR, LDWRKU, 00246 $ MAXWRK, MINMN, MINWRK, MNTHR, NCU, NCVT, NRU, 00247 $ NRVT, WRKBL 00248 INTEGER LWORK_CGEQRF, LWORK_CUNGQR_N, LWORK_CUNGQR_M, 00249 $ LWORK_CGEBRD, LWORK_CUNGBR_P, LWORK_CUNGBR_Q, 00250 $ LWORK_CGELQF, LWORK_CUNGLQ_N, LWORK_CUNGLQ_M 00251 REAL ANRM, BIGNUM, EPS, SMLNUM 00252 * .. 00253 * .. Local Arrays .. 00254 REAL DUM( 1 ) 00255 COMPLEX CDUM( 1 ) 00256 * .. 00257 * .. External Subroutines .. 00258 EXTERNAL CBDSQR, CGEBRD, CGELQF, CGEMM, CGEQRF, CLACPY, 00259 $ CLASCL, CLASET, CUNGBR, CUNGLQ, CUNGQR, CUNMBR, 00260 $ SLASCL, XERBLA 00261 * .. 00262 * .. External Functions .. 00263 LOGICAL LSAME 00264 INTEGER ILAENV 00265 REAL CLANGE, SLAMCH 00266 EXTERNAL LSAME, ILAENV, CLANGE, SLAMCH 00267 * .. 00268 * .. Intrinsic Functions .. 00269 INTRINSIC MAX, MIN, SQRT 00270 * .. 00271 * .. Executable Statements .. 00272 * 00273 * Test the input arguments 00274 * 00275 INFO = 0 00276 MINMN = MIN( M, N ) 00277 WNTUA = LSAME( JOBU, 'A' ) 00278 WNTUS = LSAME( JOBU, 'S' ) 00279 WNTUAS = WNTUA .OR. WNTUS 00280 WNTUO = LSAME( JOBU, 'O' ) 00281 WNTUN = LSAME( JOBU, 'N' ) 00282 WNTVA = LSAME( JOBVT, 'A' ) 00283 WNTVS = LSAME( JOBVT, 'S' ) 00284 WNTVAS = WNTVA .OR. WNTVS 00285 WNTVO = LSAME( JOBVT, 'O' ) 00286 WNTVN = LSAME( JOBVT, 'N' ) 00287 LQUERY = ( LWORK.EQ.-1 ) 00288 * 00289 IF( .NOT.( WNTUA .OR. WNTUS .OR. WNTUO .OR. WNTUN ) ) THEN 00290 INFO = -1 00291 ELSE IF( .NOT.( WNTVA .OR. WNTVS .OR. WNTVO .OR. WNTVN ) .OR. 00292 $ ( WNTVO .AND. WNTUO ) ) THEN 00293 INFO = -2 00294 ELSE IF( M.LT.0 ) THEN 00295 INFO = -3 00296 ELSE IF( N.LT.0 ) THEN 00297 INFO = -4 00298 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00299 INFO = -6 00300 ELSE IF( LDU.LT.1 .OR. ( WNTUAS .AND. LDU.LT.M ) ) THEN 00301 INFO = -9 00302 ELSE IF( LDVT.LT.1 .OR. ( WNTVA .AND. LDVT.LT.N ) .OR. 00303 $ ( WNTVS .AND. LDVT.LT.MINMN ) ) THEN 00304 INFO = -11 00305 END IF 00306 * 00307 * Compute workspace 00308 * (Note: Comments in the code beginning "Workspace:" describe the 00309 * minimal amount of workspace needed at that point in the code, 00310 * as well as the preferred amount for good performance. 00311 * CWorkspace refers to complex workspace, and RWorkspace to 00312 * real workspace. NB refers to the optimal block size for the 00313 * immediately following subroutine, as returned by ILAENV.) 00314 * 00315 IF( INFO.EQ.0 ) THEN 00316 MINWRK = 1 00317 MAXWRK = 1 00318 IF( M.GE.N .AND. MINMN.GT.0 ) THEN 00319 * 00320 * Space needed for ZBDSQR is BDSPAC = 5*N 00321 * 00322 MNTHR = ILAENV( 6, 'CGESVD', JOBU // JOBVT, M, N, 0, 0 ) 00323 * Compute space needed for CGEQRF 00324 CALL CGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR ) 00325 LWORK_CGEQRF=DUM(1) 00326 * Compute space needed for CUNGQR 00327 CALL CUNGQR( M, N, N, A, LDA, DUM(1), DUM(1), -1, IERR ) 00328 LWORK_CUNGQR_N=DUM(1) 00329 CALL CUNGQR( M, M, N, A, LDA, DUM(1), DUM(1), -1, IERR ) 00330 LWORK_CUNGQR_M=DUM(1) 00331 * Compute space needed for CGEBRD 00332 CALL CGEBRD( N, N, A, LDA, S, DUM(1), DUM(1), 00333 $ DUM(1), DUM(1), -1, IERR ) 00334 LWORK_CGEBRD=DUM(1) 00335 * Compute space needed for CUNGBR 00336 CALL CUNGBR( 'P', N, N, N, A, LDA, DUM(1), 00337 $ DUM(1), -1, IERR ) 00338 LWORK_CUNGBR_P=DUM(1) 00339 CALL CUNGBR( 'Q', N, N, N, A, LDA, DUM(1), 00340 $ DUM(1), -1, IERR ) 00341 LWORK_CUNGBR_Q=DUM(1) 00342 * 00343 MNTHR = ILAENV( 6, 'CGESVD', JOBU // JOBVT, M, N, 0, 0 ) 00344 IF( M.GE.MNTHR ) THEN 00345 IF( WNTUN ) THEN 00346 * 00347 * Path 1 (M much larger than N, JOBU='N') 00348 * 00349 MAXWRK = N + LWORK_CGEQRF 00350 MAXWRK = MAX( MAXWRK, 2*N+LWORK_CGEBRD ) 00351 IF( WNTVO .OR. WNTVAS ) 00352 $ MAXWRK = MAX( MAXWRK, 2*N+LWORK_CUNGBR_P ) 00353 MINWRK = 3*N 00354 ELSE IF( WNTUO .AND. WNTVN ) THEN 00355 * 00356 * Path 2 (M much larger than N, JOBU='O', JOBVT='N') 00357 * 00358 WRKBL = N + LWORK_CGEQRF 00359 WRKBL = MAX( WRKBL, N+LWORK_CUNGQR_N ) 00360 WRKBL = MAX( WRKBL, 2*N+LWORK_CGEBRD ) 00361 WRKBL = MAX( WRKBL, 2*N+LWORK_CUNGBR_Q ) 00362 MAXWRK = MAX( N*N+WRKBL, N*N+M*N ) 00363 MINWRK = 2*N + M 00364 ELSE IF( WNTUO .AND. WNTVAS ) THEN 00365 * 00366 * Path 3 (M much larger than N, JOBU='O', JOBVT='S' or 00367 * 'A') 00368 * 00369 WRKBL = N + LWORK_CGEQRF 00370 WRKBL = MAX( WRKBL, N+LWORK_CUNGQR_N ) 00371 WRKBL = MAX( WRKBL, 2*N+LWORK_CGEBRD ) 00372 WRKBL = MAX( WRKBL, 2*N+LWORK_CUNGBR_Q ) 00373 WRKBL = MAX( WRKBL, 2*N+LWORK_CUNGBR_P ) 00374 MAXWRK = MAX( N*N+WRKBL, N*N+M*N ) 00375 MINWRK = 2*N + M 00376 ELSE IF( WNTUS .AND. WNTVN ) THEN 00377 * 00378 * Path 4 (M much larger than N, JOBU='S', JOBVT='N') 00379 * 00380 WRKBL = N + LWORK_CGEQRF 00381 WRKBL = MAX( WRKBL, N+LWORK_CUNGQR_N ) 00382 WRKBL = MAX( WRKBL, 2*N+LWORK_CGEBRD ) 00383 WRKBL = MAX( WRKBL, 2*N+LWORK_CUNGBR_Q ) 00384 MAXWRK = N*N + WRKBL 00385 MINWRK = 2*N + M 00386 ELSE IF( WNTUS .AND. WNTVO ) THEN 00387 * 00388 * Path 5 (M much larger than N, JOBU='S', JOBVT='O') 00389 * 00390 WRKBL = N + LWORK_CGEQRF 00391 WRKBL = MAX( WRKBL, N+LWORK_CUNGQR_N ) 00392 WRKBL = MAX( WRKBL, 2*N+LWORK_CGEBRD ) 00393 WRKBL = MAX( WRKBL, 2*N+LWORK_CUNGBR_Q ) 00394 WRKBL = MAX( WRKBL, 2*N+LWORK_CUNGBR_P ) 00395 MAXWRK = 2*N*N + WRKBL 00396 MINWRK = 2*N + M 00397 ELSE IF( WNTUS .AND. WNTVAS ) THEN 00398 * 00399 * Path 6 (M much larger than N, JOBU='S', JOBVT='S' or 00400 * 'A') 00401 * 00402 WRKBL = N + LWORK_CGEQRF 00403 WRKBL = MAX( WRKBL, N+LWORK_CUNGQR_N ) 00404 WRKBL = MAX( WRKBL, 2*N+LWORK_CGEBRD ) 00405 WRKBL = MAX( WRKBL, 2*N+LWORK_CUNGBR_Q ) 00406 WRKBL = MAX( WRKBL, 2*N+LWORK_CUNGBR_P ) 00407 MAXWRK = N*N + WRKBL 00408 MINWRK = 2*N + M 00409 ELSE IF( WNTUA .AND. WNTVN ) THEN 00410 * 00411 * Path 7 (M much larger than N, JOBU='A', JOBVT='N') 00412 * 00413 WRKBL = N + LWORK_CGEQRF 00414 WRKBL = MAX( WRKBL, N+LWORK_CUNGQR_M ) 00415 WRKBL = MAX( WRKBL, 2*N+LWORK_CGEBRD ) 00416 WRKBL = MAX( WRKBL, 2*N+LWORK_CUNGBR_Q ) 00417 MAXWRK = N*N + WRKBL 00418 MINWRK = 2*N + M 00419 ELSE IF( WNTUA .AND. WNTVO ) THEN 00420 * 00421 * Path 8 (M much larger than N, JOBU='A', JOBVT='O') 00422 * 00423 WRKBL = N + LWORK_CGEQRF 00424 WRKBL = MAX( WRKBL, N+LWORK_CUNGQR_M ) 00425 WRKBL = MAX( WRKBL, 2*N+LWORK_CGEBRD ) 00426 WRKBL = MAX( WRKBL, 2*N+LWORK_CUNGBR_Q ) 00427 WRKBL = MAX( WRKBL, 2*N+LWORK_CUNGBR_P ) 00428 MAXWRK = 2*N*N + WRKBL 00429 MINWRK = 2*N + M 00430 ELSE IF( WNTUA .AND. WNTVAS ) THEN 00431 * 00432 * Path 9 (M much larger than N, JOBU='A', JOBVT='S' or 00433 * 'A') 00434 * 00435 WRKBL = N + LWORK_CGEQRF 00436 WRKBL = MAX( WRKBL, N+LWORK_CUNGQR_M ) 00437 WRKBL = MAX( WRKBL, 2*N+LWORK_CGEBRD ) 00438 WRKBL = MAX( WRKBL, 2*N+LWORK_CUNGBR_Q ) 00439 WRKBL = MAX( WRKBL, 2*N+LWORK_CUNGBR_P ) 00440 MAXWRK = N*N + WRKBL 00441 MINWRK = 2*N + M 00442 END IF 00443 ELSE 00444 * 00445 * Path 10 (M at least N, but not much larger) 00446 * 00447 CALL CGEBRD( M, N, A, LDA, S, DUM(1), DUM(1), 00448 $ DUM(1), DUM(1), -1, IERR ) 00449 LWORK_CGEBRD=DUM(1) 00450 MAXWRK = 2*N + LWORK_CGEBRD 00451 IF( WNTUS .OR. WNTUO ) THEN 00452 CALL CUNGBR( 'Q', M, N, N, A, LDA, DUM(1), 00453 $ DUM(1), -1, IERR ) 00454 LWORK_CUNGBR_Q=DUM(1) 00455 MAXWRK = MAX( MAXWRK, 2*N+LWORK_CUNGBR_Q ) 00456 END IF 00457 IF( WNTUA ) THEN 00458 CALL CUNGBR( 'Q', M, M, N, A, LDA, DUM(1), 00459 $ DUM(1), -1, IERR ) 00460 LWORK_CUNGBR_Q=DUM(1) 00461 MAXWRK = MAX( MAXWRK, 2*N+LWORK_CUNGBR_Q ) 00462 END IF 00463 IF( .NOT.WNTVN ) THEN 00464 MAXWRK = MAX( MAXWRK, 2*N+LWORK_CUNGBR_P ) 00465 MINWRK = 2*N + M 00466 END IF 00467 END IF 00468 ELSE IF( MINMN.GT.0 ) THEN 00469 * 00470 * Space needed for CBDSQR is BDSPAC = 5*M 00471 * 00472 MNTHR = ILAENV( 6, 'CGESVD', JOBU // JOBVT, M, N, 0, 0 ) 00473 * Compute space needed for CGELQF 00474 CALL CGELQF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR ) 00475 LWORK_CGELQF=DUM(1) 00476 * Compute space needed for CUNGLQ 00477 CALL CUNGLQ( N, N, M, DUM(1), N, DUM(1), DUM(1), -1, IERR ) 00478 LWORK_CUNGLQ_N=DUM(1) 00479 CALL CUNGLQ( M, N, M, A, LDA, DUM(1), DUM(1), -1, IERR ) 00480 LWORK_CUNGLQ_M=DUM(1) 00481 * Compute space needed for CGEBRD 00482 CALL CGEBRD( M, M, A, LDA, S, DUM(1), DUM(1), 00483 $ DUM(1), DUM(1), -1, IERR ) 00484 LWORK_CGEBRD=DUM(1) 00485 * Compute space needed for CUNGBR P 00486 CALL CUNGBR( 'P', M, M, M, A, N, DUM(1), 00487 $ DUM(1), -1, IERR ) 00488 LWORK_CUNGBR_P=DUM(1) 00489 * Compute space needed for CUNGBR Q 00490 CALL CUNGBR( 'Q', M, M, M, A, N, DUM(1), 00491 $ DUM(1), -1, IERR ) 00492 LWORK_CUNGBR_Q=DUM(1) 00493 IF( N.GE.MNTHR ) THEN 00494 IF( WNTVN ) THEN 00495 * 00496 * Path 1t(N much larger than M, JOBVT='N') 00497 * 00498 MAXWRK = M + LWORK_CGELQF 00499 MAXWRK = MAX( MAXWRK, 2*M+LWORK_CGEBRD ) 00500 IF( WNTUO .OR. WNTUAS ) 00501 $ MAXWRK = MAX( MAXWRK, 2*M+LWORK_CUNGBR_Q ) 00502 MINWRK = 3*M 00503 ELSE IF( WNTVO .AND. WNTUN ) THEN 00504 * 00505 * Path 2t(N much larger than M, JOBU='N', JOBVT='O') 00506 * 00507 WRKBL = M + LWORK_CGELQF 00508 WRKBL = MAX( WRKBL, M+LWORK_CUNGLQ_M ) 00509 WRKBL = MAX( WRKBL, 2*M+LWORK_CGEBRD ) 00510 WRKBL = MAX( WRKBL, 2*M+LWORK_CUNGBR_P ) 00511 MAXWRK = MAX( M*M+WRKBL, M*M+M*N ) 00512 MINWRK = 2*M + N 00513 ELSE IF( WNTVO .AND. WNTUAS ) THEN 00514 * 00515 * Path 3t(N much larger than M, JOBU='S' or 'A', 00516 * JOBVT='O') 00517 * 00518 WRKBL = M + LWORK_CGELQF 00519 WRKBL = MAX( WRKBL, M+LWORK_CUNGLQ_M ) 00520 WRKBL = MAX( WRKBL, 2*M+LWORK_CGEBRD ) 00521 WRKBL = MAX( WRKBL, 2*M+LWORK_CUNGBR_P ) 00522 WRKBL = MAX( WRKBL, 2*M+LWORK_CUNGBR_Q ) 00523 MAXWRK = MAX( M*M+WRKBL, M*M+M*N ) 00524 MINWRK = 2*M + N 00525 ELSE IF( WNTVS .AND. WNTUN ) THEN 00526 * 00527 * Path 4t(N much larger than M, JOBU='N', JOBVT='S') 00528 * 00529 WRKBL = M + LWORK_CGELQF 00530 WRKBL = MAX( WRKBL, M+LWORK_CUNGLQ_M ) 00531 WRKBL = MAX( WRKBL, 2*M+LWORK_CGEBRD ) 00532 WRKBL = MAX( WRKBL, 2*M+LWORK_CUNGBR_P ) 00533 MAXWRK = M*M + WRKBL 00534 MINWRK = 2*M + N 00535 ELSE IF( WNTVS .AND. WNTUO ) THEN 00536 * 00537 * Path 5t(N much larger than M, JOBU='O', JOBVT='S') 00538 * 00539 WRKBL = M + LWORK_CGELQF 00540 WRKBL = MAX( WRKBL, M+LWORK_CUNGLQ_M ) 00541 WRKBL = MAX( WRKBL, 2*M+LWORK_CGEBRD ) 00542 WRKBL = MAX( WRKBL, 2*M+LWORK_CUNGBR_P ) 00543 WRKBL = MAX( WRKBL, 2*M+LWORK_CUNGBR_Q ) 00544 MAXWRK = 2*M*M + WRKBL 00545 MINWRK = 2*M + N 00546 ELSE IF( WNTVS .AND. WNTUAS ) THEN 00547 * 00548 * Path 6t(N much larger than M, JOBU='S' or 'A', 00549 * JOBVT='S') 00550 * 00551 WRKBL = M + LWORK_CGELQF 00552 WRKBL = MAX( WRKBL, M+LWORK_CUNGLQ_M ) 00553 WRKBL = MAX( WRKBL, 2*M+LWORK_CGEBRD ) 00554 WRKBL = MAX( WRKBL, 2*M+LWORK_CUNGBR_P ) 00555 WRKBL = MAX( WRKBL, 2*M+LWORK_CUNGBR_Q ) 00556 MAXWRK = M*M + WRKBL 00557 MINWRK = 2*M + N 00558 ELSE IF( WNTVA .AND. WNTUN ) THEN 00559 * 00560 * Path 7t(N much larger than M, JOBU='N', JOBVT='A') 00561 * 00562 WRKBL = M + LWORK_CGELQF 00563 WRKBL = MAX( WRKBL, M+LWORK_CUNGLQ_N ) 00564 WRKBL = MAX( WRKBL, 2*M+LWORK_CGEBRD ) 00565 WRKBL = MAX( WRKBL, 2*M+LWORK_CUNGBR_P ) 00566 MAXWRK = M*M + WRKBL 00567 MINWRK = 2*M + N 00568 ELSE IF( WNTVA .AND. WNTUO ) THEN 00569 * 00570 * Path 8t(N much larger than M, JOBU='O', JOBVT='A') 00571 * 00572 WRKBL = M + LWORK_CGELQF 00573 WRKBL = MAX( WRKBL, M+LWORK_CUNGLQ_N ) 00574 WRKBL = MAX( WRKBL, 2*M+LWORK_CGEBRD ) 00575 WRKBL = MAX( WRKBL, 2*M+LWORK_CUNGBR_P ) 00576 WRKBL = MAX( WRKBL, 2*M+LWORK_CUNGBR_Q ) 00577 MAXWRK = 2*M*M + WRKBL 00578 MINWRK = 2*M + N 00579 ELSE IF( WNTVA .AND. WNTUAS ) THEN 00580 * 00581 * Path 9t(N much larger than M, JOBU='S' or 'A', 00582 * JOBVT='A') 00583 * 00584 WRKBL = M + LWORK_CGELQF 00585 WRKBL = MAX( WRKBL, M+LWORK_CUNGLQ_N ) 00586 WRKBL = MAX( WRKBL, 2*M+LWORK_CGEBRD ) 00587 WRKBL = MAX( WRKBL, 2*M+LWORK_CUNGBR_P ) 00588 WRKBL = MAX( WRKBL, 2*M+LWORK_CUNGBR_Q ) 00589 MAXWRK = M*M + WRKBL 00590 MINWRK = 2*M + N 00591 END IF 00592 ELSE 00593 * 00594 * Path 10t(N greater than M, but not much larger) 00595 * 00596 CALL CGEBRD( M, N, A, LDA, S, DUM(1), DUM(1), 00597 $ DUM(1), DUM(1), -1, IERR ) 00598 LWORK_CGEBRD=DUM(1) 00599 MAXWRK = 2*M + LWORK_CGEBRD 00600 IF( WNTVS .OR. WNTVO ) THEN 00601 * Compute space needed for CUNGBR P 00602 CALL CUNGBR( 'P', M, N, M, A, N, DUM(1), 00603 $ DUM(1), -1, IERR ) 00604 LWORK_CUNGBR_P=DUM(1) 00605 MAXWRK = MAX( MAXWRK, 2*M+LWORK_CUNGBR_P ) 00606 END IF 00607 IF( WNTVA ) THEN 00608 CALL CUNGBR( 'P', N, N, M, A, N, DUM(1), 00609 $ DUM(1), -1, IERR ) 00610 LWORK_CUNGBR_P=DUM(1) 00611 MAXWRK = MAX( MAXWRK, 2*M+LWORK_CUNGBR_P ) 00612 END IF 00613 IF( .NOT.WNTUN ) THEN 00614 MAXWRK = MAX( MAXWRK, 2*M+LWORK_CUNGBR_Q ) 00615 MINWRK = 2*M + N 00616 END IF 00617 END IF 00618 END IF 00619 MAXWRK = MAX( MINWRK, MAXWRK ) 00620 WORK( 1 ) = MAXWRK 00621 * 00622 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN 00623 INFO = -13 00624 END IF 00625 END IF 00626 * 00627 IF( INFO.NE.0 ) THEN 00628 CALL XERBLA( 'CGESVD', -INFO ) 00629 RETURN 00630 ELSE IF( LQUERY ) THEN 00631 RETURN 00632 END IF 00633 * 00634 * Quick return if possible 00635 * 00636 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 00637 RETURN 00638 END IF 00639 * 00640 * Get machine constants 00641 * 00642 EPS = SLAMCH( 'P' ) 00643 SMLNUM = SQRT( SLAMCH( 'S' ) ) / EPS 00644 BIGNUM = ONE / SMLNUM 00645 * 00646 * Scale A if max element outside range [SMLNUM,BIGNUM] 00647 * 00648 ANRM = CLANGE( 'M', M, N, A, LDA, DUM ) 00649 ISCL = 0 00650 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00651 ISCL = 1 00652 CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR ) 00653 ELSE IF( ANRM.GT.BIGNUM ) THEN 00654 ISCL = 1 00655 CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR ) 00656 END IF 00657 * 00658 IF( M.GE.N ) THEN 00659 * 00660 * A has at least as many rows as columns. If A has sufficiently 00661 * more rows than columns, first reduce using the QR 00662 * decomposition (if sufficient workspace available) 00663 * 00664 IF( M.GE.MNTHR ) THEN 00665 * 00666 IF( WNTUN ) THEN 00667 * 00668 * Path 1 (M much larger than N, JOBU='N') 00669 * No left singular vectors to be computed 00670 * 00671 ITAU = 1 00672 IWORK = ITAU + N 00673 * 00674 * Compute A=Q*R 00675 * (CWorkspace: need 2*N, prefer N+N*NB) 00676 * (RWorkspace: need 0) 00677 * 00678 CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ), 00679 $ LWORK-IWORK+1, IERR ) 00680 * 00681 * Zero out below R 00682 * 00683 CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ), 00684 $ LDA ) 00685 IE = 1 00686 ITAUQ = 1 00687 ITAUP = ITAUQ + N 00688 IWORK = ITAUP + N 00689 * 00690 * Bidiagonalize R in A 00691 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 00692 * (RWorkspace: need N) 00693 * 00694 CALL CGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 00695 $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, 00696 $ IERR ) 00697 NCVT = 0 00698 IF( WNTVO .OR. WNTVAS ) THEN 00699 * 00700 * If right singular vectors desired, generate P'. 00701 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) 00702 * (RWorkspace: 0) 00703 * 00704 CALL CUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), 00705 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 00706 NCVT = N 00707 END IF 00708 IRWORK = IE + N 00709 * 00710 * Perform bidiagonal QR iteration, computing right 00711 * singular vectors of A in A if desired 00712 * (CWorkspace: 0) 00713 * (RWorkspace: need BDSPAC) 00714 * 00715 CALL CBDSQR( 'U', N, NCVT, 0, 0, S, RWORK( IE ), A, LDA, 00716 $ CDUM, 1, CDUM, 1, RWORK( IRWORK ), INFO ) 00717 * 00718 * If right singular vectors desired in VT, copy them there 00719 * 00720 IF( WNTVAS ) 00721 $ CALL CLACPY( 'F', N, N, A, LDA, VT, LDVT ) 00722 * 00723 ELSE IF( WNTUO .AND. WNTVN ) THEN 00724 * 00725 * Path 2 (M much larger than N, JOBU='O', JOBVT='N') 00726 * N left singular vectors to be overwritten on A and 00727 * no right singular vectors to be computed 00728 * 00729 IF( LWORK.GE.N*N+3*N ) THEN 00730 * 00731 * Sufficient workspace for a fast algorithm 00732 * 00733 IR = 1 00734 IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*N ) THEN 00735 * 00736 * WORK(IU) is LDA by N, WORK(IR) is LDA by N 00737 * 00738 LDWRKU = LDA 00739 LDWRKR = LDA 00740 ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+N*N ) THEN 00741 * 00742 * WORK(IU) is LDA by N, WORK(IR) is N by N 00743 * 00744 LDWRKU = LDA 00745 LDWRKR = N 00746 ELSE 00747 * 00748 * WORK(IU) is LDWRKU by N, WORK(IR) is N by N 00749 * 00750 LDWRKU = ( LWORK-N*N ) / N 00751 LDWRKR = N 00752 END IF 00753 ITAU = IR + LDWRKR*N 00754 IWORK = 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 ), 00761 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 00762 * 00763 * Copy R to WORK(IR) and zero 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, 00767 $ WORK( IR+1 ), LDWRKR ) 00768 * 00769 * Generate Q in A 00770 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 00771 * (RWorkspace: 0) 00772 * 00773 CALL CUNGQR( M, N, N, A, LDA, WORK( ITAU ), 00774 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 00775 IE = 1 00776 ITAUQ = ITAU 00777 ITAUP = ITAUQ + N 00778 IWORK = 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 ), 00786 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 00787 * 00788 * Generate left vectors bidiagonalizing R 00789 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) 00790 * (RWorkspace: need 0) 00791 * 00792 CALL CUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, 00793 $ WORK( ITAUQ ), WORK( IWORK ), 00794 $ LWORK-IWORK+1, IERR ) 00795 IRWORK = IE + N 00796 * 00797 * Perform bidiagonal QR iteration, computing left 00798 * singular vectors of R in WORK(IR) 00799 * (CWorkspace: need N*N) 00800 * (RWorkspace: need BDSPAC) 00801 * 00802 CALL CBDSQR( 'U', N, 0, N, 0, S, RWORK( IE ), CDUM, 1, 00803 $ WORK( IR ), LDWRKR, CDUM, 1, 00804 $ RWORK( IRWORK ), INFO ) 00805 IU = ITAUQ 00806 * 00807 * Multiply Q in A by left singular vectors of R in 00808 * WORK(IR), storing result in WORK(IU) and copying to A 00809 * (CWorkspace: need N*N+N, prefer N*N+M*N) 00810 * (RWorkspace: 0) 00811 * 00812 DO 10 I = 1, M, LDWRKU 00813 CHUNK = MIN( M-I+1, LDWRKU ) 00814 CALL CGEMM( 'N', 'N', CHUNK, N, N, CONE, A( I, 1 ), 00815 $ LDA, WORK( IR ), LDWRKR, CZERO, 00816 $ WORK( IU ), LDWRKU ) 00817 CALL CLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU, 00818 $ A( I, 1 ), LDA ) 00819 10 CONTINUE 00820 * 00821 ELSE 00822 * 00823 * Insufficient workspace for a fast algorithm 00824 * 00825 IE = 1 00826 ITAUQ = 1 00827 ITAUP = ITAUQ + N 00828 IWORK = ITAUP + N 00829 * 00830 * Bidiagonalize A 00831 * (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB) 00832 * (RWorkspace: N) 00833 * 00834 CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), 00835 $ WORK( ITAUQ ), WORK( ITAUP ), 00836 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 00837 * 00838 * Generate left vectors bidiagonalizing A 00839 * (CWorkspace: need 3*N, prefer 2*N+N*NB) 00840 * (RWorkspace: 0) 00841 * 00842 CALL CUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ), 00843 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 00844 IRWORK = IE + N 00845 * 00846 * Perform bidiagonal QR iteration, computing left 00847 * singular vectors of A in A 00848 * (CWorkspace: need 0) 00849 * (RWorkspace: need BDSPAC) 00850 * 00851 CALL CBDSQR( 'U', N, 0, M, 0, S, RWORK( IE ), CDUM, 1, 00852 $ A, LDA, CDUM, 1, RWORK( IRWORK ), INFO ) 00853 * 00854 END IF 00855 * 00856 ELSE IF( WNTUO .AND. WNTVAS ) THEN 00857 * 00858 * Path 3 (M much larger than N, JOBU='O', JOBVT='S' or 'A') 00859 * N left singular vectors to be overwritten on A and 00860 * N right singular vectors to be computed in VT 00861 * 00862 IF( LWORK.GE.N*N+3*N ) THEN 00863 * 00864 * Sufficient workspace for a fast algorithm 00865 * 00866 IR = 1 00867 IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*N ) THEN 00868 * 00869 * WORK(IU) is LDA by N and WORK(IR) is LDA by N 00870 * 00871 LDWRKU = LDA 00872 LDWRKR = LDA 00873 ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+N*N ) THEN 00874 * 00875 * WORK(IU) is LDA by N and WORK(IR) is N by N 00876 * 00877 LDWRKU = LDA 00878 LDWRKR = N 00879 ELSE 00880 * 00881 * WORK(IU) is LDWRKU by N and WORK(IR) is N by N 00882 * 00883 LDWRKU = ( LWORK-N*N ) / N 00884 LDWRKR = N 00885 END IF 00886 ITAU = IR + LDWRKR*N 00887 IWORK = ITAU + N 00888 * 00889 * Compute A=Q*R 00890 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 00891 * (RWorkspace: 0) 00892 * 00893 CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), 00894 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 00895 * 00896 * Copy R to VT, zeroing out below it 00897 * 00898 CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT ) 00899 IF( N.GT.1 ) 00900 $ CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, 00901 $ VT( 2, 1 ), LDVT ) 00902 * 00903 * Generate Q in A 00904 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 00905 * (RWorkspace: 0) 00906 * 00907 CALL CUNGQR( M, N, N, A, LDA, WORK( ITAU ), 00908 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 00909 IE = 1 00910 ITAUQ = ITAU 00911 ITAUP = ITAUQ + N 00912 IWORK = ITAUP + N 00913 * 00914 * Bidiagonalize R in VT, copying result to WORK(IR) 00915 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) 00916 * (RWorkspace: need N) 00917 * 00918 CALL CGEBRD( N, N, VT, LDVT, S, RWORK( IE ), 00919 $ WORK( ITAUQ ), WORK( ITAUP ), 00920 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 00921 CALL CLACPY( 'L', N, N, VT, LDVT, WORK( IR ), LDWRKR ) 00922 * 00923 * Generate left vectors bidiagonalizing R in WORK(IR) 00924 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) 00925 * (RWorkspace: 0) 00926 * 00927 CALL CUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, 00928 $ WORK( ITAUQ ), WORK( IWORK ), 00929 $ LWORK-IWORK+1, IERR ) 00930 * 00931 * Generate right vectors bidiagonalizing R in VT 00932 * (CWorkspace: need N*N+3*N-1, prefer N*N+2*N+(N-1)*NB) 00933 * (RWorkspace: 0) 00934 * 00935 CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 00936 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 00937 IRWORK = IE + N 00938 * 00939 * Perform bidiagonal QR iteration, computing left 00940 * singular vectors of R in WORK(IR) and computing right 00941 * singular vectors of R in VT 00942 * (CWorkspace: need N*N) 00943 * (RWorkspace: need BDSPAC) 00944 * 00945 CALL CBDSQR( 'U', N, N, N, 0, S, RWORK( IE ), VT, 00946 $ LDVT, WORK( IR ), LDWRKR, CDUM, 1, 00947 $ RWORK( IRWORK ), INFO ) 00948 IU = ITAUQ 00949 * 00950 * Multiply Q in A by left singular vectors of R in 00951 * WORK(IR), storing result in WORK(IU) and copying to A 00952 * (CWorkspace: need N*N+N, prefer N*N+M*N) 00953 * (RWorkspace: 0) 00954 * 00955 DO 20 I = 1, M, LDWRKU 00956 CHUNK = MIN( M-I+1, LDWRKU ) 00957 CALL CGEMM( 'N', 'N', CHUNK, N, N, CONE, A( I, 1 ), 00958 $ LDA, WORK( IR ), LDWRKR, CZERO, 00959 $ WORK( IU ), LDWRKU ) 00960 CALL CLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU, 00961 $ A( I, 1 ), LDA ) 00962 20 CONTINUE 00963 * 00964 ELSE 00965 * 00966 * Insufficient workspace for a fast algorithm 00967 * 00968 ITAU = 1 00969 IWORK = ITAU + N 00970 * 00971 * Compute A=Q*R 00972 * (CWorkspace: need 2*N, prefer N+N*NB) 00973 * (RWorkspace: 0) 00974 * 00975 CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), 00976 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 00977 * 00978 * Copy R to VT, zeroing out below it 00979 * 00980 CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT ) 00981 IF( N.GT.1 ) 00982 $ CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, 00983 $ VT( 2, 1 ), LDVT ) 00984 * 00985 * Generate Q in A 00986 * (CWorkspace: need 2*N, prefer N+N*NB) 00987 * (RWorkspace: 0) 00988 * 00989 CALL CUNGQR( M, N, N, A, LDA, WORK( ITAU ), 00990 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 00991 IE = 1 00992 ITAUQ = ITAU 00993 ITAUP = ITAUQ + N 00994 IWORK = ITAUP + N 00995 * 00996 * Bidiagonalize R in VT 00997 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 00998 * (RWorkspace: N) 00999 * 01000 CALL CGEBRD( N, N, VT, LDVT, S, RWORK( IE ), 01001 $ WORK( ITAUQ ), WORK( ITAUP ), 01002 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01003 * 01004 * Multiply Q in A by left vectors bidiagonalizing R 01005 * (CWorkspace: need 2*N+M, prefer 2*N+M*NB) 01006 * (RWorkspace: 0) 01007 * 01008 CALL CUNMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT, 01009 $ WORK( ITAUQ ), A, LDA, WORK( IWORK ), 01010 $ LWORK-IWORK+1, IERR ) 01011 * 01012 * Generate right vectors bidiagonalizing R in VT 01013 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) 01014 * (RWorkspace: 0) 01015 * 01016 CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 01017 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01018 IRWORK = IE + N 01019 * 01020 * Perform bidiagonal QR iteration, computing left 01021 * singular vectors of A in A and computing right 01022 * singular vectors of A in VT 01023 * (CWorkspace: 0) 01024 * (RWorkspace: need BDSPAC) 01025 * 01026 CALL CBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), VT, 01027 $ LDVT, A, LDA, CDUM, 1, RWORK( IRWORK ), 01028 $ INFO ) 01029 * 01030 END IF 01031 * 01032 ELSE IF( WNTUS ) THEN 01033 * 01034 IF( WNTVN ) THEN 01035 * 01036 * Path 4 (M much larger than N, JOBU='S', JOBVT='N') 01037 * N left singular vectors to be computed in U and 01038 * no right singular vectors to be computed 01039 * 01040 IF( LWORK.GE.N*N+3*N ) THEN 01041 * 01042 * Sufficient workspace for a fast algorithm 01043 * 01044 IR = 1 01045 IF( LWORK.GE.WRKBL+LDA*N ) THEN 01046 * 01047 * WORK(IR) is LDA by N 01048 * 01049 LDWRKR = LDA 01050 ELSE 01051 * 01052 * WORK(IR) is N by N 01053 * 01054 LDWRKR = N 01055 END IF 01056 ITAU = IR + LDWRKR*N 01057 IWORK = ITAU + N 01058 * 01059 * Compute A=Q*R 01060 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 01061 * (RWorkspace: 0) 01062 * 01063 CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), 01064 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01065 * 01066 * Copy R to WORK(IR), zeroing out below it 01067 * 01068 CALL CLACPY( 'U', N, N, A, LDA, WORK( IR ), 01069 $ LDWRKR ) 01070 CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, 01071 $ WORK( IR+1 ), LDWRKR ) 01072 * 01073 * Generate Q in A 01074 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 01075 * (RWorkspace: 0) 01076 * 01077 CALL CUNGQR( M, N, N, A, LDA, WORK( ITAU ), 01078 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01079 IE = 1 01080 ITAUQ = ITAU 01081 ITAUP = ITAUQ + N 01082 IWORK = ITAUP + N 01083 * 01084 * Bidiagonalize R in WORK(IR) 01085 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) 01086 * (RWorkspace: need N) 01087 * 01088 CALL CGEBRD( N, N, WORK( IR ), LDWRKR, S, 01089 $ RWORK( IE ), WORK( ITAUQ ), 01090 $ WORK( ITAUP ), WORK( IWORK ), 01091 $ LWORK-IWORK+1, IERR ) 01092 * 01093 * Generate left vectors bidiagonalizing R in WORK(IR) 01094 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) 01095 * (RWorkspace: 0) 01096 * 01097 CALL CUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, 01098 $ WORK( ITAUQ ), WORK( IWORK ), 01099 $ LWORK-IWORK+1, IERR ) 01100 IRWORK = IE + N 01101 * 01102 * Perform bidiagonal QR iteration, computing left 01103 * singular vectors of R in WORK(IR) 01104 * (CWorkspace: need N*N) 01105 * (RWorkspace: need BDSPAC) 01106 * 01107 CALL CBDSQR( 'U', N, 0, N, 0, S, RWORK( IE ), CDUM, 01108 $ 1, WORK( IR ), LDWRKR, CDUM, 1, 01109 $ RWORK( IRWORK ), INFO ) 01110 * 01111 * Multiply Q in A by left singular vectors of R in 01112 * WORK(IR), storing result in U 01113 * (CWorkspace: need N*N) 01114 * (RWorkspace: 0) 01115 * 01116 CALL CGEMM( 'N', 'N', M, N, N, CONE, A, LDA, 01117 $ WORK( IR ), LDWRKR, CZERO, U, LDU ) 01118 * 01119 ELSE 01120 * 01121 * Insufficient workspace for a fast algorithm 01122 * 01123 ITAU = 1 01124 IWORK = ITAU + N 01125 * 01126 * Compute A=Q*R, copying result to U 01127 * (CWorkspace: need 2*N, prefer N+N*NB) 01128 * (RWorkspace: 0) 01129 * 01130 CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), 01131 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01132 CALL CLACPY( 'L', M, N, A, LDA, U, LDU ) 01133 * 01134 * Generate Q in U 01135 * (CWorkspace: need 2*N, prefer N+N*NB) 01136 * (RWorkspace: 0) 01137 * 01138 CALL CUNGQR( M, N, N, U, LDU, WORK( ITAU ), 01139 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01140 IE = 1 01141 ITAUQ = ITAU 01142 ITAUP = ITAUQ + N 01143 IWORK = ITAUP + N 01144 * 01145 * Zero out below R in A 01146 * 01147 CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, 01148 $ A( 2, 1 ), LDA ) 01149 * 01150 * Bidiagonalize R in A 01151 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 01152 * (RWorkspace: need N) 01153 * 01154 CALL CGEBRD( N, N, A, LDA, S, RWORK( IE ), 01155 $ WORK( ITAUQ ), WORK( ITAUP ), 01156 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01157 * 01158 * Multiply Q in U by left vectors bidiagonalizing R 01159 * (CWorkspace: need 2*N+M, prefer 2*N+M*NB) 01160 * (RWorkspace: 0) 01161 * 01162 CALL CUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA, 01163 $ WORK( ITAUQ ), U, LDU, WORK( IWORK ), 01164 $ LWORK-IWORK+1, IERR ) 01165 IRWORK = IE + N 01166 * 01167 * Perform bidiagonal QR iteration, computing left 01168 * singular vectors of A in U 01169 * (CWorkspace: 0) 01170 * (RWorkspace: need BDSPAC) 01171 * 01172 CALL CBDSQR( 'U', N, 0, M, 0, S, RWORK( IE ), CDUM, 01173 $ 1, U, LDU, CDUM, 1, RWORK( IRWORK ), 01174 $ INFO ) 01175 * 01176 END IF 01177 * 01178 ELSE IF( WNTVO ) THEN 01179 * 01180 * Path 5 (M much larger than N, JOBU='S', JOBVT='O') 01181 * N left singular vectors to be computed in U and 01182 * N right singular vectors to be overwritten on A 01183 * 01184 IF( LWORK.GE.2*N*N+3*N ) THEN 01185 * 01186 * Sufficient workspace for a fast algorithm 01187 * 01188 IU = 1 01189 IF( LWORK.GE.WRKBL+2*LDA*N ) THEN 01190 * 01191 * WORK(IU) is LDA by N and WORK(IR) is LDA by N 01192 * 01193 LDWRKU = LDA 01194 IR = IU + LDWRKU*N 01195 LDWRKR = LDA 01196 ELSE IF( LWORK.GE.WRKBL+( LDA+N )*N ) THEN 01197 * 01198 * WORK(IU) is LDA by N and WORK(IR) is N by N 01199 * 01200 LDWRKU = LDA 01201 IR = IU + LDWRKU*N 01202 LDWRKR = N 01203 ELSE 01204 * 01205 * WORK(IU) is N by N and WORK(IR) is N by N 01206 * 01207 LDWRKU = N 01208 IR = IU + LDWRKU*N 01209 LDWRKR = N 01210 END IF 01211 ITAU = IR + LDWRKR*N 01212 IWORK = ITAU + N 01213 * 01214 * Compute A=Q*R 01215 * (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB) 01216 * (RWorkspace: 0) 01217 * 01218 CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), 01219 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01220 * 01221 * Copy R to WORK(IU), zeroing out below it 01222 * 01223 CALL CLACPY( 'U', N, N, A, LDA, WORK( IU ), 01224 $ LDWRKU ) 01225 CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, 01226 $ WORK( IU+1 ), LDWRKU ) 01227 * 01228 * Generate Q in A 01229 * (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB) 01230 * (RWorkspace: 0) 01231 * 01232 CALL CUNGQR( M, N, N, A, LDA, WORK( ITAU ), 01233 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01234 IE = 1 01235 ITAUQ = ITAU 01236 ITAUP = ITAUQ + N 01237 IWORK = ITAUP + N 01238 * 01239 * Bidiagonalize R in WORK(IU), copying result to 01240 * WORK(IR) 01241 * (CWorkspace: need 2*N*N+3*N, 01242 * prefer 2*N*N+2*N+2*N*NB) 01243 * (RWorkspace: need N) 01244 * 01245 CALL CGEBRD( N, N, WORK( IU ), LDWRKU, S, 01246 $ RWORK( IE ), WORK( ITAUQ ), 01247 $ WORK( ITAUP ), WORK( IWORK ), 01248 $ LWORK-IWORK+1, IERR ) 01249 CALL CLACPY( 'U', N, N, WORK( IU ), LDWRKU, 01250 $ WORK( IR ), LDWRKR ) 01251 * 01252 * Generate left bidiagonalizing vectors in WORK(IU) 01253 * (CWorkspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB) 01254 * (RWorkspace: 0) 01255 * 01256 CALL CUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, 01257 $ WORK( ITAUQ ), WORK( IWORK ), 01258 $ LWORK-IWORK+1, IERR ) 01259 * 01260 * Generate right bidiagonalizing vectors in WORK(IR) 01261 * (CWorkspace: need 2*N*N+3*N-1, 01262 * prefer 2*N*N+2*N+(N-1)*NB) 01263 * (RWorkspace: 0) 01264 * 01265 CALL CUNGBR( 'P', N, N, N, WORK( IR ), LDWRKR, 01266 $ WORK( ITAUP ), WORK( IWORK ), 01267 $ LWORK-IWORK+1, IERR ) 01268 IRWORK = IE + N 01269 * 01270 * Perform bidiagonal QR iteration, computing left 01271 * singular vectors of R in WORK(IU) and computing 01272 * right singular vectors of R in WORK(IR) 01273 * (CWorkspace: need 2*N*N) 01274 * (RWorkspace: need BDSPAC) 01275 * 01276 CALL CBDSQR( 'U', N, N, N, 0, S, RWORK( IE ), 01277 $ WORK( IR ), LDWRKR, WORK( IU ), 01278 $ LDWRKU, CDUM, 1, RWORK( IRWORK ), 01279 $ INFO ) 01280 * 01281 * Multiply Q in A by left singular vectors of R in 01282 * WORK(IU), storing result in U 01283 * (CWorkspace: need N*N) 01284 * (RWorkspace: 0) 01285 * 01286 CALL CGEMM( 'N', 'N', M, N, N, CONE, A, LDA, 01287 $ WORK( IU ), LDWRKU, CZERO, U, LDU ) 01288 * 01289 * Copy right singular vectors of R to A 01290 * (CWorkspace: need N*N) 01291 * (RWorkspace: 0) 01292 * 01293 CALL CLACPY( 'F', N, N, WORK( IR ), LDWRKR, A, 01294 $ LDA ) 01295 * 01296 ELSE 01297 * 01298 * Insufficient workspace for a fast algorithm 01299 * 01300 ITAU = 1 01301 IWORK = ITAU + N 01302 * 01303 * Compute A=Q*R, copying result to U 01304 * (CWorkspace: need 2*N, prefer N+N*NB) 01305 * (RWorkspace: 0) 01306 * 01307 CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), 01308 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01309 CALL CLACPY( 'L', M, N, A, LDA, U, LDU ) 01310 * 01311 * Generate Q in U 01312 * (CWorkspace: need 2*N, prefer N+N*NB) 01313 * (RWorkspace: 0) 01314 * 01315 CALL CUNGQR( M, N, N, U, LDU, WORK( ITAU ), 01316 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01317 IE = 1 01318 ITAUQ = ITAU 01319 ITAUP = ITAUQ + N 01320 IWORK = ITAUP + N 01321 * 01322 * Zero out below R in A 01323 * 01324 CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, 01325 $ A( 2, 1 ), LDA ) 01326 * 01327 * Bidiagonalize R in A 01328 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 01329 * (RWorkspace: need N) 01330 * 01331 CALL CGEBRD( N, N, A, LDA, S, RWORK( IE ), 01332 $ WORK( ITAUQ ), WORK( ITAUP ), 01333 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01334 * 01335 * Multiply Q in U by left vectors bidiagonalizing R 01336 * (CWorkspace: need 2*N+M, prefer 2*N+M*NB) 01337 * (RWorkspace: 0) 01338 * 01339 CALL CUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA, 01340 $ WORK( ITAUQ ), U, LDU, WORK( IWORK ), 01341 $ LWORK-IWORK+1, IERR ) 01342 * 01343 * Generate right vectors bidiagonalizing R in A 01344 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) 01345 * (RWorkspace: 0) 01346 * 01347 CALL CUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), 01348 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01349 IRWORK = IE + N 01350 * 01351 * Perform bidiagonal QR iteration, computing left 01352 * singular vectors of A in U and computing right 01353 * singular vectors of A in A 01354 * (CWorkspace: 0) 01355 * (RWorkspace: need BDSPAC) 01356 * 01357 CALL CBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), A, 01358 $ LDA, U, LDU, CDUM, 1, RWORK( IRWORK ), 01359 $ INFO ) 01360 * 01361 END IF 01362 * 01363 ELSE IF( WNTVAS ) THEN 01364 * 01365 * Path 6 (M much larger than N, JOBU='S', JOBVT='S' 01366 * or 'A') 01367 * N left singular vectors to be computed in U and 01368 * N right singular vectors to be computed in VT 01369 * 01370 IF( LWORK.GE.N*N+3*N ) THEN 01371 * 01372 * Sufficient workspace for a fast algorithm 01373 * 01374 IU = 1 01375 IF( LWORK.GE.WRKBL+LDA*N ) THEN 01376 * 01377 * WORK(IU) is LDA by N 01378 * 01379 LDWRKU = LDA 01380 ELSE 01381 * 01382 * WORK(IU) is N by N 01383 * 01384 LDWRKU = N 01385 END IF 01386 ITAU = IU + LDWRKU*N 01387 IWORK = ITAU + N 01388 * 01389 * Compute A=Q*R 01390 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 01391 * (RWorkspace: 0) 01392 * 01393 CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), 01394 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01395 * 01396 * Copy R to WORK(IU), zeroing out below it 01397 * 01398 CALL CLACPY( 'U', N, N, A, LDA, WORK( IU ), 01399 $ LDWRKU ) 01400 CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, 01401 $ WORK( IU+1 ), LDWRKU ) 01402 * 01403 * Generate Q in A 01404 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 01405 * (RWorkspace: 0) 01406 * 01407 CALL CUNGQR( M, N, N, A, LDA, WORK( ITAU ), 01408 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01409 IE = 1 01410 ITAUQ = ITAU 01411 ITAUP = ITAUQ + N 01412 IWORK = ITAUP + N 01413 * 01414 * Bidiagonalize R in WORK(IU), copying result to VT 01415 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) 01416 * (RWorkspace: need N) 01417 * 01418 CALL CGEBRD( N, N, WORK( IU ), LDWRKU, S, 01419 $ RWORK( IE ), WORK( ITAUQ ), 01420 $ WORK( ITAUP ), WORK( IWORK ), 01421 $ LWORK-IWORK+1, IERR ) 01422 CALL CLACPY( 'U', N, N, WORK( IU ), LDWRKU, VT, 01423 $ LDVT ) 01424 * 01425 * Generate left bidiagonalizing vectors in WORK(IU) 01426 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) 01427 * (RWorkspace: 0) 01428 * 01429 CALL CUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, 01430 $ WORK( ITAUQ ), WORK( IWORK ), 01431 $ LWORK-IWORK+1, IERR ) 01432 * 01433 * Generate right bidiagonalizing vectors in VT 01434 * (CWorkspace: need N*N+3*N-1, 01435 * prefer N*N+2*N+(N-1)*NB) 01436 * (RWorkspace: 0) 01437 * 01438 CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 01439 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01440 IRWORK = IE + N 01441 * 01442 * Perform bidiagonal QR iteration, computing left 01443 * singular vectors of R in WORK(IU) and computing 01444 * right singular vectors of R in VT 01445 * (CWorkspace: need N*N) 01446 * (RWorkspace: need BDSPAC) 01447 * 01448 CALL CBDSQR( 'U', N, N, N, 0, S, RWORK( IE ), VT, 01449 $ LDVT, WORK( IU ), LDWRKU, CDUM, 1, 01450 $ RWORK( IRWORK ), INFO ) 01451 * 01452 * Multiply Q in A by left singular vectors of R in 01453 * WORK(IU), storing result in U 01454 * (CWorkspace: need N*N) 01455 * (RWorkspace: 0) 01456 * 01457 CALL CGEMM( 'N', 'N', M, N, N, CONE, A, LDA, 01458 $ WORK( IU ), LDWRKU, CZERO, U, LDU ) 01459 * 01460 ELSE 01461 * 01462 * Insufficient workspace for a fast algorithm 01463 * 01464 ITAU = 1 01465 IWORK = ITAU + N 01466 * 01467 * Compute A=Q*R, copying result to U 01468 * (CWorkspace: need 2*N, prefer N+N*NB) 01469 * (RWorkspace: 0) 01470 * 01471 CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), 01472 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01473 CALL CLACPY( 'L', M, N, A, LDA, U, LDU ) 01474 * 01475 * Generate Q in U 01476 * (CWorkspace: need 2*N, prefer N+N*NB) 01477 * (RWorkspace: 0) 01478 * 01479 CALL CUNGQR( M, N, N, U, LDU, WORK( ITAU ), 01480 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01481 * 01482 * Copy R to VT, zeroing out below it 01483 * 01484 CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT ) 01485 IF( N.GT.1 ) 01486 $ CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, 01487 $ VT( 2, 1 ), LDVT ) 01488 IE = 1 01489 ITAUQ = ITAU 01490 ITAUP = ITAUQ + N 01491 IWORK = ITAUP + N 01492 * 01493 * Bidiagonalize R in VT 01494 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 01495 * (RWorkspace: need N) 01496 * 01497 CALL CGEBRD( N, N, VT, LDVT, S, RWORK( IE ), 01498 $ WORK( ITAUQ ), WORK( ITAUP ), 01499 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01500 * 01501 * Multiply Q in U by left bidiagonalizing vectors 01502 * in VT 01503 * (CWorkspace: need 2*N+M, prefer 2*N+M*NB) 01504 * (RWorkspace: 0) 01505 * 01506 CALL CUNMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT, 01507 $ WORK( ITAUQ ), U, LDU, WORK( IWORK ), 01508 $ LWORK-IWORK+1, IERR ) 01509 * 01510 * Generate right bidiagonalizing vectors in VT 01511 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) 01512 * (RWorkspace: 0) 01513 * 01514 CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 01515 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01516 IRWORK = IE + N 01517 * 01518 * Perform bidiagonal QR iteration, computing left 01519 * singular vectors of A in U and computing right 01520 * singular vectors of A in VT 01521 * (CWorkspace: 0) 01522 * (RWorkspace: need BDSPAC) 01523 * 01524 CALL CBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), VT, 01525 $ LDVT, U, LDU, CDUM, 1, 01526 $ RWORK( IRWORK ), INFO ) 01527 * 01528 END IF 01529 * 01530 END IF 01531 * 01532 ELSE IF( WNTUA ) THEN 01533 * 01534 IF( WNTVN ) THEN 01535 * 01536 * Path 7 (M much larger than N, JOBU='A', JOBVT='N') 01537 * M left singular vectors to be computed in U and 01538 * no right singular vectors to be computed 01539 * 01540 IF( LWORK.GE.N*N+MAX( N+M, 3*N ) ) THEN 01541 * 01542 * Sufficient workspace for a fast algorithm 01543 * 01544 IR = 1 01545 IF( LWORK.GE.WRKBL+LDA*N ) THEN 01546 * 01547 * WORK(IR) is LDA by N 01548 * 01549 LDWRKR = LDA 01550 ELSE 01551 * 01552 * WORK(IR) is N by N 01553 * 01554 LDWRKR = N 01555 END IF 01556 ITAU = IR + LDWRKR*N 01557 IWORK = ITAU + N 01558 * 01559 * Compute A=Q*R, copying result to U 01560 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 01561 * (RWorkspace: 0) 01562 * 01563 CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), 01564 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01565 CALL CLACPY( 'L', M, N, A, LDA, U, LDU ) 01566 * 01567 * Copy R to WORK(IR), zeroing out below it 01568 * 01569 CALL CLACPY( 'U', N, N, A, LDA, WORK( IR ), 01570 $ LDWRKR ) 01571 CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, 01572 $ WORK( IR+1 ), LDWRKR ) 01573 * 01574 * Generate Q in U 01575 * (CWorkspace: need N*N+N+M, prefer N*N+N+M*NB) 01576 * (RWorkspace: 0) 01577 * 01578 CALL CUNGQR( M, M, N, U, LDU, WORK( ITAU ), 01579 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01580 IE = 1 01581 ITAUQ = ITAU 01582 ITAUP = ITAUQ + N 01583 IWORK = ITAUP + N 01584 * 01585 * Bidiagonalize R in WORK(IR) 01586 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) 01587 * (RWorkspace: need N) 01588 * 01589 CALL CGEBRD( N, N, WORK( IR ), LDWRKR, S, 01590 $ RWORK( IE ), WORK( ITAUQ ), 01591 $ WORK( ITAUP ), WORK( IWORK ), 01592 $ LWORK-IWORK+1, IERR ) 01593 * 01594 * Generate left bidiagonalizing vectors in WORK(IR) 01595 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) 01596 * (RWorkspace: 0) 01597 * 01598 CALL CUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, 01599 $ WORK( ITAUQ ), WORK( IWORK ), 01600 $ LWORK-IWORK+1, IERR ) 01601 IRWORK = IE + N 01602 * 01603 * Perform bidiagonal QR iteration, computing left 01604 * singular vectors of R in WORK(IR) 01605 * (CWorkspace: need N*N) 01606 * (RWorkspace: need BDSPAC) 01607 * 01608 CALL CBDSQR( 'U', N, 0, N, 0, S, RWORK( IE ), CDUM, 01609 $ 1, WORK( IR ), LDWRKR, CDUM, 1, 01610 $ RWORK( IRWORK ), INFO ) 01611 * 01612 * Multiply Q in U by left singular vectors of R in 01613 * WORK(IR), storing result in A 01614 * (CWorkspace: need N*N) 01615 * (RWorkspace: 0) 01616 * 01617 CALL CGEMM( 'N', 'N', M, N, N, CONE, U, LDU, 01618 $ WORK( IR ), LDWRKR, CZERO, A, LDA ) 01619 * 01620 * Copy left singular vectors of A from A to U 01621 * 01622 CALL CLACPY( 'F', M, N, A, LDA, U, LDU ) 01623 * 01624 ELSE 01625 * 01626 * Insufficient workspace for a fast algorithm 01627 * 01628 ITAU = 1 01629 IWORK = ITAU + N 01630 * 01631 * Compute A=Q*R, copying result to U 01632 * (CWorkspace: need 2*N, prefer N+N*NB) 01633 * (RWorkspace: 0) 01634 * 01635 CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), 01636 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01637 CALL CLACPY( 'L', M, N, A, LDA, U, LDU ) 01638 * 01639 * Generate Q in U 01640 * (CWorkspace: need N+M, prefer N+M*NB) 01641 * (RWorkspace: 0) 01642 * 01643 CALL CUNGQR( M, M, N, U, LDU, WORK( ITAU ), 01644 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01645 IE = 1 01646 ITAUQ = ITAU 01647 ITAUP = ITAUQ + N 01648 IWORK = ITAUP + N 01649 * 01650 * Zero out below R in A 01651 * 01652 CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, 01653 $ A( 2, 1 ), LDA ) 01654 * 01655 * Bidiagonalize R in A 01656 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 01657 * (RWorkspace: need N) 01658 * 01659 CALL CGEBRD( N, N, A, LDA, S, RWORK( IE ), 01660 $ WORK( ITAUQ ), WORK( ITAUP ), 01661 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01662 * 01663 * Multiply Q in U by left bidiagonalizing vectors 01664 * in A 01665 * (CWorkspace: need 2*N+M, prefer 2*N+M*NB) 01666 * (RWorkspace: 0) 01667 * 01668 CALL CUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA, 01669 $ WORK( ITAUQ ), U, LDU, WORK( IWORK ), 01670 $ LWORK-IWORK+1, IERR ) 01671 IRWORK = IE + N 01672 * 01673 * Perform bidiagonal QR iteration, computing left 01674 * singular vectors of A in U 01675 * (CWorkspace: 0) 01676 * (RWorkspace: need BDSPAC) 01677 * 01678 CALL CBDSQR( 'U', N, 0, M, 0, S, RWORK( IE ), CDUM, 01679 $ 1, U, LDU, CDUM, 1, RWORK( IRWORK ), 01680 $ INFO ) 01681 * 01682 END IF 01683 * 01684 ELSE IF( WNTVO ) THEN 01685 * 01686 * Path 8 (M much larger than N, JOBU='A', JOBVT='O') 01687 * M left singular vectors to be computed in U and 01688 * N right singular vectors to be overwritten on A 01689 * 01690 IF( LWORK.GE.2*N*N+MAX( N+M, 3*N ) ) THEN 01691 * 01692 * Sufficient workspace for a fast algorithm 01693 * 01694 IU = 1 01695 IF( LWORK.GE.WRKBL+2*LDA*N ) THEN 01696 * 01697 * WORK(IU) is LDA by N and WORK(IR) is LDA by N 01698 * 01699 LDWRKU = LDA 01700 IR = IU + LDWRKU*N 01701 LDWRKR = LDA 01702 ELSE IF( LWORK.GE.WRKBL+( LDA+N )*N ) THEN 01703 * 01704 * WORK(IU) is LDA by N and WORK(IR) is N by N 01705 * 01706 LDWRKU = LDA 01707 IR = IU + LDWRKU*N 01708 LDWRKR = N 01709 ELSE 01710 * 01711 * WORK(IU) is N by N and WORK(IR) is N by N 01712 * 01713 LDWRKU = N 01714 IR = IU + LDWRKU*N 01715 LDWRKR = N 01716 END IF 01717 ITAU = IR + LDWRKR*N 01718 IWORK = ITAU + N 01719 * 01720 * Compute A=Q*R, copying result to U 01721 * (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB) 01722 * (RWorkspace: 0) 01723 * 01724 CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), 01725 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01726 CALL CLACPY( 'L', M, N, A, LDA, U, LDU ) 01727 * 01728 * Generate Q in U 01729 * (CWorkspace: need 2*N*N+N+M, prefer 2*N*N+N+M*NB) 01730 * (RWorkspace: 0) 01731 * 01732 CALL CUNGQR( M, M, N, U, LDU, WORK( ITAU ), 01733 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01734 * 01735 * Copy R to WORK(IU), zeroing out below it 01736 * 01737 CALL CLACPY( 'U', N, N, A, LDA, WORK( IU ), 01738 $ LDWRKU ) 01739 CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, 01740 $ WORK( IU+1 ), LDWRKU ) 01741 IE = 1 01742 ITAUQ = ITAU 01743 ITAUP = ITAUQ + N 01744 IWORK = ITAUP + N 01745 * 01746 * Bidiagonalize R in WORK(IU), copying result to 01747 * WORK(IR) 01748 * (CWorkspace: need 2*N*N+3*N, 01749 * prefer 2*N*N+2*N+2*N*NB) 01750 * (RWorkspace: need N) 01751 * 01752 CALL CGEBRD( N, N, WORK( IU ), LDWRKU, S, 01753 $ RWORK( IE ), WORK( ITAUQ ), 01754 $ WORK( ITAUP ), WORK( IWORK ), 01755 $ LWORK-IWORK+1, IERR ) 01756 CALL CLACPY( 'U', N, N, WORK( IU ), LDWRKU, 01757 $ WORK( IR ), LDWRKR ) 01758 * 01759 * Generate left bidiagonalizing vectors in WORK(IU) 01760 * (CWorkspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB) 01761 * (RWorkspace: 0) 01762 * 01763 CALL CUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, 01764 $ WORK( ITAUQ ), WORK( IWORK ), 01765 $ LWORK-IWORK+1, IERR ) 01766 * 01767 * Generate right bidiagonalizing vectors in WORK(IR) 01768 * (CWorkspace: need 2*N*N+3*N-1, 01769 * prefer 2*N*N+2*N+(N-1)*NB) 01770 * (RWorkspace: 0) 01771 * 01772 CALL CUNGBR( 'P', N, N, N, WORK( IR ), LDWRKR, 01773 $ WORK( ITAUP ), WORK( IWORK ), 01774 $ LWORK-IWORK+1, IERR ) 01775 IRWORK = IE + N 01776 * 01777 * Perform bidiagonal QR iteration, computing left 01778 * singular vectors of R in WORK(IU) and computing 01779 * right singular vectors of R in WORK(IR) 01780 * (CWorkspace: need 2*N*N) 01781 * (RWorkspace: need BDSPAC) 01782 * 01783 CALL CBDSQR( 'U', N, N, N, 0, S, RWORK( IE ), 01784 $ WORK( IR ), LDWRKR, WORK( IU ), 01785 $ LDWRKU, CDUM, 1, RWORK( IRWORK ), 01786 $ INFO ) 01787 * 01788 * Multiply Q in U by left singular vectors of R in 01789 * WORK(IU), storing result in A 01790 * (CWorkspace: need N*N) 01791 * (RWorkspace: 0) 01792 * 01793 CALL CGEMM( 'N', 'N', M, N, N, CONE, U, LDU, 01794 $ WORK( IU ), LDWRKU, CZERO, A, LDA ) 01795 * 01796 * Copy left singular vectors of A from A to U 01797 * 01798 CALL CLACPY( 'F', M, N, A, LDA, U, LDU ) 01799 * 01800 * Copy right singular vectors of R from WORK(IR) to A 01801 * 01802 CALL CLACPY( 'F', N, N, WORK( IR ), LDWRKR, A, 01803 $ LDA ) 01804 * 01805 ELSE 01806 * 01807 * Insufficient workspace for a fast algorithm 01808 * 01809 ITAU = 1 01810 IWORK = ITAU + N 01811 * 01812 * Compute A=Q*R, copying result to U 01813 * (CWorkspace: need 2*N, prefer N+N*NB) 01814 * (RWorkspace: 0) 01815 * 01816 CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), 01817 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01818 CALL CLACPY( 'L', M, N, A, LDA, U, LDU ) 01819 * 01820 * Generate Q in U 01821 * (CWorkspace: need N+M, prefer N+M*NB) 01822 * (RWorkspace: 0) 01823 * 01824 CALL CUNGQR( M, M, N, U, LDU, WORK( ITAU ), 01825 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01826 IE = 1 01827 ITAUQ = ITAU 01828 ITAUP = ITAUQ + N 01829 IWORK = ITAUP + N 01830 * 01831 * Zero out below R in A 01832 * 01833 CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, 01834 $ A( 2, 1 ), LDA ) 01835 * 01836 * Bidiagonalize R in A 01837 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 01838 * (RWorkspace: need N) 01839 * 01840 CALL CGEBRD( N, N, A, LDA, S, RWORK( IE ), 01841 $ WORK( ITAUQ ), WORK( ITAUP ), 01842 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01843 * 01844 * Multiply Q in U by left bidiagonalizing vectors 01845 * in A 01846 * (CWorkspace: need 2*N+M, prefer 2*N+M*NB) 01847 * (RWorkspace: 0) 01848 * 01849 CALL CUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA, 01850 $ WORK( ITAUQ ), U, LDU, WORK( IWORK ), 01851 $ LWORK-IWORK+1, IERR ) 01852 * 01853 * Generate right bidiagonalizing vectors in A 01854 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) 01855 * (RWorkspace: 0) 01856 * 01857 CALL CUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), 01858 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01859 IRWORK = IE + N 01860 * 01861 * Perform bidiagonal QR iteration, computing left 01862 * singular vectors of A in U and computing right 01863 * singular vectors of A in A 01864 * (CWorkspace: 0) 01865 * (RWorkspace: need BDSPAC) 01866 * 01867 CALL CBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), A, 01868 $ LDA, U, LDU, CDUM, 1, RWORK( IRWORK ), 01869 $ INFO ) 01870 * 01871 END IF 01872 * 01873 ELSE IF( WNTVAS ) THEN 01874 * 01875 * Path 9 (M much larger than N, JOBU='A', JOBVT='S' 01876 * or 'A') 01877 * M left singular vectors to be computed in U and 01878 * N right singular vectors to be computed in VT 01879 * 01880 IF( LWORK.GE.N*N+MAX( N+M, 3*N ) ) THEN 01881 * 01882 * Sufficient workspace for a fast algorithm 01883 * 01884 IU = 1 01885 IF( LWORK.GE.WRKBL+LDA*N ) THEN 01886 * 01887 * WORK(IU) is LDA by N 01888 * 01889 LDWRKU = LDA 01890 ELSE 01891 * 01892 * WORK(IU) is N by N 01893 * 01894 LDWRKU = N 01895 END IF 01896 ITAU = IU + LDWRKU*N 01897 IWORK = ITAU + N 01898 * 01899 * Compute A=Q*R, copying result to U 01900 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 01901 * (RWorkspace: 0) 01902 * 01903 CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), 01904 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01905 CALL CLACPY( 'L', M, N, A, LDA, U, LDU ) 01906 * 01907 * Generate Q in U 01908 * (CWorkspace: need N*N+N+M, prefer N*N+N+M*NB) 01909 * (RWorkspace: 0) 01910 * 01911 CALL CUNGQR( M, M, N, U, LDU, WORK( ITAU ), 01912 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01913 * 01914 * Copy R to WORK(IU), zeroing out below it 01915 * 01916 CALL CLACPY( 'U', N, N, A, LDA, WORK( IU ), 01917 $ LDWRKU ) 01918 CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, 01919 $ WORK( IU+1 ), LDWRKU ) 01920 IE = 1 01921 ITAUQ = ITAU 01922 ITAUP = ITAUQ + N 01923 IWORK = ITAUP + N 01924 * 01925 * Bidiagonalize R in WORK(IU), copying result to VT 01926 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) 01927 * (RWorkspace: need N) 01928 * 01929 CALL CGEBRD( N, N, WORK( IU ), LDWRKU, S, 01930 $ RWORK( IE ), WORK( ITAUQ ), 01931 $ WORK( ITAUP ), WORK( IWORK ), 01932 $ LWORK-IWORK+1, IERR ) 01933 CALL CLACPY( 'U', N, N, WORK( IU ), LDWRKU, VT, 01934 $ LDVT ) 01935 * 01936 * Generate left bidiagonalizing vectors in WORK(IU) 01937 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) 01938 * (RWorkspace: 0) 01939 * 01940 CALL CUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, 01941 $ WORK( ITAUQ ), WORK( IWORK ), 01942 $ LWORK-IWORK+1, IERR ) 01943 * 01944 * Generate right bidiagonalizing vectors in VT 01945 * (CWorkspace: need N*N+3*N-1, 01946 * prefer N*N+2*N+(N-1)*NB) 01947 * (RWorkspace: need 0) 01948 * 01949 CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 01950 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01951 IRWORK = IE + N 01952 * 01953 * Perform bidiagonal QR iteration, computing left 01954 * singular vectors of R in WORK(IU) and computing 01955 * right singular vectors of R in VT 01956 * (CWorkspace: need N*N) 01957 * (RWorkspace: need BDSPAC) 01958 * 01959 CALL CBDSQR( 'U', N, N, N, 0, S, RWORK( IE ), VT, 01960 $ LDVT, WORK( IU ), LDWRKU, CDUM, 1, 01961 $ RWORK( IRWORK ), INFO ) 01962 * 01963 * Multiply Q in U by left singular vectors of R in 01964 * WORK(IU), storing result in A 01965 * (CWorkspace: need N*N) 01966 * (RWorkspace: 0) 01967 * 01968 CALL CGEMM( 'N', 'N', M, N, N, CONE, U, LDU, 01969 $ WORK( IU ), LDWRKU, CZERO, A, LDA ) 01970 * 01971 * Copy left singular vectors of A from A to U 01972 * 01973 CALL CLACPY( 'F', M, N, A, LDA, U, LDU ) 01974 * 01975 ELSE 01976 * 01977 * Insufficient workspace for a fast algorithm 01978 * 01979 ITAU = 1 01980 IWORK = ITAU + N 01981 * 01982 * Compute A=Q*R, copying result to U 01983 * (CWorkspace: need 2*N, prefer N+N*NB) 01984 * (RWorkspace: 0) 01985 * 01986 CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), 01987 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01988 CALL CLACPY( 'L', M, N, A, LDA, U, LDU ) 01989 * 01990 * Generate Q in U 01991 * (CWorkspace: need N+M, prefer N+M*NB) 01992 * (RWorkspace: 0) 01993 * 01994 CALL CUNGQR( M, M, N, U, LDU, WORK( ITAU ), 01995 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01996 * 01997 * Copy R from A to VT, zeroing out below it 01998 * 01999 CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT ) 02000 IF( N.GT.1 ) 02001 $ CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, 02002 $ VT( 2, 1 ), LDVT ) 02003 IE = 1 02004 ITAUQ = ITAU 02005 ITAUP = ITAUQ + N 02006 IWORK = ITAUP + N 02007 * 02008 * Bidiagonalize R in VT 02009 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 02010 * (RWorkspace: need N) 02011 * 02012 CALL CGEBRD( N, N, VT, LDVT, S, RWORK( IE ), 02013 $ WORK( ITAUQ ), WORK( ITAUP ), 02014 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02015 * 02016 * Multiply Q in U by left bidiagonalizing vectors 02017 * in VT 02018 * (CWorkspace: need 2*N+M, prefer 2*N+M*NB) 02019 * (RWorkspace: 0) 02020 * 02021 CALL CUNMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT, 02022 $ WORK( ITAUQ ), U, LDU, WORK( IWORK ), 02023 $ LWORK-IWORK+1, IERR ) 02024 * 02025 * Generate right bidiagonalizing vectors in VT 02026 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) 02027 * (RWorkspace: 0) 02028 * 02029 CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 02030 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02031 IRWORK = IE + N 02032 * 02033 * Perform bidiagonal QR iteration, computing left 02034 * singular vectors of A in U and computing right 02035 * singular vectors of A in VT 02036 * (CWorkspace: 0) 02037 * (RWorkspace: need BDSPAC) 02038 * 02039 CALL CBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), VT, 02040 $ LDVT, U, LDU, CDUM, 1, 02041 $ RWORK( IRWORK ), INFO ) 02042 * 02043 END IF 02044 * 02045 END IF 02046 * 02047 END IF 02048 * 02049 ELSE 02050 * 02051 * M .LT. MNTHR 02052 * 02053 * Path 10 (M at least N, but not much larger) 02054 * Reduce to bidiagonal form without QR decomposition 02055 * 02056 IE = 1 02057 ITAUQ = 1 02058 ITAUP = ITAUQ + N 02059 IWORK = ITAUP + N 02060 * 02061 * Bidiagonalize A 02062 * (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB) 02063 * (RWorkspace: need N) 02064 * 02065 CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 02066 $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, 02067 $ IERR ) 02068 IF( WNTUAS ) THEN 02069 * 02070 * If left singular vectors desired in U, copy result to U 02071 * and generate left bidiagonalizing vectors in U 02072 * (CWorkspace: need 2*N+NCU, prefer 2*N+NCU*NB) 02073 * (RWorkspace: 0) 02074 * 02075 CALL CLACPY( 'L', M, N, A, LDA, U, LDU ) 02076 IF( WNTUS ) 02077 $ NCU = N 02078 IF( WNTUA ) 02079 $ NCU = M 02080 CALL CUNGBR( 'Q', M, NCU, N, U, LDU, WORK( ITAUQ ), 02081 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02082 END IF 02083 IF( WNTVAS ) THEN 02084 * 02085 * If right singular vectors desired in VT, copy result to 02086 * VT and generate right bidiagonalizing vectors in VT 02087 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) 02088 * (RWorkspace: 0) 02089 * 02090 CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT ) 02091 CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 02092 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02093 END IF 02094 IF( WNTUO ) THEN 02095 * 02096 * If left singular vectors desired in A, generate left 02097 * bidiagonalizing vectors in A 02098 * (CWorkspace: need 3*N, prefer 2*N+N*NB) 02099 * (RWorkspace: 0) 02100 * 02101 CALL CUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ), 02102 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02103 END IF 02104 IF( WNTVO ) THEN 02105 * 02106 * If right singular vectors desired in A, generate right 02107 * bidiagonalizing vectors in A 02108 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) 02109 * (RWorkspace: 0) 02110 * 02111 CALL CUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), 02112 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02113 END IF 02114 IRWORK = IE + N 02115 IF( WNTUAS .OR. WNTUO ) 02116 $ NRU = M 02117 IF( WNTUN ) 02118 $ NRU = 0 02119 IF( WNTVAS .OR. WNTVO ) 02120 $ NCVT = N 02121 IF( WNTVN ) 02122 $ NCVT = 0 02123 IF( ( .NOT.WNTUO ) .AND. ( .NOT.WNTVO ) ) THEN 02124 * 02125 * Perform bidiagonal QR iteration, if desired, computing 02126 * left singular vectors in U and computing right singular 02127 * vectors in VT 02128 * (CWorkspace: 0) 02129 * (RWorkspace: need BDSPAC) 02130 * 02131 CALL CBDSQR( 'U', N, NCVT, NRU, 0, S, RWORK( IE ), VT, 02132 $ LDVT, U, LDU, CDUM, 1, RWORK( IRWORK ), 02133 $ INFO ) 02134 ELSE IF( ( .NOT.WNTUO ) .AND. WNTVO ) THEN 02135 * 02136 * Perform bidiagonal QR iteration, if desired, computing 02137 * left singular vectors in U and computing right singular 02138 * vectors in A 02139 * (CWorkspace: 0) 02140 * (RWorkspace: need BDSPAC) 02141 * 02142 CALL CBDSQR( 'U', N, NCVT, NRU, 0, S, RWORK( IE ), A, 02143 $ LDA, U, LDU, CDUM, 1, RWORK( IRWORK ), 02144 $ INFO ) 02145 ELSE 02146 * 02147 * Perform bidiagonal QR iteration, if desired, computing 02148 * left singular vectors in A and computing right singular 02149 * vectors in VT 02150 * (CWorkspace: 0) 02151 * (RWorkspace: need BDSPAC) 02152 * 02153 CALL CBDSQR( 'U', N, NCVT, NRU, 0, S, RWORK( IE ), VT, 02154 $ LDVT, A, LDA, CDUM, 1, RWORK( IRWORK ), 02155 $ INFO ) 02156 END IF 02157 * 02158 END IF 02159 * 02160 ELSE 02161 * 02162 * A has more columns than rows. If A has sufficiently more 02163 * columns than rows, first reduce using the LQ decomposition (if 02164 * sufficient workspace available) 02165 * 02166 IF( N.GE.MNTHR ) THEN 02167 * 02168 IF( WNTVN ) THEN 02169 * 02170 * Path 1t(N much larger than M, JOBVT='N') 02171 * No right singular vectors to be computed 02172 * 02173 ITAU = 1 02174 IWORK = ITAU + M 02175 * 02176 * Compute A=L*Q 02177 * (CWorkspace: need 2*M, prefer M+M*NB) 02178 * (RWorkspace: 0) 02179 * 02180 CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ), 02181 $ LWORK-IWORK+1, IERR ) 02182 * 02183 * Zero out above L 02184 * 02185 CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ), 02186 $ LDA ) 02187 IE = 1 02188 ITAUQ = 1 02189 ITAUP = ITAUQ + M 02190 IWORK = ITAUP + M 02191 * 02192 * Bidiagonalize L in A 02193 * (CWorkspace: need 3*M, prefer 2*M+2*M*NB) 02194 * (RWorkspace: need M) 02195 * 02196 CALL CGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 02197 $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, 02198 $ IERR ) 02199 IF( WNTUO .OR. WNTUAS ) THEN 02200 * 02201 * If left singular vectors desired, generate Q 02202 * (CWorkspace: need 3*M, prefer 2*M+M*NB) 02203 * (RWorkspace: 0) 02204 * 02205 CALL CUNGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ), 02206 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02207 END IF 02208 IRWORK = IE + M 02209 NRU = 0 02210 IF( WNTUO .OR. WNTUAS ) 02211 $ NRU = M 02212 * 02213 * Perform bidiagonal QR iteration, computing left singular 02214 * vectors of A in A if desired 02215 * (CWorkspace: 0) 02216 * (RWorkspace: need BDSPAC) 02217 * 02218 CALL CBDSQR( 'U', M, 0, NRU, 0, S, RWORK( IE ), CDUM, 1, 02219 $ A, LDA, CDUM, 1, RWORK( IRWORK ), INFO ) 02220 * 02221 * If left singular vectors desired in U, copy them there 02222 * 02223 IF( WNTUAS ) 02224 $ CALL CLACPY( 'F', M, M, A, LDA, U, LDU ) 02225 * 02226 ELSE IF( WNTVO .AND. WNTUN ) THEN 02227 * 02228 * Path 2t(N much larger than M, JOBU='N', JOBVT='O') 02229 * M right singular vectors to be overwritten on A and 02230 * no left singular vectors to be computed 02231 * 02232 IF( LWORK.GE.M*M+3*M ) THEN 02233 * 02234 * Sufficient workspace for a fast algorithm 02235 * 02236 IR = 1 02237 IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*M ) THEN 02238 * 02239 * WORK(IU) is LDA by N and WORK(IR) is LDA by M 02240 * 02241 LDWRKU = LDA 02242 CHUNK = N 02243 LDWRKR = LDA 02244 ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+M*M ) THEN 02245 * 02246 * WORK(IU) is LDA by N and WORK(IR) is M by M 02247 * 02248 LDWRKU = LDA 02249 CHUNK = N 02250 LDWRKR = M 02251 ELSE 02252 * 02253 * WORK(IU) is M by CHUNK and WORK(IR) is M by M 02254 * 02255 LDWRKU = M 02256 CHUNK = ( LWORK-M*M ) / M 02257 LDWRKR = M 02258 END IF 02259 ITAU = IR + LDWRKR*M 02260 IWORK = ITAU + M 02261 * 02262 * Compute A=L*Q 02263 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 02264 * (RWorkspace: 0) 02265 * 02266 CALL CGELQF( M, N, A, LDA, WORK( ITAU ), 02267 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02268 * 02269 * Copy L to WORK(IR) and zero out above it 02270 * 02271 CALL CLACPY( 'L', M, M, A, LDA, WORK( IR ), LDWRKR ) 02272 CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, 02273 $ WORK( IR+LDWRKR ), LDWRKR ) 02274 * 02275 * Generate Q in A 02276 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 02277 * (RWorkspace: 0) 02278 * 02279 CALL CUNGLQ( M, N, M, A, LDA, WORK( ITAU ), 02280 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02281 IE = 1 02282 ITAUQ = ITAU 02283 ITAUP = ITAUQ + M 02284 IWORK = ITAUP + M 02285 * 02286 * Bidiagonalize L in WORK(IR) 02287 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) 02288 * (RWorkspace: need M) 02289 * 02290 CALL CGEBRD( M, M, WORK( IR ), LDWRKR, S, RWORK( IE ), 02291 $ WORK( ITAUQ ), WORK( ITAUP ), 02292 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02293 * 02294 * Generate right vectors bidiagonalizing L 02295 * (CWorkspace: need M*M+3*M-1, prefer M*M+2*M+(M-1)*NB) 02296 * (RWorkspace: 0) 02297 * 02298 CALL CUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR, 02299 $ WORK( ITAUP ), WORK( IWORK ), 02300 $ LWORK-IWORK+1, IERR ) 02301 IRWORK = IE + M 02302 * 02303 * Perform bidiagonal QR iteration, computing right 02304 * singular vectors of L in WORK(IR) 02305 * (CWorkspace: need M*M) 02306 * (RWorkspace: need BDSPAC) 02307 * 02308 CALL CBDSQR( 'U', M, M, 0, 0, S, RWORK( IE ), 02309 $ WORK( IR ), LDWRKR, CDUM, 1, CDUM, 1, 02310 $ RWORK( IRWORK ), INFO ) 02311 IU = ITAUQ 02312 * 02313 * Multiply right singular vectors of L in WORK(IR) by Q 02314 * in A, storing result in WORK(IU) and copying to A 02315 * (CWorkspace: need M*M+M, prefer M*M+M*N) 02316 * (RWorkspace: 0) 02317 * 02318 DO 30 I = 1, N, CHUNK 02319 BLK = MIN( N-I+1, CHUNK ) 02320 CALL CGEMM( 'N', 'N', M, BLK, M, CONE, WORK( IR ), 02321 $ LDWRKR, A( 1, I ), LDA, CZERO, 02322 $ WORK( IU ), LDWRKU ) 02323 CALL CLACPY( 'F', M, BLK, WORK( IU ), LDWRKU, 02324 $ A( 1, I ), LDA ) 02325 30 CONTINUE 02326 * 02327 ELSE 02328 * 02329 * Insufficient workspace for a fast algorithm 02330 * 02331 IE = 1 02332 ITAUQ = 1 02333 ITAUP = ITAUQ + M 02334 IWORK = ITAUP + M 02335 * 02336 * Bidiagonalize A 02337 * (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB) 02338 * (RWorkspace: need M) 02339 * 02340 CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), 02341 $ WORK( ITAUQ ), WORK( ITAUP ), 02342 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02343 * 02344 * Generate right vectors bidiagonalizing A 02345 * (CWorkspace: need 3*M, prefer 2*M+M*NB) 02346 * (RWorkspace: 0) 02347 * 02348 CALL CUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), 02349 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02350 IRWORK = IE + M 02351 * 02352 * Perform bidiagonal QR iteration, computing right 02353 * singular vectors of A in A 02354 * (CWorkspace: 0) 02355 * (RWorkspace: need BDSPAC) 02356 * 02357 CALL CBDSQR( 'L', M, N, 0, 0, S, RWORK( IE ), A, LDA, 02358 $ CDUM, 1, CDUM, 1, RWORK( IRWORK ), INFO ) 02359 * 02360 END IF 02361 * 02362 ELSE IF( WNTVO .AND. WNTUAS ) THEN 02363 * 02364 * Path 3t(N much larger than M, JOBU='S' or 'A', JOBVT='O') 02365 * M right singular vectors to be overwritten on A and 02366 * M left singular vectors to be computed in U 02367 * 02368 IF( LWORK.GE.M*M+3*M ) THEN 02369 * 02370 * Sufficient workspace for a fast algorithm 02371 * 02372 IR = 1 02373 IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*M ) THEN 02374 * 02375 * WORK(IU) is LDA by N and WORK(IR) is LDA by M 02376 * 02377 LDWRKU = LDA 02378 CHUNK = N 02379 LDWRKR = LDA 02380 ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+M*M ) THEN 02381 * 02382 * WORK(IU) is LDA by N and WORK(IR) is M by M 02383 * 02384 LDWRKU = LDA 02385 CHUNK = N 02386 LDWRKR = M 02387 ELSE 02388 * 02389 * WORK(IU) is M by CHUNK and WORK(IR) is M by M 02390 * 02391 LDWRKU = M 02392 CHUNK = ( LWORK-M*M ) / M 02393 LDWRKR = M 02394 END IF 02395 ITAU = IR + LDWRKR*M 02396 IWORK = ITAU + M 02397 * 02398 * Compute A=L*Q 02399 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 02400 * (RWorkspace: 0) 02401 * 02402 CALL CGELQF( M, N, A, LDA, WORK( ITAU ), 02403 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02404 * 02405 * Copy L to U, zeroing about above it 02406 * 02407 CALL CLACPY( 'L', M, M, A, LDA, U, LDU ) 02408 CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, U( 1, 2 ), 02409 $ LDU ) 02410 * 02411 * Generate Q in A 02412 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 02413 * (RWorkspace: 0) 02414 * 02415 CALL CUNGLQ( M, N, M, A, LDA, WORK( ITAU ), 02416 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02417 IE = 1 02418 ITAUQ = ITAU 02419 ITAUP = ITAUQ + M 02420 IWORK = ITAUP + M 02421 * 02422 * Bidiagonalize L in U, copying result to WORK(IR) 02423 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) 02424 * (RWorkspace: need M) 02425 * 02426 CALL CGEBRD( M, M, U, LDU, S, RWORK( IE ), 02427 $ WORK( ITAUQ ), WORK( ITAUP ), 02428 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02429 CALL CLACPY( 'U', M, M, U, LDU, WORK( IR ), LDWRKR ) 02430 * 02431 * Generate right vectors bidiagonalizing L in WORK(IR) 02432 * (CWorkspace: need M*M+3*M-1, prefer M*M+2*M+(M-1)*NB) 02433 * (RWorkspace: 0) 02434 * 02435 CALL CUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR, 02436 $ WORK( ITAUP ), WORK( IWORK ), 02437 $ LWORK-IWORK+1, IERR ) 02438 * 02439 * Generate left vectors bidiagonalizing L in U 02440 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB) 02441 * (RWorkspace: 0) 02442 * 02443 CALL CUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), 02444 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02445 IRWORK = IE + M 02446 * 02447 * Perform bidiagonal QR iteration, computing left 02448 * singular vectors of L in U, and computing right 02449 * singular vectors of L in WORK(IR) 02450 * (CWorkspace: need M*M) 02451 * (RWorkspace: need BDSPAC) 02452 * 02453 CALL CBDSQR( 'U', M, M, M, 0, S, RWORK( IE ), 02454 $ WORK( IR ), LDWRKR, U, LDU, CDUM, 1, 02455 $ RWORK( IRWORK ), INFO ) 02456 IU = ITAUQ 02457 * 02458 * Multiply right singular vectors of L in WORK(IR) by Q 02459 * in A, storing result in WORK(IU) and copying to A 02460 * (CWorkspace: need M*M+M, prefer M*M+M*N)) 02461 * (RWorkspace: 0) 02462 * 02463 DO 40 I = 1, N, CHUNK 02464 BLK = MIN( N-I+1, CHUNK ) 02465 CALL CGEMM( 'N', 'N', M, BLK, M, CONE, WORK( IR ), 02466 $ LDWRKR, A( 1, I ), LDA, CZERO, 02467 $ WORK( IU ), LDWRKU ) 02468 CALL CLACPY( 'F', M, BLK, WORK( IU ), LDWRKU, 02469 $ A( 1, I ), LDA ) 02470 40 CONTINUE 02471 * 02472 ELSE 02473 * 02474 * Insufficient workspace for a fast algorithm 02475 * 02476 ITAU = 1 02477 IWORK = ITAU + M 02478 * 02479 * Compute A=L*Q 02480 * (CWorkspace: need 2*M, prefer M+M*NB) 02481 * (RWorkspace: 0) 02482 * 02483 CALL CGELQF( M, N, A, LDA, WORK( ITAU ), 02484 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02485 * 02486 * Copy L to U, zeroing out above it 02487 * 02488 CALL CLACPY( 'L', M, M, A, LDA, U, LDU ) 02489 CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, U( 1, 2 ), 02490 $ LDU ) 02491 * 02492 * Generate Q in A 02493 * (CWorkspace: need 2*M, prefer M+M*NB) 02494 * (RWorkspace: 0) 02495 * 02496 CALL CUNGLQ( M, N, M, A, LDA, WORK( ITAU ), 02497 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02498 IE = 1 02499 ITAUQ = ITAU 02500 ITAUP = ITAUQ + M 02501 IWORK = ITAUP + M 02502 * 02503 * Bidiagonalize L in U 02504 * (CWorkspace: need 3*M, prefer 2*M+2*M*NB) 02505 * (RWorkspace: need M) 02506 * 02507 CALL CGEBRD( M, M, U, LDU, S, RWORK( IE ), 02508 $ WORK( ITAUQ ), WORK( ITAUP ), 02509 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02510 * 02511 * Multiply right vectors bidiagonalizing L by Q in A 02512 * (CWorkspace: need 2*M+N, prefer 2*M+N*NB) 02513 * (RWorkspace: 0) 02514 * 02515 CALL CUNMBR( 'P', 'L', 'C', M, N, M, U, LDU, 02516 $ WORK( ITAUP ), A, LDA, WORK( IWORK ), 02517 $ LWORK-IWORK+1, IERR ) 02518 * 02519 * Generate left vectors bidiagonalizing L in U 02520 * (CWorkspace: need 3*M, prefer 2*M+M*NB) 02521 * (RWorkspace: 0) 02522 * 02523 CALL CUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), 02524 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02525 IRWORK = IE + M 02526 * 02527 * Perform bidiagonal QR iteration, computing left 02528 * singular vectors of A in U and computing right 02529 * singular vectors of A in A 02530 * (CWorkspace: 0) 02531 * (RWorkspace: need BDSPAC) 02532 * 02533 CALL CBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), A, LDA, 02534 $ U, LDU, CDUM, 1, RWORK( IRWORK ), INFO ) 02535 * 02536 END IF 02537 * 02538 ELSE IF( WNTVS ) THEN 02539 * 02540 IF( WNTUN ) THEN 02541 * 02542 * Path 4t(N much larger than M, JOBU='N', JOBVT='S') 02543 * M right singular vectors to be computed in VT and 02544 * no left singular vectors to be computed 02545 * 02546 IF( LWORK.GE.M*M+3*M ) THEN 02547 * 02548 * Sufficient workspace for a fast algorithm 02549 * 02550 IR = 1 02551 IF( LWORK.GE.WRKBL+LDA*M ) THEN 02552 * 02553 * WORK(IR) is LDA by M 02554 * 02555 LDWRKR = LDA 02556 ELSE 02557 * 02558 * WORK(IR) is M by M 02559 * 02560 LDWRKR = M 02561 END IF 02562 ITAU = IR + LDWRKR*M 02563 IWORK = ITAU + M 02564 * 02565 * Compute A=L*Q 02566 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 02567 * (RWorkspace: 0) 02568 * 02569 CALL CGELQF( M, N, A, LDA, WORK( ITAU ), 02570 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02571 * 02572 * Copy L to WORK(IR), zeroing out above it 02573 * 02574 CALL CLACPY( 'L', M, M, A, LDA, WORK( IR ), 02575 $ LDWRKR ) 02576 CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, 02577 $ WORK( IR+LDWRKR ), LDWRKR ) 02578 * 02579 * Generate Q in A 02580 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 02581 * (RWorkspace: 0) 02582 * 02583 CALL CUNGLQ( M, N, M, A, LDA, WORK( ITAU ), 02584 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02585 IE = 1 02586 ITAUQ = ITAU 02587 ITAUP = ITAUQ + M 02588 IWORK = ITAUP + M 02589 * 02590 * Bidiagonalize L in WORK(IR) 02591 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) 02592 * (RWorkspace: need M) 02593 * 02594 CALL CGEBRD( M, M, WORK( IR ), LDWRKR, S, 02595 $ RWORK( IE ), WORK( ITAUQ ), 02596 $ WORK( ITAUP ), WORK( IWORK ), 02597 $ LWORK-IWORK+1, IERR ) 02598 * 02599 * Generate right vectors bidiagonalizing L in 02600 * WORK(IR) 02601 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+(M-1)*NB) 02602 * (RWorkspace: 0) 02603 * 02604 CALL CUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR, 02605 $ WORK( ITAUP ), WORK( IWORK ), 02606 $ LWORK-IWORK+1, IERR ) 02607 IRWORK = IE + M 02608 * 02609 * Perform bidiagonal QR iteration, computing right 02610 * singular vectors of L in WORK(IR) 02611 * (CWorkspace: need M*M) 02612 * (RWorkspace: need BDSPAC) 02613 * 02614 CALL CBDSQR( 'U', M, M, 0, 0, S, RWORK( IE ), 02615 $ WORK( IR ), LDWRKR, CDUM, 1, CDUM, 1, 02616 $ RWORK( IRWORK ), INFO ) 02617 * 02618 * Multiply right singular vectors of L in WORK(IR) by 02619 * Q in A, storing result in VT 02620 * (CWorkspace: need M*M) 02621 * (RWorkspace: 0) 02622 * 02623 CALL CGEMM( 'N', 'N', M, N, M, CONE, WORK( IR ), 02624 $ LDWRKR, A, LDA, CZERO, VT, LDVT ) 02625 * 02626 ELSE 02627 * 02628 * Insufficient workspace for a fast algorithm 02629 * 02630 ITAU = 1 02631 IWORK = ITAU + M 02632 * 02633 * Compute A=L*Q 02634 * (CWorkspace: need 2*M, prefer M+M*NB) 02635 * (RWorkspace: 0) 02636 * 02637 CALL CGELQF( M, N, A, LDA, WORK( ITAU ), 02638 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02639 * 02640 * Copy result to VT 02641 * 02642 CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT ) 02643 * 02644 * Generate Q in VT 02645 * (CWorkspace: need 2*M, prefer M+M*NB) 02646 * (RWorkspace: 0) 02647 * 02648 CALL CUNGLQ( M, N, M, VT, LDVT, WORK( ITAU ), 02649 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02650 IE = 1 02651 ITAUQ = ITAU 02652 ITAUP = ITAUQ + M 02653 IWORK = ITAUP + M 02654 * 02655 * Zero out above L in A 02656 * 02657 CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, 02658 $ A( 1, 2 ), LDA ) 02659 * 02660 * Bidiagonalize L in A 02661 * (CWorkspace: need 3*M, prefer 2*M+2*M*NB) 02662 * (RWorkspace: need M) 02663 * 02664 CALL CGEBRD( M, M, A, LDA, S, RWORK( IE ), 02665 $ WORK( ITAUQ ), WORK( ITAUP ), 02666 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02667 * 02668 * Multiply right vectors bidiagonalizing L by Q in VT 02669 * (CWorkspace: need 2*M+N, prefer 2*M+N*NB) 02670 * (RWorkspace: 0) 02671 * 02672 CALL CUNMBR( 'P', 'L', 'C', M, N, M, A, LDA, 02673 $ WORK( ITAUP ), VT, LDVT, 02674 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02675 IRWORK = IE + M 02676 * 02677 * Perform bidiagonal QR iteration, computing right 02678 * singular vectors of A in VT 02679 * (CWorkspace: 0) 02680 * (RWorkspace: need BDSPAC) 02681 * 02682 CALL CBDSQR( 'U', M, N, 0, 0, S, RWORK( IE ), VT, 02683 $ LDVT, CDUM, 1, CDUM, 1, 02684 $ RWORK( IRWORK ), INFO ) 02685 * 02686 END IF 02687 * 02688 ELSE IF( WNTUO ) THEN 02689 * 02690 * Path 5t(N much larger than M, JOBU='O', JOBVT='S') 02691 * M right singular vectors to be computed in VT and 02692 * M left singular vectors to be overwritten on A 02693 * 02694 IF( LWORK.GE.2*M*M+3*M ) THEN 02695 * 02696 * Sufficient workspace for a fast algorithm 02697 * 02698 IU = 1 02699 IF( LWORK.GE.WRKBL+2*LDA*M ) THEN 02700 * 02701 * WORK(IU) is LDA by M and WORK(IR) is LDA by M 02702 * 02703 LDWRKU = LDA 02704 IR = IU + LDWRKU*M 02705 LDWRKR = LDA 02706 ELSE IF( LWORK.GE.WRKBL+( LDA+M )*M ) THEN 02707 * 02708 * WORK(IU) is LDA by M and WORK(IR) is M by M 02709 * 02710 LDWRKU = LDA 02711 IR = IU + LDWRKU*M 02712 LDWRKR = M 02713 ELSE 02714 * 02715 * WORK(IU) is M by M and WORK(IR) is M by M 02716 * 02717 LDWRKU = M 02718 IR = IU + LDWRKU*M 02719 LDWRKR = M 02720 END IF 02721 ITAU = IR + LDWRKR*M 02722 IWORK = ITAU + M 02723 * 02724 * Compute A=L*Q 02725 * (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB) 02726 * (RWorkspace: 0) 02727 * 02728 CALL CGELQF( M, N, A, LDA, WORK( ITAU ), 02729 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02730 * 02731 * Copy L to WORK(IU), zeroing out below it 02732 * 02733 CALL CLACPY( 'L', M, M, A, LDA, WORK( IU ), 02734 $ LDWRKU ) 02735 CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, 02736 $ WORK( IU+LDWRKU ), LDWRKU ) 02737 * 02738 * Generate Q in A 02739 * (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB) 02740 * (RWorkspace: 0) 02741 * 02742 CALL CUNGLQ( M, N, M, A, LDA, WORK( ITAU ), 02743 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02744 IE = 1 02745 ITAUQ = ITAU 02746 ITAUP = ITAUQ + M 02747 IWORK = ITAUP + M 02748 * 02749 * Bidiagonalize L in WORK(IU), copying result to 02750 * WORK(IR) 02751 * (CWorkspace: need 2*M*M+3*M, 02752 * prefer 2*M*M+2*M+2*M*NB) 02753 * (RWorkspace: need M) 02754 * 02755 CALL CGEBRD( M, M, WORK( IU ), LDWRKU, S, 02756 $ RWORK( IE ), WORK( ITAUQ ), 02757 $ WORK( ITAUP ), WORK( IWORK ), 02758 $ LWORK-IWORK+1, IERR ) 02759 CALL CLACPY( 'L', M, M, WORK( IU ), LDWRKU, 02760 $ WORK( IR ), LDWRKR ) 02761 * 02762 * Generate right bidiagonalizing vectors in WORK(IU) 02763 * (CWorkspace: need 2*M*M+3*M-1, 02764 * prefer 2*M*M+2*M+(M-1)*NB) 02765 * (RWorkspace: 0) 02766 * 02767 CALL CUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU, 02768 $ WORK( ITAUP ), WORK( IWORK ), 02769 $ LWORK-IWORK+1, IERR ) 02770 * 02771 * Generate left bidiagonalizing vectors in WORK(IR) 02772 * (CWorkspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB) 02773 * (RWorkspace: 0) 02774 * 02775 CALL CUNGBR( 'Q', M, M, M, WORK( IR ), LDWRKR, 02776 $ WORK( ITAUQ ), WORK( IWORK ), 02777 $ LWORK-IWORK+1, IERR ) 02778 IRWORK = IE + M 02779 * 02780 * Perform bidiagonal QR iteration, computing left 02781 * singular vectors of L in WORK(IR) and computing 02782 * right singular vectors of L in WORK(IU) 02783 * (CWorkspace: need 2*M*M) 02784 * (RWorkspace: need BDSPAC) 02785 * 02786 CALL CBDSQR( 'U', M, M, M, 0, S, RWORK( IE ), 02787 $ WORK( IU ), LDWRKU, WORK( IR ), 02788 $ LDWRKR, CDUM, 1, RWORK( IRWORK ), 02789 $ INFO ) 02790 * 02791 * Multiply right singular vectors of L in WORK(IU) by 02792 * Q in A, storing result in VT 02793 * (CWorkspace: need M*M) 02794 * (RWorkspace: 0) 02795 * 02796 CALL CGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ), 02797 $ LDWRKU, A, LDA, CZERO, VT, LDVT ) 02798 * 02799 * Copy left singular vectors of L to A 02800 * (CWorkspace: need M*M) 02801 * (RWorkspace: 0) 02802 * 02803 CALL CLACPY( 'F', M, M, WORK( IR ), LDWRKR, A, 02804 $ LDA ) 02805 * 02806 ELSE 02807 * 02808 * Insufficient workspace for a fast algorithm 02809 * 02810 ITAU = 1 02811 IWORK = ITAU + M 02812 * 02813 * Compute A=L*Q, copying result to VT 02814 * (CWorkspace: need 2*M, prefer M+M*NB) 02815 * (RWorkspace: 0) 02816 * 02817 CALL CGELQF( M, N, A, LDA, WORK( ITAU ), 02818 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02819 CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT ) 02820 * 02821 * Generate Q in VT 02822 * (CWorkspace: need 2*M, prefer M+M*NB) 02823 * (RWorkspace: 0) 02824 * 02825 CALL CUNGLQ( M, N, M, VT, LDVT, WORK( ITAU ), 02826 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02827 IE = 1 02828 ITAUQ = ITAU 02829 ITAUP = ITAUQ + M 02830 IWORK = ITAUP + M 02831 * 02832 * Zero out above L in A 02833 * 02834 CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, 02835 $ A( 1, 2 ), LDA ) 02836 * 02837 * Bidiagonalize L in A 02838 * (CWorkspace: need 3*M, prefer 2*M+2*M*NB) 02839 * (RWorkspace: need M) 02840 * 02841 CALL CGEBRD( M, M, A, LDA, S, RWORK( IE ), 02842 $ WORK( ITAUQ ), WORK( ITAUP ), 02843 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02844 * 02845 * Multiply right vectors bidiagonalizing L by Q in VT 02846 * (CWorkspace: need 2*M+N, prefer 2*M+N*NB) 02847 * (RWorkspace: 0) 02848 * 02849 CALL CUNMBR( 'P', 'L', 'C', M, N, M, A, LDA, 02850 $ WORK( ITAUP ), VT, LDVT, 02851 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02852 * 02853 * Generate left bidiagonalizing vectors of L in A 02854 * (CWorkspace: need 3*M, prefer 2*M+M*NB) 02855 * (RWorkspace: 0) 02856 * 02857 CALL CUNGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ), 02858 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02859 IRWORK = IE + M 02860 * 02861 * Perform bidiagonal QR iteration, computing left 02862 * singular vectors of A in A and computing right 02863 * singular vectors of A in VT 02864 * (CWorkspace: 0) 02865 * (RWorkspace: need BDSPAC) 02866 * 02867 CALL CBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT, 02868 $ LDVT, A, LDA, CDUM, 1, 02869 $ RWORK( IRWORK ), INFO ) 02870 * 02871 END IF 02872 * 02873 ELSE IF( WNTUAS ) THEN 02874 * 02875 * Path 6t(N much larger than M, JOBU='S' or 'A', 02876 * JOBVT='S') 02877 * M right singular vectors to be computed in VT and 02878 * M left singular vectors to be computed in U 02879 * 02880 IF( LWORK.GE.M*M+3*M ) THEN 02881 * 02882 * Sufficient workspace for a fast algorithm 02883 * 02884 IU = 1 02885 IF( LWORK.GE.WRKBL+LDA*M ) THEN 02886 * 02887 * WORK(IU) is LDA by N 02888 * 02889 LDWRKU = LDA 02890 ELSE 02891 * 02892 * WORK(IU) is LDA by M 02893 * 02894 LDWRKU = M 02895 END IF 02896 ITAU = IU + LDWRKU*M 02897 IWORK = ITAU + M 02898 * 02899 * Compute A=L*Q 02900 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 02901 * (RWorkspace: 0) 02902 * 02903 CALL CGELQF( M, N, A, LDA, WORK( ITAU ), 02904 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02905 * 02906 * Copy L to WORK(IU), zeroing out above it 02907 * 02908 CALL CLACPY( 'L', M, M, A, LDA, WORK( IU ), 02909 $ LDWRKU ) 02910 CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, 02911 $ WORK( IU+LDWRKU ), LDWRKU ) 02912 * 02913 * Generate Q in A 02914 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 02915 * (RWorkspace: 0) 02916 * 02917 CALL CUNGLQ( M, N, M, A, LDA, WORK( ITAU ), 02918 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02919 IE = 1 02920 ITAUQ = ITAU 02921 ITAUP = ITAUQ + M 02922 IWORK = ITAUP + M 02923 * 02924 * Bidiagonalize L in WORK(IU), copying result to U 02925 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) 02926 * (RWorkspace: need M) 02927 * 02928 CALL CGEBRD( M, M, WORK( IU ), LDWRKU, S, 02929 $ RWORK( IE ), WORK( ITAUQ ), 02930 $ WORK( ITAUP ), WORK( IWORK ), 02931 $ LWORK-IWORK+1, IERR ) 02932 CALL CLACPY( 'L', M, M, WORK( IU ), LDWRKU, U, 02933 $ LDU ) 02934 * 02935 * Generate right bidiagonalizing vectors in WORK(IU) 02936 * (CWorkspace: need M*M+3*M-1, 02937 * prefer M*M+2*M+(M-1)*NB) 02938 * (RWorkspace: 0) 02939 * 02940 CALL CUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU, 02941 $ WORK( ITAUP ), WORK( IWORK ), 02942 $ LWORK-IWORK+1, IERR ) 02943 * 02944 * Generate left bidiagonalizing vectors in U 02945 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB) 02946 * (RWorkspace: 0) 02947 * 02948 CALL CUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), 02949 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02950 IRWORK = IE + M 02951 * 02952 * Perform bidiagonal QR iteration, computing left 02953 * singular vectors of L in U and computing right 02954 * singular vectors of L in WORK(IU) 02955 * (CWorkspace: need M*M) 02956 * (RWorkspace: need BDSPAC) 02957 * 02958 CALL CBDSQR( 'U', M, M, M, 0, S, RWORK( IE ), 02959 $ WORK( IU ), LDWRKU, U, LDU, CDUM, 1, 02960 $ RWORK( IRWORK ), INFO ) 02961 * 02962 * Multiply right singular vectors of L in WORK(IU) by 02963 * Q in A, storing result in VT 02964 * (CWorkspace: need M*M) 02965 * (RWorkspace: 0) 02966 * 02967 CALL CGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ), 02968 $ LDWRKU, A, LDA, CZERO, VT, LDVT ) 02969 * 02970 ELSE 02971 * 02972 * Insufficient workspace for a fast algorithm 02973 * 02974 ITAU = 1 02975 IWORK = ITAU + M 02976 * 02977 * Compute A=L*Q, copying result to VT 02978 * (CWorkspace: need 2*M, prefer M+M*NB) 02979 * (RWorkspace: 0) 02980 * 02981 CALL CGELQF( M, N, A, LDA, WORK( ITAU ), 02982 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02983 CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT ) 02984 * 02985 * Generate Q in VT 02986 * (CWorkspace: need 2*M, prefer M+M*NB) 02987 * (RWorkspace: 0) 02988 * 02989 CALL CUNGLQ( M, N, M, VT, LDVT, WORK( ITAU ), 02990 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02991 * 02992 * Copy L to U, zeroing out above it 02993 * 02994 CALL CLACPY( 'L', M, M, A, LDA, U, LDU ) 02995 CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, 02996 $ U( 1, 2 ), LDU ) 02997 IE = 1 02998 ITAUQ = ITAU 02999 ITAUP = ITAUQ + M 03000 IWORK = ITAUP + M 03001 * 03002 * Bidiagonalize L in U 03003 * (CWorkspace: need 3*M, prefer 2*M+2*M*NB) 03004 * (RWorkspace: need M) 03005 * 03006 CALL CGEBRD( M, M, U, LDU, S, RWORK( IE ), 03007 $ WORK( ITAUQ ), WORK( ITAUP ), 03008 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03009 * 03010 * Multiply right bidiagonalizing vectors in U by Q 03011 * in VT 03012 * (CWorkspace: need 2*M+N, prefer 2*M+N*NB) 03013 * (RWorkspace: 0) 03014 * 03015 CALL CUNMBR( 'P', 'L', 'C', M, N, M, U, LDU, 03016 $ WORK( ITAUP ), VT, LDVT, 03017 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03018 * 03019 * Generate left bidiagonalizing vectors in U 03020 * (CWorkspace: need 3*M, prefer 2*M+M*NB) 03021 * (RWorkspace: 0) 03022 * 03023 CALL CUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), 03024 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03025 IRWORK = IE + M 03026 * 03027 * Perform bidiagonal QR iteration, computing left 03028 * singular vectors of A in U and computing right 03029 * singular vectors of A in VT 03030 * (CWorkspace: 0) 03031 * (RWorkspace: need BDSPAC) 03032 * 03033 CALL CBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT, 03034 $ LDVT, U, LDU, CDUM, 1, 03035 $ RWORK( IRWORK ), INFO ) 03036 * 03037 END IF 03038 * 03039 END IF 03040 * 03041 ELSE IF( WNTVA ) THEN 03042 * 03043 IF( WNTUN ) THEN 03044 * 03045 * Path 7t(N much larger than M, JOBU='N', JOBVT='A') 03046 * N right singular vectors to be computed in VT and 03047 * no left singular vectors to be computed 03048 * 03049 IF( LWORK.GE.M*M+MAX( N+M, 3*M ) ) THEN 03050 * 03051 * Sufficient workspace for a fast algorithm 03052 * 03053 IR = 1 03054 IF( LWORK.GE.WRKBL+LDA*M ) THEN 03055 * 03056 * WORK(IR) is LDA by M 03057 * 03058 LDWRKR = LDA 03059 ELSE 03060 * 03061 * WORK(IR) is M by M 03062 * 03063 LDWRKR = M 03064 END IF 03065 ITAU = IR + LDWRKR*M 03066 IWORK = ITAU + M 03067 * 03068 * Compute A=L*Q, copying result to VT 03069 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 03070 * (RWorkspace: 0) 03071 * 03072 CALL CGELQF( M, N, A, LDA, WORK( ITAU ), 03073 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03074 CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT ) 03075 * 03076 * Copy L to WORK(IR), zeroing out above it 03077 * 03078 CALL CLACPY( 'L', M, M, A, LDA, WORK( IR ), 03079 $ LDWRKR ) 03080 CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, 03081 $ WORK( IR+LDWRKR ), LDWRKR ) 03082 * 03083 * Generate Q in VT 03084 * (CWorkspace: need M*M+M+N, prefer M*M+M+N*NB) 03085 * (RWorkspace: 0) 03086 * 03087 CALL CUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ), 03088 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03089 IE = 1 03090 ITAUQ = ITAU 03091 ITAUP = ITAUQ + M 03092 IWORK = ITAUP + M 03093 * 03094 * Bidiagonalize L in WORK(IR) 03095 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) 03096 * (RWorkspace: need M) 03097 * 03098 CALL CGEBRD( M, M, WORK( IR ), LDWRKR, S, 03099 $ RWORK( IE ), WORK( ITAUQ ), 03100 $ WORK( ITAUP ), WORK( IWORK ), 03101 $ LWORK-IWORK+1, IERR ) 03102 * 03103 * Generate right bidiagonalizing vectors in WORK(IR) 03104 * (CWorkspace: need M*M+3*M-1, 03105 * prefer M*M+2*M+(M-1)*NB) 03106 * (RWorkspace: 0) 03107 * 03108 CALL CUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR, 03109 $ WORK( ITAUP ), WORK( IWORK ), 03110 $ LWORK-IWORK+1, IERR ) 03111 IRWORK = IE + M 03112 * 03113 * Perform bidiagonal QR iteration, computing right 03114 * singular vectors of L in WORK(IR) 03115 * (CWorkspace: need M*M) 03116 * (RWorkspace: need BDSPAC) 03117 * 03118 CALL CBDSQR( 'U', M, M, 0, 0, S, RWORK( IE ), 03119 $ WORK( IR ), LDWRKR, CDUM, 1, CDUM, 1, 03120 $ RWORK( IRWORK ), INFO ) 03121 * 03122 * Multiply right singular vectors of L in WORK(IR) by 03123 * Q in VT, storing result in A 03124 * (CWorkspace: need M*M) 03125 * (RWorkspace: 0) 03126 * 03127 CALL CGEMM( 'N', 'N', M, N, M, CONE, WORK( IR ), 03128 $ LDWRKR, VT, LDVT, CZERO, A, LDA ) 03129 * 03130 * Copy right singular vectors of A from A to VT 03131 * 03132 CALL CLACPY( 'F', M, N, A, LDA, VT, LDVT ) 03133 * 03134 ELSE 03135 * 03136 * Insufficient workspace for a fast algorithm 03137 * 03138 ITAU = 1 03139 IWORK = ITAU + M 03140 * 03141 * Compute A=L*Q, copying result to VT 03142 * (CWorkspace: need 2*M, prefer M+M*NB) 03143 * (RWorkspace: 0) 03144 * 03145 CALL CGELQF( M, N, A, LDA, WORK( ITAU ), 03146 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03147 CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT ) 03148 * 03149 * Generate Q in VT 03150 * (CWorkspace: need M+N, prefer M+N*NB) 03151 * (RWorkspace: 0) 03152 * 03153 CALL CUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ), 03154 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03155 IE = 1 03156 ITAUQ = ITAU 03157 ITAUP = ITAUQ + M 03158 IWORK = ITAUP + M 03159 * 03160 * Zero out above L in A 03161 * 03162 CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, 03163 $ A( 1, 2 ), LDA ) 03164 * 03165 * Bidiagonalize L in A 03166 * (CWorkspace: need 3*M, prefer 2*M+2*M*NB) 03167 * (RWorkspace: need M) 03168 * 03169 CALL CGEBRD( M, M, A, LDA, S, RWORK( IE ), 03170 $ WORK( ITAUQ ), WORK( ITAUP ), 03171 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03172 * 03173 * Multiply right bidiagonalizing vectors in A by Q 03174 * in VT 03175 * (CWorkspace: need 2*M+N, prefer 2*M+N*NB) 03176 * (RWorkspace: 0) 03177 * 03178 CALL CUNMBR( 'P', 'L', 'C', M, N, M, A, LDA, 03179 $ WORK( ITAUP ), VT, LDVT, 03180 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03181 IRWORK = IE + M 03182 * 03183 * Perform bidiagonal QR iteration, computing right 03184 * singular vectors of A in VT 03185 * (CWorkspace: 0) 03186 * (RWorkspace: need BDSPAC) 03187 * 03188 CALL CBDSQR( 'U', M, N, 0, 0, S, RWORK( IE ), VT, 03189 $ LDVT, CDUM, 1, CDUM, 1, 03190 $ RWORK( IRWORK ), INFO ) 03191 * 03192 END IF 03193 * 03194 ELSE IF( WNTUO ) THEN 03195 * 03196 * Path 8t(N much larger than M, JOBU='O', JOBVT='A') 03197 * N right singular vectors to be computed in VT and 03198 * M left singular vectors to be overwritten on A 03199 * 03200 IF( LWORK.GE.2*M*M+MAX( N+M, 3*M ) ) THEN 03201 * 03202 * Sufficient workspace for a fast algorithm 03203 * 03204 IU = 1 03205 IF( LWORK.GE.WRKBL+2*LDA*M ) THEN 03206 * 03207 * WORK(IU) is LDA by M and WORK(IR) is LDA by M 03208 * 03209 LDWRKU = LDA 03210 IR = IU + LDWRKU*M 03211 LDWRKR = LDA 03212 ELSE IF( LWORK.GE.WRKBL+( LDA+M )*M ) THEN 03213 * 03214 * WORK(IU) is LDA by M and WORK(IR) is M by M 03215 * 03216 LDWRKU = LDA 03217 IR = IU + LDWRKU*M 03218 LDWRKR = M 03219 ELSE 03220 * 03221 * WORK(IU) is M by M and WORK(IR) is M by M 03222 * 03223 LDWRKU = M 03224 IR = IU + LDWRKU*M 03225 LDWRKR = M 03226 END IF 03227 ITAU = IR + LDWRKR*M 03228 IWORK = ITAU + M 03229 * 03230 * Compute A=L*Q, copying result to VT 03231 * (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB) 03232 * (RWorkspace: 0) 03233 * 03234 CALL CGELQF( M, N, A, LDA, WORK( ITAU ), 03235 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03236 CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT ) 03237 * 03238 * Generate Q in VT 03239 * (CWorkspace: need 2*M*M+M+N, prefer 2*M*M+M+N*NB) 03240 * (RWorkspace: 0) 03241 * 03242 CALL CUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ), 03243 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03244 * 03245 * Copy L to WORK(IU), zeroing out above it 03246 * 03247 CALL CLACPY( 'L', M, M, A, LDA, WORK( IU ), 03248 $ LDWRKU ) 03249 CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, 03250 $ WORK( IU+LDWRKU ), LDWRKU ) 03251 IE = 1 03252 ITAUQ = ITAU 03253 ITAUP = ITAUQ + M 03254 IWORK = ITAUP + M 03255 * 03256 * Bidiagonalize L in WORK(IU), copying result to 03257 * WORK(IR) 03258 * (CWorkspace: need 2*M*M+3*M, 03259 * prefer 2*M*M+2*M+2*M*NB) 03260 * (RWorkspace: need M) 03261 * 03262 CALL CGEBRD( M, M, WORK( IU ), LDWRKU, S, 03263 $ RWORK( IE ), WORK( ITAUQ ), 03264 $ WORK( ITAUP ), WORK( IWORK ), 03265 $ LWORK-IWORK+1, IERR ) 03266 CALL CLACPY( 'L', M, M, WORK( IU ), LDWRKU, 03267 $ WORK( IR ), LDWRKR ) 03268 * 03269 * Generate right bidiagonalizing vectors in WORK(IU) 03270 * (CWorkspace: need 2*M*M+3*M-1, 03271 * prefer 2*M*M+2*M+(M-1)*NB) 03272 * (RWorkspace: 0) 03273 * 03274 CALL CUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU, 03275 $ WORK( ITAUP ), WORK( IWORK ), 03276 $ LWORK-IWORK+1, IERR ) 03277 * 03278 * Generate left bidiagonalizing vectors in WORK(IR) 03279 * (CWorkspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB) 03280 * (RWorkspace: 0) 03281 * 03282 CALL CUNGBR( 'Q', M, M, M, WORK( IR ), LDWRKR, 03283 $ WORK( ITAUQ ), WORK( IWORK ), 03284 $ LWORK-IWORK+1, IERR ) 03285 IRWORK = IE + M 03286 * 03287 * Perform bidiagonal QR iteration, computing left 03288 * singular vectors of L in WORK(IR) and computing 03289 * right singular vectors of L in WORK(IU) 03290 * (CWorkspace: need 2*M*M) 03291 * (RWorkspace: need BDSPAC) 03292 * 03293 CALL CBDSQR( 'U', M, M, M, 0, S, RWORK( IE ), 03294 $ WORK( IU ), LDWRKU, WORK( IR ), 03295 $ LDWRKR, CDUM, 1, RWORK( IRWORK ), 03296 $ INFO ) 03297 * 03298 * Multiply right singular vectors of L in WORK(IU) by 03299 * Q in VT, storing result in A 03300 * (CWorkspace: need M*M) 03301 * (RWorkspace: 0) 03302 * 03303 CALL CGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ), 03304 $ LDWRKU, VT, LDVT, CZERO, A, LDA ) 03305 * 03306 * Copy right singular vectors of A from A to VT 03307 * 03308 CALL CLACPY( 'F', M, N, A, LDA, VT, LDVT ) 03309 * 03310 * Copy left singular vectors of A from WORK(IR) to A 03311 * 03312 CALL CLACPY( 'F', M, M, WORK( IR ), LDWRKR, A, 03313 $ LDA ) 03314 * 03315 ELSE 03316 * 03317 * Insufficient workspace for a fast algorithm 03318 * 03319 ITAU = 1 03320 IWORK = ITAU + M 03321 * 03322 * Compute A=L*Q, copying result to VT 03323 * (CWorkspace: need 2*M, prefer M+M*NB) 03324 * (RWorkspace: 0) 03325 * 03326 CALL CGELQF( M, N, A, LDA, WORK( ITAU ), 03327 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03328 CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT ) 03329 * 03330 * Generate Q in VT 03331 * (CWorkspace: need M+N, prefer M+N*NB) 03332 * (RWorkspace: 0) 03333 * 03334 CALL CUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ), 03335 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03336 IE = 1 03337 ITAUQ = ITAU 03338 ITAUP = ITAUQ + M 03339 IWORK = ITAUP + M 03340 * 03341 * Zero out above L in A 03342 * 03343 CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, 03344 $ A( 1, 2 ), LDA ) 03345 * 03346 * Bidiagonalize L in A 03347 * (CWorkspace: need 3*M, prefer 2*M+2*M*NB) 03348 * (RWorkspace: need M) 03349 * 03350 CALL CGEBRD( M, M, A, LDA, S, RWORK( IE ), 03351 $ WORK( ITAUQ ), WORK( ITAUP ), 03352 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03353 * 03354 * Multiply right bidiagonalizing vectors in A by Q 03355 * in VT 03356 * (CWorkspace: need 2*M+N, prefer 2*M+N*NB) 03357 * (RWorkspace: 0) 03358 * 03359 CALL CUNMBR( 'P', 'L', 'C', M, N, M, A, LDA, 03360 $ WORK( ITAUP ), VT, LDVT, 03361 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03362 * 03363 * Generate left bidiagonalizing vectors in A 03364 * (CWorkspace: need 3*M, prefer 2*M+M*NB) 03365 * (RWorkspace: 0) 03366 * 03367 CALL CUNGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ), 03368 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03369 IRWORK = IE + M 03370 * 03371 * Perform bidiagonal QR iteration, computing left 03372 * singular vectors of A in A and computing right 03373 * singular vectors of A in VT 03374 * (CWorkspace: 0) 03375 * (RWorkspace: need BDSPAC) 03376 * 03377 CALL CBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT, 03378 $ LDVT, A, LDA, CDUM, 1, 03379 $ RWORK( IRWORK ), INFO ) 03380 * 03381 END IF 03382 * 03383 ELSE IF( WNTUAS ) THEN 03384 * 03385 * Path 9t(N much larger than M, JOBU='S' or 'A', 03386 * JOBVT='A') 03387 * N right singular vectors to be computed in VT and 03388 * M left singular vectors to be computed in U 03389 * 03390 IF( LWORK.GE.M*M+MAX( N+M, 3*M ) ) THEN 03391 * 03392 * Sufficient workspace for a fast algorithm 03393 * 03394 IU = 1 03395 IF( LWORK.GE.WRKBL+LDA*M ) THEN 03396 * 03397 * WORK(IU) is LDA by M 03398 * 03399 LDWRKU = LDA 03400 ELSE 03401 * 03402 * WORK(IU) is M by M 03403 * 03404 LDWRKU = M 03405 END IF 03406 ITAU = IU + LDWRKU*M 03407 IWORK = ITAU + M 03408 * 03409 * Compute A=L*Q, copying result to VT 03410 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 03411 * (RWorkspace: 0) 03412 * 03413 CALL CGELQF( M, N, A, LDA, WORK( ITAU ), 03414 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03415 CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT ) 03416 * 03417 * Generate Q in VT 03418 * (CWorkspace: need M*M+M+N, prefer M*M+M+N*NB) 03419 * (RWorkspace: 0) 03420 * 03421 CALL CUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ), 03422 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03423 * 03424 * Copy L to WORK(IU), zeroing out above it 03425 * 03426 CALL CLACPY( 'L', M, M, A, LDA, WORK( IU ), 03427 $ LDWRKU ) 03428 CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, 03429 $ WORK( IU+LDWRKU ), LDWRKU ) 03430 IE = 1 03431 ITAUQ = ITAU 03432 ITAUP = ITAUQ + M 03433 IWORK = ITAUP + M 03434 * 03435 * Bidiagonalize L in WORK(IU), copying result to U 03436 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) 03437 * (RWorkspace: need M) 03438 * 03439 CALL CGEBRD( M, M, WORK( IU ), LDWRKU, S, 03440 $ RWORK( IE ), WORK( ITAUQ ), 03441 $ WORK( ITAUP ), WORK( IWORK ), 03442 $ LWORK-IWORK+1, IERR ) 03443 CALL CLACPY( 'L', M, M, WORK( IU ), LDWRKU, U, 03444 $ LDU ) 03445 * 03446 * Generate right bidiagonalizing vectors in WORK(IU) 03447 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+(M-1)*NB) 03448 * (RWorkspace: 0) 03449 * 03450 CALL CUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU, 03451 $ WORK( ITAUP ), WORK( IWORK ), 03452 $ LWORK-IWORK+1, IERR ) 03453 * 03454 * Generate left bidiagonalizing vectors in U 03455 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB) 03456 * (RWorkspace: 0) 03457 * 03458 CALL CUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), 03459 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03460 IRWORK = IE + M 03461 * 03462 * Perform bidiagonal QR iteration, computing left 03463 * singular vectors of L in U and computing right 03464 * singular vectors of L in WORK(IU) 03465 * (CWorkspace: need M*M) 03466 * (RWorkspace: need BDSPAC) 03467 * 03468 CALL CBDSQR( 'U', M, M, M, 0, S, RWORK( IE ), 03469 $ WORK( IU ), LDWRKU, U, LDU, CDUM, 1, 03470 $ RWORK( IRWORK ), INFO ) 03471 * 03472 * Multiply right singular vectors of L in WORK(IU) by 03473 * Q in VT, storing result in A 03474 * (CWorkspace: need M*M) 03475 * (RWorkspace: 0) 03476 * 03477 CALL CGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ), 03478 $ LDWRKU, VT, LDVT, CZERO, A, LDA ) 03479 * 03480 * Copy right singular vectors of A from A to VT 03481 * 03482 CALL CLACPY( 'F', M, N, A, LDA, VT, LDVT ) 03483 * 03484 ELSE 03485 * 03486 * Insufficient workspace for a fast algorithm 03487 * 03488 ITAU = 1 03489 IWORK = ITAU + M 03490 * 03491 * Compute A=L*Q, copying result to VT 03492 * (CWorkspace: need 2*M, prefer M+M*NB) 03493 * (RWorkspace: 0) 03494 * 03495 CALL CGELQF( M, N, A, LDA, WORK( ITAU ), 03496 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03497 CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT ) 03498 * 03499 * Generate Q in VT 03500 * (CWorkspace: need M+N, prefer M+N*NB) 03501 * (RWorkspace: 0) 03502 * 03503 CALL CUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ), 03504 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03505 * 03506 * Copy L to U, zeroing out above it 03507 * 03508 CALL CLACPY( 'L', M, M, A, LDA, U, LDU ) 03509 CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, 03510 $ U( 1, 2 ), LDU ) 03511 IE = 1 03512 ITAUQ = ITAU 03513 ITAUP = ITAUQ + M 03514 IWORK = ITAUP + M 03515 * 03516 * Bidiagonalize L in U 03517 * (CWorkspace: need 3*M, prefer 2*M+2*M*NB) 03518 * (RWorkspace: need M) 03519 * 03520 CALL CGEBRD( M, M, U, LDU, S, RWORK( IE ), 03521 $ WORK( ITAUQ ), WORK( ITAUP ), 03522 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03523 * 03524 * Multiply right bidiagonalizing vectors in U by Q 03525 * in VT 03526 * (CWorkspace: need 2*M+N, prefer 2*M+N*NB) 03527 * (RWorkspace: 0) 03528 * 03529 CALL CUNMBR( 'P', 'L', 'C', M, N, M, U, LDU, 03530 $ WORK( ITAUP ), VT, LDVT, 03531 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03532 * 03533 * Generate left bidiagonalizing vectors in U 03534 * (CWorkspace: need 3*M, prefer 2*M+M*NB) 03535 * (RWorkspace: 0) 03536 * 03537 CALL CUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), 03538 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03539 IRWORK = IE + M 03540 * 03541 * Perform bidiagonal QR iteration, computing left 03542 * singular vectors of A in U and computing right 03543 * singular vectors of A in VT 03544 * (CWorkspace: 0) 03545 * (RWorkspace: need BDSPAC) 03546 * 03547 CALL CBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT, 03548 $ LDVT, U, LDU, CDUM, 1, 03549 $ RWORK( IRWORK ), INFO ) 03550 * 03551 END IF 03552 * 03553 END IF 03554 * 03555 END IF 03556 * 03557 ELSE 03558 * 03559 * N .LT. MNTHR 03560 * 03561 * Path 10t(N greater than M, but not much larger) 03562 * Reduce to bidiagonal form without LQ decomposition 03563 * 03564 IE = 1 03565 ITAUQ = 1 03566 ITAUP = ITAUQ + M 03567 IWORK = ITAUP + M 03568 * 03569 * Bidiagonalize A 03570 * (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB) 03571 * (RWorkspace: M) 03572 * 03573 CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 03574 $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, 03575 $ IERR ) 03576 IF( WNTUAS ) THEN 03577 * 03578 * If left singular vectors desired in U, copy result to U 03579 * and generate left bidiagonalizing vectors in U 03580 * (CWorkspace: need 3*M-1, prefer 2*M+(M-1)*NB) 03581 * (RWorkspace: 0) 03582 * 03583 CALL CLACPY( 'L', M, M, A, LDA, U, LDU ) 03584 CALL CUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ), 03585 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03586 END IF 03587 IF( WNTVAS ) THEN 03588 * 03589 * If right singular vectors desired in VT, copy result to 03590 * VT and generate right bidiagonalizing vectors in VT 03591 * (CWorkspace: need 2*M+NRVT, prefer 2*M+NRVT*NB) 03592 * (RWorkspace: 0) 03593 * 03594 CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT ) 03595 IF( WNTVA ) 03596 $ NRVT = N 03597 IF( WNTVS ) 03598 $ NRVT = M 03599 CALL CUNGBR( 'P', NRVT, N, M, VT, LDVT, WORK( ITAUP ), 03600 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03601 END IF 03602 IF( WNTUO ) THEN 03603 * 03604 * If left singular vectors desired in A, generate left 03605 * bidiagonalizing vectors in A 03606 * (CWorkspace: need 3*M-1, prefer 2*M+(M-1)*NB) 03607 * (RWorkspace: 0) 03608 * 03609 CALL CUNGBR( 'Q', M, M, N, A, LDA, WORK( ITAUQ ), 03610 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03611 END IF 03612 IF( WNTVO ) THEN 03613 * 03614 * If right singular vectors desired in A, generate right 03615 * bidiagonalizing vectors in A 03616 * (CWorkspace: need 3*M, prefer 2*M+M*NB) 03617 * (RWorkspace: 0) 03618 * 03619 CALL CUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), 03620 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03621 END IF 03622 IRWORK = IE + M 03623 IF( WNTUAS .OR. WNTUO ) 03624 $ NRU = M 03625 IF( WNTUN ) 03626 $ NRU = 0 03627 IF( WNTVAS .OR. WNTVO ) 03628 $ NCVT = N 03629 IF( WNTVN ) 03630 $ NCVT = 0 03631 IF( ( .NOT.WNTUO ) .AND. ( .NOT.WNTVO ) ) THEN 03632 * 03633 * Perform bidiagonal QR iteration, if desired, computing 03634 * left singular vectors in U and computing right singular 03635 * vectors in VT 03636 * (CWorkspace: 0) 03637 * (RWorkspace: need BDSPAC) 03638 * 03639 CALL CBDSQR( 'L', M, NCVT, NRU, 0, S, RWORK( IE ), VT, 03640 $ LDVT, U, LDU, CDUM, 1, RWORK( IRWORK ), 03641 $ INFO ) 03642 ELSE IF( ( .NOT.WNTUO ) .AND. WNTVO ) THEN 03643 * 03644 * Perform bidiagonal QR iteration, if desired, computing 03645 * left singular vectors in U and computing right singular 03646 * vectors in A 03647 * (CWorkspace: 0) 03648 * (RWorkspace: need BDSPAC) 03649 * 03650 CALL CBDSQR( 'L', M, NCVT, NRU, 0, S, RWORK( IE ), A, 03651 $ LDA, U, LDU, CDUM, 1, RWORK( IRWORK ), 03652 $ INFO ) 03653 ELSE 03654 * 03655 * Perform bidiagonal QR iteration, if desired, computing 03656 * left singular vectors in A and computing right singular 03657 * vectors in VT 03658 * (CWorkspace: 0) 03659 * (RWorkspace: need BDSPAC) 03660 * 03661 CALL CBDSQR( 'L', M, NCVT, NRU, 0, S, RWORK( IE ), VT, 03662 $ LDVT, A, LDA, CDUM, 1, RWORK( IRWORK ), 03663 $ INFO ) 03664 END IF 03665 * 03666 END IF 03667 * 03668 END IF 03669 * 03670 * Undo scaling if necessary 03671 * 03672 IF( ISCL.EQ.1 ) THEN 03673 IF( ANRM.GT.BIGNUM ) 03674 $ CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN, 03675 $ IERR ) 03676 IF( INFO.NE.0 .AND. ANRM.GT.BIGNUM ) 03677 $ CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN-1, 1, 03678 $ RWORK( IE ), MINMN, IERR ) 03679 IF( ANRM.LT.SMLNUM ) 03680 $ CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN, 03681 $ IERR ) 03682 IF( INFO.NE.0 .AND. ANRM.LT.SMLNUM ) 03683 $ CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN-1, 1, 03684 $ RWORK( IE ), MINMN, IERR ) 03685 END IF 03686 * 03687 * Return optimal workspace in WORK(1) 03688 * 03689 WORK( 1 ) = MAXWRK 03690 * 03691 RETURN 03692 * 03693 * End of CGESVD 03694 * 03695 END