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