![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> DGELSD 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 DGELSD + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgelsd.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgelsd.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgelsd.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, 00022 * WORK, LWORK, IWORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK 00026 * DOUBLE PRECISION RCOND 00027 * .. 00028 * .. Array Arguments .. 00029 * INTEGER IWORK( * ) 00030 * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> DGELSD 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 WORK array, returns 00167 *> this value as the first entry of the WORK array, and no error 00168 *> message related to LWORK is issued by XERBLA. 00169 *> \endverbatim 00170 *> 00171 *> \param[out] IWORK 00172 *> \verbatim 00173 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK)) 00174 *> LIWORK >= max(1, 3 * MINMN * NLVL + 11 * MINMN), 00175 *> where MINMN = MIN( M,N ). 00176 *> On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK. 00177 *> \endverbatim 00178 *> 00179 *> \param[out] INFO 00180 *> \verbatim 00181 *> INFO is INTEGER 00182 *> = 0: successful exit 00183 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00184 *> > 0: the algorithm for computing the SVD failed to converge; 00185 *> if INFO = i, i off-diagonal elements of an intermediate 00186 *> bidiagonal form did not converge to zero. 00187 *> \endverbatim 00188 * 00189 * Authors: 00190 * ======== 00191 * 00192 *> \author Univ. of Tennessee 00193 *> \author Univ. of California Berkeley 00194 *> \author Univ. of Colorado Denver 00195 *> \author NAG Ltd. 00196 * 00197 *> \date November 2011 00198 * 00199 *> \ingroup doubleGEsolve 00200 * 00201 *> \par Contributors: 00202 * ================== 00203 *> 00204 *> Ming Gu and Ren-Cang Li, Computer Science Division, University of 00205 *> California at Berkeley, USA \n 00206 *> Osni Marques, LBNL/NERSC, USA \n 00207 * 00208 * ===================================================================== 00209 SUBROUTINE DGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, 00210 $ WORK, LWORK, IWORK, INFO ) 00211 * 00212 * -- LAPACK driver routine (version 3.4.0) -- 00213 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00214 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00215 * November 2011 00216 * 00217 * .. Scalar Arguments .. 00218 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK 00219 DOUBLE PRECISION RCOND 00220 * .. 00221 * .. Array Arguments .. 00222 INTEGER IWORK( * ) 00223 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * ) 00224 * .. 00225 * 00226 * ===================================================================== 00227 * 00228 * .. Parameters .. 00229 DOUBLE PRECISION ZERO, ONE, TWO 00230 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 ) 00231 * .. 00232 * .. Local Scalars .. 00233 LOGICAL LQUERY 00234 INTEGER IASCL, IBSCL, IE, IL, ITAU, ITAUP, ITAUQ, 00235 $ LDWORK, LIWORK, MAXMN, MAXWRK, MINMN, MINWRK, 00236 $ MM, MNTHR, NLVL, NWORK, SMLSIZ, WLALSD 00237 DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM 00238 * .. 00239 * .. External Subroutines .. 00240 EXTERNAL DGEBRD, DGELQF, DGEQRF, DLABAD, DLACPY, DLALSD, 00241 $ DLASCL, DLASET, DORMBR, DORMLQ, DORMQR, XERBLA 00242 * .. 00243 * .. External Functions .. 00244 INTEGER ILAENV 00245 DOUBLE PRECISION DLAMCH, DLANGE 00246 EXTERNAL ILAENV, DLAMCH, DLANGE 00247 * .. 00248 * .. Intrinsic Functions .. 00249 INTRINSIC DBLE, INT, LOG, MAX, MIN 00250 * .. 00251 * .. Executable Statements .. 00252 * 00253 * Test the input arguments. 00254 * 00255 INFO = 0 00256 MINMN = MIN( M, N ) 00257 MAXMN = MAX( M, N ) 00258 MNTHR = ILAENV( 6, 'DGELSD', ' ', M, N, NRHS, -1 ) 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 SMLSIZ = ILAENV( 9, 'DGELSD', ' ', 0, 0, 0, 0 ) 00273 * 00274 * Compute workspace. 00275 * (Note: Comments in the code beginning "Workspace:" describe the 00276 * minimal amount of workspace needed at that point in the code, 00277 * as well as the preferred amount for good performance. 00278 * NB refers to the optimal block size for the immediately 00279 * following subroutine, as returned by ILAENV.) 00280 * 00281 MINWRK = 1 00282 LIWORK = 1 00283 MINMN = MAX( 1, MINMN ) 00284 NLVL = MAX( INT( LOG( DBLE( MINMN ) / DBLE( SMLSIZ+1 ) ) / 00285 $ LOG( TWO ) ) + 1, 0 ) 00286 * 00287 IF( INFO.EQ.0 ) THEN 00288 MAXWRK = 0 00289 LIWORK = 3*MINMN*NLVL + 11*MINMN 00290 MM = M 00291 IF( M.GE.N .AND. M.GE.MNTHR ) THEN 00292 * 00293 * Path 1a - overdetermined, with many more rows than columns. 00294 * 00295 MM = N 00296 MAXWRK = MAX( MAXWRK, N+N*ILAENV( 1, 'DGEQRF', ' ', M, N, 00297 $ -1, -1 ) ) 00298 MAXWRK = MAX( MAXWRK, N+NRHS* 00299 $ ILAENV( 1, 'DORMQR', 'LT', 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 )* 00306 $ ILAENV( 1, 'DGEBRD', ' ', MM, N, -1, -1 ) ) 00307 MAXWRK = MAX( MAXWRK, 3*N+NRHS* 00308 $ ILAENV( 1, 'DORMBR', 'QLT', MM, NRHS, N, -1 ) ) 00309 MAXWRK = MAX( MAXWRK, 3*N+( N-1 )* 00310 $ ILAENV( 1, 'DORMBR', 'PLN', N, NRHS, N, -1 ) ) 00311 WLALSD = 9*N+2*N*SMLSIZ+8*N*NLVL+N*NRHS+(SMLSIZ+1)**2 00312 MAXWRK = MAX( MAXWRK, 3*N+WLALSD ) 00313 MINWRK = MAX( 3*N+MM, 3*N+NRHS, 3*N+WLALSD ) 00314 END IF 00315 IF( N.GT.M ) THEN 00316 WLALSD = 9*M+2*M*SMLSIZ+8*M*NLVL+M*NRHS+(SMLSIZ+1)**2 00317 IF( N.GE.MNTHR ) THEN 00318 * 00319 * Path 2a - underdetermined, with many more columns 00320 * than rows. 00321 * 00322 MAXWRK = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) 00323 MAXWRK = MAX( MAXWRK, M*M+4*M+2*M* 00324 $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) 00325 MAXWRK = MAX( MAXWRK, M*M+4*M+NRHS* 00326 $ ILAENV( 1, 'DORMBR', 'QLT', M, NRHS, M, -1 ) ) 00327 MAXWRK = MAX( MAXWRK, M*M+4*M+( M-1 )* 00328 $ ILAENV( 1, 'DORMBR', 'PLN', M, NRHS, M, -1 ) ) 00329 IF( NRHS.GT.1 ) THEN 00330 MAXWRK = MAX( MAXWRK, M*M+M+M*NRHS ) 00331 ELSE 00332 MAXWRK = MAX( MAXWRK, M*M+2*M ) 00333 END IF 00334 MAXWRK = MAX( MAXWRK, M+NRHS* 00335 $ ILAENV( 1, 'DORMLQ', 'LT', N, NRHS, M, -1 ) ) 00336 MAXWRK = MAX( MAXWRK, M*M+4*M+WLALSD ) 00337 ! XXX: Ensure the Path 2a case below is triggered. The workspace 00338 ! calculation should use queries for all routines eventually. 00339 MAXWRK = MAX( MAXWRK, 00340 $ 4*M+M*M+MAX( M, 2*M-4, NRHS, N-3*M ) ) 00341 ELSE 00342 * 00343 * Path 2 - remaining underdetermined cases. 00344 * 00345 MAXWRK = 3*M + ( N+M )*ILAENV( 1, 'DGEBRD', ' ', M, N, 00346 $ -1, -1 ) 00347 MAXWRK = MAX( MAXWRK, 3*M+NRHS* 00348 $ ILAENV( 1, 'DORMBR', 'QLT', M, NRHS, N, -1 ) ) 00349 MAXWRK = MAX( MAXWRK, 3*M+M* 00350 $ ILAENV( 1, 'DORMBR', 'PLN', N, NRHS, M, -1 ) ) 00351 MAXWRK = MAX( MAXWRK, 3*M+WLALSD ) 00352 END IF 00353 MINWRK = MAX( 3*M+NRHS, 3*M+M, 3*M+WLALSD ) 00354 END IF 00355 MINWRK = MIN( MINWRK, MAXWRK ) 00356 WORK( 1 ) = MAXWRK 00357 IWORK( 1 ) = LIWORK 00358 00359 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN 00360 INFO = -12 00361 END IF 00362 END IF 00363 * 00364 IF( INFO.NE.0 ) THEN 00365 CALL XERBLA( 'DGELSD', -INFO ) 00366 RETURN 00367 ELSE IF( LQUERY ) THEN 00368 GO TO 10 00369 END IF 00370 * 00371 * Quick return if possible. 00372 * 00373 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 00374 RANK = 0 00375 RETURN 00376 END IF 00377 * 00378 * Get machine parameters. 00379 * 00380 EPS = DLAMCH( 'P' ) 00381 SFMIN = DLAMCH( 'S' ) 00382 SMLNUM = SFMIN / EPS 00383 BIGNUM = ONE / SMLNUM 00384 CALL DLABAD( SMLNUM, BIGNUM ) 00385 * 00386 * Scale A if max entry outside range [SMLNUM,BIGNUM]. 00387 * 00388 ANRM = DLANGE( 'M', M, N, A, LDA, WORK ) 00389 IASCL = 0 00390 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00391 * 00392 * Scale matrix norm up to SMLNUM. 00393 * 00394 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO ) 00395 IASCL = 1 00396 ELSE IF( ANRM.GT.BIGNUM ) THEN 00397 * 00398 * Scale matrix norm down to BIGNUM. 00399 * 00400 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO ) 00401 IASCL = 2 00402 ELSE IF( ANRM.EQ.ZERO ) THEN 00403 * 00404 * Matrix all zero. Return zero solution. 00405 * 00406 CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB ) 00407 CALL DLASET( 'F', MINMN, 1, ZERO, ZERO, S, 1 ) 00408 RANK = 0 00409 GO TO 10 00410 END IF 00411 * 00412 * Scale B if max entry outside range [SMLNUM,BIGNUM]. 00413 * 00414 BNRM = DLANGE( 'M', M, NRHS, B, LDB, WORK ) 00415 IBSCL = 0 00416 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 00417 * 00418 * Scale matrix norm up to SMLNUM. 00419 * 00420 CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO ) 00421 IBSCL = 1 00422 ELSE IF( BNRM.GT.BIGNUM ) THEN 00423 * 00424 * Scale matrix norm down to BIGNUM. 00425 * 00426 CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO ) 00427 IBSCL = 2 00428 END IF 00429 * 00430 * If M < N make sure certain entries of B are zero. 00431 * 00432 IF( M.LT.N ) 00433 $ CALL DLASET( 'F', N-M, NRHS, ZERO, ZERO, B( M+1, 1 ), LDB ) 00434 * 00435 * Overdetermined case. 00436 * 00437 IF( M.GE.N ) THEN 00438 * 00439 * Path 1 - overdetermined or exactly determined. 00440 * 00441 MM = M 00442 IF( M.GE.MNTHR ) THEN 00443 * 00444 * Path 1a - overdetermined, with many more rows than columns. 00445 * 00446 MM = N 00447 ITAU = 1 00448 NWORK = ITAU + N 00449 * 00450 * Compute A=Q*R. 00451 * (Workspace: need 2*N, prefer N+N*NB) 00452 * 00453 CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 00454 $ LWORK-NWORK+1, INFO ) 00455 * 00456 * Multiply B by transpose(Q). 00457 * (Workspace: need N+NRHS, prefer N+NRHS*NB) 00458 * 00459 CALL DORMQR( 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAU ), B, 00460 $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00461 * 00462 * Zero out below R. 00463 * 00464 IF( N.GT.1 ) THEN 00465 CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA ) 00466 END IF 00467 END IF 00468 * 00469 IE = 1 00470 ITAUQ = IE + N 00471 ITAUP = ITAUQ + N 00472 NWORK = ITAUP + N 00473 * 00474 * Bidiagonalize R in A. 00475 * (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB) 00476 * 00477 CALL DGEBRD( MM, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 00478 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 00479 $ INFO ) 00480 * 00481 * Multiply B by transpose of left bidiagonalizing vectors of R. 00482 * (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB) 00483 * 00484 CALL DORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, WORK( ITAUQ ), 00485 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00486 * 00487 * Solve the bidiagonal least squares problem. 00488 * 00489 CALL DLALSD( 'U', SMLSIZ, N, NRHS, S, WORK( IE ), B, LDB, 00490 $ RCOND, RANK, WORK( NWORK ), IWORK, INFO ) 00491 IF( INFO.NE.0 ) THEN 00492 GO TO 10 00493 END IF 00494 * 00495 * Multiply B by right bidiagonalizing vectors of R. 00496 * 00497 CALL DORMBR( 'P', 'L', 'N', N, NRHS, N, A, LDA, WORK( ITAUP ), 00498 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00499 * 00500 ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+ 00501 $ MAX( M, 2*M-4, NRHS, N-3*M, WLALSD ) ) THEN 00502 * 00503 * Path 2a - underdetermined, with many more columns than rows 00504 * and sufficient workspace for an efficient algorithm. 00505 * 00506 LDWORK = M 00507 IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ), 00508 $ M*LDA+M+M*NRHS, 4*M+M*LDA+WLALSD ) )LDWORK = LDA 00509 ITAU = 1 00510 NWORK = M + 1 00511 * 00512 * Compute A=L*Q. 00513 * (Workspace: need 2*M, prefer M+M*NB) 00514 * 00515 CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 00516 $ LWORK-NWORK+1, INFO ) 00517 IL = NWORK 00518 * 00519 * Copy L to WORK(IL), zeroing out above its diagonal. 00520 * 00521 CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK ) 00522 CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, WORK( IL+LDWORK ), 00523 $ LDWORK ) 00524 IE = IL + LDWORK*M 00525 ITAUQ = IE + M 00526 ITAUP = ITAUQ + M 00527 NWORK = ITAUP + M 00528 * 00529 * Bidiagonalize L in WORK(IL). 00530 * (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB) 00531 * 00532 CALL DGEBRD( M, M, WORK( IL ), LDWORK, S, WORK( IE ), 00533 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), 00534 $ LWORK-NWORK+1, INFO ) 00535 * 00536 * Multiply B by transpose of left bidiagonalizing vectors of L. 00537 * (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB) 00538 * 00539 CALL DORMBR( 'Q', 'L', 'T', M, NRHS, M, WORK( IL ), LDWORK, 00540 $ WORK( ITAUQ ), B, LDB, WORK( NWORK ), 00541 $ LWORK-NWORK+1, INFO ) 00542 * 00543 * Solve the bidiagonal least squares problem. 00544 * 00545 CALL DLALSD( 'U', SMLSIZ, M, NRHS, S, WORK( IE ), B, LDB, 00546 $ RCOND, RANK, WORK( NWORK ), IWORK, INFO ) 00547 IF( INFO.NE.0 ) THEN 00548 GO TO 10 00549 END IF 00550 * 00551 * Multiply B by right bidiagonalizing vectors of L. 00552 * 00553 CALL DORMBR( 'P', 'L', 'N', M, NRHS, M, WORK( IL ), LDWORK, 00554 $ WORK( ITAUP ), B, LDB, WORK( NWORK ), 00555 $ LWORK-NWORK+1, INFO ) 00556 * 00557 * Zero out below first M rows of B. 00558 * 00559 CALL DLASET( 'F', N-M, NRHS, ZERO, ZERO, B( M+1, 1 ), LDB ) 00560 NWORK = ITAU + M 00561 * 00562 * Multiply transpose(Q) by B. 00563 * (Workspace: need M+NRHS, prefer M+NRHS*NB) 00564 * 00565 CALL DORMLQ( 'L', 'T', N, NRHS, M, A, LDA, WORK( ITAU ), B, 00566 $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00567 * 00568 ELSE 00569 * 00570 * Path 2 - remaining underdetermined cases. 00571 * 00572 IE = 1 00573 ITAUQ = IE + M 00574 ITAUP = ITAUQ + M 00575 NWORK = ITAUP + M 00576 * 00577 * Bidiagonalize A. 00578 * (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB) 00579 * 00580 CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 00581 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 00582 $ INFO ) 00583 * 00584 * Multiply B by transpose of left bidiagonalizing vectors. 00585 * (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB) 00586 * 00587 CALL DORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAUQ ), 00588 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00589 * 00590 * Solve the bidiagonal least squares problem. 00591 * 00592 CALL DLALSD( 'L', SMLSIZ, M, NRHS, S, WORK( IE ), B, LDB, 00593 $ RCOND, RANK, WORK( NWORK ), IWORK, INFO ) 00594 IF( INFO.NE.0 ) THEN 00595 GO TO 10 00596 END IF 00597 * 00598 * Multiply B by right bidiagonalizing vectors of A. 00599 * 00600 CALL DORMBR( 'P', 'L', 'N', N, NRHS, M, A, LDA, WORK( ITAUP ), 00601 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00602 * 00603 END IF 00604 * 00605 * Undo scaling. 00606 * 00607 IF( IASCL.EQ.1 ) THEN 00608 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO ) 00609 CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN, 00610 $ INFO ) 00611 ELSE IF( IASCL.EQ.2 ) THEN 00612 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO ) 00613 CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN, 00614 $ INFO ) 00615 END IF 00616 IF( IBSCL.EQ.1 ) THEN 00617 CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO ) 00618 ELSE IF( IBSCL.EQ.2 ) THEN 00619 CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO ) 00620 END IF 00621 * 00622 10 CONTINUE 00623 WORK( 1 ) = MAXWRK 00624 IWORK( 1 ) = LIWORK 00625 RETURN 00626 * 00627 * End of DGELSD 00628 * 00629 END