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