![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DGESDD 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DGESDD + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesdd.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesdd.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesdd.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, 00022 * LWORK, IWORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER JOBZ 00026 * INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N 00027 * .. 00028 * .. Array Arguments .. 00029 * INTEGER IWORK( * ) 00030 * DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ), 00031 * $ VT( LDVT, * ), WORK( * ) 00032 * .. 00033 * 00034 * 00035 *> \par Purpose: 00036 * ============= 00037 *> 00038 *> \verbatim 00039 *> 00040 *> DGESDD computes the singular value decomposition (SVD) of a real 00041 *> M-by-N matrix A, optionally computing the left and right singular 00042 *> vectors. If singular vectors are desired, it uses a 00043 *> divide-and-conquer algorithm. 00044 *> 00045 *> The SVD is written 00046 *> 00047 *> A = U * SIGMA * transpose(V) 00048 *> 00049 *> where SIGMA is an M-by-N matrix which is zero except for its 00050 *> min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and 00051 *> V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA 00052 *> are the singular values of A; they are real and non-negative, and 00053 *> are returned in descending order. The first min(m,n) columns of 00054 *> U and V are the left and right singular vectors of A. 00055 *> 00056 *> Note that the routine returns VT = V**T, not V. 00057 *> 00058 *> The divide and conquer algorithm makes very mild assumptions about 00059 *> floating point arithmetic. It will work on machines with a guard 00060 *> digit in add/subtract, or on those binary machines without guard 00061 *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or 00062 *> Cray-2. It could conceivably fail on hexadecimal or decimal machines 00063 *> without guard digits, but we know of none. 00064 *> \endverbatim 00065 * 00066 * Arguments: 00067 * ========== 00068 * 00069 *> \param[in] JOBZ 00070 *> \verbatim 00071 *> JOBZ is CHARACTER*1 00072 *> Specifies options for computing all or part of the matrix U: 00073 *> = 'A': all M columns of U and all N rows of V**T are 00074 *> returned in the arrays U and VT; 00075 *> = 'S': the first min(M,N) columns of U and the first 00076 *> min(M,N) rows of V**T are returned in the arrays U 00077 *> and VT; 00078 *> = 'O': If M >= N, the first N columns of U are overwritten 00079 *> on the array A and all rows of V**T are returned in 00080 *> the array VT; 00081 *> otherwise, all columns of U are returned in the 00082 *> array U and the first M rows of V**T are overwritten 00083 *> in the array A; 00084 *> = 'N': no columns of U or rows of V**T are computed. 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 JOBZ = 'O', A is overwritten with the first N columns 00105 *> of U (the left singular vectors, stored 00106 *> columnwise) if M >= N; 00107 *> A is overwritten with the first M rows 00108 *> of V**T (the right singular vectors, stored 00109 *> rowwise) otherwise. 00110 *> if JOBZ .ne. 'O', the contents of A are destroyed. 00111 *> \endverbatim 00112 *> 00113 *> \param[in] LDA 00114 *> \verbatim 00115 *> LDA is INTEGER 00116 *> The leading dimension of the array A. LDA >= max(1,M). 00117 *> \endverbatim 00118 *> 00119 *> \param[out] S 00120 *> \verbatim 00121 *> S is DOUBLE PRECISION array, dimension (min(M,N)) 00122 *> The singular values of A, sorted so that S(i) >= S(i+1). 00123 *> \endverbatim 00124 *> 00125 *> \param[out] U 00126 *> \verbatim 00127 *> U is DOUBLE PRECISION array, dimension (LDU,UCOL) 00128 *> UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N; 00129 *> UCOL = min(M,N) if JOBZ = 'S'. 00130 *> If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M 00131 *> orthogonal matrix U; 00132 *> if JOBZ = 'S', U contains the first min(M,N) columns of U 00133 *> (the left singular vectors, stored columnwise); 00134 *> if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced. 00135 *> \endverbatim 00136 *> 00137 *> \param[in] LDU 00138 *> \verbatim 00139 *> LDU is INTEGER 00140 *> The leading dimension of the array U. LDU >= 1; if 00141 *> JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M. 00142 *> \endverbatim 00143 *> 00144 *> \param[out] VT 00145 *> \verbatim 00146 *> VT is DOUBLE PRECISION array, dimension (LDVT,N) 00147 *> If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the 00148 *> N-by-N orthogonal matrix V**T; 00149 *> if JOBZ = 'S', VT contains the first min(M,N) rows of 00150 *> V**T (the right singular vectors, stored rowwise); 00151 *> if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced. 00152 *> \endverbatim 00153 *> 00154 *> \param[in] LDVT 00155 *> \verbatim 00156 *> LDVT is INTEGER 00157 *> The leading dimension of the array VT. LDVT >= 1; if 00158 *> JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N; 00159 *> if JOBZ = 'S', LDVT >= min(M,N). 00160 *> \endverbatim 00161 *> 00162 *> \param[out] WORK 00163 *> \verbatim 00164 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 00165 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK; 00166 *> \endverbatim 00167 *> 00168 *> \param[in] LWORK 00169 *> \verbatim 00170 *> LWORK is INTEGER 00171 *> The dimension of the array WORK. LWORK >= 1. 00172 *> If JOBZ = 'N', 00173 *> LWORK >= 3*min(M,N) + max(max(M,N),7*min(M,N)). 00174 *> If JOBZ = 'O', 00175 *> LWORK >= 3*min(M,N) + 00176 *> max(max(M,N),5*min(M,N)*min(M,N)+4*min(M,N)). 00177 *> If JOBZ = 'S' or 'A' 00178 *> LWORK >= 3*min(M,N) + 00179 *> max(max(M,N),4*min(M,N)*min(M,N)+4*min(M,N)). 00180 *> For good performance, LWORK should generally be larger. 00181 *> If LWORK = -1 but other input arguments are legal, WORK(1) 00182 *> returns the optimal LWORK. 00183 *> \endverbatim 00184 *> 00185 *> \param[out] IWORK 00186 *> \verbatim 00187 *> IWORK is INTEGER array, dimension (8*min(M,N)) 00188 *> \endverbatim 00189 *> 00190 *> \param[out] INFO 00191 *> \verbatim 00192 *> INFO is INTEGER 00193 *> = 0: successful exit. 00194 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00195 *> > 0: DBDSDC did not converge, updating process failed. 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 November 2011 00207 * 00208 *> \ingroup doubleGEsing 00209 * 00210 *> \par Contributors: 00211 * ================== 00212 *> 00213 *> Ming Gu and Huan Ren, Computer Science Division, University of 00214 *> California at Berkeley, USA 00215 *> 00216 * ===================================================================== 00217 SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, 00218 $ LWORK, IWORK, INFO ) 00219 * 00220 * -- LAPACK driver routine (version 3.4.0) -- 00221 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00222 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00223 * November 2011 00224 * 00225 * .. Scalar Arguments .. 00226 CHARACTER JOBZ 00227 INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N 00228 * .. 00229 * .. Array Arguments .. 00230 INTEGER IWORK( * ) 00231 DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ), 00232 $ VT( LDVT, * ), WORK( * ) 00233 * .. 00234 * 00235 * ===================================================================== 00236 * 00237 * .. Parameters .. 00238 DOUBLE PRECISION ZERO, ONE 00239 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00240 * .. 00241 * .. Local Scalars .. 00242 LOGICAL LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS 00243 INTEGER BDSPAC, BLK, CHUNK, I, IE, IERR, IL, 00244 $ IR, ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT, 00245 $ LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK, 00246 $ MNTHR, NWORK, WRKBL 00247 DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM 00248 * .. 00249 * .. Local Arrays .. 00250 INTEGER IDUM( 1 ) 00251 DOUBLE PRECISION DUM( 1 ) 00252 * .. 00253 * .. External Subroutines .. 00254 EXTERNAL DBDSDC, DGEBRD, DGELQF, DGEMM, DGEQRF, DLACPY, 00255 $ DLASCL, DLASET, DORGBR, DORGLQ, DORGQR, DORMBR, 00256 $ XERBLA 00257 * .. 00258 * .. External Functions .. 00259 LOGICAL LSAME 00260 INTEGER ILAENV 00261 DOUBLE PRECISION DLAMCH, DLANGE 00262 EXTERNAL DLAMCH, DLANGE, ILAENV, LSAME 00263 * .. 00264 * .. Intrinsic Functions .. 00265 INTRINSIC INT, MAX, MIN, SQRT 00266 * .. 00267 * .. Executable Statements .. 00268 * 00269 * Test the input arguments 00270 * 00271 INFO = 0 00272 MINMN = MIN( M, N ) 00273 WNTQA = LSAME( JOBZ, 'A' ) 00274 WNTQS = LSAME( JOBZ, 'S' ) 00275 WNTQAS = WNTQA .OR. WNTQS 00276 WNTQO = LSAME( JOBZ, 'O' ) 00277 WNTQN = LSAME( JOBZ, 'N' ) 00278 LQUERY = ( LWORK.EQ.-1 ) 00279 * 00280 IF( .NOT.( WNTQA .OR. WNTQS .OR. WNTQO .OR. WNTQN ) ) THEN 00281 INFO = -1 00282 ELSE IF( M.LT.0 ) THEN 00283 INFO = -2 00284 ELSE IF( N.LT.0 ) THEN 00285 INFO = -3 00286 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00287 INFO = -5 00288 ELSE IF( LDU.LT.1 .OR. ( WNTQAS .AND. LDU.LT.M ) .OR. 00289 $ ( WNTQO .AND. M.LT.N .AND. LDU.LT.M ) ) THEN 00290 INFO = -8 00291 ELSE IF( LDVT.LT.1 .OR. ( WNTQA .AND. LDVT.LT.N ) .OR. 00292 $ ( WNTQS .AND. LDVT.LT.MINMN ) .OR. 00293 $ ( WNTQO .AND. M.GE.N .AND. LDVT.LT.N ) ) THEN 00294 INFO = -10 00295 END IF 00296 * 00297 * Compute workspace 00298 * (Note: Comments in the code beginning "Workspace:" describe the 00299 * minimal amount of workspace needed at that point in the code, 00300 * as well as the preferred amount for good performance. 00301 * NB refers to the optimal block size for the immediately 00302 * following subroutine, as returned by ILAENV.) 00303 * 00304 IF( INFO.EQ.0 ) THEN 00305 MINWRK = 1 00306 MAXWRK = 1 00307 IF( M.GE.N .AND. MINMN.GT.0 ) THEN 00308 * 00309 * Compute space needed for DBDSDC 00310 * 00311 MNTHR = INT( MINMN*11.0D0 / 6.0D0 ) 00312 IF( WNTQN ) THEN 00313 BDSPAC = 7*N 00314 ELSE 00315 BDSPAC = 3*N*N + 4*N 00316 END IF 00317 IF( M.GE.MNTHR ) THEN 00318 IF( WNTQN ) THEN 00319 * 00320 * Path 1 (M much larger than N, JOBZ='N') 00321 * 00322 WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, 00323 $ -1 ) 00324 WRKBL = MAX( WRKBL, 3*N+2*N* 00325 $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) 00326 MAXWRK = MAX( WRKBL, BDSPAC+N ) 00327 MINWRK = BDSPAC + N 00328 ELSE IF( WNTQO ) THEN 00329 * 00330 * Path 2 (M much larger than N, JOBZ='O') 00331 * 00332 WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) 00333 WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M, 00334 $ N, N, -1 ) ) 00335 WRKBL = MAX( WRKBL, 3*N+2*N* 00336 $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) 00337 WRKBL = MAX( WRKBL, 3*N+N* 00338 $ ILAENV( 1, 'DORMBR', 'QLN', N, N, N, -1 ) ) 00339 WRKBL = MAX( WRKBL, 3*N+N* 00340 $ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) ) 00341 WRKBL = MAX( WRKBL, BDSPAC+3*N ) 00342 MAXWRK = WRKBL + 2*N*N 00343 MINWRK = BDSPAC + 2*N*N + 3*N 00344 ELSE IF( WNTQS ) THEN 00345 * 00346 * Path 3 (M much larger than N, JOBZ='S') 00347 * 00348 WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) 00349 WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M, 00350 $ N, N, -1 ) ) 00351 WRKBL = MAX( WRKBL, 3*N+2*N* 00352 $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) 00353 WRKBL = MAX( WRKBL, 3*N+N* 00354 $ ILAENV( 1, 'DORMBR', 'QLN', N, N, N, -1 ) ) 00355 WRKBL = MAX( WRKBL, 3*N+N* 00356 $ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) ) 00357 WRKBL = MAX( WRKBL, BDSPAC+3*N ) 00358 MAXWRK = WRKBL + N*N 00359 MINWRK = BDSPAC + N*N + 3*N 00360 ELSE IF( WNTQA ) THEN 00361 * 00362 * Path 4 (M much larger than N, JOBZ='A') 00363 * 00364 WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) 00365 WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'DORGQR', ' ', M, 00366 $ M, N, -1 ) ) 00367 WRKBL = MAX( WRKBL, 3*N+2*N* 00368 $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) 00369 WRKBL = MAX( WRKBL, 3*N+N* 00370 $ ILAENV( 1, 'DORMBR', 'QLN', N, N, N, -1 ) ) 00371 WRKBL = MAX( WRKBL, 3*N+N* 00372 $ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) ) 00373 WRKBL = MAX( WRKBL, BDSPAC+3*N ) 00374 MAXWRK = WRKBL + N*N 00375 MINWRK = BDSPAC + N*N + 3*N 00376 END IF 00377 ELSE 00378 * 00379 * Path 5 (M at least N, but not much larger) 00380 * 00381 WRKBL = 3*N + ( M+N )*ILAENV( 1, 'DGEBRD', ' ', M, N, -1, 00382 $ -1 ) 00383 IF( WNTQN ) THEN 00384 MAXWRK = MAX( WRKBL, BDSPAC+3*N ) 00385 MINWRK = 3*N + MAX( M, BDSPAC ) 00386 ELSE IF( WNTQO ) THEN 00387 WRKBL = MAX( WRKBL, 3*N+N* 00388 $ ILAENV( 1, 'DORMBR', 'QLN', M, N, N, -1 ) ) 00389 WRKBL = MAX( WRKBL, 3*N+N* 00390 $ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) ) 00391 WRKBL = MAX( WRKBL, BDSPAC+3*N ) 00392 MAXWRK = WRKBL + M*N 00393 MINWRK = 3*N + MAX( M, N*N+BDSPAC ) 00394 ELSE IF( WNTQS ) THEN 00395 WRKBL = MAX( WRKBL, 3*N+N* 00396 $ ILAENV( 1, 'DORMBR', 'QLN', M, N, N, -1 ) ) 00397 WRKBL = MAX( WRKBL, 3*N+N* 00398 $ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) ) 00399 MAXWRK = MAX( WRKBL, BDSPAC+3*N ) 00400 MINWRK = 3*N + MAX( M, BDSPAC ) 00401 ELSE IF( WNTQA ) THEN 00402 WRKBL = MAX( WRKBL, 3*N+M* 00403 $ ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) ) 00404 WRKBL = MAX( WRKBL, 3*N+N* 00405 $ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) ) 00406 MAXWRK = MAX( MAXWRK, BDSPAC+3*N ) 00407 MINWRK = 3*N + MAX( M, BDSPAC ) 00408 END IF 00409 END IF 00410 ELSE IF( MINMN.GT.0 ) THEN 00411 * 00412 * Compute space needed for DBDSDC 00413 * 00414 MNTHR = INT( MINMN*11.0D0 / 6.0D0 ) 00415 IF( WNTQN ) THEN 00416 BDSPAC = 7*M 00417 ELSE 00418 BDSPAC = 3*M*M + 4*M 00419 END IF 00420 IF( N.GE.MNTHR ) THEN 00421 IF( WNTQN ) THEN 00422 * 00423 * Path 1t (N much larger than M, JOBZ='N') 00424 * 00425 WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, 00426 $ -1 ) 00427 WRKBL = MAX( WRKBL, 3*M+2*M* 00428 $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) 00429 MAXWRK = MAX( WRKBL, BDSPAC+M ) 00430 MINWRK = BDSPAC + M 00431 ELSE IF( WNTQO ) THEN 00432 * 00433 * Path 2t (N much larger than M, JOBZ='O') 00434 * 00435 WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) 00436 WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M, 00437 $ N, M, -1 ) ) 00438 WRKBL = MAX( WRKBL, 3*M+2*M* 00439 $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) 00440 WRKBL = MAX( WRKBL, 3*M+M* 00441 $ ILAENV( 1, 'DORMBR', 'QLN', M, M, M, -1 ) ) 00442 WRKBL = MAX( WRKBL, 3*M+M* 00443 $ ILAENV( 1, 'DORMBR', 'PRT', M, M, M, -1 ) ) 00444 WRKBL = MAX( WRKBL, BDSPAC+3*M ) 00445 MAXWRK = WRKBL + 2*M*M 00446 MINWRK = BDSPAC + 2*M*M + 3*M 00447 ELSE IF( WNTQS ) THEN 00448 * 00449 * Path 3t (N much larger than M, JOBZ='S') 00450 * 00451 WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) 00452 WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M, 00453 $ N, M, -1 ) ) 00454 WRKBL = MAX( WRKBL, 3*M+2*M* 00455 $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) 00456 WRKBL = MAX( WRKBL, 3*M+M* 00457 $ ILAENV( 1, 'DORMBR', 'QLN', M, M, M, -1 ) ) 00458 WRKBL = MAX( WRKBL, 3*M+M* 00459 $ ILAENV( 1, 'DORMBR', 'PRT', M, M, M, -1 ) ) 00460 WRKBL = MAX( WRKBL, BDSPAC+3*M ) 00461 MAXWRK = WRKBL + M*M 00462 MINWRK = BDSPAC + M*M + 3*M 00463 ELSE IF( WNTQA ) THEN 00464 * 00465 * Path 4t (N much larger than M, JOBZ='A') 00466 * 00467 WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) 00468 WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'DORGLQ', ' ', N, 00469 $ N, M, -1 ) ) 00470 WRKBL = MAX( WRKBL, 3*M+2*M* 00471 $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) 00472 WRKBL = MAX( WRKBL, 3*M+M* 00473 $ ILAENV( 1, 'DORMBR', 'QLN', M, M, M, -1 ) ) 00474 WRKBL = MAX( WRKBL, 3*M+M* 00475 $ ILAENV( 1, 'DORMBR', 'PRT', M, M, M, -1 ) ) 00476 WRKBL = MAX( WRKBL, BDSPAC+3*M ) 00477 MAXWRK = WRKBL + M*M 00478 MINWRK = BDSPAC + M*M + 3*M 00479 END IF 00480 ELSE 00481 * 00482 * Path 5t (N greater than M, but not much larger) 00483 * 00484 WRKBL = 3*M + ( M+N )*ILAENV( 1, 'DGEBRD', ' ', M, N, -1, 00485 $ -1 ) 00486 IF( WNTQN ) THEN 00487 MAXWRK = MAX( WRKBL, BDSPAC+3*M ) 00488 MINWRK = 3*M + MAX( N, BDSPAC ) 00489 ELSE IF( WNTQO ) THEN 00490 WRKBL = MAX( WRKBL, 3*M+M* 00491 $ ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) ) 00492 WRKBL = MAX( WRKBL, 3*M+M* 00493 $ ILAENV( 1, 'DORMBR', 'PRT', M, N, M, -1 ) ) 00494 WRKBL = MAX( WRKBL, BDSPAC+3*M ) 00495 MAXWRK = WRKBL + M*N 00496 MINWRK = 3*M + MAX( N, M*M+BDSPAC ) 00497 ELSE IF( WNTQS ) THEN 00498 WRKBL = MAX( WRKBL, 3*M+M* 00499 $ ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) ) 00500 WRKBL = MAX( WRKBL, 3*M+M* 00501 $ ILAENV( 1, 'DORMBR', 'PRT', M, N, M, -1 ) ) 00502 MAXWRK = MAX( WRKBL, BDSPAC+3*M ) 00503 MINWRK = 3*M + MAX( N, BDSPAC ) 00504 ELSE IF( WNTQA ) THEN 00505 WRKBL = MAX( WRKBL, 3*M+M* 00506 $ ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) ) 00507 WRKBL = MAX( WRKBL, 3*M+M* 00508 $ ILAENV( 1, 'DORMBR', 'PRT', N, N, M, -1 ) ) 00509 MAXWRK = MAX( WRKBL, BDSPAC+3*M ) 00510 MINWRK = 3*M + MAX( N, BDSPAC ) 00511 END IF 00512 END IF 00513 END IF 00514 MAXWRK = MAX( MAXWRK, MINWRK ) 00515 WORK( 1 ) = MAXWRK 00516 * 00517 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN 00518 INFO = -12 00519 END IF 00520 END IF 00521 * 00522 IF( INFO.NE.0 ) THEN 00523 CALL XERBLA( 'DGESDD', -INFO ) 00524 RETURN 00525 ELSE IF( LQUERY ) THEN 00526 RETURN 00527 END IF 00528 * 00529 * Quick return if possible 00530 * 00531 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 00532 RETURN 00533 END IF 00534 * 00535 * Get machine constants 00536 * 00537 EPS = DLAMCH( 'P' ) 00538 SMLNUM = SQRT( DLAMCH( 'S' ) ) / EPS 00539 BIGNUM = ONE / SMLNUM 00540 * 00541 * Scale A if max element outside range [SMLNUM,BIGNUM] 00542 * 00543 ANRM = DLANGE( 'M', M, N, A, LDA, DUM ) 00544 ISCL = 0 00545 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00546 ISCL = 1 00547 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR ) 00548 ELSE IF( ANRM.GT.BIGNUM ) THEN 00549 ISCL = 1 00550 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR ) 00551 END IF 00552 * 00553 IF( M.GE.N ) THEN 00554 * 00555 * A has at least as many rows as columns. If A has sufficiently 00556 * more rows than columns, first reduce using the QR 00557 * decomposition (if sufficient workspace available) 00558 * 00559 IF( M.GE.MNTHR ) THEN 00560 * 00561 IF( WNTQN ) THEN 00562 * 00563 * Path 1 (M much larger than N, JOBZ='N') 00564 * No singular vectors to be computed 00565 * 00566 ITAU = 1 00567 NWORK = ITAU + N 00568 * 00569 * Compute A=Q*R 00570 * (Workspace: need 2*N, prefer N+N*NB) 00571 * 00572 CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 00573 $ LWORK-NWORK+1, IERR ) 00574 * 00575 * Zero out below R 00576 * 00577 CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA ) 00578 IE = 1 00579 ITAUQ = IE + N 00580 ITAUP = ITAUQ + N 00581 NWORK = ITAUP + N 00582 * 00583 * Bidiagonalize R in A 00584 * (Workspace: need 4*N, prefer 3*N+2*N*NB) 00585 * 00586 CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 00587 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 00588 $ IERR ) 00589 NWORK = IE + N 00590 * 00591 * Perform bidiagonal SVD, computing singular values only 00592 * (Workspace: need N+BDSPAC) 00593 * 00594 CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1, 00595 $ DUM, IDUM, WORK( NWORK ), IWORK, INFO ) 00596 * 00597 ELSE IF( WNTQO ) THEN 00598 * 00599 * Path 2 (M much larger than N, JOBZ = 'O') 00600 * N left singular vectors to be overwritten on A and 00601 * N right singular vectors to be computed in VT 00602 * 00603 IR = 1 00604 * 00605 * WORK(IR) is LDWRKR by N 00606 * 00607 IF( LWORK.GE.LDA*N+N*N+3*N+BDSPAC ) THEN 00608 LDWRKR = LDA 00609 ELSE 00610 LDWRKR = ( LWORK-N*N-3*N-BDSPAC ) / N 00611 END IF 00612 ITAU = IR + LDWRKR*N 00613 NWORK = ITAU + N 00614 * 00615 * Compute A=Q*R 00616 * (Workspace: need N*N+2*N, prefer N*N+N+N*NB) 00617 * 00618 CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 00619 $ LWORK-NWORK+1, IERR ) 00620 * 00621 * Copy R to WORK(IR), zeroing out below it 00622 * 00623 CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR ) 00624 CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IR+1 ), 00625 $ LDWRKR ) 00626 * 00627 * Generate Q in A 00628 * (Workspace: need N*N+2*N, prefer N*N+N+N*NB) 00629 * 00630 CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), 00631 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 00632 IE = ITAU 00633 ITAUQ = IE + N 00634 ITAUP = ITAUQ + N 00635 NWORK = ITAUP + N 00636 * 00637 * Bidiagonalize R in VT, copying result to WORK(IR) 00638 * (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) 00639 * 00640 CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ), 00641 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), 00642 $ LWORK-NWORK+1, IERR ) 00643 * 00644 * WORK(IU) is N by N 00645 * 00646 IU = NWORK 00647 NWORK = IU + N*N 00648 * 00649 * Perform bidiagonal SVD, computing left singular vectors 00650 * of bidiagonal matrix in WORK(IU) and computing right 00651 * singular vectors of bidiagonal matrix in VT 00652 * (Workspace: need N+N*N+BDSPAC) 00653 * 00654 CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N, 00655 $ VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK, 00656 $ INFO ) 00657 * 00658 * Overwrite WORK(IU) by left singular vectors of R 00659 * and VT by right singular vectors of R 00660 * (Workspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB) 00661 * 00662 CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR, 00663 $ WORK( ITAUQ ), WORK( IU ), N, WORK( NWORK ), 00664 $ LWORK-NWORK+1, IERR ) 00665 CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR, 00666 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 00667 $ LWORK-NWORK+1, IERR ) 00668 * 00669 * Multiply Q in A by left singular vectors of R in 00670 * WORK(IU), storing result in WORK(IR) and copying to A 00671 * (Workspace: need 2*N*N, prefer N*N+M*N) 00672 * 00673 DO 10 I = 1, M, LDWRKR 00674 CHUNK = MIN( M-I+1, LDWRKR ) 00675 CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ), 00676 $ LDA, WORK( IU ), N, ZERO, WORK( IR ), 00677 $ LDWRKR ) 00678 CALL DLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR, 00679 $ A( I, 1 ), LDA ) 00680 10 CONTINUE 00681 * 00682 ELSE IF( WNTQS ) THEN 00683 * 00684 * Path 3 (M much larger than N, JOBZ='S') 00685 * N left singular vectors to be computed in U and 00686 * N right singular vectors to be computed in VT 00687 * 00688 IR = 1 00689 * 00690 * WORK(IR) is N by N 00691 * 00692 LDWRKR = N 00693 ITAU = IR + LDWRKR*N 00694 NWORK = ITAU + N 00695 * 00696 * Compute A=Q*R 00697 * (Workspace: need N*N+2*N, prefer N*N+N+N*NB) 00698 * 00699 CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 00700 $ LWORK-NWORK+1, IERR ) 00701 * 00702 * Copy R to WORK(IR), zeroing out below it 00703 * 00704 CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR ) 00705 CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IR+1 ), 00706 $ LDWRKR ) 00707 * 00708 * Generate Q in A 00709 * (Workspace: need N*N+2*N, prefer N*N+N+N*NB) 00710 * 00711 CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), 00712 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 00713 IE = ITAU 00714 ITAUQ = IE + N 00715 ITAUP = ITAUQ + N 00716 NWORK = ITAUP + N 00717 * 00718 * Bidiagonalize R in WORK(IR) 00719 * (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) 00720 * 00721 CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ), 00722 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), 00723 $ LWORK-NWORK+1, IERR ) 00724 * 00725 * Perform bidiagonal SVD, computing left singular vectors 00726 * of bidiagoal matrix in U and computing right singular 00727 * vectors of bidiagonal matrix in VT 00728 * (Workspace: need N+BDSPAC) 00729 * 00730 CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT, 00731 $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK, 00732 $ INFO ) 00733 * 00734 * Overwrite U by left singular vectors of R and VT 00735 * by right singular vectors of R 00736 * (Workspace: need N*N+3*N, prefer N*N+2*N+N*NB) 00737 * 00738 CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR, 00739 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 00740 $ LWORK-NWORK+1, IERR ) 00741 * 00742 CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR, 00743 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 00744 $ LWORK-NWORK+1, IERR ) 00745 * 00746 * Multiply Q in A by left singular vectors of R in 00747 * WORK(IR), storing result in U 00748 * (Workspace: need N*N) 00749 * 00750 CALL DLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR ) 00751 CALL DGEMM( 'N', 'N', M, N, N, ONE, A, LDA, WORK( IR ), 00752 $ LDWRKR, ZERO, U, LDU ) 00753 * 00754 ELSE IF( WNTQA ) THEN 00755 * 00756 * Path 4 (M much larger than N, JOBZ='A') 00757 * M left singular vectors to be computed in U and 00758 * N right singular vectors to be computed in VT 00759 * 00760 IU = 1 00761 * 00762 * WORK(IU) is N by N 00763 * 00764 LDWRKU = N 00765 ITAU = IU + LDWRKU*N 00766 NWORK = ITAU + N 00767 * 00768 * Compute A=Q*R, copying result to U 00769 * (Workspace: need N*N+2*N, prefer N*N+N+N*NB) 00770 * 00771 CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 00772 $ LWORK-NWORK+1, IERR ) 00773 CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) 00774 * 00775 * Generate Q in U 00776 * (Workspace: need N*N+2*N, prefer N*N+N+N*NB) 00777 CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ), 00778 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 00779 * 00780 * Produce R in A, zeroing out other entries 00781 * 00782 CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA ) 00783 IE = ITAU 00784 ITAUQ = IE + N 00785 ITAUP = ITAUQ + N 00786 NWORK = ITAUP + N 00787 * 00788 * Bidiagonalize R in A 00789 * (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) 00790 * 00791 CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 00792 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 00793 $ IERR ) 00794 * 00795 * Perform bidiagonal SVD, computing left singular vectors 00796 * of bidiagonal matrix in WORK(IU) and computing right 00797 * singular vectors of bidiagonal matrix in VT 00798 * (Workspace: need N+N*N+BDSPAC) 00799 * 00800 CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N, 00801 $ VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK, 00802 $ INFO ) 00803 * 00804 * Overwrite WORK(IU) by left singular vectors of R and VT 00805 * by right singular vectors of R 00806 * (Workspace: need N*N+3*N, prefer N*N+2*N+N*NB) 00807 * 00808 CALL DORMBR( 'Q', 'L', 'N', N, N, N, A, LDA, 00809 $ WORK( ITAUQ ), WORK( IU ), LDWRKU, 00810 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 00811 CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA, 00812 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 00813 $ LWORK-NWORK+1, IERR ) 00814 * 00815 * Multiply Q in U by left singular vectors of R in 00816 * WORK(IU), storing result in A 00817 * (Workspace: need N*N) 00818 * 00819 CALL DGEMM( 'N', 'N', M, N, N, ONE, U, LDU, WORK( IU ), 00820 $ LDWRKU, ZERO, A, LDA ) 00821 * 00822 * Copy left singular vectors of A from A to U 00823 * 00824 CALL DLACPY( 'F', M, N, A, LDA, U, LDU ) 00825 * 00826 END IF 00827 * 00828 ELSE 00829 * 00830 * M .LT. MNTHR 00831 * 00832 * Path 5 (M at least N, but not much larger) 00833 * Reduce to bidiagonal form without QR decomposition 00834 * 00835 IE = 1 00836 ITAUQ = IE + N 00837 ITAUP = ITAUQ + N 00838 NWORK = ITAUP + N 00839 * 00840 * Bidiagonalize A 00841 * (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB) 00842 * 00843 CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 00844 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 00845 $ IERR ) 00846 IF( WNTQN ) THEN 00847 * 00848 * Perform bidiagonal SVD, only computing singular values 00849 * (Workspace: need N+BDSPAC) 00850 * 00851 CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1, 00852 $ DUM, IDUM, WORK( NWORK ), IWORK, INFO ) 00853 ELSE IF( WNTQO ) THEN 00854 IU = NWORK 00855 IF( LWORK.GE.M*N+3*N+BDSPAC ) THEN 00856 * 00857 * WORK( IU ) is M by N 00858 * 00859 LDWRKU = M 00860 NWORK = IU + LDWRKU*N 00861 CALL DLASET( 'F', M, N, ZERO, ZERO, WORK( IU ), 00862 $ LDWRKU ) 00863 ELSE 00864 * 00865 * WORK( IU ) is N by N 00866 * 00867 LDWRKU = N 00868 NWORK = IU + LDWRKU*N 00869 * 00870 * WORK(IR) is LDWRKR by N 00871 * 00872 IR = NWORK 00873 LDWRKR = ( LWORK-N*N-3*N ) / N 00874 END IF 00875 NWORK = IU + LDWRKU*N 00876 * 00877 * Perform bidiagonal SVD, computing left singular vectors 00878 * of bidiagonal matrix in WORK(IU) and computing right 00879 * singular vectors of bidiagonal matrix in VT 00880 * (Workspace: need N+N*N+BDSPAC) 00881 * 00882 CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), 00883 $ LDWRKU, VT, LDVT, DUM, IDUM, WORK( NWORK ), 00884 $ IWORK, INFO ) 00885 * 00886 * Overwrite VT by right singular vectors of A 00887 * (Workspace: need N*N+2*N, prefer N*N+N+N*NB) 00888 * 00889 CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA, 00890 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 00891 $ LWORK-NWORK+1, IERR ) 00892 * 00893 IF( LWORK.GE.M*N+3*N+BDSPAC ) THEN 00894 * 00895 * Overwrite WORK(IU) by left singular vectors of A 00896 * (Workspace: need N*N+2*N, prefer N*N+N+N*NB) 00897 * 00898 CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA, 00899 $ WORK( ITAUQ ), WORK( IU ), LDWRKU, 00900 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 00901 * 00902 * Copy left singular vectors of A from WORK(IU) to A 00903 * 00904 CALL DLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA ) 00905 ELSE 00906 * 00907 * Generate Q in A 00908 * (Workspace: need N*N+2*N, prefer N*N+N+N*NB) 00909 * 00910 CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ), 00911 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 00912 * 00913 * Multiply Q in A by left singular vectors of 00914 * bidiagonal matrix in WORK(IU), storing result in 00915 * WORK(IR) and copying to A 00916 * (Workspace: need 2*N*N, prefer N*N+M*N) 00917 * 00918 DO 20 I = 1, M, LDWRKR 00919 CHUNK = MIN( M-I+1, LDWRKR ) 00920 CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ), 00921 $ LDA, WORK( IU ), LDWRKU, ZERO, 00922 $ WORK( IR ), LDWRKR ) 00923 CALL DLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR, 00924 $ A( I, 1 ), LDA ) 00925 20 CONTINUE 00926 END IF 00927 * 00928 ELSE IF( WNTQS ) THEN 00929 * 00930 * Perform bidiagonal SVD, computing left singular vectors 00931 * of bidiagonal matrix in U and computing right singular 00932 * vectors of bidiagonal matrix in VT 00933 * (Workspace: need N+BDSPAC) 00934 * 00935 CALL DLASET( 'F', M, N, ZERO, ZERO, U, LDU ) 00936 CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT, 00937 $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK, 00938 $ INFO ) 00939 * 00940 * Overwrite U by left singular vectors of A and VT 00941 * by right singular vectors of A 00942 * (Workspace: need 3*N, prefer 2*N+N*NB) 00943 * 00944 CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA, 00945 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 00946 $ LWORK-NWORK+1, IERR ) 00947 CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA, 00948 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 00949 $ LWORK-NWORK+1, IERR ) 00950 ELSE IF( WNTQA ) THEN 00951 * 00952 * Perform bidiagonal SVD, computing left singular vectors 00953 * of bidiagonal matrix in U and computing right singular 00954 * vectors of bidiagonal matrix in VT 00955 * (Workspace: need N+BDSPAC) 00956 * 00957 CALL DLASET( 'F', M, M, ZERO, ZERO, U, LDU ) 00958 CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT, 00959 $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK, 00960 $ INFO ) 00961 * 00962 * Set the right corner of U to identity matrix 00963 * 00964 IF( M.GT.N ) THEN 00965 CALL DLASET( 'F', M-N, M-N, ZERO, ONE, U( N+1, N+1 ), 00966 $ LDU ) 00967 END IF 00968 * 00969 * Overwrite U by left singular vectors of A and VT 00970 * by right singular vectors of A 00971 * (Workspace: need N*N+2*N+M, prefer N*N+2*N+M*NB) 00972 * 00973 CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA, 00974 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 00975 $ LWORK-NWORK+1, IERR ) 00976 CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA, 00977 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 00978 $ LWORK-NWORK+1, IERR ) 00979 END IF 00980 * 00981 END IF 00982 * 00983 ELSE 00984 * 00985 * A has more columns than rows. If A has sufficiently more 00986 * columns than rows, first reduce using the LQ decomposition (if 00987 * sufficient workspace available) 00988 * 00989 IF( N.GE.MNTHR ) THEN 00990 * 00991 IF( WNTQN ) THEN 00992 * 00993 * Path 1t (N much larger than M, JOBZ='N') 00994 * No singular vectors to be computed 00995 * 00996 ITAU = 1 00997 NWORK = ITAU + M 00998 * 00999 * Compute A=L*Q 01000 * (Workspace: need 2*M, prefer M+M*NB) 01001 * 01002 CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 01003 $ LWORK-NWORK+1, IERR ) 01004 * 01005 * Zero out above L 01006 * 01007 CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA ) 01008 IE = 1 01009 ITAUQ = IE + M 01010 ITAUP = ITAUQ + M 01011 NWORK = ITAUP + M 01012 * 01013 * Bidiagonalize L in A 01014 * (Workspace: need 4*M, prefer 3*M+2*M*NB) 01015 * 01016 CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 01017 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 01018 $ IERR ) 01019 NWORK = IE + M 01020 * 01021 * Perform bidiagonal SVD, computing singular values only 01022 * (Workspace: need M+BDSPAC) 01023 * 01024 CALL DBDSDC( 'U', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1, 01025 $ DUM, IDUM, WORK( NWORK ), IWORK, INFO ) 01026 * 01027 ELSE IF( WNTQO ) THEN 01028 * 01029 * Path 2t (N much larger than M, JOBZ='O') 01030 * M right singular vectors to be overwritten on A and 01031 * M left singular vectors to be computed in U 01032 * 01033 IVT = 1 01034 * 01035 * IVT is M by M 01036 * 01037 IL = IVT + M*M 01038 IF( LWORK.GE.M*N+M*M+3*M+BDSPAC ) THEN 01039 * 01040 * WORK(IL) is M by N 01041 * 01042 LDWRKL = M 01043 CHUNK = N 01044 ELSE 01045 LDWRKL = M 01046 CHUNK = ( LWORK-M*M ) / M 01047 END IF 01048 ITAU = IL + LDWRKL*M 01049 NWORK = ITAU + M 01050 * 01051 * Compute A=L*Q 01052 * (Workspace: need M*M+2*M, prefer M*M+M+M*NB) 01053 * 01054 CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 01055 $ LWORK-NWORK+1, IERR ) 01056 * 01057 * Copy L to WORK(IL), zeroing about above it 01058 * 01059 CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL ) 01060 CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, 01061 $ WORK( IL+LDWRKL ), LDWRKL ) 01062 * 01063 * Generate Q in A 01064 * (Workspace: need M*M+2*M, prefer M*M+M+M*NB) 01065 * 01066 CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), 01067 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01068 IE = ITAU 01069 ITAUQ = IE + M 01070 ITAUP = ITAUQ + M 01071 NWORK = ITAUP + M 01072 * 01073 * Bidiagonalize L in WORK(IL) 01074 * (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) 01075 * 01076 CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ), 01077 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), 01078 $ LWORK-NWORK+1, IERR ) 01079 * 01080 * Perform bidiagonal SVD, computing left singular vectors 01081 * of bidiagonal matrix in U, and computing right singular 01082 * vectors of bidiagonal matrix in WORK(IVT) 01083 * (Workspace: need M+M*M+BDSPAC) 01084 * 01085 CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU, 01086 $ WORK( IVT ), M, DUM, IDUM, WORK( NWORK ), 01087 $ IWORK, INFO ) 01088 * 01089 * Overwrite U by left singular vectors of L and WORK(IVT) 01090 * by right singular vectors of L 01091 * (Workspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB) 01092 * 01093 CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL, 01094 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 01095 $ LWORK-NWORK+1, IERR ) 01096 CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL, 01097 $ WORK( ITAUP ), WORK( IVT ), M, 01098 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01099 * 01100 * Multiply right singular vectors of L in WORK(IVT) by Q 01101 * in A, storing result in WORK(IL) and copying to A 01102 * (Workspace: need 2*M*M, prefer M*M+M*N) 01103 * 01104 DO 30 I = 1, N, CHUNK 01105 BLK = MIN( N-I+1, CHUNK ) 01106 CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ), M, 01107 $ A( 1, I ), LDA, ZERO, WORK( IL ), LDWRKL ) 01108 CALL DLACPY( 'F', M, BLK, WORK( IL ), LDWRKL, 01109 $ A( 1, I ), LDA ) 01110 30 CONTINUE 01111 * 01112 ELSE IF( WNTQS ) THEN 01113 * 01114 * Path 3t (N much larger than M, JOBZ='S') 01115 * M right singular vectors to be computed in VT and 01116 * M left singular vectors to be computed in U 01117 * 01118 IL = 1 01119 * 01120 * WORK(IL) is M by M 01121 * 01122 LDWRKL = M 01123 ITAU = IL + LDWRKL*M 01124 NWORK = ITAU + M 01125 * 01126 * Compute A=L*Q 01127 * (Workspace: need M*M+2*M, prefer M*M+M+M*NB) 01128 * 01129 CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 01130 $ LWORK-NWORK+1, IERR ) 01131 * 01132 * Copy L to WORK(IL), zeroing out above it 01133 * 01134 CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL ) 01135 CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, 01136 $ WORK( IL+LDWRKL ), LDWRKL ) 01137 * 01138 * Generate Q in A 01139 * (Workspace: need M*M+2*M, prefer M*M+M+M*NB) 01140 * 01141 CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), 01142 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01143 IE = ITAU 01144 ITAUQ = IE + M 01145 ITAUP = ITAUQ + M 01146 NWORK = ITAUP + M 01147 * 01148 * Bidiagonalize L in WORK(IU), copying result to U 01149 * (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) 01150 * 01151 CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ), 01152 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), 01153 $ LWORK-NWORK+1, IERR ) 01154 * 01155 * Perform bidiagonal SVD, computing left singular vectors 01156 * of bidiagonal matrix in U and computing right singular 01157 * vectors of bidiagonal matrix in VT 01158 * (Workspace: need M+BDSPAC) 01159 * 01160 CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU, VT, 01161 $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK, 01162 $ INFO ) 01163 * 01164 * Overwrite U by left singular vectors of L and VT 01165 * by right singular vectors of L 01166 * (Workspace: need M*M+3*M, prefer M*M+2*M+M*NB) 01167 * 01168 CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL, 01169 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 01170 $ LWORK-NWORK+1, IERR ) 01171 CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL, 01172 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 01173 $ LWORK-NWORK+1, IERR ) 01174 * 01175 * Multiply right singular vectors of L in WORK(IL) by 01176 * Q in A, storing result in VT 01177 * (Workspace: need M*M) 01178 * 01179 CALL DLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL ) 01180 CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IL ), LDWRKL, 01181 $ A, LDA, ZERO, VT, LDVT ) 01182 * 01183 ELSE IF( WNTQA ) THEN 01184 * 01185 * Path 4t (N much larger than M, JOBZ='A') 01186 * N right singular vectors to be computed in VT and 01187 * M left singular vectors to be computed in U 01188 * 01189 IVT = 1 01190 * 01191 * WORK(IVT) is M by M 01192 * 01193 LDWKVT = M 01194 ITAU = IVT + LDWKVT*M 01195 NWORK = ITAU + M 01196 * 01197 * Compute A=L*Q, copying result to VT 01198 * (Workspace: need M*M+2*M, prefer M*M+M+M*NB) 01199 * 01200 CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 01201 $ LWORK-NWORK+1, IERR ) 01202 CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) 01203 * 01204 * Generate Q in VT 01205 * (Workspace: need M*M+2*M, prefer M*M+M+M*NB) 01206 * 01207 CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ), 01208 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01209 * 01210 * Produce L in A, zeroing out other entries 01211 * 01212 CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA ) 01213 IE = ITAU 01214 ITAUQ = IE + M 01215 ITAUP = ITAUQ + M 01216 NWORK = ITAUP + M 01217 * 01218 * Bidiagonalize L in A 01219 * (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) 01220 * 01221 CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 01222 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 01223 $ IERR ) 01224 * 01225 * Perform bidiagonal SVD, computing left singular vectors 01226 * of bidiagonal matrix in U and computing right singular 01227 * vectors of bidiagonal matrix in WORK(IVT) 01228 * (Workspace: need M+M*M+BDSPAC) 01229 * 01230 CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU, 01231 $ WORK( IVT ), LDWKVT, DUM, IDUM, 01232 $ WORK( NWORK ), IWORK, INFO ) 01233 * 01234 * Overwrite U by left singular vectors of L and WORK(IVT) 01235 * by right singular vectors of L 01236 * (Workspace: need M*M+3*M, prefer M*M+2*M+M*NB) 01237 * 01238 CALL DORMBR( 'Q', 'L', 'N', M, M, M, A, LDA, 01239 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 01240 $ LWORK-NWORK+1, IERR ) 01241 CALL DORMBR( 'P', 'R', 'T', M, M, M, A, LDA, 01242 $ WORK( ITAUP ), WORK( IVT ), LDWKVT, 01243 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01244 * 01245 * Multiply right singular vectors of L in WORK(IVT) by 01246 * Q in VT, storing result in A 01247 * (Workspace: need M*M) 01248 * 01249 CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IVT ), LDWKVT, 01250 $ VT, LDVT, ZERO, A, LDA ) 01251 * 01252 * Copy right singular vectors of A from A to VT 01253 * 01254 CALL DLACPY( 'F', M, N, A, LDA, VT, LDVT ) 01255 * 01256 END IF 01257 * 01258 ELSE 01259 * 01260 * N .LT. MNTHR 01261 * 01262 * Path 5t (N greater than M, but not much larger) 01263 * Reduce to bidiagonal form without LQ decomposition 01264 * 01265 IE = 1 01266 ITAUQ = IE + M 01267 ITAUP = ITAUQ + M 01268 NWORK = ITAUP + M 01269 * 01270 * Bidiagonalize A 01271 * (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB) 01272 * 01273 CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 01274 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 01275 $ IERR ) 01276 IF( WNTQN ) THEN 01277 * 01278 * Perform bidiagonal SVD, only computing singular values 01279 * (Workspace: need M+BDSPAC) 01280 * 01281 CALL DBDSDC( 'L', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1, 01282 $ DUM, IDUM, WORK( NWORK ), IWORK, INFO ) 01283 ELSE IF( WNTQO ) THEN 01284 LDWKVT = M 01285 IVT = NWORK 01286 IF( LWORK.GE.M*N+3*M+BDSPAC ) THEN 01287 * 01288 * WORK( IVT ) is M by N 01289 * 01290 CALL DLASET( 'F', M, N, ZERO, ZERO, WORK( IVT ), 01291 $ LDWKVT ) 01292 NWORK = IVT + LDWKVT*N 01293 ELSE 01294 * 01295 * WORK( IVT ) is M by M 01296 * 01297 NWORK = IVT + LDWKVT*M 01298 IL = NWORK 01299 * 01300 * WORK(IL) is M by CHUNK 01301 * 01302 CHUNK = ( LWORK-M*M-3*M ) / M 01303 END IF 01304 * 01305 * Perform bidiagonal SVD, computing left singular vectors 01306 * of bidiagonal matrix in U and computing right singular 01307 * vectors of bidiagonal matrix in WORK(IVT) 01308 * (Workspace: need M*M+BDSPAC) 01309 * 01310 CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, 01311 $ WORK( IVT ), LDWKVT, DUM, IDUM, 01312 $ WORK( NWORK ), IWORK, INFO ) 01313 * 01314 * Overwrite U by left singular vectors of A 01315 * (Workspace: need M*M+2*M, prefer M*M+M+M*NB) 01316 * 01317 CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA, 01318 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 01319 $ LWORK-NWORK+1, IERR ) 01320 * 01321 IF( LWORK.GE.M*N+3*M+BDSPAC ) THEN 01322 * 01323 * Overwrite WORK(IVT) by left singular vectors of A 01324 * (Workspace: need M*M+2*M, prefer M*M+M+M*NB) 01325 * 01326 CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA, 01327 $ WORK( ITAUP ), WORK( IVT ), LDWKVT, 01328 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01329 * 01330 * Copy right singular vectors of A from WORK(IVT) to A 01331 * 01332 CALL DLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA ) 01333 ELSE 01334 * 01335 * Generate P**T in A 01336 * (Workspace: need M*M+2*M, prefer M*M+M+M*NB) 01337 * 01338 CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), 01339 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01340 * 01341 * Multiply Q in A by right singular vectors of 01342 * bidiagonal matrix in WORK(IVT), storing result in 01343 * WORK(IL) and copying to A 01344 * (Workspace: need 2*M*M, prefer M*M+M*N) 01345 * 01346 DO 40 I = 1, N, CHUNK 01347 BLK = MIN( N-I+1, CHUNK ) 01348 CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ), 01349 $ LDWKVT, A( 1, I ), LDA, ZERO, 01350 $ WORK( IL ), M ) 01351 CALL DLACPY( 'F', M, BLK, WORK( IL ), M, A( 1, I ), 01352 $ LDA ) 01353 40 CONTINUE 01354 END IF 01355 ELSE IF( WNTQS ) THEN 01356 * 01357 * Perform bidiagonal SVD, computing left singular vectors 01358 * of bidiagonal matrix in U and computing right singular 01359 * vectors of bidiagonal matrix in VT 01360 * (Workspace: need M+BDSPAC) 01361 * 01362 CALL DLASET( 'F', M, N, ZERO, ZERO, VT, LDVT ) 01363 CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT, 01364 $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK, 01365 $ INFO ) 01366 * 01367 * Overwrite U by left singular vectors of A and VT 01368 * by right singular vectors of A 01369 * (Workspace: need 3*M, prefer 2*M+M*NB) 01370 * 01371 CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA, 01372 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 01373 $ LWORK-NWORK+1, IERR ) 01374 CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA, 01375 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 01376 $ LWORK-NWORK+1, IERR ) 01377 ELSE IF( WNTQA ) THEN 01378 * 01379 * Perform bidiagonal SVD, computing left singular vectors 01380 * of bidiagonal matrix in U and computing right singular 01381 * vectors of bidiagonal matrix in VT 01382 * (Workspace: need M+BDSPAC) 01383 * 01384 CALL DLASET( 'F', N, N, ZERO, ZERO, VT, LDVT ) 01385 CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT, 01386 $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK, 01387 $ INFO ) 01388 * 01389 * Set the right corner of VT to identity matrix 01390 * 01391 IF( N.GT.M ) THEN 01392 CALL DLASET( 'F', N-M, N-M, ZERO, ONE, VT( M+1, M+1 ), 01393 $ LDVT ) 01394 END IF 01395 * 01396 * Overwrite U by left singular vectors of A and VT 01397 * by right singular vectors of A 01398 * (Workspace: need 2*M+N, prefer 2*M+N*NB) 01399 * 01400 CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA, 01401 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 01402 $ LWORK-NWORK+1, IERR ) 01403 CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA, 01404 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 01405 $ LWORK-NWORK+1, IERR ) 01406 END IF 01407 * 01408 END IF 01409 * 01410 END IF 01411 * 01412 * Undo scaling if necessary 01413 * 01414 IF( ISCL.EQ.1 ) THEN 01415 IF( ANRM.GT.BIGNUM ) 01416 $ CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN, 01417 $ IERR ) 01418 IF( ANRM.LT.SMLNUM ) 01419 $ CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN, 01420 $ IERR ) 01421 END IF 01422 * 01423 * Return optimal workspace in WORK(1) 01424 * 01425 WORK( 1 ) = MAXWRK 01426 * 01427 RETURN 01428 * 01429 * End of DGESDD 01430 * 01431 END