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