![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> SGELSS solves overdetermined or underdetermined systems 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 SGELSS + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgelss.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgelss.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgelss.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, 00022 * WORK, LWORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK 00026 * REAL RCOND 00027 * .. 00028 * .. Array Arguments .. 00029 * REAL A( LDA, * ), B( LDB, * ), S( * ), WORK( * ) 00030 * .. 00031 * 00032 * 00033 *> \par Purpose: 00034 * ============= 00035 *> 00036 *> \verbatim 00037 *> 00038 *> SGELSS computes the minimum norm solution to a real linear least 00039 *> squares problem: 00040 *> 00041 *> Minimize 2-norm(| b - A*x |). 00042 *> 00043 *> using the singular value decomposition (SVD) of A. A is an M-by-N 00044 *> matrix which may be rank-deficient. 00045 *> 00046 *> Several right hand side vectors b and solution vectors x can be 00047 *> handled in a single call; they are stored as the columns of the 00048 *> M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix 00049 *> X. 00050 *> 00051 *> The effective rank of A is determined by treating as zero those 00052 *> singular values which are less than RCOND times the largest singular 00053 *> value. 00054 *> \endverbatim 00055 * 00056 * Arguments: 00057 * ========== 00058 * 00059 *> \param[in] M 00060 *> \verbatim 00061 *> M is INTEGER 00062 *> The number of rows of the matrix A. M >= 0. 00063 *> \endverbatim 00064 *> 00065 *> \param[in] N 00066 *> \verbatim 00067 *> N is INTEGER 00068 *> The number of columns of the matrix A. N >= 0. 00069 *> \endverbatim 00070 *> 00071 *> \param[in] NRHS 00072 *> \verbatim 00073 *> NRHS is INTEGER 00074 *> The number of right hand sides, i.e., the number of columns 00075 *> of the matrices B and X. NRHS >= 0. 00076 *> \endverbatim 00077 *> 00078 *> \param[in,out] A 00079 *> \verbatim 00080 *> A is REAL array, dimension (LDA,N) 00081 *> On entry, the M-by-N matrix A. 00082 *> On exit, the first min(m,n) rows of A are overwritten with 00083 *> its right singular vectors, stored rowwise. 00084 *> \endverbatim 00085 *> 00086 *> \param[in] LDA 00087 *> \verbatim 00088 *> LDA is INTEGER 00089 *> The leading dimension of the array A. LDA >= max(1,M). 00090 *> \endverbatim 00091 *> 00092 *> \param[in,out] B 00093 *> \verbatim 00094 *> B is REAL array, dimension (LDB,NRHS) 00095 *> On entry, the M-by-NRHS right hand side matrix B. 00096 *> On exit, B is overwritten by the N-by-NRHS solution 00097 *> matrix X. If m >= n and RANK = n, the residual 00098 *> sum-of-squares for the solution in the i-th column is given 00099 *> by the sum of squares of elements n+1:m in that column. 00100 *> \endverbatim 00101 *> 00102 *> \param[in] LDB 00103 *> \verbatim 00104 *> LDB is INTEGER 00105 *> The leading dimension of the array B. LDB >= max(1,max(M,N)). 00106 *> \endverbatim 00107 *> 00108 *> \param[out] S 00109 *> \verbatim 00110 *> S is REAL array, dimension (min(M,N)) 00111 *> The singular values of A in decreasing order. 00112 *> The condition number of A in the 2-norm = S(1)/S(min(m,n)). 00113 *> \endverbatim 00114 *> 00115 *> \param[in] RCOND 00116 *> \verbatim 00117 *> RCOND is REAL 00118 *> RCOND is used to determine the effective rank of A. 00119 *> Singular values S(i) <= RCOND*S(1) are treated as zero. 00120 *> If RCOND < 0, machine precision is used instead. 00121 *> \endverbatim 00122 *> 00123 *> \param[out] RANK 00124 *> \verbatim 00125 *> RANK is INTEGER 00126 *> The effective rank of A, i.e., the number of singular values 00127 *> which are greater than RCOND*S(1). 00128 *> \endverbatim 00129 *> 00130 *> \param[out] WORK 00131 *> \verbatim 00132 *> WORK is REAL array, dimension (MAX(1,LWORK)) 00133 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00134 *> \endverbatim 00135 *> 00136 *> \param[in] LWORK 00137 *> \verbatim 00138 *> LWORK is INTEGER 00139 *> The dimension of the array WORK. LWORK >= 1, and also: 00140 *> LWORK >= 3*min(M,N) + max( 2*min(M,N), max(M,N), NRHS ) 00141 *> For good performance, LWORK should generally be larger. 00142 *> 00143 *> If LWORK = -1, then a workspace query is assumed; the routine 00144 *> only calculates the optimal size of the WORK array, returns 00145 *> this value as the first entry of the WORK array, and no error 00146 *> message related to LWORK is issued by XERBLA. 00147 *> \endverbatim 00148 *> 00149 *> \param[out] INFO 00150 *> \verbatim 00151 *> INFO is INTEGER 00152 *> = 0: successful exit 00153 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00154 *> > 0: the algorithm for computing the SVD failed to converge; 00155 *> if INFO = i, i off-diagonal elements of an intermediate 00156 *> bidiagonal form did not converge to zero. 00157 *> \endverbatim 00158 * 00159 * Authors: 00160 * ======== 00161 * 00162 *> \author Univ. of Tennessee 00163 *> \author Univ. of California Berkeley 00164 *> \author Univ. of Colorado Denver 00165 *> \author NAG Ltd. 00166 * 00167 *> \date November 2011 00168 * 00169 *> \ingroup realGEsolve 00170 * 00171 * ===================================================================== 00172 SUBROUTINE SGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, 00173 $ WORK, LWORK, INFO ) 00174 * 00175 * -- LAPACK driver routine (version 3.4.0) -- 00176 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00177 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00178 * November 2011 00179 * 00180 * .. Scalar Arguments .. 00181 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK 00182 REAL RCOND 00183 * .. 00184 * .. Array Arguments .. 00185 REAL A( LDA, * ), B( LDB, * ), S( * ), WORK( * ) 00186 * .. 00187 * 00188 * ===================================================================== 00189 * 00190 * .. Parameters .. 00191 REAL ZERO, ONE 00192 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00193 * .. 00194 * .. Local Scalars .. 00195 LOGICAL LQUERY 00196 INTEGER BDSPAC, BL, CHUNK, I, IASCL, IBSCL, IE, IL, 00197 $ ITAU, ITAUP, ITAUQ, IWORK, LDWORK, MAXMN, 00198 $ MAXWRK, MINMN, MINWRK, MM, MNTHR 00199 INTEGER LWORK_SGEQRF, LWORK_SORMQR, LWORK_SGEBRD, 00200 $ LWORK_SORMBR, LWORK_SORGBR, LWORK_SORMLQ 00201 REAL ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR 00202 * .. 00203 * .. Local Arrays .. 00204 REAL DUM( 1 ) 00205 * .. 00206 * .. External Subroutines .. 00207 EXTERNAL SBDSQR, SCOPY, SGEBRD, SGELQF, SGEMM, SGEMV, 00208 $ SGEQRF, SLABAD, SLACPY, SLASCL, SLASET, SORGBR, 00209 $ SORMBR, SORMLQ, SORMQR, SRSCL, XERBLA 00210 * .. 00211 * .. External Functions .. 00212 INTEGER ILAENV 00213 REAL SLAMCH, SLANGE 00214 EXTERNAL ILAENV, SLAMCH, SLANGE 00215 * .. 00216 * .. Intrinsic Functions .. 00217 INTRINSIC MAX, MIN 00218 * .. 00219 * .. Executable Statements .. 00220 * 00221 * Test the input arguments 00222 * 00223 INFO = 0 00224 MINMN = MIN( M, N ) 00225 MAXMN = MAX( M, N ) 00226 LQUERY = ( LWORK.EQ.-1 ) 00227 IF( M.LT.0 ) THEN 00228 INFO = -1 00229 ELSE IF( N.LT.0 ) THEN 00230 INFO = -2 00231 ELSE IF( NRHS.LT.0 ) THEN 00232 INFO = -3 00233 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00234 INFO = -5 00235 ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN 00236 INFO = -7 00237 END IF 00238 * 00239 * Compute workspace 00240 * (Note: Comments in the code beginning "Workspace:" describe the 00241 * minimal amount of workspace needed at that point in the code, 00242 * as well as the preferred amount for good performance. 00243 * NB refers to the optimal block size for the immediately 00244 * following subroutine, as returned by ILAENV.) 00245 * 00246 IF( INFO.EQ.0 ) THEN 00247 MINWRK = 1 00248 MAXWRK = 1 00249 IF( MINMN.GT.0 ) THEN 00250 MM = M 00251 MNTHR = ILAENV( 6, 'SGELSS', ' ', M, N, NRHS, -1 ) 00252 IF( M.GE.N .AND. M.GE.MNTHR ) THEN 00253 * 00254 * Path 1a - overdetermined, with many more rows than 00255 * columns 00256 * 00257 * Compute space needed for SGEQRF 00258 CALL SGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, INFO ) 00259 LWORK_SGEQRF=DUM(1) 00260 * Compute space needed for SORMQR 00261 CALL SORMQR( 'L', 'T', M, NRHS, N, A, LDA, DUM(1), B, 00262 $ LDB, DUM(1), -1, INFO ) 00263 LWORK_SORMQR=DUM(1) 00264 MM = N 00265 MAXWRK = MAX( MAXWRK, N + LWORK_SGEQRF ) 00266 MAXWRK = MAX( MAXWRK, N + LWORK_SORMQR ) 00267 END IF 00268 IF( M.GE.N ) THEN 00269 * 00270 * Path 1 - overdetermined or exactly determined 00271 * 00272 * Compute workspace needed for SBDSQR 00273 * 00274 BDSPAC = MAX( 1, 5*N ) 00275 * Compute space needed for SGEBRD 00276 CALL SGEBRD( MM, N, A, LDA, S, DUM(1), DUM(1), 00277 $ DUM(1), DUM(1), -1, INFO ) 00278 LWORK_SGEBRD=DUM(1) 00279 * Compute space needed for SORMBR 00280 CALL SORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, DUM(1), 00281 $ B, LDB, DUM(1), -1, INFO ) 00282 LWORK_SORMBR=DUM(1) 00283 * Compute space needed for SORGBR 00284 CALL SORGBR( 'P', N, N, N, A, LDA, DUM(1), 00285 $ DUM(1), -1, INFO ) 00286 LWORK_SORGBR=DUM(1) 00287 * Compute total workspace needed 00288 MAXWRK = MAX( MAXWRK, 3*N + LWORK_SGEBRD ) 00289 MAXWRK = MAX( MAXWRK, 3*N + LWORK_SORMBR ) 00290 MAXWRK = MAX( MAXWRK, 3*N + LWORK_SORGBR ) 00291 MAXWRK = MAX( MAXWRK, BDSPAC ) 00292 MAXWRK = MAX( MAXWRK, N*NRHS ) 00293 MINWRK = MAX( 3*N + MM, 3*N + NRHS, BDSPAC ) 00294 MAXWRK = MAX( MINWRK, MAXWRK ) 00295 END IF 00296 IF( N.GT.M ) THEN 00297 * 00298 * Compute workspace needed for SBDSQR 00299 * 00300 BDSPAC = MAX( 1, 5*M ) 00301 MINWRK = MAX( 3*M+NRHS, 3*M+N, BDSPAC ) 00302 IF( N.GE.MNTHR ) THEN 00303 * 00304 * Path 2a - underdetermined, with many more columns 00305 * than rows 00306 * 00307 * Compute space needed for SGEBRD 00308 CALL SGEBRD( M, M, A, LDA, S, DUM(1), DUM(1), 00309 $ DUM(1), DUM(1), -1, INFO ) 00310 LWORK_SGEBRD=DUM(1) 00311 * Compute space needed for SORMBR 00312 CALL SORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA, 00313 $ DUM(1), B, LDB, DUM(1), -1, INFO ) 00314 LWORK_SORMBR=DUM(1) 00315 * Compute space needed for SORGBR 00316 CALL SORGBR( 'P', M, M, M, A, LDA, DUM(1), 00317 $ DUM(1), -1, INFO ) 00318 LWORK_SORGBR=DUM(1) 00319 * Compute space needed for SORMLQ 00320 CALL SORMLQ( 'L', 'T', N, NRHS, M, A, LDA, DUM(1), 00321 $ B, LDB, DUM(1), -1, INFO ) 00322 LWORK_SORMLQ=DUM(1) 00323 * Compute total workspace needed 00324 MAXWRK = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, 00325 $ -1 ) 00326 MAXWRK = MAX( MAXWRK, M*M + 4*M + LWORK_SGEBRD ) 00327 MAXWRK = MAX( MAXWRK, M*M + 4*M + LWORK_SORMBR ) 00328 MAXWRK = MAX( MAXWRK, M*M + 4*M + LWORK_SORGBR ) 00329 MAXWRK = MAX( MAXWRK, M*M + M + BDSPAC ) 00330 IF( NRHS.GT.1 ) THEN 00331 MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS ) 00332 ELSE 00333 MAXWRK = MAX( MAXWRK, M*M + 2*M ) 00334 END IF 00335 MAXWRK = MAX( MAXWRK, M + LWORK_SORMLQ ) 00336 ELSE 00337 * 00338 * Path 2 - underdetermined 00339 * 00340 * Compute space needed for SGEBRD 00341 CALL SGEBRD( M, N, A, LDA, S, DUM(1), DUM(1), 00342 $ DUM(1), DUM(1), -1, INFO ) 00343 LWORK_SGEBRD=DUM(1) 00344 * Compute space needed for SORMBR 00345 CALL SORMBR( 'Q', 'L', 'T', M, NRHS, M, A, LDA, 00346 $ DUM(1), B, LDB, DUM(1), -1, INFO ) 00347 LWORK_SORMBR=DUM(1) 00348 * Compute space needed for SORGBR 00349 CALL SORGBR( 'P', M, N, M, A, LDA, DUM(1), 00350 $ DUM(1), -1, INFO ) 00351 LWORK_SORGBR=DUM(1) 00352 MAXWRK = 3*M + LWORK_SGEBRD 00353 MAXWRK = MAX( MAXWRK, 3*M + LWORK_SORMBR ) 00354 MAXWRK = MAX( MAXWRK, 3*M + LWORK_SORGBR ) 00355 MAXWRK = MAX( MAXWRK, BDSPAC ) 00356 MAXWRK = MAX( MAXWRK, N*NRHS ) 00357 END IF 00358 END IF 00359 MAXWRK = MAX( MINWRK, MAXWRK ) 00360 END IF 00361 WORK( 1 ) = MAXWRK 00362 * 00363 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) 00364 $ INFO = -12 00365 END IF 00366 * 00367 IF( INFO.NE.0 ) THEN 00368 CALL XERBLA( 'SGELSS', -INFO ) 00369 RETURN 00370 ELSE IF( LQUERY ) THEN 00371 RETURN 00372 END IF 00373 * 00374 * Quick return if possible 00375 * 00376 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 00377 RANK = 0 00378 RETURN 00379 END IF 00380 * 00381 * Get machine parameters 00382 * 00383 EPS = SLAMCH( 'P' ) 00384 SFMIN = SLAMCH( 'S' ) 00385 SMLNUM = SFMIN / EPS 00386 BIGNUM = ONE / SMLNUM 00387 CALL SLABAD( SMLNUM, BIGNUM ) 00388 * 00389 * Scale A if max element outside range [SMLNUM,BIGNUM] 00390 * 00391 ANRM = SLANGE( 'M', M, N, A, LDA, WORK ) 00392 IASCL = 0 00393 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00394 * 00395 * Scale matrix norm up to SMLNUM 00396 * 00397 CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO ) 00398 IASCL = 1 00399 ELSE IF( ANRM.GT.BIGNUM ) THEN 00400 * 00401 * Scale matrix norm down to BIGNUM 00402 * 00403 CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO ) 00404 IASCL = 2 00405 ELSE IF( ANRM.EQ.ZERO ) THEN 00406 * 00407 * Matrix all zero. Return zero solution. 00408 * 00409 CALL SLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB ) 00410 CALL SLASET( 'F', MINMN, 1, ZERO, ZERO, S, 1 ) 00411 RANK = 0 00412 GO TO 70 00413 END IF 00414 * 00415 * Scale B if max element outside range [SMLNUM,BIGNUM] 00416 * 00417 BNRM = SLANGE( 'M', M, NRHS, B, LDB, WORK ) 00418 IBSCL = 0 00419 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 00420 * 00421 * Scale matrix norm up to SMLNUM 00422 * 00423 CALL SLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO ) 00424 IBSCL = 1 00425 ELSE IF( BNRM.GT.BIGNUM ) THEN 00426 * 00427 * Scale matrix norm down to BIGNUM 00428 * 00429 CALL SLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO ) 00430 IBSCL = 2 00431 END IF 00432 * 00433 * Overdetermined case 00434 * 00435 IF( M.GE.N ) THEN 00436 * 00437 * Path 1 - overdetermined or exactly determined 00438 * 00439 MM = M 00440 IF( M.GE.MNTHR ) THEN 00441 * 00442 * Path 1a - overdetermined, with many more rows than columns 00443 * 00444 MM = N 00445 ITAU = 1 00446 IWORK = ITAU + N 00447 * 00448 * Compute A=Q*R 00449 * (Workspace: need 2*N, prefer N+N*NB) 00450 * 00451 CALL SGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ), 00452 $ LWORK-IWORK+1, INFO ) 00453 * 00454 * Multiply B by transpose(Q) 00455 * (Workspace: need N+NRHS, prefer N+NRHS*NB) 00456 * 00457 CALL SORMQR( 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAU ), B, 00458 $ LDB, WORK( IWORK ), LWORK-IWORK+1, INFO ) 00459 * 00460 * Zero out below R 00461 * 00462 IF( N.GT.1 ) 00463 $ CALL SLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA ) 00464 END IF 00465 * 00466 IE = 1 00467 ITAUQ = IE + N 00468 ITAUP = ITAUQ + N 00469 IWORK = ITAUP + N 00470 * 00471 * Bidiagonalize R in A 00472 * (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB) 00473 * 00474 CALL SGEBRD( MM, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 00475 $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, 00476 $ INFO ) 00477 * 00478 * Multiply B by transpose of left bidiagonalizing vectors of R 00479 * (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB) 00480 * 00481 CALL SORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, WORK( ITAUQ ), 00482 $ B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO ) 00483 * 00484 * Generate right bidiagonalizing vectors of R in A 00485 * (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) 00486 * 00487 CALL SORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), 00488 $ WORK( IWORK ), LWORK-IWORK+1, INFO ) 00489 IWORK = IE + N 00490 * 00491 * Perform bidiagonal QR iteration 00492 * multiply B by transpose of left singular vectors 00493 * compute right singular vectors in A 00494 * (Workspace: need BDSPAC) 00495 * 00496 CALL SBDSQR( 'U', N, N, 0, NRHS, S, WORK( IE ), A, LDA, DUM, 00497 $ 1, B, LDB, WORK( IWORK ), INFO ) 00498 IF( INFO.NE.0 ) 00499 $ GO TO 70 00500 * 00501 * Multiply B by reciprocals of singular values 00502 * 00503 THR = MAX( RCOND*S( 1 ), SFMIN ) 00504 IF( RCOND.LT.ZERO ) 00505 $ THR = MAX( EPS*S( 1 ), SFMIN ) 00506 RANK = 0 00507 DO 10 I = 1, N 00508 IF( S( I ).GT.THR ) THEN 00509 CALL SRSCL( NRHS, S( I ), B( I, 1 ), LDB ) 00510 RANK = RANK + 1 00511 ELSE 00512 CALL SLASET( 'F', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB ) 00513 END IF 00514 10 CONTINUE 00515 * 00516 * Multiply B by right singular vectors 00517 * (Workspace: need N, prefer N*NRHS) 00518 * 00519 IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN 00520 CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, A, LDA, B, LDB, ZERO, 00521 $ WORK, LDB ) 00522 CALL SLACPY( 'G', N, NRHS, WORK, LDB, B, LDB ) 00523 ELSE IF( NRHS.GT.1 ) THEN 00524 CHUNK = LWORK / N 00525 DO 20 I = 1, NRHS, CHUNK 00526 BL = MIN( NRHS-I+1, CHUNK ) 00527 CALL SGEMM( 'T', 'N', N, BL, N, ONE, A, LDA, B( 1, I ), 00528 $ LDB, ZERO, WORK, N ) 00529 CALL SLACPY( 'G', N, BL, WORK, N, B( 1, I ), LDB ) 00530 20 CONTINUE 00531 ELSE 00532 CALL SGEMV( 'T', N, N, ONE, A, LDA, B, 1, ZERO, WORK, 1 ) 00533 CALL SCOPY( N, WORK, 1, B, 1 ) 00534 END IF 00535 * 00536 ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+ 00537 $ MAX( M, 2*M-4, NRHS, N-3*M ) ) THEN 00538 * 00539 * Path 2a - underdetermined, with many more columns than rows 00540 * and sufficient workspace for an efficient algorithm 00541 * 00542 LDWORK = M 00543 IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ), 00544 $ M*LDA+M+M*NRHS ) )LDWORK = LDA 00545 ITAU = 1 00546 IWORK = M + 1 00547 * 00548 * Compute A=L*Q 00549 * (Workspace: need 2*M, prefer M+M*NB) 00550 * 00551 CALL SGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ), 00552 $ LWORK-IWORK+1, INFO ) 00553 IL = IWORK 00554 * 00555 * Copy L to WORK(IL), zeroing out above it 00556 * 00557 CALL SLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK ) 00558 CALL SLASET( 'U', M-1, M-1, ZERO, ZERO, WORK( IL+LDWORK ), 00559 $ LDWORK ) 00560 IE = IL + LDWORK*M 00561 ITAUQ = IE + M 00562 ITAUP = ITAUQ + M 00563 IWORK = ITAUP + M 00564 * 00565 * Bidiagonalize L in WORK(IL) 00566 * (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB) 00567 * 00568 CALL SGEBRD( M, M, WORK( IL ), LDWORK, S, WORK( IE ), 00569 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( IWORK ), 00570 $ LWORK-IWORK+1, INFO ) 00571 * 00572 * Multiply B by transpose of left bidiagonalizing vectors of L 00573 * (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB) 00574 * 00575 CALL SORMBR( 'Q', 'L', 'T', M, NRHS, M, WORK( IL ), LDWORK, 00576 $ WORK( ITAUQ ), B, LDB, WORK( IWORK ), 00577 $ LWORK-IWORK+1, INFO ) 00578 * 00579 * Generate right bidiagonalizing vectors of R in WORK(IL) 00580 * (Workspace: need M*M+5*M-1, prefer M*M+4*M+(M-1)*NB) 00581 * 00582 CALL SORGBR( 'P', M, M, M, WORK( IL ), LDWORK, WORK( ITAUP ), 00583 $ WORK( IWORK ), LWORK-IWORK+1, INFO ) 00584 IWORK = IE + M 00585 * 00586 * Perform bidiagonal QR iteration, 00587 * computing right singular vectors of L in WORK(IL) and 00588 * multiplying B by transpose of left singular vectors 00589 * (Workspace: need M*M+M+BDSPAC) 00590 * 00591 CALL SBDSQR( 'U', M, M, 0, NRHS, S, WORK( IE ), WORK( IL ), 00592 $ LDWORK, A, LDA, B, LDB, WORK( IWORK ), INFO ) 00593 IF( INFO.NE.0 ) 00594 $ GO TO 70 00595 * 00596 * Multiply B by reciprocals of singular values 00597 * 00598 THR = MAX( RCOND*S( 1 ), SFMIN ) 00599 IF( RCOND.LT.ZERO ) 00600 $ THR = MAX( EPS*S( 1 ), SFMIN ) 00601 RANK = 0 00602 DO 30 I = 1, M 00603 IF( S( I ).GT.THR ) THEN 00604 CALL SRSCL( NRHS, S( I ), B( I, 1 ), LDB ) 00605 RANK = RANK + 1 00606 ELSE 00607 CALL SLASET( 'F', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB ) 00608 END IF 00609 30 CONTINUE 00610 IWORK = IE 00611 * 00612 * Multiply B by right singular vectors of L in WORK(IL) 00613 * (Workspace: need M*M+2*M, prefer M*M+M+M*NRHS) 00614 * 00615 IF( LWORK.GE.LDB*NRHS+IWORK-1 .AND. NRHS.GT.1 ) THEN 00616 CALL SGEMM( 'T', 'N', M, NRHS, M, ONE, WORK( IL ), LDWORK, 00617 $ B, LDB, ZERO, WORK( IWORK ), LDB ) 00618 CALL SLACPY( 'G', M, NRHS, WORK( IWORK ), LDB, B, LDB ) 00619 ELSE IF( NRHS.GT.1 ) THEN 00620 CHUNK = ( LWORK-IWORK+1 ) / M 00621 DO 40 I = 1, NRHS, CHUNK 00622 BL = MIN( NRHS-I+1, CHUNK ) 00623 CALL SGEMM( 'T', 'N', M, BL, M, ONE, WORK( IL ), LDWORK, 00624 $ B( 1, I ), LDB, ZERO, WORK( IWORK ), M ) 00625 CALL SLACPY( 'G', M, BL, WORK( IWORK ), M, B( 1, I ), 00626 $ LDB ) 00627 40 CONTINUE 00628 ELSE 00629 CALL SGEMV( 'T', M, M, ONE, WORK( IL ), LDWORK, B( 1, 1 ), 00630 $ 1, ZERO, WORK( IWORK ), 1 ) 00631 CALL SCOPY( M, WORK( IWORK ), 1, B( 1, 1 ), 1 ) 00632 END IF 00633 * 00634 * Zero out below first M rows of B 00635 * 00636 CALL SLASET( 'F', N-M, NRHS, ZERO, ZERO, B( M+1, 1 ), LDB ) 00637 IWORK = ITAU + M 00638 * 00639 * Multiply transpose(Q) by B 00640 * (Workspace: need M+NRHS, prefer M+NRHS*NB) 00641 * 00642 CALL SORMLQ( 'L', 'T', N, NRHS, M, A, LDA, WORK( ITAU ), B, 00643 $ LDB, WORK( IWORK ), LWORK-IWORK+1, INFO ) 00644 * 00645 ELSE 00646 * 00647 * Path 2 - remaining underdetermined cases 00648 * 00649 IE = 1 00650 ITAUQ = IE + M 00651 ITAUP = ITAUQ + M 00652 IWORK = ITAUP + M 00653 * 00654 * Bidiagonalize A 00655 * (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB) 00656 * 00657 CALL SGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 00658 $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, 00659 $ INFO ) 00660 * 00661 * Multiply B by transpose of left bidiagonalizing vectors 00662 * (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB) 00663 * 00664 CALL SORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAUQ ), 00665 $ B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO ) 00666 * 00667 * Generate right bidiagonalizing vectors in A 00668 * (Workspace: need 4*M, prefer 3*M+M*NB) 00669 * 00670 CALL SORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), 00671 $ WORK( IWORK ), LWORK-IWORK+1, INFO ) 00672 IWORK = IE + M 00673 * 00674 * Perform bidiagonal QR iteration, 00675 * computing right singular vectors of A in A and 00676 * multiplying B by transpose of left singular vectors 00677 * (Workspace: need BDSPAC) 00678 * 00679 CALL SBDSQR( 'L', M, N, 0, NRHS, S, WORK( IE ), A, LDA, DUM, 00680 $ 1, B, LDB, WORK( IWORK ), INFO ) 00681 IF( INFO.NE.0 ) 00682 $ GO TO 70 00683 * 00684 * Multiply B by reciprocals of singular values 00685 * 00686 THR = MAX( RCOND*S( 1 ), SFMIN ) 00687 IF( RCOND.LT.ZERO ) 00688 $ THR = MAX( EPS*S( 1 ), SFMIN ) 00689 RANK = 0 00690 DO 50 I = 1, M 00691 IF( S( I ).GT.THR ) THEN 00692 CALL SRSCL( NRHS, S( I ), B( I, 1 ), LDB ) 00693 RANK = RANK + 1 00694 ELSE 00695 CALL SLASET( 'F', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB ) 00696 END IF 00697 50 CONTINUE 00698 * 00699 * Multiply B by right singular vectors of A 00700 * (Workspace: need N, prefer N*NRHS) 00701 * 00702 IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN 00703 CALL SGEMM( 'T', 'N', N, NRHS, M, ONE, A, LDA, B, LDB, ZERO, 00704 $ WORK, LDB ) 00705 CALL SLACPY( 'F', N, NRHS, WORK, LDB, B, LDB ) 00706 ELSE IF( NRHS.GT.1 ) THEN 00707 CHUNK = LWORK / N 00708 DO 60 I = 1, NRHS, CHUNK 00709 BL = MIN( NRHS-I+1, CHUNK ) 00710 CALL SGEMM( 'T', 'N', N, BL, M, ONE, A, LDA, B( 1, I ), 00711 $ LDB, ZERO, WORK, N ) 00712 CALL SLACPY( 'F', N, BL, WORK, N, B( 1, I ), LDB ) 00713 60 CONTINUE 00714 ELSE 00715 CALL SGEMV( 'T', M, N, ONE, A, LDA, B, 1, ZERO, WORK, 1 ) 00716 CALL SCOPY( N, WORK, 1, B, 1 ) 00717 END IF 00718 END IF 00719 * 00720 * Undo scaling 00721 * 00722 IF( IASCL.EQ.1 ) THEN 00723 CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO ) 00724 CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN, 00725 $ INFO ) 00726 ELSE IF( IASCL.EQ.2 ) THEN 00727 CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO ) 00728 CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN, 00729 $ INFO ) 00730 END IF 00731 IF( IBSCL.EQ.1 ) THEN 00732 CALL SLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO ) 00733 ELSE IF( IBSCL.EQ.2 ) THEN 00734 CALL SLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO ) 00735 END IF 00736 * 00737 70 CONTINUE 00738 WORK( 1 ) = MAXWRK 00739 RETURN 00740 * 00741 * End of SGELSS 00742 * 00743 END