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