 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, 00002 $ WORK, LWORK, IWORK, INFO ) 00003 * 00004 * -- LAPACK driver routine (version 3.2.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 * June 2010 00008 * 00009 * .. Scalar Arguments .. 00010 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK 00011 DOUBLE PRECISION RCOND 00012 * .. 00013 * .. Array Arguments .. 00014 INTEGER IWORK( * ) 00015 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * DGELSD 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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 WORK array, returns 00113 * this value as the first entry of the WORK array, and no error 00114 * message related to LWORK is issued by XERBLA. 00115 * 00116 * IWORK (workspace) INTEGER array, dimension (MAX(1,LIWORK)) 00117 * LIWORK >= max(1, 3 * MINMN * NLVL + 11 * MINMN), 00118 * where MINMN = MIN( M,N ). 00119 * On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK. 00120 * 00121 * INFO (output) INTEGER 00122 * = 0: successful exit 00123 * < 0: if INFO = -i, the i-th argument had an illegal value. 00124 * > 0: the algorithm for computing the SVD failed to converge; 00125 * if INFO = i, i off-diagonal elements of an intermediate 00126 * bidiagonal form did not converge to zero. 00127 * 00128 * Further Details 00129 * =============== 00130 * 00131 * Based on contributions by 00132 * Ming Gu and Ren-Cang Li, Computer Science Division, University of 00133 * California at Berkeley, USA 00134 * Osni Marques, LBNL/NERSC, USA 00135 * 00136 * ===================================================================== 00137 * 00138 * .. Parameters .. 00139 DOUBLE PRECISION ZERO, ONE, TWO 00140 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 ) 00141 * .. 00142 * .. Local Scalars .. 00143 LOGICAL LQUERY 00144 INTEGER IASCL, IBSCL, IE, IL, ITAU, ITAUP, ITAUQ, 00145 $ LDWORK, LIWORK, MAXMN, MAXWRK, MINMN, MINWRK, 00146 $ MM, MNTHR, NLVL, NWORK, SMLSIZ, WLALSD 00147 DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM 00148 * .. 00149 * .. External Subroutines .. 00150 EXTERNAL DGEBRD, DGELQF, DGEQRF, DLABAD, DLACPY, DLALSD, 00151 $ DLASCL, DLASET, DORMBR, DORMLQ, DORMQR, XERBLA 00152 * .. 00153 * .. External Functions .. 00154 INTEGER ILAENV 00155 DOUBLE PRECISION DLAMCH, DLANGE 00156 EXTERNAL ILAENV, DLAMCH, DLANGE 00157 * .. 00158 * .. Intrinsic Functions .. 00159 INTRINSIC DBLE, INT, LOG, MAX, MIN 00160 * .. 00161 * .. Executable Statements .. 00162 * 00163 * Test the input arguments. 00164 * 00165 INFO = 0 00166 MINMN = MIN( M, N ) 00167 MAXMN = MAX( M, N ) 00168 MNTHR = ILAENV( 6, 'DGELSD', ' ', M, N, NRHS, -1 ) 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 SMLSIZ = ILAENV( 9, 'DGELSD', ' ', 0, 0, 0, 0 ) 00183 * 00184 * Compute workspace. 00185 * (Note: Comments in the code beginning "Workspace:" describe the 00186 * minimal amount of workspace needed at that point in the code, 00187 * as well as the preferred amount for good performance. 00188 * NB refers to the optimal block size for the immediately 00189 * following subroutine, as returned by ILAENV.) 00190 * 00191 MINWRK = 1 00192 LIWORK = 1 00193 MINMN = MAX( 1, MINMN ) 00194 NLVL = MAX( INT( LOG( DBLE( MINMN ) / DBLE( SMLSIZ+1 ) ) / 00195 $ LOG( TWO ) ) + 1, 0 ) 00196 * 00197 IF( INFO.EQ.0 ) THEN 00198 MAXWRK = 0 00199 LIWORK = 3*MINMN*NLVL + 11*MINMN 00200 MM = M 00201 IF( M.GE.N .AND. M.GE.MNTHR ) THEN 00202 * 00203 * Path 1a - overdetermined, with many more rows than columns. 00204 * 00205 MM = N 00206 MAXWRK = MAX( MAXWRK, N+N*ILAENV( 1, 'DGEQRF', ' ', M, N, 00207 $ -1, -1 ) ) 00208 MAXWRK = MAX( MAXWRK, N+NRHS* 00209 $ ILAENV( 1, 'DORMQR', 'LT', 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 )* 00216 $ ILAENV( 1, 'DGEBRD', ' ', MM, N, -1, -1 ) ) 00217 MAXWRK = MAX( MAXWRK, 3*N+NRHS* 00218 $ ILAENV( 1, 'DORMBR', 'QLT', MM, NRHS, N, -1 ) ) 00219 MAXWRK = MAX( MAXWRK, 3*N+( N-1 )* 00220 $ ILAENV( 1, 'DORMBR', 'PLN', N, NRHS, N, -1 ) ) 00221 WLALSD = 9*N+2*N*SMLSIZ+8*N*NLVL+N*NRHS+(SMLSIZ+1)**2 00222 MAXWRK = MAX( MAXWRK, 3*N+WLALSD ) 00223 MINWRK = MAX( 3*N+MM, 3*N+NRHS, 3*N+WLALSD ) 00224 END IF 00225 IF( N.GT.M ) THEN 00226 WLALSD = 9*M+2*M*SMLSIZ+8*M*NLVL+M*NRHS+(SMLSIZ+1)**2 00227 IF( N.GE.MNTHR ) THEN 00228 * 00229 * Path 2a - underdetermined, with many more columns 00230 * than rows. 00231 * 00232 MAXWRK = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) 00233 MAXWRK = MAX( MAXWRK, M*M+4*M+2*M* 00234 $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) 00235 MAXWRK = MAX( MAXWRK, M*M+4*M+NRHS* 00236 $ ILAENV( 1, 'DORMBR', 'QLT', M, NRHS, M, -1 ) ) 00237 MAXWRK = MAX( MAXWRK, M*M+4*M+( M-1 )* 00238 $ ILAENV( 1, 'DORMBR', 'PLN', M, NRHS, M, -1 ) ) 00239 IF( NRHS.GT.1 ) THEN 00240 MAXWRK = MAX( MAXWRK, M*M+M+M*NRHS ) 00241 ELSE 00242 MAXWRK = MAX( MAXWRK, M*M+2*M ) 00243 END IF 00244 MAXWRK = MAX( MAXWRK, M+NRHS* 00245 $ ILAENV( 1, 'DORMLQ', 'LT', N, NRHS, M, -1 ) ) 00246 MAXWRK = MAX( MAXWRK, M*M+4*M+WLALSD ) 00247 ! XXX: Ensure the Path 2a case below is triggered. The workspace 00248 ! calculation should use queries for all routines eventually. 00249 MAXWRK = MAX( MAXWRK, 00250 $ 4*M+M*M+MAX( M, 2*M-4, NRHS, N-3*M ) ) 00251 ELSE 00252 * 00253 * Path 2 - remaining underdetermined cases. 00254 * 00255 MAXWRK = 3*M + ( N+M )*ILAENV( 1, 'DGEBRD', ' ', M, N, 00256 $ -1, -1 ) 00257 MAXWRK = MAX( MAXWRK, 3*M+NRHS* 00258 $ ILAENV( 1, 'DORMBR', 'QLT', M, NRHS, N, -1 ) ) 00259 MAXWRK = MAX( MAXWRK, 3*M+M* 00260 $ ILAENV( 1, 'DORMBR', 'PLN', N, NRHS, M, -1 ) ) 00261 MAXWRK = MAX( MAXWRK, 3*M+WLALSD ) 00262 END IF 00263 MINWRK = MAX( 3*M+NRHS, 3*M+M, 3*M+WLALSD ) 00264 END IF 00265 MINWRK = MIN( MINWRK, MAXWRK ) 00266 WORK( 1 ) = MAXWRK 00267 IWORK( 1 ) = LIWORK 00268 00269 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN 00270 INFO = -12 00271 END IF 00272 END IF 00273 * 00274 IF( INFO.NE.0 ) THEN 00275 CALL XERBLA( 'DGELSD', -INFO ) 00276 RETURN 00277 ELSE IF( LQUERY ) THEN 00278 GO TO 10 00279 END IF 00280 * 00281 * Quick return if possible. 00282 * 00283 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 00284 RANK = 0 00285 RETURN 00286 END IF 00287 * 00288 * Get machine parameters. 00289 * 00290 EPS = DLAMCH( 'P' ) 00291 SFMIN = DLAMCH( 'S' ) 00292 SMLNUM = SFMIN / EPS 00293 BIGNUM = ONE / SMLNUM 00294 CALL DLABAD( SMLNUM, BIGNUM ) 00295 * 00296 * Scale A if max entry outside range [SMLNUM,BIGNUM]. 00297 * 00298 ANRM = DLANGE( 'M', M, N, A, LDA, WORK ) 00299 IASCL = 0 00300 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00301 * 00302 * Scale matrix norm up to SMLNUM. 00303 * 00304 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO ) 00305 IASCL = 1 00306 ELSE IF( ANRM.GT.BIGNUM ) THEN 00307 * 00308 * Scale matrix norm down to BIGNUM. 00309 * 00310 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO ) 00311 IASCL = 2 00312 ELSE IF( ANRM.EQ.ZERO ) THEN 00313 * 00314 * Matrix all zero. Return zero solution. 00315 * 00316 CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB ) 00317 CALL DLASET( 'F', MINMN, 1, ZERO, ZERO, S, 1 ) 00318 RANK = 0 00319 GO TO 10 00320 END IF 00321 * 00322 * Scale B if max entry outside range [SMLNUM,BIGNUM]. 00323 * 00324 BNRM = DLANGE( 'M', M, NRHS, B, LDB, WORK ) 00325 IBSCL = 0 00326 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 00327 * 00328 * Scale matrix norm up to SMLNUM. 00329 * 00330 CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO ) 00331 IBSCL = 1 00332 ELSE IF( BNRM.GT.BIGNUM ) THEN 00333 * 00334 * Scale matrix norm down to BIGNUM. 00335 * 00336 CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO ) 00337 IBSCL = 2 00338 END IF 00339 * 00340 * If M < N make sure certain entries of B are zero. 00341 * 00342 IF( M.LT.N ) 00343 $ CALL DLASET( 'F', N-M, NRHS, ZERO, ZERO, B( M+1, 1 ), LDB ) 00344 * 00345 * Overdetermined case. 00346 * 00347 IF( M.GE.N ) THEN 00348 * 00349 * Path 1 - overdetermined or exactly determined. 00350 * 00351 MM = M 00352 IF( M.GE.MNTHR ) THEN 00353 * 00354 * Path 1a - overdetermined, with many more rows than columns. 00355 * 00356 MM = N 00357 ITAU = 1 00358 NWORK = ITAU + N 00359 * 00360 * Compute A=Q*R. 00361 * (Workspace: need 2*N, prefer N+N*NB) 00362 * 00363 CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 00364 $ LWORK-NWORK+1, INFO ) 00365 * 00366 * Multiply B by transpose(Q). 00367 * (Workspace: need N+NRHS, prefer N+NRHS*NB) 00368 * 00369 CALL DORMQR( 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAU ), B, 00370 $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00371 * 00372 * Zero out below R. 00373 * 00374 IF( N.GT.1 ) THEN 00375 CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA ) 00376 END IF 00377 END IF 00378 * 00379 IE = 1 00380 ITAUQ = IE + N 00381 ITAUP = ITAUQ + N 00382 NWORK = ITAUP + N 00383 * 00384 * Bidiagonalize R in A. 00385 * (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB) 00386 * 00387 CALL DGEBRD( MM, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 00388 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 00389 $ INFO ) 00390 * 00391 * Multiply B by transpose of left bidiagonalizing vectors of R. 00392 * (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB) 00393 * 00394 CALL DORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, WORK( ITAUQ ), 00395 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00396 * 00397 * Solve the bidiagonal least squares problem. 00398 * 00399 CALL DLALSD( 'U', SMLSIZ, N, NRHS, S, WORK( IE ), B, LDB, 00400 $ RCOND, RANK, WORK( NWORK ), IWORK, INFO ) 00401 IF( INFO.NE.0 ) THEN 00402 GO TO 10 00403 END IF 00404 * 00405 * Multiply B by right bidiagonalizing vectors of R. 00406 * 00407 CALL DORMBR( 'P', 'L', 'N', N, NRHS, N, A, LDA, WORK( ITAUP ), 00408 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00409 * 00410 ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+ 00411 $ MAX( M, 2*M-4, NRHS, N-3*M, WLALSD ) ) THEN 00412 * 00413 * Path 2a - underdetermined, with many more columns than rows 00414 * and sufficient workspace for an efficient algorithm. 00415 * 00416 LDWORK = M 00417 IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ), 00418 $ M*LDA+M+M*NRHS, 4*M+M*LDA+WLALSD ) )LDWORK = LDA 00419 ITAU = 1 00420 NWORK = M + 1 00421 * 00422 * Compute A=L*Q. 00423 * (Workspace: need 2*M, prefer M+M*NB) 00424 * 00425 CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 00426 $ LWORK-NWORK+1, INFO ) 00427 IL = NWORK 00428 * 00429 * Copy L to WORK(IL), zeroing out above its diagonal. 00430 * 00431 CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK ) 00432 CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, WORK( IL+LDWORK ), 00433 $ LDWORK ) 00434 IE = IL + LDWORK*M 00435 ITAUQ = IE + M 00436 ITAUP = ITAUQ + M 00437 NWORK = ITAUP + M 00438 * 00439 * Bidiagonalize L in WORK(IL). 00440 * (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB) 00441 * 00442 CALL DGEBRD( M, M, WORK( IL ), LDWORK, S, WORK( IE ), 00443 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), 00444 $ LWORK-NWORK+1, INFO ) 00445 * 00446 * Multiply B by transpose of left bidiagonalizing vectors of L. 00447 * (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB) 00448 * 00449 CALL DORMBR( 'Q', 'L', 'T', M, NRHS, M, WORK( IL ), LDWORK, 00450 $ WORK( ITAUQ ), B, LDB, WORK( NWORK ), 00451 $ LWORK-NWORK+1, INFO ) 00452 * 00453 * Solve the bidiagonal least squares problem. 00454 * 00455 CALL DLALSD( 'U', SMLSIZ, M, NRHS, S, WORK( IE ), B, LDB, 00456 $ RCOND, RANK, WORK( NWORK ), IWORK, INFO ) 00457 IF( INFO.NE.0 ) THEN 00458 GO TO 10 00459 END IF 00460 * 00461 * Multiply B by right bidiagonalizing vectors of L. 00462 * 00463 CALL DORMBR( 'P', 'L', 'N', M, NRHS, M, WORK( IL ), LDWORK, 00464 $ WORK( ITAUP ), B, LDB, WORK( NWORK ), 00465 $ LWORK-NWORK+1, INFO ) 00466 * 00467 * Zero out below first M rows of B. 00468 * 00469 CALL DLASET( 'F', N-M, NRHS, ZERO, ZERO, B( M+1, 1 ), LDB ) 00470 NWORK = ITAU + M 00471 * 00472 * Multiply transpose(Q) by B. 00473 * (Workspace: need M+NRHS, prefer M+NRHS*NB) 00474 * 00475 CALL DORMLQ( 'L', 'T', N, NRHS, M, A, LDA, WORK( ITAU ), B, 00476 $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00477 * 00478 ELSE 00479 * 00480 * Path 2 - remaining underdetermined cases. 00481 * 00482 IE = 1 00483 ITAUQ = IE + M 00484 ITAUP = ITAUQ + M 00485 NWORK = ITAUP + M 00486 * 00487 * Bidiagonalize A. 00488 * (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB) 00489 * 00490 CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 00491 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 00492 $ INFO ) 00493 * 00494 * Multiply B by transpose of left bidiagonalizing vectors. 00495 * (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB) 00496 * 00497 CALL DORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAUQ ), 00498 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00499 * 00500 * Solve the bidiagonal least squares problem. 00501 * 00502 CALL DLALSD( 'L', SMLSIZ, M, NRHS, S, WORK( IE ), B, LDB, 00503 $ RCOND, RANK, WORK( NWORK ), IWORK, INFO ) 00504 IF( INFO.NE.0 ) THEN 00505 GO TO 10 00506 END IF 00507 * 00508 * Multiply B by right bidiagonalizing vectors of A. 00509 * 00510 CALL DORMBR( 'P', 'L', 'N', N, NRHS, M, A, LDA, WORK( ITAUP ), 00511 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00512 * 00513 END IF 00514 * 00515 * Undo scaling. 00516 * 00517 IF( IASCL.EQ.1 ) THEN 00518 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO ) 00519 CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN, 00520 $ INFO ) 00521 ELSE IF( IASCL.EQ.2 ) THEN 00522 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO ) 00523 CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN, 00524 $ INFO ) 00525 END IF 00526 IF( IBSCL.EQ.1 ) THEN 00527 CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO ) 00528 ELSE IF( IBSCL.EQ.2 ) THEN 00529 CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO ) 00530 END IF 00531 * 00532 10 CONTINUE 00533 WORK( 1 ) = MAXWRK 00534 IWORK( 1 ) = LIWORK 00535 RETURN 00536 * 00537 * End of DGELSD 00538 * 00539 END
1.7.3 
