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