|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DGELSX( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, 00002 $ WORK, 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, M, N, NRHS, RANK 00011 DOUBLE PRECISION RCOND 00012 * .. 00013 * .. Array Arguments .. 00014 INTEGER JPVT( * ) 00015 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * This routine is deprecated and has been replaced by routine DGELSY. 00022 * 00023 * DGELSX computes the minimum-norm solution to a real linear least 00024 * squares problem: 00025 * minimize || A * X - B || 00026 * using a complete orthogonal factorization of A. A is an M-by-N 00027 * matrix which may be rank-deficient. 00028 * 00029 * Several right hand side vectors b and solution vectors x can be 00030 * handled in a single call; they are stored as the columns of the 00031 * M-by-NRHS right hand side matrix B and the N-by-NRHS solution 00032 * matrix X. 00033 * 00034 * The routine first computes a QR factorization with column pivoting: 00035 * A * P = Q * [ R11 R12 ] 00036 * [ 0 R22 ] 00037 * with R11 defined as the largest leading submatrix whose estimated 00038 * condition number is less than 1/RCOND. The order of R11, RANK, 00039 * is the effective rank of A. 00040 * 00041 * Then, R22 is considered to be negligible, and R12 is annihilated 00042 * by orthogonal transformations from the right, arriving at the 00043 * complete orthogonal factorization: 00044 * A * P = Q * [ T11 0 ] * Z 00045 * [ 0 0 ] 00046 * The minimum-norm solution is then 00047 * X = P * Z' [ inv(T11)*Q1'*B ] 00048 * [ 0 ] 00049 * where Q1 consists of the first RANK columns of Q. 00050 * 00051 * Arguments 00052 * ========= 00053 * 00054 * M (input) INTEGER 00055 * The number of rows of the matrix A. M >= 0. 00056 * 00057 * N (input) INTEGER 00058 * The number of columns of the matrix A. N >= 0. 00059 * 00060 * NRHS (input) INTEGER 00061 * The number of right hand sides, i.e., the number of 00062 * columns of matrices B and X. NRHS >= 0. 00063 * 00064 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 00065 * On entry, the M-by-N matrix A. 00066 * On exit, A has been overwritten by details of its 00067 * complete orthogonal factorization. 00068 * 00069 * LDA (input) INTEGER 00070 * The leading dimension of the array A. LDA >= max(1,M). 00071 * 00072 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) 00073 * On entry, the M-by-NRHS right hand side matrix B. 00074 * On exit, the N-by-NRHS solution matrix X. 00075 * If m >= n and RANK = n, the residual sum-of-squares for 00076 * the solution in the i-th column is given by the sum of 00077 * squares of elements N+1:M in that column. 00078 * 00079 * LDB (input) INTEGER 00080 * The leading dimension of the array B. LDB >= max(1,M,N). 00081 * 00082 * JPVT (input/output) INTEGER array, dimension (N) 00083 * On entry, if JPVT(i) .ne. 0, the i-th column of A is an 00084 * initial column, otherwise it is a free column. Before 00085 * the QR factorization of A, all initial columns are 00086 * permuted to the leading positions; only the remaining 00087 * free columns are moved as a result of column pivoting 00088 * during the factorization. 00089 * On exit, if JPVT(i) = k, then the i-th column of A*P 00090 * was the k-th column of A. 00091 * 00092 * RCOND (input) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension 00104 * (max( min(M,N)+3*N, 2*min(M,N)+NRHS )), 00105 * 00106 * INFO (output) INTEGER 00107 * = 0: successful exit 00108 * < 0: if INFO = -i, the i-th argument had an illegal value 00109 * 00110 * ===================================================================== 00111 * 00112 * .. Parameters .. 00113 INTEGER IMAX, IMIN 00114 PARAMETER ( IMAX = 1, IMIN = 2 ) 00115 DOUBLE PRECISION ZERO, ONE, DONE, NTDONE 00116 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, DONE = ZERO, 00117 $ NTDONE = ONE ) 00118 * .. 00119 * .. Local Scalars .. 00120 INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, K, MN 00121 DOUBLE PRECISION ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX, 00122 $ SMAXPR, SMIN, SMINPR, SMLNUM, T1, T2 00123 * .. 00124 * .. External Functions .. 00125 DOUBLE PRECISION DLAMCH, DLANGE 00126 EXTERNAL DLAMCH, DLANGE 00127 * .. 00128 * .. External Subroutines .. 00129 EXTERNAL DGEQPF, DLAIC1, DLASCL, DLASET, DLATZM, DORM2R, 00130 $ DTRSM, DTZRQF, XERBLA 00131 * .. 00132 * .. Intrinsic Functions .. 00133 INTRINSIC ABS, MAX, MIN 00134 * .. 00135 * .. Executable Statements .. 00136 * 00137 MN = MIN( M, N ) 00138 ISMIN = MN + 1 00139 ISMAX = 2*MN + 1 00140 * 00141 * Test the input arguments. 00142 * 00143 INFO = 0 00144 IF( M.LT.0 ) THEN 00145 INFO = -1 00146 ELSE IF( N.LT.0 ) THEN 00147 INFO = -2 00148 ELSE IF( NRHS.LT.0 ) THEN 00149 INFO = -3 00150 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00151 INFO = -5 00152 ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN 00153 INFO = -7 00154 END IF 00155 * 00156 IF( INFO.NE.0 ) THEN 00157 CALL XERBLA( 'DGELSX', -INFO ) 00158 RETURN 00159 END IF 00160 * 00161 * Quick return if possible 00162 * 00163 IF( MIN( M, N, NRHS ).EQ.0 ) THEN 00164 RANK = 0 00165 RETURN 00166 END IF 00167 * 00168 * Get machine parameters 00169 * 00170 SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' ) 00171 BIGNUM = ONE / SMLNUM 00172 CALL DLABAD( SMLNUM, BIGNUM ) 00173 * 00174 * Scale A, B if max elements outside range [SMLNUM,BIGNUM] 00175 * 00176 ANRM = DLANGE( 'M', M, N, A, LDA, WORK ) 00177 IASCL = 0 00178 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00179 * 00180 * Scale matrix norm up to SMLNUM 00181 * 00182 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO ) 00183 IASCL = 1 00184 ELSE IF( ANRM.GT.BIGNUM ) THEN 00185 * 00186 * Scale matrix norm down to BIGNUM 00187 * 00188 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO ) 00189 IASCL = 2 00190 ELSE IF( ANRM.EQ.ZERO ) THEN 00191 * 00192 * Matrix all zero. Return zero solution. 00193 * 00194 CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB ) 00195 RANK = 0 00196 GO TO 100 00197 END IF 00198 * 00199 BNRM = DLANGE( 'M', M, NRHS, B, LDB, WORK ) 00200 IBSCL = 0 00201 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 00202 * 00203 * Scale matrix norm up to SMLNUM 00204 * 00205 CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO ) 00206 IBSCL = 1 00207 ELSE IF( BNRM.GT.BIGNUM ) THEN 00208 * 00209 * Scale matrix norm down to BIGNUM 00210 * 00211 CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO ) 00212 IBSCL = 2 00213 END IF 00214 * 00215 * Compute QR factorization with column pivoting of A: 00216 * A * P = Q * R 00217 * 00218 CALL DGEQPF( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ), INFO ) 00219 * 00220 * workspace 3*N. Details of Householder rotations stored 00221 * in WORK(1:MN). 00222 * 00223 * Determine RANK using incremental condition estimation 00224 * 00225 WORK( ISMIN ) = ONE 00226 WORK( ISMAX ) = ONE 00227 SMAX = ABS( A( 1, 1 ) ) 00228 SMIN = SMAX 00229 IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN 00230 RANK = 0 00231 CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB ) 00232 GO TO 100 00233 ELSE 00234 RANK = 1 00235 END IF 00236 * 00237 10 CONTINUE 00238 IF( RANK.LT.MN ) THEN 00239 I = RANK + 1 00240 CALL DLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ), 00241 $ A( I, I ), SMINPR, S1, C1 ) 00242 CALL DLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ), 00243 $ A( I, I ), SMAXPR, S2, C2 ) 00244 * 00245 IF( SMAXPR*RCOND.LE.SMINPR ) THEN 00246 DO 20 I = 1, RANK 00247 WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 ) 00248 WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 ) 00249 20 CONTINUE 00250 WORK( ISMIN+RANK ) = C1 00251 WORK( ISMAX+RANK ) = C2 00252 SMIN = SMINPR 00253 SMAX = SMAXPR 00254 RANK = RANK + 1 00255 GO TO 10 00256 END IF 00257 END IF 00258 * 00259 * Logically partition R = [ R11 R12 ] 00260 * [ 0 R22 ] 00261 * where R11 = R(1:RANK,1:RANK) 00262 * 00263 * [R11,R12] = [ T11, 0 ] * Y 00264 * 00265 IF( RANK.LT.N ) 00266 $ CALL DTZRQF( RANK, N, A, LDA, WORK( MN+1 ), INFO ) 00267 * 00268 * Details of Householder rotations stored in WORK(MN+1:2*MN) 00269 * 00270 * B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS) 00271 * 00272 CALL DORM2R( 'Left', 'Transpose', M, NRHS, MN, A, LDA, WORK( 1 ), 00273 $ B, LDB, WORK( 2*MN+1 ), INFO ) 00274 * 00275 * workspace NRHS 00276 * 00277 * B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS) 00278 * 00279 CALL DTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK, 00280 $ NRHS, ONE, A, LDA, B, LDB ) 00281 * 00282 DO 40 I = RANK + 1, N 00283 DO 30 J = 1, NRHS 00284 B( I, J ) = ZERO 00285 30 CONTINUE 00286 40 CONTINUE 00287 * 00288 * B(1:N,1:NRHS) := Y' * B(1:N,1:NRHS) 00289 * 00290 IF( RANK.LT.N ) THEN 00291 DO 50 I = 1, RANK 00292 CALL DLATZM( 'Left', N-RANK+1, NRHS, A( I, RANK+1 ), LDA, 00293 $ WORK( MN+I ), B( I, 1 ), B( RANK+1, 1 ), LDB, 00294 $ WORK( 2*MN+1 ) ) 00295 50 CONTINUE 00296 END IF 00297 * 00298 * workspace NRHS 00299 * 00300 * B(1:N,1:NRHS) := P * B(1:N,1:NRHS) 00301 * 00302 DO 90 J = 1, NRHS 00303 DO 60 I = 1, N 00304 WORK( 2*MN+I ) = NTDONE 00305 60 CONTINUE 00306 DO 80 I = 1, N 00307 IF( WORK( 2*MN+I ).EQ.NTDONE ) THEN 00308 IF( JPVT( I ).NE.I ) THEN 00309 K = I 00310 T1 = B( K, J ) 00311 T2 = B( JPVT( K ), J ) 00312 70 CONTINUE 00313 B( JPVT( K ), J ) = T1 00314 WORK( 2*MN+K ) = DONE 00315 T1 = T2 00316 K = JPVT( K ) 00317 T2 = B( JPVT( K ), J ) 00318 IF( JPVT( K ).NE.I ) 00319 $ GO TO 70 00320 B( I, J ) = T1 00321 WORK( 2*MN+K ) = DONE 00322 END IF 00323 END IF 00324 80 CONTINUE 00325 90 CONTINUE 00326 * 00327 * Undo scaling 00328 * 00329 IF( IASCL.EQ.1 ) THEN 00330 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO ) 00331 CALL DLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA, 00332 $ INFO ) 00333 ELSE IF( IASCL.EQ.2 ) THEN 00334 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO ) 00335 CALL DLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA, 00336 $ INFO ) 00337 END IF 00338 IF( IBSCL.EQ.1 ) THEN 00339 CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO ) 00340 ELSE IF( IBSCL.EQ.2 ) THEN 00341 CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO ) 00342 END IF 00343 * 00344 100 CONTINUE 00345 * 00346 RETURN 00347 * 00348 * End of DGELSX 00349 * 00350 END
1.7.3 
