![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> DGELS 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 DGELS + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgels.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgels.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgels.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DGELS( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, 00022 * INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER TRANS 00026 * INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS 00027 * .. 00028 * .. Array Arguments .. 00029 * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * ) 00030 * .. 00031 * 00032 * 00033 *> \par Purpose: 00034 * ============= 00035 *> 00036 *> \verbatim 00037 *> 00038 *> DGELS solves overdetermined or underdetermined real linear systems 00039 *> involving an M-by-N matrix A, or its transpose, using a QR or LQ 00040 *> factorization of A. It is assumed that A has full rank. 00041 *> 00042 *> The following options are provided: 00043 *> 00044 *> 1. If TRANS = 'N' and m >= n: find the least squares solution of 00045 *> an overdetermined system, i.e., solve the least squares problem 00046 *> minimize || B - A*X ||. 00047 *> 00048 *> 2. If TRANS = 'N' and m < n: find the minimum norm solution of 00049 *> an underdetermined system A * X = B. 00050 *> 00051 *> 3. If TRANS = 'T' and m >= n: find the minimum norm solution of 00052 *> an undetermined system A**T * X = B. 00053 *> 00054 *> 4. If TRANS = 'T' and m < n: find the least squares solution of 00055 *> an overdetermined system, i.e., solve the least squares problem 00056 *> minimize || B - A**T * X ||. 00057 *> 00058 *> Several right hand side vectors b and solution vectors x can be 00059 *> handled in a single call; they are stored as the columns of the 00060 *> M-by-NRHS right hand side matrix B and the N-by-NRHS solution 00061 *> matrix X. 00062 *> \endverbatim 00063 * 00064 * Arguments: 00065 * ========== 00066 * 00067 *> \param[in] TRANS 00068 *> \verbatim 00069 *> TRANS is CHARACTER*1 00070 *> = 'N': the linear system involves A; 00071 *> = 'T': the linear system involves A**T. 00072 *> \endverbatim 00073 *> 00074 *> \param[in] M 00075 *> \verbatim 00076 *> M is INTEGER 00077 *> The number of rows of the matrix A. M >= 0. 00078 *> \endverbatim 00079 *> 00080 *> \param[in] N 00081 *> \verbatim 00082 *> N is INTEGER 00083 *> The number of columns of the matrix A. N >= 0. 00084 *> \endverbatim 00085 *> 00086 *> \param[in] NRHS 00087 *> \verbatim 00088 *> NRHS is INTEGER 00089 *> The number of right hand sides, i.e., the number of 00090 *> columns of the matrices B and X. NRHS >=0. 00091 *> \endverbatim 00092 *> 00093 *> \param[in,out] A 00094 *> \verbatim 00095 *> A is DOUBLE PRECISION array, dimension (LDA,N) 00096 *> On entry, the M-by-N matrix A. 00097 *> On exit, 00098 *> if M >= N, A is overwritten by details of its QR 00099 *> factorization as returned by DGEQRF; 00100 *> if M < N, A is overwritten by details of its LQ 00101 *> factorization as returned by DGELQF. 00102 *> \endverbatim 00103 *> 00104 *> \param[in] LDA 00105 *> \verbatim 00106 *> LDA is INTEGER 00107 *> The leading dimension of the array A. LDA >= max(1,M). 00108 *> \endverbatim 00109 *> 00110 *> \param[in,out] B 00111 *> \verbatim 00112 *> B is DOUBLE PRECISION array, dimension (LDB,NRHS) 00113 *> On entry, the matrix B of right hand side vectors, stored 00114 *> columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS 00115 *> if TRANS = 'T'. 00116 *> On exit, if INFO = 0, B is overwritten by the solution 00117 *> vectors, stored columnwise: 00118 *> if TRANS = 'N' and m >= n, rows 1 to n of B contain the least 00119 *> squares solution vectors; the residual sum of squares for the 00120 *> solution in each column is given by the sum of squares of 00121 *> elements N+1 to M in that column; 00122 *> if TRANS = 'N' and m < n, rows 1 to N of B contain the 00123 *> minimum norm solution vectors; 00124 *> if TRANS = 'T' and m >= n, rows 1 to M of B contain the 00125 *> minimum norm solution vectors; 00126 *> if TRANS = 'T' and m < n, rows 1 to M of B contain the 00127 *> least squares solution vectors; the residual sum of squares 00128 *> for the solution in each column is given by the sum of 00129 *> squares of elements M+1 to N in that column. 00130 *> \endverbatim 00131 *> 00132 *> \param[in] LDB 00133 *> \verbatim 00134 *> LDB is INTEGER 00135 *> The leading dimension of the array B. LDB >= MAX(1,M,N). 00136 *> \endverbatim 00137 *> 00138 *> \param[out] WORK 00139 *> \verbatim 00140 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 00141 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00142 *> \endverbatim 00143 *> 00144 *> \param[in] LWORK 00145 *> \verbatim 00146 *> LWORK is INTEGER 00147 *> The dimension of the array WORK. 00148 *> LWORK >= max( 1, MN + max( MN, NRHS ) ). 00149 *> For optimal performance, 00150 *> LWORK >= max( 1, MN + max( MN, NRHS )*NB ). 00151 *> where MN = min(M,N) and NB is the optimum block size. 00152 *> 00153 *> If LWORK = -1, then a workspace query is assumed; the routine 00154 *> only calculates the optimal size of the WORK array, returns 00155 *> this value as the first entry of the WORK array, and no error 00156 *> message related to LWORK is issued by XERBLA. 00157 *> \endverbatim 00158 *> 00159 *> \param[out] INFO 00160 *> \verbatim 00161 *> INFO is INTEGER 00162 *> = 0: successful exit 00163 *> < 0: if INFO = -i, the i-th argument had an illegal value 00164 *> > 0: if INFO = i, the i-th diagonal element of the 00165 *> triangular factor of A is zero, so that A does not have 00166 *> full rank; the least squares solution could not be 00167 *> computed. 00168 *> \endverbatim 00169 * 00170 * Authors: 00171 * ======== 00172 * 00173 *> \author Univ. of Tennessee 00174 *> \author Univ. of California Berkeley 00175 *> \author Univ. of Colorado Denver 00176 *> \author NAG Ltd. 00177 * 00178 *> \date November 2011 00179 * 00180 *> \ingroup doubleGEsolve 00181 * 00182 * ===================================================================== 00183 SUBROUTINE DGELS( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, 00184 $ INFO ) 00185 * 00186 * -- LAPACK driver routine (version 3.4.0) -- 00187 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00188 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00189 * November 2011 00190 * 00191 * .. Scalar Arguments .. 00192 CHARACTER TRANS 00193 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS 00194 * .. 00195 * .. Array Arguments .. 00196 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * ) 00197 * .. 00198 * 00199 * ===================================================================== 00200 * 00201 * .. Parameters .. 00202 DOUBLE PRECISION ZERO, ONE 00203 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00204 * .. 00205 * .. Local Scalars .. 00206 LOGICAL LQUERY, TPSD 00207 INTEGER BROW, I, IASCL, IBSCL, J, MN, NB, SCLLEN, WSIZE 00208 DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM 00209 * .. 00210 * .. Local Arrays .. 00211 DOUBLE PRECISION RWORK( 1 ) 00212 * .. 00213 * .. External Functions .. 00214 LOGICAL LSAME 00215 INTEGER ILAENV 00216 DOUBLE PRECISION DLAMCH, DLANGE 00217 EXTERNAL LSAME, ILAENV, DLABAD, DLAMCH, DLANGE 00218 * .. 00219 * .. External Subroutines .. 00220 EXTERNAL DGELQF, DGEQRF, DLASCL, DLASET, DORMLQ, DORMQR, 00221 $ DTRTRS, XERBLA 00222 * .. 00223 * .. Intrinsic Functions .. 00224 INTRINSIC DBLE, MAX, MIN 00225 * .. 00226 * .. Executable Statements .. 00227 * 00228 * Test the input arguments. 00229 * 00230 INFO = 0 00231 MN = MIN( M, N ) 00232 LQUERY = ( LWORK.EQ.-1 ) 00233 IF( .NOT.( LSAME( TRANS, 'N' ) .OR. LSAME( TRANS, 'T' ) ) ) THEN 00234 INFO = -1 00235 ELSE IF( M.LT.0 ) THEN 00236 INFO = -2 00237 ELSE IF( N.LT.0 ) THEN 00238 INFO = -3 00239 ELSE IF( NRHS.LT.0 ) THEN 00240 INFO = -4 00241 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00242 INFO = -6 00243 ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN 00244 INFO = -8 00245 ELSE IF( LWORK.LT.MAX( 1, MN+MAX( MN, NRHS ) ) .AND. .NOT.LQUERY ) 00246 $ THEN 00247 INFO = -10 00248 END IF 00249 * 00250 * Figure out optimal block size 00251 * 00252 IF( INFO.EQ.0 .OR. INFO.EQ.-10 ) THEN 00253 * 00254 TPSD = .TRUE. 00255 IF( LSAME( TRANS, 'N' ) ) 00256 $ TPSD = .FALSE. 00257 * 00258 IF( M.GE.N ) THEN 00259 NB = ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) 00260 IF( TPSD ) THEN 00261 NB = MAX( NB, ILAENV( 1, 'DORMQR', 'LN', M, NRHS, N, 00262 $ -1 ) ) 00263 ELSE 00264 NB = MAX( NB, ILAENV( 1, 'DORMQR', 'LT', M, NRHS, N, 00265 $ -1 ) ) 00266 END IF 00267 ELSE 00268 NB = ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) 00269 IF( TPSD ) THEN 00270 NB = MAX( NB, ILAENV( 1, 'DORMLQ', 'LT', N, NRHS, M, 00271 $ -1 ) ) 00272 ELSE 00273 NB = MAX( NB, ILAENV( 1, 'DORMLQ', 'LN', N, NRHS, M, 00274 $ -1 ) ) 00275 END IF 00276 END IF 00277 * 00278 WSIZE = MAX( 1, MN+MAX( MN, NRHS )*NB ) 00279 WORK( 1 ) = DBLE( WSIZE ) 00280 * 00281 END IF 00282 * 00283 IF( INFO.NE.0 ) THEN 00284 CALL XERBLA( 'DGELS ', -INFO ) 00285 RETURN 00286 ELSE IF( LQUERY ) THEN 00287 RETURN 00288 END IF 00289 * 00290 * Quick return if possible 00291 * 00292 IF( MIN( M, N, NRHS ).EQ.0 ) THEN 00293 CALL DLASET( 'Full', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB ) 00294 RETURN 00295 END IF 00296 * 00297 * Get machine parameters 00298 * 00299 SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' ) 00300 BIGNUM = ONE / SMLNUM 00301 CALL DLABAD( SMLNUM, BIGNUM ) 00302 * 00303 * Scale A, B if max element outside range [SMLNUM,BIGNUM] 00304 * 00305 ANRM = DLANGE( 'M', M, N, A, LDA, RWORK ) 00306 IASCL = 0 00307 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00308 * 00309 * Scale matrix norm up to SMLNUM 00310 * 00311 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO ) 00312 IASCL = 1 00313 ELSE IF( ANRM.GT.BIGNUM ) THEN 00314 * 00315 * Scale matrix norm down to BIGNUM 00316 * 00317 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO ) 00318 IASCL = 2 00319 ELSE IF( ANRM.EQ.ZERO ) THEN 00320 * 00321 * Matrix all zero. Return zero solution. 00322 * 00323 CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB ) 00324 GO TO 50 00325 END IF 00326 * 00327 BROW = M 00328 IF( TPSD ) 00329 $ BROW = N 00330 BNRM = DLANGE( 'M', BROW, NRHS, B, LDB, RWORK ) 00331 IBSCL = 0 00332 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 00333 * 00334 * Scale matrix norm up to SMLNUM 00335 * 00336 CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, BROW, NRHS, B, LDB, 00337 $ INFO ) 00338 IBSCL = 1 00339 ELSE IF( BNRM.GT.BIGNUM ) THEN 00340 * 00341 * Scale matrix norm down to BIGNUM 00342 * 00343 CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, BROW, NRHS, B, LDB, 00344 $ INFO ) 00345 IBSCL = 2 00346 END IF 00347 * 00348 IF( M.GE.N ) THEN 00349 * 00350 * compute QR factorization of A 00351 * 00352 CALL DGEQRF( M, N, A, LDA, WORK( 1 ), WORK( MN+1 ), LWORK-MN, 00353 $ INFO ) 00354 * 00355 * workspace at least N, optimally N*NB 00356 * 00357 IF( .NOT.TPSD ) THEN 00358 * 00359 * Least-Squares Problem min || A * X - B || 00360 * 00361 * B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS) 00362 * 00363 CALL DORMQR( 'Left', 'Transpose', M, NRHS, N, A, LDA, 00364 $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN, 00365 $ INFO ) 00366 * 00367 * workspace at least NRHS, optimally NRHS*NB 00368 * 00369 * B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS) 00370 * 00371 CALL DTRTRS( 'Upper', 'No transpose', 'Non-unit', N, NRHS, 00372 $ A, LDA, B, LDB, INFO ) 00373 * 00374 IF( INFO.GT.0 ) THEN 00375 RETURN 00376 END IF 00377 * 00378 SCLLEN = N 00379 * 00380 ELSE 00381 * 00382 * Overdetermined system of equations A**T * X = B 00383 * 00384 * B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS) 00385 * 00386 CALL DTRTRS( 'Upper', 'Transpose', 'Non-unit', N, NRHS, 00387 $ A, LDA, B, LDB, INFO ) 00388 * 00389 IF( INFO.GT.0 ) THEN 00390 RETURN 00391 END IF 00392 * 00393 * B(N+1:M,1:NRHS) = ZERO 00394 * 00395 DO 20 J = 1, NRHS 00396 DO 10 I = N + 1, M 00397 B( I, J ) = ZERO 00398 10 CONTINUE 00399 20 CONTINUE 00400 * 00401 * B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS) 00402 * 00403 CALL DORMQR( 'Left', 'No transpose', M, NRHS, N, A, LDA, 00404 $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN, 00405 $ INFO ) 00406 * 00407 * workspace at least NRHS, optimally NRHS*NB 00408 * 00409 SCLLEN = M 00410 * 00411 END IF 00412 * 00413 ELSE 00414 * 00415 * Compute LQ factorization of A 00416 * 00417 CALL DGELQF( M, N, A, LDA, WORK( 1 ), WORK( MN+1 ), LWORK-MN, 00418 $ INFO ) 00419 * 00420 * workspace at least M, optimally M*NB. 00421 * 00422 IF( .NOT.TPSD ) THEN 00423 * 00424 * underdetermined system of equations A * X = B 00425 * 00426 * B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS) 00427 * 00428 CALL DTRTRS( 'Lower', 'No transpose', 'Non-unit', M, NRHS, 00429 $ A, LDA, B, LDB, INFO ) 00430 * 00431 IF( INFO.GT.0 ) THEN 00432 RETURN 00433 END IF 00434 * 00435 * B(M+1:N,1:NRHS) = 0 00436 * 00437 DO 40 J = 1, NRHS 00438 DO 30 I = M + 1, N 00439 B( I, J ) = ZERO 00440 30 CONTINUE 00441 40 CONTINUE 00442 * 00443 * B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS) 00444 * 00445 CALL DORMLQ( 'Left', 'Transpose', N, NRHS, M, A, LDA, 00446 $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN, 00447 $ INFO ) 00448 * 00449 * workspace at least NRHS, optimally NRHS*NB 00450 * 00451 SCLLEN = N 00452 * 00453 ELSE 00454 * 00455 * overdetermined system min || A**T * X - B || 00456 * 00457 * B(1:N,1:NRHS) := Q * B(1:N,1:NRHS) 00458 * 00459 CALL DORMLQ( 'Left', 'No transpose', N, NRHS, M, A, LDA, 00460 $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN, 00461 $ INFO ) 00462 * 00463 * workspace at least NRHS, optimally NRHS*NB 00464 * 00465 * B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS) 00466 * 00467 CALL DTRTRS( 'Lower', 'Transpose', 'Non-unit', M, NRHS, 00468 $ A, LDA, B, LDB, INFO ) 00469 * 00470 IF( INFO.GT.0 ) THEN 00471 RETURN 00472 END IF 00473 * 00474 SCLLEN = M 00475 * 00476 END IF 00477 * 00478 END IF 00479 * 00480 * Undo scaling 00481 * 00482 IF( IASCL.EQ.1 ) THEN 00483 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB, 00484 $ INFO ) 00485 ELSE IF( IASCL.EQ.2 ) THEN 00486 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB, 00487 $ INFO ) 00488 END IF 00489 IF( IBSCL.EQ.1 ) THEN 00490 CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB, 00491 $ INFO ) 00492 ELSE IF( IBSCL.EQ.2 ) THEN 00493 CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB, 00494 $ INFO ) 00495 END IF 00496 * 00497 50 CONTINUE 00498 WORK( 1 ) = DBLE( WSIZE ) 00499 * 00500 RETURN 00501 * 00502 * End of DGELS 00503 * 00504 END