![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> DGELSS 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 DGELSS + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgelss.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgelss.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgelss.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DGELSS( 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 * DOUBLE PRECISION RCOND 00027 * .. 00028 * .. Array Arguments .. 00029 * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * ) 00030 * .. 00031 * 00032 * 00033 *> \par Purpose: 00034 * ============= 00035 *> 00036 *> \verbatim 00037 *> 00038 *> DGELSS 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 doubleGEsolve 00170 * 00171 * ===================================================================== 00172 SUBROUTINE DGELSS( 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 DOUBLE PRECISION RCOND 00183 * .. 00184 * .. Array Arguments .. 00185 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * ) 00186 * .. 00187 * 00188 * ===================================================================== 00189 * 00190 * .. Parameters .. 00191 DOUBLE PRECISION ZERO, ONE 00192 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+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_DGEQRF, LWORK_DORMQR, LWORK_DGEBRD, 00200 $ LWORK_DORMBR, LWORK_DORGBR, LWORK_DORMLQ, 00201 $ LWORK_DGELQF 00202 DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR 00203 * .. 00204 * .. Local Arrays .. 00205 DOUBLE PRECISION DUM( 1 ) 00206 * .. 00207 * .. External Subroutines .. 00208 EXTERNAL DBDSQR, DCOPY, DGEBRD, DGELQF, DGEMM, DGEMV, 00209 $ DGEQRF, DLABAD, DLACPY, DLASCL, DLASET, DORGBR, 00210 $ DORMBR, DORMLQ, DORMQR, DRSCL, XERBLA 00211 * .. 00212 * .. External Functions .. 00213 INTEGER ILAENV 00214 DOUBLE PRECISION DLAMCH, DLANGE 00215 EXTERNAL ILAENV, DLAMCH, DLANGE 00216 * .. 00217 * .. Intrinsic Functions .. 00218 INTRINSIC MAX, MIN 00219 * .. 00220 * .. Executable Statements .. 00221 * 00222 * Test the input arguments 00223 * 00224 INFO = 0 00225 MINMN = MIN( M, N ) 00226 MAXMN = MAX( M, N ) 00227 LQUERY = ( LWORK.EQ.-1 ) 00228 IF( M.LT.0 ) THEN 00229 INFO = -1 00230 ELSE IF( N.LT.0 ) THEN 00231 INFO = -2 00232 ELSE IF( NRHS.LT.0 ) THEN 00233 INFO = -3 00234 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00235 INFO = -5 00236 ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN 00237 INFO = -7 00238 END IF 00239 * 00240 * Compute workspace 00241 * (Note: Comments in the code beginning "Workspace:" describe the 00242 * minimal amount of workspace needed at that point in the code, 00243 * as well as the preferred amount for good performance. 00244 * NB refers to the optimal block size for the immediately 00245 * following subroutine, as returned by ILAENV.) 00246 * 00247 IF( INFO.EQ.0 ) THEN 00248 MINWRK = 1 00249 MAXWRK = 1 00250 IF( MINMN.GT.0 ) THEN 00251 MM = M 00252 MNTHR = ILAENV( 6, 'DGELSS', ' ', M, N, NRHS, -1 ) 00253 IF( M.GE.N .AND. M.GE.MNTHR ) THEN 00254 * 00255 * Path 1a - overdetermined, with many more rows than 00256 * columns 00257 * 00258 * Compute space needed for DGEQRF 00259 CALL DGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, INFO ) 00260 LWORK_DGEQRF=DUM(1) 00261 * Compute space needed for DORMQR 00262 CALL DORMQR( 'L', 'T', M, NRHS, N, A, LDA, DUM(1), B, 00263 $ LDB, DUM(1), -1, INFO ) 00264 LWORK_DORMQR=DUM(1) 00265 MM = N 00266 MAXWRK = MAX( MAXWRK, N + LWORK_DGEQRF ) 00267 MAXWRK = MAX( MAXWRK, N + LWORK_DORMQR ) 00268 END IF 00269 IF( M.GE.N ) THEN 00270 * 00271 * Path 1 - overdetermined or exactly determined 00272 * 00273 * Compute workspace needed for DBDSQR 00274 * 00275 BDSPAC = MAX( 1, 5*N ) 00276 * Compute space needed for DGEBRD 00277 CALL DGEBRD( MM, N, A, LDA, S, DUM(1), DUM(1), 00278 $ DUM(1), DUM(1), -1, INFO ) 00279 LWORK_DGEBRD=DUM(1) 00280 * Compute space needed for DORMBR 00281 CALL DORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, DUM(1), 00282 $ B, LDB, DUM(1), -1, INFO ) 00283 LWORK_DORMBR=DUM(1) 00284 * Compute space needed for DORGBR 00285 CALL DORGBR( 'P', N, N, N, A, LDA, DUM(1), 00286 $ DUM(1), -1, INFO ) 00287 LWORK_DORGBR=DUM(1) 00288 * Compute total workspace needed 00289 MAXWRK = MAX( MAXWRK, 3*N + LWORK_DGEBRD ) 00290 MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORMBR ) 00291 MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR ) 00292 MAXWRK = MAX( MAXWRK, BDSPAC ) 00293 MAXWRK = MAX( MAXWRK, N*NRHS ) 00294 MINWRK = MAX( 3*N + MM, 3*N + NRHS, BDSPAC ) 00295 MAXWRK = MAX( MINWRK, MAXWRK ) 00296 END IF 00297 IF( N.GT.M ) THEN 00298 * 00299 * Compute workspace needed for DBDSQR 00300 * 00301 BDSPAC = MAX( 1, 5*M ) 00302 MINWRK = MAX( 3*M+NRHS, 3*M+N, BDSPAC ) 00303 IF( N.GE.MNTHR ) THEN 00304 * 00305 * Path 2a - underdetermined, with many more columns 00306 * than rows 00307 * 00308 * Compute space needed for DGELQF 00309 CALL DGELQF( M, N, A, LDA, DUM(1), DUM(1), 00310 $ -1, INFO ) 00311 LWORK_DGELQF=DUM(1) 00312 * Compute space needed for DGEBRD 00313 CALL DGEBRD( M, M, A, LDA, S, DUM(1), DUM(1), 00314 $ DUM(1), DUM(1), -1, INFO ) 00315 LWORK_DGEBRD=DUM(1) 00316 * Compute space needed for DORMBR 00317 CALL DORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA, 00318 $ DUM(1), B, LDB, DUM(1), -1, INFO ) 00319 LWORK_DORMBR=DUM(1) 00320 * Compute space needed for DORGBR 00321 CALL DORGBR( 'P', M, M, M, A, LDA, DUM(1), 00322 $ DUM(1), -1, INFO ) 00323 LWORK_DORGBR=DUM(1) 00324 * Compute space needed for DORMLQ 00325 CALL DORMLQ( 'L', 'T', N, NRHS, M, A, LDA, DUM(1), 00326 $ B, LDB, DUM(1), -1, INFO ) 00327 LWORK_DORMLQ=DUM(1) 00328 * Compute total workspace needed 00329 MAXWRK = M + LWORK_DGELQF 00330 MAXWRK = MAX( MAXWRK, M*M + 4*M + LWORK_DGEBRD ) 00331 MAXWRK = MAX( MAXWRK, M*M + 4*M + LWORK_DORMBR ) 00332 MAXWRK = MAX( MAXWRK, M*M + 4*M + LWORK_DORGBR ) 00333 MAXWRK = MAX( MAXWRK, M*M + M + BDSPAC ) 00334 IF( NRHS.GT.1 ) THEN 00335 MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS ) 00336 ELSE 00337 MAXWRK = MAX( MAXWRK, M*M + 2*M ) 00338 END IF 00339 MAXWRK = MAX( MAXWRK, M + LWORK_DORMLQ ) 00340 ELSE 00341 * 00342 * Path 2 - underdetermined 00343 * 00344 * Compute space needed for DGEBRD 00345 CALL DGEBRD( M, N, A, LDA, S, DUM(1), DUM(1), 00346 $ DUM(1), DUM(1), -1, INFO ) 00347 LWORK_DGEBRD=DUM(1) 00348 * Compute space needed for DORMBR 00349 CALL DORMBR( 'Q', 'L', 'T', M, NRHS, M, A, LDA, 00350 $ DUM(1), B, LDB, DUM(1), -1, INFO ) 00351 LWORK_DORMBR=DUM(1) 00352 * Compute space needed for DORGBR 00353 CALL DORGBR( 'P', M, N, M, A, LDA, DUM(1), 00354 $ DUM(1), -1, INFO ) 00355 LWORK_DORGBR=DUM(1) 00356 MAXWRK = 3*M + LWORK_DGEBRD 00357 MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORMBR ) 00358 MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR ) 00359 MAXWRK = MAX( MAXWRK, BDSPAC ) 00360 MAXWRK = MAX( MAXWRK, N*NRHS ) 00361 END IF 00362 END IF 00363 MAXWRK = MAX( MINWRK, MAXWRK ) 00364 END IF 00365 WORK( 1 ) = MAXWRK 00366 * 00367 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) 00368 $ INFO = -12 00369 END IF 00370 * 00371 IF( INFO.NE.0 ) THEN 00372 CALL XERBLA( 'DGELSS', -INFO ) 00373 RETURN 00374 ELSE IF( LQUERY ) THEN 00375 RETURN 00376 END IF 00377 * 00378 * Quick return if possible 00379 * 00380 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 00381 RANK = 0 00382 RETURN 00383 END IF 00384 * 00385 * Get machine parameters 00386 * 00387 EPS = DLAMCH( 'P' ) 00388 SFMIN = DLAMCH( 'S' ) 00389 SMLNUM = SFMIN / EPS 00390 BIGNUM = ONE / SMLNUM 00391 CALL DLABAD( SMLNUM, BIGNUM ) 00392 * 00393 * Scale A if max element outside range [SMLNUM,BIGNUM] 00394 * 00395 ANRM = DLANGE( 'M', M, N, A, LDA, WORK ) 00396 IASCL = 0 00397 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00398 * 00399 * Scale matrix norm up to SMLNUM 00400 * 00401 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO ) 00402 IASCL = 1 00403 ELSE IF( ANRM.GT.BIGNUM ) THEN 00404 * 00405 * Scale matrix norm down to BIGNUM 00406 * 00407 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO ) 00408 IASCL = 2 00409 ELSE IF( ANRM.EQ.ZERO ) THEN 00410 * 00411 * Matrix all zero. Return zero solution. 00412 * 00413 CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB ) 00414 CALL DLASET( 'F', MINMN, 1, ZERO, ZERO, S, MINMN ) 00415 RANK = 0 00416 GO TO 70 00417 END IF 00418 * 00419 * Scale B if max element outside range [SMLNUM,BIGNUM] 00420 * 00421 BNRM = DLANGE( 'M', M, NRHS, B, LDB, WORK ) 00422 IBSCL = 0 00423 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 00424 * 00425 * Scale matrix norm up to SMLNUM 00426 * 00427 CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO ) 00428 IBSCL = 1 00429 ELSE IF( BNRM.GT.BIGNUM ) THEN 00430 * 00431 * Scale matrix norm down to BIGNUM 00432 * 00433 CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO ) 00434 IBSCL = 2 00435 END IF 00436 * 00437 * Overdetermined case 00438 * 00439 IF( M.GE.N ) THEN 00440 * 00441 * Path 1 - overdetermined or exactly determined 00442 * 00443 MM = M 00444 IF( M.GE.MNTHR ) THEN 00445 * 00446 * Path 1a - overdetermined, with many more rows than columns 00447 * 00448 MM = N 00449 ITAU = 1 00450 IWORK = ITAU + N 00451 * 00452 * Compute A=Q*R 00453 * (Workspace: need 2*N, prefer N+N*NB) 00454 * 00455 CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ), 00456 $ LWORK-IWORK+1, INFO ) 00457 * 00458 * Multiply B by transpose(Q) 00459 * (Workspace: need N+NRHS, prefer N+NRHS*NB) 00460 * 00461 CALL DORMQR( 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAU ), B, 00462 $ LDB, WORK( IWORK ), LWORK-IWORK+1, INFO ) 00463 * 00464 * Zero out below R 00465 * 00466 IF( N.GT.1 ) 00467 $ CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA ) 00468 END IF 00469 * 00470 IE = 1 00471 ITAUQ = IE + N 00472 ITAUP = ITAUQ + N 00473 IWORK = ITAUP + N 00474 * 00475 * Bidiagonalize R in A 00476 * (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB) 00477 * 00478 CALL DGEBRD( MM, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 00479 $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, 00480 $ INFO ) 00481 * 00482 * Multiply B by transpose of left bidiagonalizing vectors of R 00483 * (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB) 00484 * 00485 CALL DORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, WORK( ITAUQ ), 00486 $ B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO ) 00487 * 00488 * Generate right bidiagonalizing vectors of R in A 00489 * (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) 00490 * 00491 CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), 00492 $ WORK( IWORK ), LWORK-IWORK+1, INFO ) 00493 IWORK = IE + N 00494 * 00495 * Perform bidiagonal QR iteration 00496 * multiply B by transpose of left singular vectors 00497 * compute right singular vectors in A 00498 * (Workspace: need BDSPAC) 00499 * 00500 CALL DBDSQR( 'U', N, N, 0, NRHS, S, WORK( IE ), A, LDA, DUM, 00501 $ 1, B, LDB, WORK( IWORK ), INFO ) 00502 IF( INFO.NE.0 ) 00503 $ GO TO 70 00504 * 00505 * Multiply B by reciprocals of singular values 00506 * 00507 THR = MAX( RCOND*S( 1 ), SFMIN ) 00508 IF( RCOND.LT.ZERO ) 00509 $ THR = MAX( EPS*S( 1 ), SFMIN ) 00510 RANK = 0 00511 DO 10 I = 1, N 00512 IF( S( I ).GT.THR ) THEN 00513 CALL DRSCL( NRHS, S( I ), B( I, 1 ), LDB ) 00514 RANK = RANK + 1 00515 ELSE 00516 CALL DLASET( 'F', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB ) 00517 END IF 00518 10 CONTINUE 00519 * 00520 * Multiply B by right singular vectors 00521 * (Workspace: need N, prefer N*NRHS) 00522 * 00523 IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN 00524 CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, A, LDA, B, LDB, ZERO, 00525 $ WORK, LDB ) 00526 CALL DLACPY( 'G', N, NRHS, WORK, LDB, B, LDB ) 00527 ELSE IF( NRHS.GT.1 ) THEN 00528 CHUNK = LWORK / N 00529 DO 20 I = 1, NRHS, CHUNK 00530 BL = MIN( NRHS-I+1, CHUNK ) 00531 CALL DGEMM( 'T', 'N', N, BL, N, ONE, A, LDA, B( 1, I ), 00532 $ LDB, ZERO, WORK, N ) 00533 CALL DLACPY( 'G', N, BL, WORK, N, B( 1, I ), LDB ) 00534 20 CONTINUE 00535 ELSE 00536 CALL DGEMV( 'T', N, N, ONE, A, LDA, B, 1, ZERO, WORK, 1 ) 00537 CALL DCOPY( N, WORK, 1, B, 1 ) 00538 END IF 00539 * 00540 ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+ 00541 $ MAX( M, 2*M-4, NRHS, N-3*M ) ) THEN 00542 * 00543 * Path 2a - underdetermined, with many more columns than rows 00544 * and sufficient workspace for an efficient algorithm 00545 * 00546 LDWORK = M 00547 IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ), 00548 $ M*LDA+M+M*NRHS ) )LDWORK = LDA 00549 ITAU = 1 00550 IWORK = M + 1 00551 * 00552 * Compute A=L*Q 00553 * (Workspace: need 2*M, prefer M+M*NB) 00554 * 00555 CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ), 00556 $ LWORK-IWORK+1, INFO ) 00557 IL = IWORK 00558 * 00559 * Copy L to WORK(IL), zeroing out above it 00560 * 00561 CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK ) 00562 CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, WORK( IL+LDWORK ), 00563 $ LDWORK ) 00564 IE = IL + LDWORK*M 00565 ITAUQ = IE + M 00566 ITAUP = ITAUQ + M 00567 IWORK = ITAUP + M 00568 * 00569 * Bidiagonalize L in WORK(IL) 00570 * (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB) 00571 * 00572 CALL DGEBRD( M, M, WORK( IL ), LDWORK, S, WORK( IE ), 00573 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( IWORK ), 00574 $ LWORK-IWORK+1, INFO ) 00575 * 00576 * Multiply B by transpose of left bidiagonalizing vectors of L 00577 * (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB) 00578 * 00579 CALL DORMBR( 'Q', 'L', 'T', M, NRHS, M, WORK( IL ), LDWORK, 00580 $ WORK( ITAUQ ), B, LDB, WORK( IWORK ), 00581 $ LWORK-IWORK+1, INFO ) 00582 * 00583 * Generate right bidiagonalizing vectors of R in WORK(IL) 00584 * (Workspace: need M*M+5*M-1, prefer M*M+4*M+(M-1)*NB) 00585 * 00586 CALL DORGBR( 'P', M, M, M, WORK( IL ), LDWORK, WORK( ITAUP ), 00587 $ WORK( IWORK ), LWORK-IWORK+1, INFO ) 00588 IWORK = IE + M 00589 * 00590 * Perform bidiagonal QR iteration, 00591 * computing right singular vectors of L in WORK(IL) and 00592 * multiplying B by transpose of left singular vectors 00593 * (Workspace: need M*M+M+BDSPAC) 00594 * 00595 CALL DBDSQR( 'U', M, M, 0, NRHS, S, WORK( IE ), WORK( IL ), 00596 $ LDWORK, A, LDA, B, LDB, WORK( IWORK ), INFO ) 00597 IF( INFO.NE.0 ) 00598 $ GO TO 70 00599 * 00600 * Multiply B by reciprocals of singular values 00601 * 00602 THR = MAX( RCOND*S( 1 ), SFMIN ) 00603 IF( RCOND.LT.ZERO ) 00604 $ THR = MAX( EPS*S( 1 ), SFMIN ) 00605 RANK = 0 00606 DO 30 I = 1, M 00607 IF( S( I ).GT.THR ) THEN 00608 CALL DRSCL( NRHS, S( I ), B( I, 1 ), LDB ) 00609 RANK = RANK + 1 00610 ELSE 00611 CALL DLASET( 'F', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB ) 00612 END IF 00613 30 CONTINUE 00614 IWORK = IE 00615 * 00616 * Multiply B by right singular vectors of L in WORK(IL) 00617 * (Workspace: need M*M+2*M, prefer M*M+M+M*NRHS) 00618 * 00619 IF( LWORK.GE.LDB*NRHS+IWORK-1 .AND. NRHS.GT.1 ) THEN 00620 CALL DGEMM( 'T', 'N', M, NRHS, M, ONE, WORK( IL ), LDWORK, 00621 $ B, LDB, ZERO, WORK( IWORK ), LDB ) 00622 CALL DLACPY( 'G', M, NRHS, WORK( IWORK ), LDB, B, LDB ) 00623 ELSE IF( NRHS.GT.1 ) THEN 00624 CHUNK = ( LWORK-IWORK+1 ) / M 00625 DO 40 I = 1, NRHS, CHUNK 00626 BL = MIN( NRHS-I+1, CHUNK ) 00627 CALL DGEMM( 'T', 'N', M, BL, M, ONE, WORK( IL ), LDWORK, 00628 $ B( 1, I ), LDB, ZERO, WORK( IWORK ), M ) 00629 CALL DLACPY( 'G', M, BL, WORK( IWORK ), M, B( 1, I ), 00630 $ LDB ) 00631 40 CONTINUE 00632 ELSE 00633 CALL DGEMV( 'T', M, M, ONE, WORK( IL ), LDWORK, B( 1, 1 ), 00634 $ 1, ZERO, WORK( IWORK ), 1 ) 00635 CALL DCOPY( M, WORK( IWORK ), 1, B( 1, 1 ), 1 ) 00636 END IF 00637 * 00638 * Zero out below first M rows of B 00639 * 00640 CALL DLASET( 'F', N-M, NRHS, ZERO, ZERO, B( M+1, 1 ), LDB ) 00641 IWORK = ITAU + M 00642 * 00643 * Multiply transpose(Q) by B 00644 * (Workspace: need M+NRHS, prefer M+NRHS*NB) 00645 * 00646 CALL DORMLQ( 'L', 'T', N, NRHS, M, A, LDA, WORK( ITAU ), B, 00647 $ LDB, WORK( IWORK ), LWORK-IWORK+1, INFO ) 00648 * 00649 ELSE 00650 * 00651 * Path 2 - remaining underdetermined cases 00652 * 00653 IE = 1 00654 ITAUQ = IE + M 00655 ITAUP = ITAUQ + M 00656 IWORK = ITAUP + M 00657 * 00658 * Bidiagonalize A 00659 * (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB) 00660 * 00661 CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 00662 $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, 00663 $ INFO ) 00664 * 00665 * Multiply B by transpose of left bidiagonalizing vectors 00666 * (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB) 00667 * 00668 CALL DORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAUQ ), 00669 $ B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO ) 00670 * 00671 * Generate right bidiagonalizing vectors in A 00672 * (Workspace: need 4*M, prefer 3*M+M*NB) 00673 * 00674 CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), 00675 $ WORK( IWORK ), LWORK-IWORK+1, INFO ) 00676 IWORK = IE + M 00677 * 00678 * Perform bidiagonal QR iteration, 00679 * computing right singular vectors of A in A and 00680 * multiplying B by transpose of left singular vectors 00681 * (Workspace: need BDSPAC) 00682 * 00683 CALL DBDSQR( 'L', M, N, 0, NRHS, S, WORK( IE ), A, LDA, DUM, 00684 $ 1, B, LDB, WORK( IWORK ), INFO ) 00685 IF( INFO.NE.0 ) 00686 $ GO TO 70 00687 * 00688 * Multiply B by reciprocals of singular values 00689 * 00690 THR = MAX( RCOND*S( 1 ), SFMIN ) 00691 IF( RCOND.LT.ZERO ) 00692 $ THR = MAX( EPS*S( 1 ), SFMIN ) 00693 RANK = 0 00694 DO 50 I = 1, M 00695 IF( S( I ).GT.THR ) THEN 00696 CALL DRSCL( NRHS, S( I ), B( I, 1 ), LDB ) 00697 RANK = RANK + 1 00698 ELSE 00699 CALL DLASET( 'F', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB ) 00700 END IF 00701 50 CONTINUE 00702 * 00703 * Multiply B by right singular vectors of A 00704 * (Workspace: need N, prefer N*NRHS) 00705 * 00706 IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN 00707 CALL DGEMM( 'T', 'N', N, NRHS, M, ONE, A, LDA, B, LDB, ZERO, 00708 $ WORK, LDB ) 00709 CALL DLACPY( 'F', N, NRHS, WORK, LDB, B, LDB ) 00710 ELSE IF( NRHS.GT.1 ) THEN 00711 CHUNK = LWORK / N 00712 DO 60 I = 1, NRHS, CHUNK 00713 BL = MIN( NRHS-I+1, CHUNK ) 00714 CALL DGEMM( 'T', 'N', N, BL, M, ONE, A, LDA, B( 1, I ), 00715 $ LDB, ZERO, WORK, N ) 00716 CALL DLACPY( 'F', N, BL, WORK, N, B( 1, I ), LDB ) 00717 60 CONTINUE 00718 ELSE 00719 CALL DGEMV( 'T', M, N, ONE, A, LDA, B, 1, ZERO, WORK, 1 ) 00720 CALL DCOPY( N, WORK, 1, B, 1 ) 00721 END IF 00722 END IF 00723 * 00724 * Undo scaling 00725 * 00726 IF( IASCL.EQ.1 ) THEN 00727 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO ) 00728 CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN, 00729 $ INFO ) 00730 ELSE IF( IASCL.EQ.2 ) THEN 00731 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO ) 00732 CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN, 00733 $ INFO ) 00734 END IF 00735 IF( IBSCL.EQ.1 ) THEN 00736 CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO ) 00737 ELSE IF( IBSCL.EQ.2 ) THEN 00738 CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO ) 00739 END IF 00740 * 00741 70 CONTINUE 00742 WORK( 1 ) = MAXWRK 00743 RETURN 00744 * 00745 * End of DGELSS 00746 * 00747 END