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