 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE SGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, 00002 $ RANK, WORK, LWORK, IWORK, 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 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK 00011 REAL RCOND 00012 * .. 00013 * .. Array Arguments .. 00014 INTEGER IWORK( * ) 00015 REAL A( LDA, * ), B( LDB, * ), S( * ), WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * SGELSD computes the minimum-norm solution to a real linear least 00022 * squares problem: 00023 * minimize 2-norm(| b - A*x |) 00024 * using the singular value decomposition (SVD) of A. A is an M-by-N 00025 * matrix which may be rank-deficient. 00026 * 00027 * Several right hand side vectors b and solution vectors x can be 00028 * handled in a single call; they are stored as the columns of the 00029 * M-by-NRHS right hand side matrix B and the N-by-NRHS solution 00030 * matrix X. 00031 * 00032 * The problem is solved in three steps: 00033 * (1) Reduce the coefficient matrix A to bidiagonal form with 00034 * Householder transformations, reducing the original problem 00035 * into a "bidiagonal least squares problem" (BLS) 00036 * (2) Solve the BLS using a divide and conquer approach. 00037 * (3) Apply back all the Householder tranformations to solve 00038 * the original least squares problem. 00039 * 00040 * The effective rank of A is determined by treating as zero those 00041 * singular values which are less than RCOND times the largest singular 00042 * value. 00043 * 00044 * The divide and conquer algorithm makes very mild assumptions about 00045 * floating point arithmetic. It will work on machines with a guard 00046 * digit in add/subtract, or on those binary machines without guard 00047 * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or 00048 * Cray-2. It could conceivably fail on hexadecimal or decimal machines 00049 * without guard digits, but we know of none. 00050 * 00051 * Arguments 00052 * ========= 00053 * 00054 * M (input) INTEGER 00055 * The number of rows of A. M >= 0. 00056 * 00057 * N (input) INTEGER 00058 * The number of columns of A. N >= 0. 00059 * 00060 * NRHS (input) INTEGER 00061 * The number of right hand sides, i.e., the number of columns 00062 * of the matrices B and X. NRHS >= 0. 00063 * 00064 * A (input) REAL array, dimension (LDA,N) 00065 * On entry, the M-by-N matrix A. 00066 * On exit, A has been destroyed. 00067 * 00068 * LDA (input) INTEGER 00069 * The leading dimension of the array A. LDA >= max(1,M). 00070 * 00071 * B (input/output) REAL array, dimension (LDB,NRHS) 00072 * On entry, the M-by-NRHS right hand side matrix B. 00073 * On exit, B is overwritten by the N-by-NRHS solution 00074 * matrix X. If m >= n and RANK = n, the residual 00075 * sum-of-squares for the solution in the i-th column is given 00076 * by the sum of squares of elements n+1:m in that column. 00077 * 00078 * LDB (input) INTEGER 00079 * The leading dimension of the array B. LDB >= max(1,max(M,N)). 00080 * 00081 * S (output) REAL array, dimension (min(M,N)) 00082 * The singular values of A in decreasing order. 00083 * The condition number of A in the 2-norm = S(1)/S(min(m,n)). 00084 * 00085 * RCOND (input) REAL 00086 * RCOND is used to determine the effective rank of A. 00087 * Singular values S(i) <= RCOND*S(1) are treated as zero. 00088 * If RCOND < 0, machine precision is used instead. 00089 * 00090 * RANK (output) INTEGER 00091 * The effective rank of A, i.e., the number of singular values 00092 * which are greater than RCOND*S(1). 00093 * 00094 * WORK (workspace/output) REAL array, dimension (MAX(1,LWORK)) 00095 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00096 * 00097 * LWORK (input) INTEGER 00098 * The dimension of the array WORK. LWORK must be at least 1. 00099 * The exact minimum amount of workspace needed depends on M, 00100 * N and NRHS. As long as LWORK is at least 00101 * 12*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2, 00102 * if M is greater than or equal to N or 00103 * 12*M + 2*M*SMLSIZ + 8*M*NLVL + M*NRHS + (SMLSIZ+1)**2, 00104 * if M is less than N, the code will execute correctly. 00105 * SMLSIZ is returned by ILAENV and is equal to the maximum 00106 * size of the subproblems at the bottom of the computation 00107 * tree (usually about 25), and 00108 * NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 ) 00109 * For good performance, LWORK should generally be larger. 00110 * 00111 * If LWORK = -1, then a workspace query is assumed; the routine 00112 * only calculates the optimal size of the array WORK and the 00113 * minimum size of the array IWORK, and returns these values as 00114 * the first entries of the WORK and IWORK arrays, and no error 00115 * message related to LWORK is issued by XERBLA. 00116 * 00117 * IWORK (workspace) INTEGER array, dimension (MAX(1,LIWORK)) 00118 * LIWORK >= max(1, 3*MINMN*NLVL + 11*MINMN), 00119 * where MINMN = MIN( M,N ). 00120 * On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK. 00121 * 00122 * INFO (output) INTEGER 00123 * = 0: successful exit 00124 * < 0: if INFO = -i, the i-th argument had an illegal value. 00125 * > 0: the algorithm for computing the SVD failed to converge; 00126 * if INFO = i, i off-diagonal elements of an intermediate 00127 * bidiagonal form did not converge to zero. 00128 * 00129 * Further Details 00130 * =============== 00131 * 00132 * Based on contributions by 00133 * Ming Gu and Ren-Cang Li, Computer Science Division, University of 00134 * California at Berkeley, USA 00135 * Osni Marques, LBNL/NERSC, USA 00136 * 00137 * ===================================================================== 00138 * 00139 * .. Parameters .. 00140 REAL ZERO, ONE, TWO 00141 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 ) 00142 * .. 00143 * .. Local Scalars .. 00144 LOGICAL LQUERY 00145 INTEGER IASCL, IBSCL, IE, IL, ITAU, ITAUP, ITAUQ, 00146 $ LDWORK, LIWORK, MAXMN, MAXWRK, MINMN, MINWRK, 00147 $ MM, MNTHR, NLVL, NWORK, SMLSIZ, WLALSD 00148 REAL ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM 00149 * .. 00150 * .. External Subroutines .. 00151 EXTERNAL SGEBRD, SGELQF, SGEQRF, SLABAD, SLACPY, SLALSD, 00152 $ SLASCL, SLASET, SORMBR, SORMLQ, SORMQR, XERBLA 00153 * .. 00154 * .. External Functions .. 00155 INTEGER ILAENV 00156 REAL SLAMCH, SLANGE 00157 EXTERNAL SLAMCH, SLANGE, ILAENV 00158 * .. 00159 * .. Intrinsic Functions .. 00160 INTRINSIC INT, LOG, MAX, MIN, REAL 00161 * .. 00162 * .. Executable Statements .. 00163 * 00164 * Test the input arguments. 00165 * 00166 INFO = 0 00167 MINMN = MIN( M, N ) 00168 MAXMN = MAX( M, N ) 00169 LQUERY = ( LWORK.EQ.-1 ) 00170 IF( M.LT.0 ) THEN 00171 INFO = -1 00172 ELSE IF( N.LT.0 ) THEN 00173 INFO = -2 00174 ELSE IF( NRHS.LT.0 ) THEN 00175 INFO = -3 00176 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00177 INFO = -5 00178 ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN 00179 INFO = -7 00180 END IF 00181 * 00182 * Compute workspace. 00183 * (Note: Comments in the code beginning "Workspace:" describe the 00184 * minimal amount of workspace needed at that point in the code, 00185 * as well as the preferred amount for good performance. 00186 * NB refers to the optimal block size for the immediately 00187 * following subroutine, as returned by ILAENV.) 00188 * 00189 IF( INFO.EQ.0 ) THEN 00190 MINWRK = 1 00191 MAXWRK = 1 00192 LIWORK = 1 00193 IF( MINMN.GT.0 ) THEN 00194 SMLSIZ = ILAENV( 9, 'SGELSD', ' ', 0, 0, 0, 0 ) 00195 MNTHR = ILAENV( 6, 'SGELSD', ' ', M, N, NRHS, -1 ) 00196 NLVL = MAX( INT( LOG( REAL( MINMN ) / REAL( SMLSIZ + 1 ) ) / 00197 $ LOG( TWO ) ) + 1, 0 ) 00198 LIWORK = 3*MINMN*NLVL + 11*MINMN 00199 MM = M 00200 IF( M.GE.N .AND. M.GE.MNTHR ) THEN 00201 * 00202 * Path 1a - overdetermined, with many more rows than 00203 * columns. 00204 * 00205 MM = N 00206 MAXWRK = MAX( MAXWRK, N + N*ILAENV( 1, 'SGEQRF', ' ', M, 00207 $ N, -1, -1 ) ) 00208 MAXWRK = MAX( MAXWRK, N + NRHS*ILAENV( 1, 'SORMQR', 'LT', 00209 $ M, NRHS, N, -1 ) ) 00210 END IF 00211 IF( M.GE.N ) THEN 00212 * 00213 * Path 1 - overdetermined or exactly determined. 00214 * 00215 MAXWRK = MAX( MAXWRK, 3*N + ( MM + N )*ILAENV( 1, 00216 $ 'SGEBRD', ' ', MM, N, -1, -1 ) ) 00217 MAXWRK = MAX( MAXWRK, 3*N + NRHS*ILAENV( 1, 'SORMBR', 00218 $ 'QLT', MM, NRHS, N, -1 ) ) 00219 MAXWRK = MAX( MAXWRK, 3*N + ( N - 1 )*ILAENV( 1, 00220 $ 'SORMBR', 'PLN', N, NRHS, N, -1 ) ) 00221 WLALSD = 9*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + 00222 $ ( SMLSIZ + 1 )**2 00223 MAXWRK = MAX( MAXWRK, 3*N + WLALSD ) 00224 MINWRK = MAX( 3*N + MM, 3*N + NRHS, 3*N + WLALSD ) 00225 END IF 00226 IF( N.GT.M ) THEN 00227 WLALSD = 9*M + 2*M*SMLSIZ + 8*M*NLVL + M*NRHS + 00228 $ ( SMLSIZ + 1 )**2 00229 IF( N.GE.MNTHR ) THEN 00230 * 00231 * Path 2a - underdetermined, with many more columns 00232 * than rows. 00233 * 00234 MAXWRK = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, 00235 $ -1 ) 00236 MAXWRK = MAX( MAXWRK, M*M + 4*M + 2*M*ILAENV( 1, 00237 $ 'SGEBRD', ' ', M, M, -1, -1 ) ) 00238 MAXWRK = MAX( MAXWRK, M*M + 4*M + NRHS*ILAENV( 1, 00239 $ 'SORMBR', 'QLT', M, NRHS, M, -1 ) ) 00240 MAXWRK = MAX( MAXWRK, M*M + 4*M + ( M - 1 )*ILAENV( 1, 00241 $ 'SORMBR', 'PLN', M, NRHS, M, -1 ) ) 00242 IF( NRHS.GT.1 ) THEN 00243 MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS ) 00244 ELSE 00245 MAXWRK = MAX( MAXWRK, M*M + 2*M ) 00246 END IF 00247 MAXWRK = MAX( MAXWRK, M + NRHS*ILAENV( 1, 'SORMLQ', 00248 $ 'LT', N, NRHS, M, -1 ) ) 00249 MAXWRK = MAX( MAXWRK, M*M + 4*M + WLALSD ) 00250 ! XXX: Ensure the Path 2a case below is triggered. The workspace 00251 ! calculation should use queries for all routines eventually. 00252 MAXWRK = MAX( MAXWRK, 00253 $ 4*M+M*M+MAX( M, 2*M-4, NRHS, N-3*M ) ) 00254 ELSE 00255 * 00256 * Path 2 - remaining underdetermined cases. 00257 * 00258 MAXWRK = 3*M + ( N + M )*ILAENV( 1, 'SGEBRD', ' ', M, 00259 $ N, -1, -1 ) 00260 MAXWRK = MAX( MAXWRK, 3*M + NRHS*ILAENV( 1, 'SORMBR', 00261 $ 'QLT', M, NRHS, N, -1 ) ) 00262 MAXWRK = MAX( MAXWRK, 3*M + M*ILAENV( 1, 'SORMBR', 00263 $ 'PLN', N, NRHS, M, -1 ) ) 00264 MAXWRK = MAX( MAXWRK, 3*M + WLALSD ) 00265 END IF 00266 MINWRK = MAX( 3*M + NRHS, 3*M + M, 3*M + WLALSD ) 00267 END IF 00268 END IF 00269 MINWRK = MIN( MINWRK, MAXWRK ) 00270 WORK( 1 ) = MAXWRK 00271 IWORK( 1 ) = LIWORK 00272 * 00273 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN 00274 INFO = -12 00275 END IF 00276 END IF 00277 * 00278 IF( INFO.NE.0 ) THEN 00279 CALL XERBLA( 'SGELSD', -INFO ) 00280 RETURN 00281 ELSE IF( LQUERY ) THEN 00282 RETURN 00283 END IF 00284 * 00285 * Quick return if possible. 00286 * 00287 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 00288 RANK = 0 00289 RETURN 00290 END IF 00291 * 00292 * Get machine parameters. 00293 * 00294 EPS = SLAMCH( 'P' ) 00295 SFMIN = SLAMCH( 'S' ) 00296 SMLNUM = SFMIN / EPS 00297 BIGNUM = ONE / SMLNUM 00298 CALL SLABAD( SMLNUM, BIGNUM ) 00299 * 00300 * Scale A if max entry outside range [SMLNUM,BIGNUM]. 00301 * 00302 ANRM = SLANGE( 'M', M, N, A, LDA, WORK ) 00303 IASCL = 0 00304 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00305 * 00306 * Scale matrix norm up to SMLNUM. 00307 * 00308 CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO ) 00309 IASCL = 1 00310 ELSE IF( ANRM.GT.BIGNUM ) THEN 00311 * 00312 * Scale matrix norm down to BIGNUM. 00313 * 00314 CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO ) 00315 IASCL = 2 00316 ELSE IF( ANRM.EQ.ZERO ) THEN 00317 * 00318 * Matrix all zero. Return zero solution. 00319 * 00320 CALL SLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB ) 00321 CALL SLASET( 'F', MINMN, 1, ZERO, ZERO, S, 1 ) 00322 RANK = 0 00323 GO TO 10 00324 END IF 00325 * 00326 * Scale B if max entry outside range [SMLNUM,BIGNUM]. 00327 * 00328 BNRM = SLANGE( 'M', M, NRHS, B, LDB, WORK ) 00329 IBSCL = 0 00330 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 00331 * 00332 * Scale matrix norm up to SMLNUM. 00333 * 00334 CALL SLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO ) 00335 IBSCL = 1 00336 ELSE IF( BNRM.GT.BIGNUM ) THEN 00337 * 00338 * Scale matrix norm down to BIGNUM. 00339 * 00340 CALL SLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO ) 00341 IBSCL = 2 00342 END IF 00343 * 00344 * If M < N make sure certain entries of B are zero. 00345 * 00346 IF( M.LT.N ) 00347 $ CALL SLASET( 'F', N-M, NRHS, ZERO, ZERO, B( M+1, 1 ), LDB ) 00348 * 00349 * Overdetermined case. 00350 * 00351 IF( M.GE.N ) THEN 00352 * 00353 * Path 1 - overdetermined or exactly determined. 00354 * 00355 MM = M 00356 IF( M.GE.MNTHR ) THEN 00357 * 00358 * Path 1a - overdetermined, with many more rows than columns. 00359 * 00360 MM = N 00361 ITAU = 1 00362 NWORK = ITAU + N 00363 * 00364 * Compute A=Q*R. 00365 * (Workspace: need 2*N, prefer N+N*NB) 00366 * 00367 CALL SGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 00368 $ LWORK-NWORK+1, INFO ) 00369 * 00370 * Multiply B by transpose(Q). 00371 * (Workspace: need N+NRHS, prefer N+NRHS*NB) 00372 * 00373 CALL SORMQR( 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAU ), B, 00374 $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00375 * 00376 * Zero out below R. 00377 * 00378 IF( N.GT.1 ) THEN 00379 CALL SLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA ) 00380 END IF 00381 END IF 00382 * 00383 IE = 1 00384 ITAUQ = IE + N 00385 ITAUP = ITAUQ + N 00386 NWORK = ITAUP + N 00387 * 00388 * Bidiagonalize R in A. 00389 * (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB) 00390 * 00391 CALL SGEBRD( MM, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 00392 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 00393 $ INFO ) 00394 * 00395 * Multiply B by transpose of left bidiagonalizing vectors of R. 00396 * (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB) 00397 * 00398 CALL SORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, WORK( ITAUQ ), 00399 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00400 * 00401 * Solve the bidiagonal least squares problem. 00402 * 00403 CALL SLALSD( 'U', SMLSIZ, N, NRHS, S, WORK( IE ), B, LDB, 00404 $ RCOND, RANK, WORK( NWORK ), IWORK, INFO ) 00405 IF( INFO.NE.0 ) THEN 00406 GO TO 10 00407 END IF 00408 * 00409 * Multiply B by right bidiagonalizing vectors of R. 00410 * 00411 CALL SORMBR( 'P', 'L', 'N', N, NRHS, N, A, LDA, WORK( ITAUP ), 00412 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00413 * 00414 ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+ 00415 $ MAX( M, 2*M-4, NRHS, N-3*M, WLALSD ) ) THEN 00416 * 00417 * Path 2a - underdetermined, with many more columns than rows 00418 * and sufficient workspace for an efficient algorithm. 00419 * 00420 LDWORK = M 00421 IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ), 00422 $ M*LDA+M+M*NRHS, 4*M+M*LDA+WLALSD ) )LDWORK = LDA 00423 ITAU = 1 00424 NWORK = M + 1 00425 * 00426 * Compute A=L*Q. 00427 * (Workspace: need 2*M, prefer M+M*NB) 00428 * 00429 CALL SGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 00430 $ LWORK-NWORK+1, INFO ) 00431 IL = NWORK 00432 * 00433 * Copy L to WORK(IL), zeroing out above its diagonal. 00434 * 00435 CALL SLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK ) 00436 CALL SLASET( 'U', M-1, M-1, ZERO, ZERO, WORK( IL+LDWORK ), 00437 $ LDWORK ) 00438 IE = IL + LDWORK*M 00439 ITAUQ = IE + M 00440 ITAUP = ITAUQ + M 00441 NWORK = ITAUP + M 00442 * 00443 * Bidiagonalize L in WORK(IL). 00444 * (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB) 00445 * 00446 CALL SGEBRD( M, M, WORK( IL ), LDWORK, S, WORK( IE ), 00447 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), 00448 $ LWORK-NWORK+1, INFO ) 00449 * 00450 * Multiply B by transpose of left bidiagonalizing vectors of L. 00451 * (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB) 00452 * 00453 CALL SORMBR( 'Q', 'L', 'T', M, NRHS, M, WORK( IL ), LDWORK, 00454 $ WORK( ITAUQ ), B, LDB, WORK( NWORK ), 00455 $ LWORK-NWORK+1, INFO ) 00456 * 00457 * Solve the bidiagonal least squares problem. 00458 * 00459 CALL SLALSD( 'U', SMLSIZ, M, NRHS, S, WORK( IE ), B, LDB, 00460 $ RCOND, RANK, WORK( NWORK ), IWORK, INFO ) 00461 IF( INFO.NE.0 ) THEN 00462 GO TO 10 00463 END IF 00464 * 00465 * Multiply B by right bidiagonalizing vectors of L. 00466 * 00467 CALL SORMBR( 'P', 'L', 'N', M, NRHS, M, WORK( IL ), LDWORK, 00468 $ WORK( ITAUP ), B, LDB, WORK( NWORK ), 00469 $ LWORK-NWORK+1, INFO ) 00470 * 00471 * Zero out below first M rows of B. 00472 * 00473 CALL SLASET( 'F', N-M, NRHS, ZERO, ZERO, B( M+1, 1 ), LDB ) 00474 NWORK = ITAU + M 00475 * 00476 * Multiply transpose(Q) by B. 00477 * (Workspace: need M+NRHS, prefer M+NRHS*NB) 00478 * 00479 CALL SORMLQ( 'L', 'T', N, NRHS, M, A, LDA, WORK( ITAU ), B, 00480 $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00481 * 00482 ELSE 00483 * 00484 * Path 2 - remaining underdetermined cases. 00485 * 00486 IE = 1 00487 ITAUQ = IE + M 00488 ITAUP = ITAUQ + M 00489 NWORK = ITAUP + M 00490 * 00491 * Bidiagonalize A. 00492 * (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB) 00493 * 00494 CALL SGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 00495 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 00496 $ INFO ) 00497 * 00498 * Multiply B by transpose of left bidiagonalizing vectors. 00499 * (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB) 00500 * 00501 CALL SORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAUQ ), 00502 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00503 * 00504 * Solve the bidiagonal least squares problem. 00505 * 00506 CALL SLALSD( 'L', SMLSIZ, M, NRHS, S, WORK( IE ), B, LDB, 00507 $ RCOND, RANK, WORK( NWORK ), IWORK, INFO ) 00508 IF( INFO.NE.0 ) THEN 00509 GO TO 10 00510 END IF 00511 * 00512 * Multiply B by right bidiagonalizing vectors of A. 00513 * 00514 CALL SORMBR( 'P', 'L', 'N', N, NRHS, M, A, LDA, WORK( ITAUP ), 00515 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00516 * 00517 END IF 00518 * 00519 * Undo scaling. 00520 * 00521 IF( IASCL.EQ.1 ) THEN 00522 CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO ) 00523 CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN, 00524 $ INFO ) 00525 ELSE IF( IASCL.EQ.2 ) THEN 00526 CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO ) 00527 CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN, 00528 $ INFO ) 00529 END IF 00530 IF( IBSCL.EQ.1 ) THEN 00531 CALL SLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO ) 00532 ELSE IF( IBSCL.EQ.2 ) THEN 00533 CALL SLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO ) 00534 END IF 00535 * 00536 10 CONTINUE 00537 WORK( 1 ) = MAXWRK 00538 IWORK( 1 ) = LIWORK 00539 RETURN 00540 * 00541 * End of SGELSD 00542 * 00543 END
1.7.3 
