|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, 00002 $ WORK, LWORK, 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 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * DGELSS computes the minimum norm solution to a real linear least 00021 * squares problem: 00022 * 00023 * Minimize 2-norm(| b - A*x |). 00024 * 00025 * using the singular value decomposition (SVD) of A. A is an M-by-N 00026 * matrix which may be rank-deficient. 00027 * 00028 * Several right hand side vectors b and solution vectors x can be 00029 * handled in a single call; they are stored as the columns of the 00030 * M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix 00031 * X. 00032 * 00033 * The effective rank of A is determined by treating as zero those 00034 * singular values which are less than RCOND times the largest singular 00035 * value. 00036 * 00037 * Arguments 00038 * ========= 00039 * 00040 * M (input) INTEGER 00041 * The number of rows of the matrix A. M >= 0. 00042 * 00043 * N (input) INTEGER 00044 * The number of columns of the matrix A. N >= 0. 00045 * 00046 * NRHS (input) INTEGER 00047 * The number of right hand sides, i.e., the number of columns 00048 * of the matrices B and X. NRHS >= 0. 00049 * 00050 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 00051 * On entry, the M-by-N matrix A. 00052 * On exit, the first min(m,n) rows of A are overwritten with 00053 * its right singular vectors, stored rowwise. 00054 * 00055 * LDA (input) INTEGER 00056 * The leading dimension of the array A. LDA >= max(1,M). 00057 * 00058 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) 00059 * On entry, the M-by-NRHS right hand side matrix B. 00060 * On exit, B is overwritten by the N-by-NRHS solution 00061 * matrix X. If m >= n and RANK = n, the residual 00062 * sum-of-squares for the solution in the i-th column is given 00063 * by the sum of squares of elements n+1:m in that column. 00064 * 00065 * LDB (input) INTEGER 00066 * The leading dimension of the array B. LDB >= max(1,max(M,N)). 00067 * 00068 * S (output) DOUBLE PRECISION array, dimension (min(M,N)) 00069 * The singular values of A in decreasing order. 00070 * The condition number of A in the 2-norm = S(1)/S(min(m,n)). 00071 * 00072 * RCOND (input) DOUBLE PRECISION 00073 * RCOND is used to determine the effective rank of A. 00074 * Singular values S(i) <= RCOND*S(1) are treated as zero. 00075 * If RCOND < 0, machine precision is used instead. 00076 * 00077 * RANK (output) INTEGER 00078 * The effective rank of A, i.e., the number of singular values 00079 * which are greater than RCOND*S(1). 00080 * 00081 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 00082 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00083 * 00084 * LWORK (input) INTEGER 00085 * The dimension of the array WORK. LWORK >= 1, and also: 00086 * LWORK >= 3*min(M,N) + max( 2*min(M,N), max(M,N), NRHS ) 00087 * For good performance, LWORK should generally be larger. 00088 * 00089 * If LWORK = -1, then a workspace query is assumed; the routine 00090 * only calculates the optimal size of the WORK array, returns 00091 * this value as the first entry of the WORK array, and no error 00092 * message related to LWORK is issued by XERBLA. 00093 * 00094 * INFO (output) INTEGER 00095 * = 0: successful exit 00096 * < 0: if INFO = -i, the i-th argument had an illegal value. 00097 * > 0: the algorithm for computing the SVD failed to converge; 00098 * if INFO = i, i off-diagonal elements of an intermediate 00099 * bidiagonal form did not converge to zero. 00100 * 00101 * ===================================================================== 00102 * 00103 * .. Parameters .. 00104 DOUBLE PRECISION ZERO, ONE 00105 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00106 * .. 00107 * .. Local Scalars .. 00108 LOGICAL LQUERY 00109 INTEGER BDSPAC, BL, CHUNK, I, IASCL, IBSCL, IE, IL, 00110 $ ITAU, ITAUP, ITAUQ, IWORK, LDWORK, MAXMN, 00111 $ MAXWRK, MINMN, MINWRK, MM, MNTHR 00112 DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR 00113 * .. 00114 * .. Local Arrays .. 00115 DOUBLE PRECISION VDUM( 1 ) 00116 * .. 00117 * .. External Subroutines .. 00118 EXTERNAL DBDSQR, DCOPY, DGEBRD, DGELQF, DGEMM, DGEMV, 00119 $ DGEQRF, DLABAD, DLACPY, DLASCL, DLASET, DORGBR, 00120 $ DORMBR, DORMLQ, DORMQR, DRSCL, XERBLA 00121 * .. 00122 * .. External Functions .. 00123 INTEGER ILAENV 00124 DOUBLE PRECISION DLAMCH, DLANGE 00125 EXTERNAL ILAENV, DLAMCH, DLANGE 00126 * .. 00127 * .. Intrinsic Functions .. 00128 INTRINSIC MAX, MIN 00129 * .. 00130 * .. Executable Statements .. 00131 * 00132 * Test the input arguments 00133 * 00134 INFO = 0 00135 MINMN = MIN( M, N ) 00136 MAXMN = MAX( M, N ) 00137 LQUERY = ( LWORK.EQ.-1 ) 00138 IF( M.LT.0 ) THEN 00139 INFO = -1 00140 ELSE IF( N.LT.0 ) THEN 00141 INFO = -2 00142 ELSE IF( NRHS.LT.0 ) THEN 00143 INFO = -3 00144 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00145 INFO = -5 00146 ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN 00147 INFO = -7 00148 END IF 00149 * 00150 * Compute workspace 00151 * (Note: Comments in the code beginning "Workspace:" describe the 00152 * minimal amount of workspace needed at that point in the code, 00153 * as well as the preferred amount for good performance. 00154 * NB refers to the optimal block size for the immediately 00155 * following subroutine, as returned by ILAENV.) 00156 * 00157 IF( INFO.EQ.0 ) THEN 00158 MINWRK = 1 00159 MAXWRK = 1 00160 IF( MINMN.GT.0 ) THEN 00161 MM = M 00162 MNTHR = ILAENV( 6, 'DGELSS', ' ', M, N, NRHS, -1 ) 00163 IF( M.GE.N .AND. M.GE.MNTHR ) THEN 00164 * 00165 * Path 1a - overdetermined, with many more rows than 00166 * columns 00167 * 00168 MM = N 00169 MAXWRK = MAX( MAXWRK, N + N*ILAENV( 1, 'DGEQRF', ' ', M, 00170 $ N, -1, -1 ) ) 00171 MAXWRK = MAX( MAXWRK, N + NRHS*ILAENV( 1, 'DORMQR', 'LT', 00172 $ M, NRHS, N, -1 ) ) 00173 END IF 00174 IF( M.GE.N ) THEN 00175 * 00176 * Path 1 - overdetermined or exactly determined 00177 * 00178 * Compute workspace needed for DBDSQR 00179 * 00180 BDSPAC = MAX( 1, 5*N ) 00181 MAXWRK = MAX( MAXWRK, 3*N + ( MM + N )*ILAENV( 1, 00182 $ 'DGEBRD', ' ', MM, N, -1, -1 ) ) 00183 MAXWRK = MAX( MAXWRK, 3*N + NRHS*ILAENV( 1, 'DORMBR', 00184 $ 'QLT', MM, NRHS, N, -1 ) ) 00185 MAXWRK = MAX( MAXWRK, 3*N + ( N - 1 )*ILAENV( 1, 00186 $ 'DORGBR', 'P', N, N, N, -1 ) ) 00187 MAXWRK = MAX( MAXWRK, BDSPAC ) 00188 MAXWRK = MAX( MAXWRK, N*NRHS ) 00189 MINWRK = MAX( 3*N + MM, 3*N + NRHS, BDSPAC ) 00190 MAXWRK = MAX( MINWRK, MAXWRK ) 00191 END IF 00192 IF( N.GT.M ) THEN 00193 * 00194 * Compute workspace needed for DBDSQR 00195 * 00196 BDSPAC = MAX( 1, 5*M ) 00197 MINWRK = MAX( 3*M+NRHS, 3*M+N, BDSPAC ) 00198 IF( N.GE.MNTHR ) THEN 00199 * 00200 * Path 2a - underdetermined, with many more columns 00201 * than rows 00202 * 00203 MAXWRK = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, 00204 $ -1 ) 00205 MAXWRK = MAX( MAXWRK, M*M + 4*M + 2*M*ILAENV( 1, 00206 $ 'DGEBRD', ' ', M, M, -1, -1 ) ) 00207 MAXWRK = MAX( MAXWRK, M*M + 4*M + NRHS*ILAENV( 1, 00208 $ 'DORMBR', 'QLT', M, NRHS, M, -1 ) ) 00209 MAXWRK = MAX( MAXWRK, M*M + 4*M + 00210 $ ( M - 1 )*ILAENV( 1, 'DORGBR', 'P', M, 00211 $ M, M, -1 ) ) 00212 MAXWRK = MAX( MAXWRK, M*M + M + BDSPAC ) 00213 IF( NRHS.GT.1 ) THEN 00214 MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS ) 00215 ELSE 00216 MAXWRK = MAX( MAXWRK, M*M + 2*M ) 00217 END IF 00218 MAXWRK = MAX( MAXWRK, M + NRHS*ILAENV( 1, 'DORMLQ', 00219 $ 'LT', N, NRHS, M, -1 ) ) 00220 ELSE 00221 * 00222 * Path 2 - underdetermined 00223 * 00224 MAXWRK = 3*M + ( N + M )*ILAENV( 1, 'DGEBRD', ' ', M, 00225 $ N, -1, -1 ) 00226 MAXWRK = MAX( MAXWRK, 3*M + NRHS*ILAENV( 1, 'DORMBR', 00227 $ 'QLT', M, NRHS, M, -1 ) ) 00228 MAXWRK = MAX( MAXWRK, 3*M + M*ILAENV( 1, 'DORGBR', 00229 $ 'P', M, N, M, -1 ) ) 00230 MAXWRK = MAX( MAXWRK, BDSPAC ) 00231 MAXWRK = MAX( MAXWRK, N*NRHS ) 00232 END IF 00233 END IF 00234 MAXWRK = MAX( MINWRK, MAXWRK ) 00235 END IF 00236 WORK( 1 ) = MAXWRK 00237 * 00238 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) 00239 $ INFO = -12 00240 END IF 00241 * 00242 IF( INFO.NE.0 ) THEN 00243 CALL XERBLA( 'DGELSS', -INFO ) 00244 RETURN 00245 ELSE IF( LQUERY ) THEN 00246 RETURN 00247 END IF 00248 * 00249 * Quick return if possible 00250 * 00251 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 00252 RANK = 0 00253 RETURN 00254 END IF 00255 * 00256 * Get machine parameters 00257 * 00258 EPS = DLAMCH( 'P' ) 00259 SFMIN = DLAMCH( 'S' ) 00260 SMLNUM = SFMIN / EPS 00261 BIGNUM = ONE / SMLNUM 00262 CALL DLABAD( SMLNUM, BIGNUM ) 00263 * 00264 * Scale A if max element outside range [SMLNUM,BIGNUM] 00265 * 00266 ANRM = DLANGE( 'M', M, N, A, LDA, WORK ) 00267 IASCL = 0 00268 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00269 * 00270 * Scale matrix norm up to SMLNUM 00271 * 00272 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO ) 00273 IASCL = 1 00274 ELSE IF( ANRM.GT.BIGNUM ) THEN 00275 * 00276 * Scale matrix norm down to BIGNUM 00277 * 00278 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO ) 00279 IASCL = 2 00280 ELSE IF( ANRM.EQ.ZERO ) THEN 00281 * 00282 * Matrix all zero. Return zero solution. 00283 * 00284 CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB ) 00285 CALL DLASET( 'F', MINMN, 1, ZERO, ZERO, S, MINMN ) 00286 RANK = 0 00287 GO TO 70 00288 END IF 00289 * 00290 * Scale B if max element outside range [SMLNUM,BIGNUM] 00291 * 00292 BNRM = DLANGE( 'M', M, NRHS, B, LDB, WORK ) 00293 IBSCL = 0 00294 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 00295 * 00296 * Scale matrix norm up to SMLNUM 00297 * 00298 CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO ) 00299 IBSCL = 1 00300 ELSE IF( BNRM.GT.BIGNUM ) THEN 00301 * 00302 * Scale matrix norm down to BIGNUM 00303 * 00304 CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO ) 00305 IBSCL = 2 00306 END IF 00307 * 00308 * Overdetermined case 00309 * 00310 IF( M.GE.N ) THEN 00311 * 00312 * Path 1 - overdetermined or exactly determined 00313 * 00314 MM = M 00315 IF( M.GE.MNTHR ) THEN 00316 * 00317 * Path 1a - overdetermined, with many more rows than columns 00318 * 00319 MM = N 00320 ITAU = 1 00321 IWORK = ITAU + N 00322 * 00323 * Compute A=Q*R 00324 * (Workspace: need 2*N, prefer N+N*NB) 00325 * 00326 CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ), 00327 $ LWORK-IWORK+1, INFO ) 00328 * 00329 * Multiply B by transpose(Q) 00330 * (Workspace: need N+NRHS, prefer N+NRHS*NB) 00331 * 00332 CALL DORMQR( 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAU ), B, 00333 $ LDB, WORK( IWORK ), LWORK-IWORK+1, INFO ) 00334 * 00335 * Zero out below R 00336 * 00337 IF( N.GT.1 ) 00338 $ CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA ) 00339 END IF 00340 * 00341 IE = 1 00342 ITAUQ = IE + N 00343 ITAUP = ITAUQ + N 00344 IWORK = ITAUP + N 00345 * 00346 * Bidiagonalize R in A 00347 * (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB) 00348 * 00349 CALL DGEBRD( MM, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 00350 $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, 00351 $ INFO ) 00352 * 00353 * Multiply B by transpose of left bidiagonalizing vectors of R 00354 * (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB) 00355 * 00356 CALL DORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, WORK( ITAUQ ), 00357 $ B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO ) 00358 * 00359 * Generate right bidiagonalizing vectors of R in A 00360 * (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) 00361 * 00362 CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), 00363 $ WORK( IWORK ), LWORK-IWORK+1, INFO ) 00364 IWORK = IE + N 00365 * 00366 * Perform bidiagonal QR iteration 00367 * multiply B by transpose of left singular vectors 00368 * compute right singular vectors in A 00369 * (Workspace: need BDSPAC) 00370 * 00371 CALL DBDSQR( 'U', N, N, 0, NRHS, S, WORK( IE ), A, LDA, VDUM, 00372 $ 1, B, LDB, WORK( IWORK ), INFO ) 00373 IF( INFO.NE.0 ) 00374 $ GO TO 70 00375 * 00376 * Multiply B by reciprocals of singular values 00377 * 00378 THR = MAX( RCOND*S( 1 ), SFMIN ) 00379 IF( RCOND.LT.ZERO ) 00380 $ THR = MAX( EPS*S( 1 ), SFMIN ) 00381 RANK = 0 00382 DO 10 I = 1, N 00383 IF( S( I ).GT.THR ) THEN 00384 CALL DRSCL( NRHS, S( I ), B( I, 1 ), LDB ) 00385 RANK = RANK + 1 00386 ELSE 00387 CALL DLASET( 'F', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB ) 00388 END IF 00389 10 CONTINUE 00390 * 00391 * Multiply B by right singular vectors 00392 * (Workspace: need N, prefer N*NRHS) 00393 * 00394 IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN 00395 CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, A, LDA, B, LDB, ZERO, 00396 $ WORK, LDB ) 00397 CALL DLACPY( 'G', N, NRHS, WORK, LDB, B, LDB ) 00398 ELSE IF( NRHS.GT.1 ) THEN 00399 CHUNK = LWORK / N 00400 DO 20 I = 1, NRHS, CHUNK 00401 BL = MIN( NRHS-I+1, CHUNK ) 00402 CALL DGEMM( 'T', 'N', N, BL, N, ONE, A, LDA, B( 1, I ), 00403 $ LDB, ZERO, WORK, N ) 00404 CALL DLACPY( 'G', N, BL, WORK, N, B( 1, I ), LDB ) 00405 20 CONTINUE 00406 ELSE 00407 CALL DGEMV( 'T', N, N, ONE, A, LDA, B, 1, ZERO, WORK, 1 ) 00408 CALL DCOPY( N, WORK, 1, B, 1 ) 00409 END IF 00410 * 00411 ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+ 00412 $ MAX( M, 2*M-4, NRHS, N-3*M ) ) THEN 00413 * 00414 * Path 2a - underdetermined, with many more columns than rows 00415 * and sufficient workspace for an efficient algorithm 00416 * 00417 LDWORK = M 00418 IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ), 00419 $ M*LDA+M+M*NRHS ) )LDWORK = LDA 00420 ITAU = 1 00421 IWORK = M + 1 00422 * 00423 * Compute A=L*Q 00424 * (Workspace: need 2*M, prefer M+M*NB) 00425 * 00426 CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ), 00427 $ LWORK-IWORK+1, INFO ) 00428 IL = IWORK 00429 * 00430 * Copy L to WORK(IL), zeroing out above it 00431 * 00432 CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK ) 00433 CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, WORK( IL+LDWORK ), 00434 $ LDWORK ) 00435 IE = IL + LDWORK*M 00436 ITAUQ = IE + M 00437 ITAUP = ITAUQ + M 00438 IWORK = ITAUP + M 00439 * 00440 * Bidiagonalize L in WORK(IL) 00441 * (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB) 00442 * 00443 CALL DGEBRD( M, M, WORK( IL ), LDWORK, S, WORK( IE ), 00444 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( IWORK ), 00445 $ LWORK-IWORK+1, INFO ) 00446 * 00447 * Multiply B by transpose of left bidiagonalizing vectors of L 00448 * (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB) 00449 * 00450 CALL DORMBR( 'Q', 'L', 'T', M, NRHS, M, WORK( IL ), LDWORK, 00451 $ WORK( ITAUQ ), B, LDB, WORK( IWORK ), 00452 $ LWORK-IWORK+1, INFO ) 00453 * 00454 * Generate right bidiagonalizing vectors of R in WORK(IL) 00455 * (Workspace: need M*M+5*M-1, prefer M*M+4*M+(M-1)*NB) 00456 * 00457 CALL DORGBR( 'P', M, M, M, WORK( IL ), LDWORK, WORK( ITAUP ), 00458 $ WORK( IWORK ), LWORK-IWORK+1, INFO ) 00459 IWORK = IE + M 00460 * 00461 * Perform bidiagonal QR iteration, 00462 * computing right singular vectors of L in WORK(IL) and 00463 * multiplying B by transpose of left singular vectors 00464 * (Workspace: need M*M+M+BDSPAC) 00465 * 00466 CALL DBDSQR( 'U', M, M, 0, NRHS, S, WORK( IE ), WORK( IL ), 00467 $ LDWORK, A, LDA, B, LDB, WORK( IWORK ), INFO ) 00468 IF( INFO.NE.0 ) 00469 $ GO TO 70 00470 * 00471 * Multiply B by reciprocals of singular values 00472 * 00473 THR = MAX( RCOND*S( 1 ), SFMIN ) 00474 IF( RCOND.LT.ZERO ) 00475 $ THR = MAX( EPS*S( 1 ), SFMIN ) 00476 RANK = 0 00477 DO 30 I = 1, M 00478 IF( S( I ).GT.THR ) THEN 00479 CALL DRSCL( NRHS, S( I ), B( I, 1 ), LDB ) 00480 RANK = RANK + 1 00481 ELSE 00482 CALL DLASET( 'F', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB ) 00483 END IF 00484 30 CONTINUE 00485 IWORK = IE 00486 * 00487 * Multiply B by right singular vectors of L in WORK(IL) 00488 * (Workspace: need M*M+2*M, prefer M*M+M+M*NRHS) 00489 * 00490 IF( LWORK.GE.LDB*NRHS+IWORK-1 .AND. NRHS.GT.1 ) THEN 00491 CALL DGEMM( 'T', 'N', M, NRHS, M, ONE, WORK( IL ), LDWORK, 00492 $ B, LDB, ZERO, WORK( IWORK ), LDB ) 00493 CALL DLACPY( 'G', M, NRHS, WORK( IWORK ), LDB, B, LDB ) 00494 ELSE IF( NRHS.GT.1 ) THEN 00495 CHUNK = ( LWORK-IWORK+1 ) / M 00496 DO 40 I = 1, NRHS, CHUNK 00497 BL = MIN( NRHS-I+1, CHUNK ) 00498 CALL DGEMM( 'T', 'N', M, BL, M, ONE, WORK( IL ), LDWORK, 00499 $ B( 1, I ), LDB, ZERO, WORK( IWORK ), M ) 00500 CALL DLACPY( 'G', M, BL, WORK( IWORK ), M, B( 1, I ), 00501 $ LDB ) 00502 40 CONTINUE 00503 ELSE 00504 CALL DGEMV( 'T', M, M, ONE, WORK( IL ), LDWORK, B( 1, 1 ), 00505 $ 1, ZERO, WORK( IWORK ), 1 ) 00506 CALL DCOPY( M, WORK( IWORK ), 1, B( 1, 1 ), 1 ) 00507 END IF 00508 * 00509 * Zero out below first M rows of B 00510 * 00511 CALL DLASET( 'F', N-M, NRHS, ZERO, ZERO, B( M+1, 1 ), LDB ) 00512 IWORK = ITAU + M 00513 * 00514 * Multiply transpose(Q) by B 00515 * (Workspace: need M+NRHS, prefer M+NRHS*NB) 00516 * 00517 CALL DORMLQ( 'L', 'T', N, NRHS, M, A, LDA, WORK( ITAU ), B, 00518 $ LDB, WORK( IWORK ), LWORK-IWORK+1, INFO ) 00519 * 00520 ELSE 00521 * 00522 * Path 2 - remaining underdetermined cases 00523 * 00524 IE = 1 00525 ITAUQ = IE + M 00526 ITAUP = ITAUQ + M 00527 IWORK = ITAUP + M 00528 * 00529 * Bidiagonalize A 00530 * (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB) 00531 * 00532 CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 00533 $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, 00534 $ INFO ) 00535 * 00536 * Multiply B by transpose of left bidiagonalizing vectors 00537 * (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB) 00538 * 00539 CALL DORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAUQ ), 00540 $ B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO ) 00541 * 00542 * Generate right bidiagonalizing vectors in A 00543 * (Workspace: need 4*M, prefer 3*M+M*NB) 00544 * 00545 CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), 00546 $ WORK( IWORK ), LWORK-IWORK+1, INFO ) 00547 IWORK = IE + M 00548 * 00549 * Perform bidiagonal QR iteration, 00550 * computing right singular vectors of A in A and 00551 * multiplying B by transpose of left singular vectors 00552 * (Workspace: need BDSPAC) 00553 * 00554 CALL DBDSQR( 'L', M, N, 0, NRHS, S, WORK( IE ), A, LDA, VDUM, 00555 $ 1, B, LDB, WORK( IWORK ), INFO ) 00556 IF( INFO.NE.0 ) 00557 $ GO TO 70 00558 * 00559 * Multiply B by reciprocals of singular values 00560 * 00561 THR = MAX( RCOND*S( 1 ), SFMIN ) 00562 IF( RCOND.LT.ZERO ) 00563 $ THR = MAX( EPS*S( 1 ), SFMIN ) 00564 RANK = 0 00565 DO 50 I = 1, M 00566 IF( S( I ).GT.THR ) THEN 00567 CALL DRSCL( NRHS, S( I ), B( I, 1 ), LDB ) 00568 RANK = RANK + 1 00569 ELSE 00570 CALL DLASET( 'F', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB ) 00571 END IF 00572 50 CONTINUE 00573 * 00574 * Multiply B by right singular vectors of A 00575 * (Workspace: need N, prefer N*NRHS) 00576 * 00577 IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN 00578 CALL DGEMM( 'T', 'N', N, NRHS, M, ONE, A, LDA, B, LDB, ZERO, 00579 $ WORK, LDB ) 00580 CALL DLACPY( 'F', N, NRHS, WORK, LDB, B, LDB ) 00581 ELSE IF( NRHS.GT.1 ) THEN 00582 CHUNK = LWORK / N 00583 DO 60 I = 1, NRHS, CHUNK 00584 BL = MIN( NRHS-I+1, CHUNK ) 00585 CALL DGEMM( 'T', 'N', N, BL, M, ONE, A, LDA, B( 1, I ), 00586 $ LDB, ZERO, WORK, N ) 00587 CALL DLACPY( 'F', N, BL, WORK, N, B( 1, I ), LDB ) 00588 60 CONTINUE 00589 ELSE 00590 CALL DGEMV( 'T', M, N, ONE, A, LDA, B, 1, ZERO, WORK, 1 ) 00591 CALL DCOPY( N, WORK, 1, B, 1 ) 00592 END IF 00593 END IF 00594 * 00595 * Undo scaling 00596 * 00597 IF( IASCL.EQ.1 ) THEN 00598 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO ) 00599 CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN, 00600 $ INFO ) 00601 ELSE IF( IASCL.EQ.2 ) THEN 00602 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO ) 00603 CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN, 00604 $ INFO ) 00605 END IF 00606 IF( IBSCL.EQ.1 ) THEN 00607 CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO ) 00608 ELSE IF( IBSCL.EQ.2 ) THEN 00609 CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO ) 00610 END IF 00611 * 00612 70 CONTINUE 00613 WORK( 1 ) = MAXWRK 00614 RETURN 00615 * 00616 * End of DGELSS 00617 * 00618 END
1.7.3 
