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