![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> SGELSX 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 SGELSX + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgelsx.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgelsx.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgelsx.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SGELSX( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, 00022 * WORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * INTEGER INFO, LDA, LDB, 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 *> This routine is deprecated and has been replaced by routine SGELSY. 00040 *> 00041 *> SGELSX computes the minimum-norm solution to a real linear least 00042 *> squares problem: 00043 *> minimize || A * X - B || 00044 *> using a complete orthogonal factorization of A. A is an M-by-N 00045 *> matrix which may be rank-deficient. 00046 *> 00047 *> Several right hand side vectors b and solution vectors x can be 00048 *> handled in a single call; they are stored as the columns of the 00049 *> M-by-NRHS right hand side matrix B and the N-by-NRHS solution 00050 *> matrix X. 00051 *> 00052 *> The routine first computes a QR factorization with column pivoting: 00053 *> A * P = Q * [ R11 R12 ] 00054 *> [ 0 R22 ] 00055 *> with R11 defined as the largest leading submatrix whose estimated 00056 *> condition number is less than 1/RCOND. The order of R11, RANK, 00057 *> is the effective rank of A. 00058 *> 00059 *> Then, R22 is considered to be negligible, and R12 is annihilated 00060 *> by orthogonal transformations from the right, arriving at the 00061 *> complete orthogonal factorization: 00062 *> A * P = Q * [ T11 0 ] * Z 00063 *> [ 0 0 ] 00064 *> The minimum-norm solution is then 00065 *> X = P * Z**T [ inv(T11)*Q1**T*B ] 00066 *> [ 0 ] 00067 *> where Q1 consists of the first RANK columns of Q. 00068 *> \endverbatim 00069 * 00070 * Arguments: 00071 * ========== 00072 * 00073 *> \param[in] M 00074 *> \verbatim 00075 *> M is INTEGER 00076 *> The number of rows of the matrix A. M >= 0. 00077 *> \endverbatim 00078 *> 00079 *> \param[in] N 00080 *> \verbatim 00081 *> N is INTEGER 00082 *> The number of columns of the matrix A. N >= 0. 00083 *> \endverbatim 00084 *> 00085 *> \param[in] NRHS 00086 *> \verbatim 00087 *> NRHS is INTEGER 00088 *> The number of right hand sides, i.e., the number of 00089 *> columns of matrices B and X. NRHS >= 0. 00090 *> \endverbatim 00091 *> 00092 *> \param[in,out] A 00093 *> \verbatim 00094 *> A is REAL array, dimension (LDA,N) 00095 *> On entry, the M-by-N matrix A. 00096 *> On exit, A has been overwritten by details of its 00097 *> complete orthogonal factorization. 00098 *> \endverbatim 00099 *> 00100 *> \param[in] LDA 00101 *> \verbatim 00102 *> LDA is INTEGER 00103 *> The leading dimension of the array A. LDA >= max(1,M). 00104 *> \endverbatim 00105 *> 00106 *> \param[in,out] B 00107 *> \verbatim 00108 *> B is REAL array, dimension (LDB,NRHS) 00109 *> On entry, the M-by-NRHS right hand side matrix B. 00110 *> On exit, the N-by-NRHS solution matrix X. 00111 *> If m >= n and RANK = n, the residual sum-of-squares for 00112 *> the solution in the i-th column is given by the sum of 00113 *> squares of elements N+1:M in that column. 00114 *> \endverbatim 00115 *> 00116 *> \param[in] LDB 00117 *> \verbatim 00118 *> LDB is INTEGER 00119 *> The leading dimension of the array B. LDB >= max(1,M,N). 00120 *> \endverbatim 00121 *> 00122 *> \param[in,out] JPVT 00123 *> \verbatim 00124 *> JPVT is INTEGER array, dimension (N) 00125 *> On entry, if JPVT(i) .ne. 0, the i-th column of A is an 00126 *> initial column, otherwise it is a free column. Before 00127 *> the QR factorization of A, all initial columns are 00128 *> permuted to the leading positions; only the remaining 00129 *> free columns are moved as a result of column pivoting 00130 *> during the factorization. 00131 *> On exit, if JPVT(i) = k, then the i-th column of A*P 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 00155 *> (max( min(M,N)+3*N, 2*min(M,N)+NRHS )), 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 *> \endverbatim 00164 * 00165 * Authors: 00166 * ======== 00167 * 00168 *> \author Univ. of Tennessee 00169 *> \author Univ. of California Berkeley 00170 *> \author Univ. of Colorado Denver 00171 *> \author NAG Ltd. 00172 * 00173 *> \date November 2011 00174 * 00175 *> \ingroup realGEsolve 00176 * 00177 * ===================================================================== 00178 SUBROUTINE SGELSX( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, 00179 $ WORK, INFO ) 00180 * 00181 * -- LAPACK driver routine (version 3.4.0) -- 00182 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00183 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00184 * November 2011 00185 * 00186 * .. Scalar Arguments .. 00187 INTEGER INFO, LDA, LDB, M, N, NRHS, RANK 00188 REAL RCOND 00189 * .. 00190 * .. Array Arguments .. 00191 INTEGER JPVT( * ) 00192 REAL A( LDA, * ), B( LDB, * ), WORK( * ) 00193 * .. 00194 * 00195 * ===================================================================== 00196 * 00197 * .. Parameters .. 00198 INTEGER IMAX, IMIN 00199 PARAMETER ( IMAX = 1, IMIN = 2 ) 00200 REAL ZERO, ONE, DONE, NTDONE 00201 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, DONE = ZERO, 00202 $ NTDONE = ONE ) 00203 * .. 00204 * .. Local Scalars .. 00205 INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, K, MN 00206 REAL ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX, 00207 $ SMAXPR, SMIN, SMINPR, SMLNUM, T1, T2 00208 * .. 00209 * .. External Functions .. 00210 REAL SLAMCH, SLANGE 00211 EXTERNAL SLAMCH, SLANGE 00212 * .. 00213 * .. External Subroutines .. 00214 EXTERNAL SGEQPF, SLABAD, SLAIC1, SLASCL, SLASET, SLATZM, 00215 $ SORM2R, STRSM, STZRQF, XERBLA 00216 * .. 00217 * .. Intrinsic Functions .. 00218 INTRINSIC ABS, MAX, MIN 00219 * .. 00220 * .. Executable Statements .. 00221 * 00222 MN = MIN( M, N ) 00223 ISMIN = MN + 1 00224 ISMAX = 2*MN + 1 00225 * 00226 * Test the input arguments. 00227 * 00228 INFO = 0 00229 IF( M.LT.0 ) THEN 00230 INFO = -1 00231 ELSE IF( N.LT.0 ) THEN 00232 INFO = -2 00233 ELSE IF( NRHS.LT.0 ) THEN 00234 INFO = -3 00235 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00236 INFO = -5 00237 ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN 00238 INFO = -7 00239 END IF 00240 * 00241 IF( INFO.NE.0 ) THEN 00242 CALL XERBLA( 'SGELSX', -INFO ) 00243 RETURN 00244 END IF 00245 * 00246 * Quick return if possible 00247 * 00248 IF( MIN( M, N, NRHS ).EQ.0 ) THEN 00249 RANK = 0 00250 RETURN 00251 END IF 00252 * 00253 * Get machine parameters 00254 * 00255 SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' ) 00256 BIGNUM = ONE / SMLNUM 00257 CALL SLABAD( SMLNUM, BIGNUM ) 00258 * 00259 * Scale A, B if max elements outside range [SMLNUM,BIGNUM] 00260 * 00261 ANRM = SLANGE( 'M', M, N, A, LDA, WORK ) 00262 IASCL = 0 00263 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00264 * 00265 * Scale matrix norm up to SMLNUM 00266 * 00267 CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO ) 00268 IASCL = 1 00269 ELSE IF( ANRM.GT.BIGNUM ) THEN 00270 * 00271 * Scale matrix norm down to BIGNUM 00272 * 00273 CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO ) 00274 IASCL = 2 00275 ELSE IF( ANRM.EQ.ZERO ) THEN 00276 * 00277 * Matrix all zero. Return zero solution. 00278 * 00279 CALL SLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB ) 00280 RANK = 0 00281 GO TO 100 00282 END IF 00283 * 00284 BNRM = SLANGE( 'M', M, NRHS, B, LDB, WORK ) 00285 IBSCL = 0 00286 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 00287 * 00288 * Scale matrix norm up to SMLNUM 00289 * 00290 CALL SLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO ) 00291 IBSCL = 1 00292 ELSE IF( BNRM.GT.BIGNUM ) THEN 00293 * 00294 * Scale matrix norm down to BIGNUM 00295 * 00296 CALL SLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO ) 00297 IBSCL = 2 00298 END IF 00299 * 00300 * Compute QR factorization with column pivoting of A: 00301 * A * P = Q * R 00302 * 00303 CALL SGEQPF( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ), INFO ) 00304 * 00305 * workspace 3*N. Details of Householder rotations stored 00306 * in WORK(1:MN). 00307 * 00308 * Determine RANK using incremental condition estimation 00309 * 00310 WORK( ISMIN ) = ONE 00311 WORK( ISMAX ) = ONE 00312 SMAX = ABS( A( 1, 1 ) ) 00313 SMIN = SMAX 00314 IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN 00315 RANK = 0 00316 CALL SLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB ) 00317 GO TO 100 00318 ELSE 00319 RANK = 1 00320 END IF 00321 * 00322 10 CONTINUE 00323 IF( RANK.LT.MN ) THEN 00324 I = RANK + 1 00325 CALL SLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ), 00326 $ A( I, I ), SMINPR, S1, C1 ) 00327 CALL SLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ), 00328 $ A( I, I ), SMAXPR, S2, C2 ) 00329 * 00330 IF( SMAXPR*RCOND.LE.SMINPR ) THEN 00331 DO 20 I = 1, RANK 00332 WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 ) 00333 WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 ) 00334 20 CONTINUE 00335 WORK( ISMIN+RANK ) = C1 00336 WORK( ISMAX+RANK ) = C2 00337 SMIN = SMINPR 00338 SMAX = SMAXPR 00339 RANK = RANK + 1 00340 GO TO 10 00341 END IF 00342 END IF 00343 * 00344 * Logically partition R = [ R11 R12 ] 00345 * [ 0 R22 ] 00346 * where R11 = R(1:RANK,1:RANK) 00347 * 00348 * [R11,R12] = [ T11, 0 ] * Y 00349 * 00350 IF( RANK.LT.N ) 00351 $ CALL STZRQF( RANK, N, A, LDA, WORK( MN+1 ), INFO ) 00352 * 00353 * Details of Householder rotations stored in WORK(MN+1:2*MN) 00354 * 00355 * B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS) 00356 * 00357 CALL SORM2R( 'Left', 'Transpose', M, NRHS, MN, A, LDA, WORK( 1 ), 00358 $ B, LDB, WORK( 2*MN+1 ), INFO ) 00359 * 00360 * workspace NRHS 00361 * 00362 * B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS) 00363 * 00364 CALL STRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK, 00365 $ NRHS, ONE, A, LDA, B, LDB ) 00366 * 00367 DO 40 I = RANK + 1, N 00368 DO 30 J = 1, NRHS 00369 B( I, J ) = ZERO 00370 30 CONTINUE 00371 40 CONTINUE 00372 * 00373 * B(1:N,1:NRHS) := Y**T * B(1:N,1:NRHS) 00374 * 00375 IF( RANK.LT.N ) THEN 00376 DO 50 I = 1, RANK 00377 CALL SLATZM( 'Left', N-RANK+1, NRHS, A( I, RANK+1 ), LDA, 00378 $ WORK( MN+I ), B( I, 1 ), B( RANK+1, 1 ), LDB, 00379 $ WORK( 2*MN+1 ) ) 00380 50 CONTINUE 00381 END IF 00382 * 00383 * workspace NRHS 00384 * 00385 * B(1:N,1:NRHS) := P * B(1:N,1:NRHS) 00386 * 00387 DO 90 J = 1, NRHS 00388 DO 60 I = 1, N 00389 WORK( 2*MN+I ) = NTDONE 00390 60 CONTINUE 00391 DO 80 I = 1, N 00392 IF( WORK( 2*MN+I ).EQ.NTDONE ) THEN 00393 IF( JPVT( I ).NE.I ) THEN 00394 K = I 00395 T1 = B( K, J ) 00396 T2 = B( JPVT( K ), J ) 00397 70 CONTINUE 00398 B( JPVT( K ), J ) = T1 00399 WORK( 2*MN+K ) = DONE 00400 T1 = T2 00401 K = JPVT( K ) 00402 T2 = B( JPVT( K ), J ) 00403 IF( JPVT( K ).NE.I ) 00404 $ GO TO 70 00405 B( I, J ) = T1 00406 WORK( 2*MN+K ) = DONE 00407 END IF 00408 END IF 00409 80 CONTINUE 00410 90 CONTINUE 00411 * 00412 * Undo scaling 00413 * 00414 IF( IASCL.EQ.1 ) THEN 00415 CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO ) 00416 CALL SLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA, 00417 $ INFO ) 00418 ELSE IF( IASCL.EQ.2 ) THEN 00419 CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO ) 00420 CALL SLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA, 00421 $ INFO ) 00422 END IF 00423 IF( IBSCL.EQ.1 ) THEN 00424 CALL SLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO ) 00425 ELSE IF( IBSCL.EQ.2 ) THEN 00426 CALL SLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO ) 00427 END IF 00428 * 00429 100 CONTINUE 00430 * 00431 RETURN 00432 * 00433 * End of SGELSX 00434 * 00435 END