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