![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> SGELSY 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 SGELSY + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgelsy.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgelsy.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgelsy.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, 00022 * WORK, LWORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK 00026 * REAL RCOND 00027 * .. 00028 * .. Array Arguments .. 00029 * INTEGER JPVT( * ) 00030 * REAL A( LDA, * ), B( LDB, * ), WORK( * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> SGELSY computes the minimum-norm solution to a real linear least 00040 *> squares problem: 00041 *> minimize || A * X - B || 00042 *> using a complete orthogonal factorization of A. A is an M-by-N 00043 *> matrix which may be rank-deficient. 00044 *> 00045 *> Several right hand side vectors b and solution vectors x can be 00046 *> handled in a single call; they are stored as the columns of the 00047 *> M-by-NRHS right hand side matrix B and the N-by-NRHS solution 00048 *> matrix X. 00049 *> 00050 *> The routine first computes a QR factorization with column pivoting: 00051 *> A * P = Q * [ R11 R12 ] 00052 *> [ 0 R22 ] 00053 *> with R11 defined as the largest leading submatrix whose estimated 00054 *> condition number is less than 1/RCOND. The order of R11, RANK, 00055 *> is the effective rank of A. 00056 *> 00057 *> Then, R22 is considered to be negligible, and R12 is annihilated 00058 *> by orthogonal transformations from the right, arriving at the 00059 *> complete orthogonal factorization: 00060 *> A * P = Q * [ T11 0 ] * Z 00061 *> [ 0 0 ] 00062 *> The minimum-norm solution is then 00063 *> X = P * Z**T [ inv(T11)*Q1**T*B ] 00064 *> [ 0 ] 00065 *> where Q1 consists of the first RANK columns of Q. 00066 *> 00067 *> This routine is basically identical to the original xGELSX except 00068 *> three differences: 00069 *> o The call to the subroutine xGEQPF has been substituted by the 00070 *> the call to the subroutine xGEQP3. This subroutine is a Blas-3 00071 *> version of the QR factorization with column pivoting. 00072 *> o Matrix B (the right hand side) is updated with Blas-3. 00073 *> o The permutation of matrix B (the right hand side) is faster and 00074 *> more simple. 00075 *> \endverbatim 00076 * 00077 * Arguments: 00078 * ========== 00079 * 00080 *> \param[in] M 00081 *> \verbatim 00082 *> M is INTEGER 00083 *> The number of rows of the matrix A. M >= 0. 00084 *> \endverbatim 00085 *> 00086 *> \param[in] N 00087 *> \verbatim 00088 *> N is INTEGER 00089 *> The number of columns of the matrix A. N >= 0. 00090 *> \endverbatim 00091 *> 00092 *> \param[in] NRHS 00093 *> \verbatim 00094 *> NRHS is INTEGER 00095 *> The number of right hand sides, i.e., the number of 00096 *> columns of matrices B and X. NRHS >= 0. 00097 *> \endverbatim 00098 *> 00099 *> \param[in,out] A 00100 *> \verbatim 00101 *> A is REAL array, dimension (LDA,N) 00102 *> On entry, the M-by-N matrix A. 00103 *> On exit, A has been overwritten by details of its 00104 *> complete orthogonal factorization. 00105 *> \endverbatim 00106 *> 00107 *> \param[in] LDA 00108 *> \verbatim 00109 *> LDA is INTEGER 00110 *> The leading dimension of the array A. LDA >= max(1,M). 00111 *> \endverbatim 00112 *> 00113 *> \param[in,out] B 00114 *> \verbatim 00115 *> B is REAL array, dimension (LDB,NRHS) 00116 *> On entry, the M-by-NRHS right hand side matrix B. 00117 *> On exit, the N-by-NRHS solution matrix X. 00118 *> \endverbatim 00119 *> 00120 *> \param[in] LDB 00121 *> \verbatim 00122 *> LDB is INTEGER 00123 *> The leading dimension of the array B. LDB >= max(1,M,N). 00124 *> \endverbatim 00125 *> 00126 *> \param[in,out] JPVT 00127 *> \verbatim 00128 *> JPVT is INTEGER array, dimension (N) 00129 *> On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted 00130 *> to the front of AP, otherwise column i is a free column. 00131 *> On exit, if JPVT(i) = k, then the i-th column of AP 00132 *> was the k-th column of A. 00133 *> \endverbatim 00134 *> 00135 *> \param[in] RCOND 00136 *> \verbatim 00137 *> RCOND is REAL 00138 *> RCOND is used to determine the effective rank of A, which 00139 *> is defined as the order of the largest leading triangular 00140 *> submatrix R11 in the QR factorization with pivoting of A, 00141 *> whose estimated condition number < 1/RCOND. 00142 *> \endverbatim 00143 *> 00144 *> \param[out] RANK 00145 *> \verbatim 00146 *> RANK is INTEGER 00147 *> The effective rank of A, i.e., the order of the submatrix 00148 *> R11. This is the same as the order of the submatrix T11 00149 *> in the complete orthogonal factorization of A. 00150 *> \endverbatim 00151 *> 00152 *> \param[out] WORK 00153 *> \verbatim 00154 *> WORK is REAL array, dimension (MAX(1,LWORK)) 00155 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00156 *> \endverbatim 00157 *> 00158 *> \param[in] LWORK 00159 *> \verbatim 00160 *> LWORK is INTEGER 00161 *> The dimension of the array WORK. 00162 *> The unblocked strategy requires that: 00163 *> LWORK >= MAX( MN+3*N+1, 2*MN+NRHS ), 00164 *> where MN = min( M, N ). 00165 *> The block algorithm requires that: 00166 *> LWORK >= MAX( MN+2*N+NB*(N+1), 2*MN+NB*NRHS ), 00167 *> where NB is an upper bound on the blocksize returned 00168 *> by ILAENV for the routines SGEQP3, STZRZF, STZRQF, SORMQR, 00169 *> and SORMRZ. 00170 *> 00171 *> If LWORK = -1, then a workspace query is assumed; the routine 00172 *> only calculates the optimal size of the WORK array, returns 00173 *> this value as the first entry of the WORK array, and no error 00174 *> message related to LWORK is issued by XERBLA. 00175 *> \endverbatim 00176 *> 00177 *> \param[out] INFO 00178 *> \verbatim 00179 *> INFO is INTEGER 00180 *> = 0: successful exit 00181 *> < 0: If INFO = -i, the i-th argument had an illegal value. 00182 *> \endverbatim 00183 * 00184 * Authors: 00185 * ======== 00186 * 00187 *> \author Univ. of Tennessee 00188 *> \author Univ. of California Berkeley 00189 *> \author Univ. of Colorado Denver 00190 *> \author NAG Ltd. 00191 * 00192 *> \date November 2011 00193 * 00194 *> \ingroup realGEsolve 00195 * 00196 *> \par Contributors: 00197 * ================== 00198 *> 00199 *> A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA \n 00200 *> E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n 00201 *> G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n 00202 *> 00203 * ===================================================================== 00204 SUBROUTINE SGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, 00205 $ WORK, LWORK, INFO ) 00206 * 00207 * -- LAPACK driver routine (version 3.4.0) -- 00208 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00209 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00210 * November 2011 00211 * 00212 * .. Scalar Arguments .. 00213 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK 00214 REAL RCOND 00215 * .. 00216 * .. Array Arguments .. 00217 INTEGER JPVT( * ) 00218 REAL A( LDA, * ), B( LDB, * ), WORK( * ) 00219 * .. 00220 * 00221 * ===================================================================== 00222 * 00223 * .. Parameters .. 00224 INTEGER IMAX, IMIN 00225 PARAMETER ( IMAX = 1, IMIN = 2 ) 00226 REAL ZERO, ONE 00227 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00228 * .. 00229 * .. Local Scalars .. 00230 LOGICAL LQUERY 00231 INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, LWKMIN, 00232 $ LWKOPT, MN, NB, NB1, NB2, NB3, NB4 00233 REAL ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX, 00234 $ SMAXPR, SMIN, SMINPR, SMLNUM, WSIZE 00235 * .. 00236 * .. External Functions .. 00237 INTEGER ILAENV 00238 REAL SLAMCH, SLANGE 00239 EXTERNAL ILAENV, SLAMCH, SLANGE 00240 * .. 00241 * .. External Subroutines .. 00242 EXTERNAL SCOPY, SGEQP3, SLABAD, SLAIC1, SLASCL, SLASET, 00243 $ SORMQR, SORMRZ, STRSM, STZRZF, XERBLA 00244 * .. 00245 * .. Intrinsic Functions .. 00246 INTRINSIC ABS, MAX, MIN 00247 * .. 00248 * .. Executable Statements .. 00249 * 00250 MN = MIN( M, N ) 00251 ISMIN = MN + 1 00252 ISMAX = 2*MN + 1 00253 * 00254 * Test the input arguments. 00255 * 00256 INFO = 0 00257 LQUERY = ( LWORK.EQ.-1 ) 00258 IF( M.LT.0 ) THEN 00259 INFO = -1 00260 ELSE IF( N.LT.0 ) THEN 00261 INFO = -2 00262 ELSE IF( NRHS.LT.0 ) THEN 00263 INFO = -3 00264 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00265 INFO = -5 00266 ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN 00267 INFO = -7 00268 END IF 00269 * 00270 * Figure out optimal block size 00271 * 00272 IF( INFO.EQ.0 ) THEN 00273 IF( MN.EQ.0 .OR. NRHS.EQ.0 ) THEN 00274 LWKMIN = 1 00275 LWKOPT = 1 00276 ELSE 00277 NB1 = ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 ) 00278 NB2 = ILAENV( 1, 'SGERQF', ' ', M, N, -1, -1 ) 00279 NB3 = ILAENV( 1, 'SORMQR', ' ', M, N, NRHS, -1 ) 00280 NB4 = ILAENV( 1, 'SORMRQ', ' ', M, N, NRHS, -1 ) 00281 NB = MAX( NB1, NB2, NB3, NB4 ) 00282 LWKMIN = MN + MAX( 2*MN, N + 1, MN + NRHS ) 00283 LWKOPT = MAX( LWKMIN, 00284 $ MN + 2*N + NB*( N + 1 ), 2*MN + NB*NRHS ) 00285 END IF 00286 WORK( 1 ) = LWKOPT 00287 * 00288 IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN 00289 INFO = -12 00290 END IF 00291 END IF 00292 * 00293 IF( INFO.NE.0 ) THEN 00294 CALL XERBLA( 'SGELSY', -INFO ) 00295 RETURN 00296 ELSE IF( LQUERY ) THEN 00297 RETURN 00298 END IF 00299 * 00300 * Quick return if possible 00301 * 00302 IF( MN.EQ.0 .OR. NRHS.EQ.0 ) THEN 00303 RANK = 0 00304 RETURN 00305 END IF 00306 * 00307 * Get machine parameters 00308 * 00309 SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' ) 00310 BIGNUM = ONE / SMLNUM 00311 CALL SLABAD( SMLNUM, BIGNUM ) 00312 * 00313 * Scale A, B if max entries outside range [SMLNUM,BIGNUM] 00314 * 00315 ANRM = SLANGE( 'M', M, N, A, LDA, WORK ) 00316 IASCL = 0 00317 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00318 * 00319 * Scale matrix norm up to SMLNUM 00320 * 00321 CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO ) 00322 IASCL = 1 00323 ELSE IF( ANRM.GT.BIGNUM ) THEN 00324 * 00325 * Scale matrix norm down to BIGNUM 00326 * 00327 CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO ) 00328 IASCL = 2 00329 ELSE IF( ANRM.EQ.ZERO ) THEN 00330 * 00331 * Matrix all zero. Return zero solution. 00332 * 00333 CALL SLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB ) 00334 RANK = 0 00335 GO TO 70 00336 END IF 00337 * 00338 BNRM = SLANGE( 'M', M, NRHS, B, LDB, WORK ) 00339 IBSCL = 0 00340 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 00341 * 00342 * Scale matrix norm up to SMLNUM 00343 * 00344 CALL SLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO ) 00345 IBSCL = 1 00346 ELSE IF( BNRM.GT.BIGNUM ) THEN 00347 * 00348 * Scale matrix norm down to BIGNUM 00349 * 00350 CALL SLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO ) 00351 IBSCL = 2 00352 END IF 00353 * 00354 * Compute QR factorization with column pivoting of A: 00355 * A * P = Q * R 00356 * 00357 CALL SGEQP3( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ), 00358 $ LWORK-MN, INFO ) 00359 WSIZE = MN + WORK( MN+1 ) 00360 * 00361 * workspace: MN+2*N+NB*(N+1). 00362 * Details of Householder rotations stored in WORK(1:MN). 00363 * 00364 * Determine RANK using incremental condition estimation 00365 * 00366 WORK( ISMIN ) = ONE 00367 WORK( ISMAX ) = ONE 00368 SMAX = ABS( A( 1, 1 ) ) 00369 SMIN = SMAX 00370 IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN 00371 RANK = 0 00372 CALL SLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB ) 00373 GO TO 70 00374 ELSE 00375 RANK = 1 00376 END IF 00377 * 00378 10 CONTINUE 00379 IF( RANK.LT.MN ) THEN 00380 I = RANK + 1 00381 CALL SLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ), 00382 $ A( I, I ), SMINPR, S1, C1 ) 00383 CALL SLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ), 00384 $ A( I, I ), SMAXPR, S2, C2 ) 00385 * 00386 IF( SMAXPR*RCOND.LE.SMINPR ) THEN 00387 DO 20 I = 1, RANK 00388 WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 ) 00389 WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 ) 00390 20 CONTINUE 00391 WORK( ISMIN+RANK ) = C1 00392 WORK( ISMAX+RANK ) = C2 00393 SMIN = SMINPR 00394 SMAX = SMAXPR 00395 RANK = RANK + 1 00396 GO TO 10 00397 END IF 00398 END IF 00399 * 00400 * workspace: 3*MN. 00401 * 00402 * Logically partition R = [ R11 R12 ] 00403 * [ 0 R22 ] 00404 * where R11 = R(1:RANK,1:RANK) 00405 * 00406 * [R11,R12] = [ T11, 0 ] * Y 00407 * 00408 IF( RANK.LT.N ) 00409 $ CALL STZRZF( RANK, N, A, LDA, WORK( MN+1 ), WORK( 2*MN+1 ), 00410 $ LWORK-2*MN, INFO ) 00411 * 00412 * workspace: 2*MN. 00413 * Details of Householder rotations stored in WORK(MN+1:2*MN) 00414 * 00415 * B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS) 00416 * 00417 CALL SORMQR( 'Left', 'Transpose', M, NRHS, MN, A, LDA, WORK( 1 ), 00418 $ B, LDB, WORK( 2*MN+1 ), LWORK-2*MN, INFO ) 00419 WSIZE = MAX( WSIZE, 2*MN+WORK( 2*MN+1 ) ) 00420 * 00421 * workspace: 2*MN+NB*NRHS. 00422 * 00423 * B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS) 00424 * 00425 CALL STRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK, 00426 $ NRHS, ONE, A, LDA, B, LDB ) 00427 * 00428 DO 40 J = 1, NRHS 00429 DO 30 I = RANK + 1, N 00430 B( I, J ) = ZERO 00431 30 CONTINUE 00432 40 CONTINUE 00433 * 00434 * B(1:N,1:NRHS) := Y**T * B(1:N,1:NRHS) 00435 * 00436 IF( RANK.LT.N ) THEN 00437 CALL SORMRZ( 'Left', 'Transpose', N, NRHS, RANK, N-RANK, A, 00438 $ LDA, WORK( MN+1 ), B, LDB, WORK( 2*MN+1 ), 00439 $ LWORK-2*MN, INFO ) 00440 END IF 00441 * 00442 * workspace: 2*MN+NRHS. 00443 * 00444 * B(1:N,1:NRHS) := P * B(1:N,1:NRHS) 00445 * 00446 DO 60 J = 1, NRHS 00447 DO 50 I = 1, N 00448 WORK( JPVT( I ) ) = B( I, J ) 00449 50 CONTINUE 00450 CALL SCOPY( N, WORK( 1 ), 1, B( 1, J ), 1 ) 00451 60 CONTINUE 00452 * 00453 * workspace: N. 00454 * 00455 * Undo scaling 00456 * 00457 IF( IASCL.EQ.1 ) THEN 00458 CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO ) 00459 CALL SLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA, 00460 $ INFO ) 00461 ELSE IF( IASCL.EQ.2 ) THEN 00462 CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO ) 00463 CALL SLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA, 00464 $ INFO ) 00465 END IF 00466 IF( IBSCL.EQ.1 ) THEN 00467 CALL SLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO ) 00468 ELSE IF( IBSCL.EQ.2 ) THEN 00469 CALL SLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO ) 00470 END IF 00471 * 00472 70 CONTINUE 00473 WORK( 1 ) = LWKOPT 00474 * 00475 RETURN 00476 * 00477 * End of SGELSY 00478 * 00479 END