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