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