|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DGESVX( FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, 00002 $ EQUED, R, C, B, LDB, X, LDX, RCOND, FERR, BERR, 00003 $ WORK, IWORK, INFO ) 00004 * 00005 * -- LAPACK driver routine (version 3.2) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * November 2006 00009 * 00010 * .. Scalar Arguments .. 00011 CHARACTER EQUED, FACT, TRANS 00012 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS 00013 DOUBLE PRECISION RCOND 00014 * .. 00015 * .. Array Arguments .. 00016 INTEGER IPIV( * ), IWORK( * ) 00017 DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 00018 $ BERR( * ), C( * ), FERR( * ), R( * ), 00019 $ WORK( * ), X( LDX, * ) 00020 * .. 00021 * 00022 * Purpose 00023 * ======= 00024 * 00025 * DGESVX uses the LU factorization to compute the solution to a real 00026 * system of linear equations 00027 * A * X = B, 00028 * where A is an N-by-N matrix and X and B are N-by-NRHS matrices. 00029 * 00030 * Error bounds on the solution and a condition estimate are also 00031 * provided. 00032 * 00033 * Description 00034 * =========== 00035 * 00036 * The following steps are performed: 00037 * 00038 * 1. If FACT = 'E', real scaling factors are computed to equilibrate 00039 * the system: 00040 * TRANS = 'N': diag(R)*A*diag(C) *inv(diag(C))*X = diag(R)*B 00041 * TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B 00042 * TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B 00043 * Whether or not the system will be equilibrated depends on the 00044 * scaling of the matrix A, but if equilibration is used, A is 00045 * overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N') 00046 * or diag(C)*B (if TRANS = 'T' or 'C'). 00047 * 00048 * 2. If FACT = 'N' or 'E', the LU decomposition is used to factor the 00049 * matrix A (after equilibration if FACT = 'E') as 00050 * A = P * L * U, 00051 * where P is a permutation matrix, L is a unit lower triangular 00052 * matrix, and U is upper triangular. 00053 * 00054 * 3. If some U(i,i)=0, so that U is exactly singular, then the routine 00055 * returns with INFO = i. Otherwise, the factored form of A is used 00056 * to estimate the condition number of the matrix A. If the 00057 * reciprocal of the condition number is less than machine precision, 00058 * INFO = N+1 is returned as a warning, but the routine still goes on 00059 * to solve for X and compute error bounds as described below. 00060 * 00061 * 4. The system of equations is solved for X using the factored form 00062 * of A. 00063 * 00064 * 5. Iterative refinement is applied to improve the computed solution 00065 * matrix and calculate error bounds and backward error estimates 00066 * for it. 00067 * 00068 * 6. If equilibration was used, the matrix X is premultiplied by 00069 * diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so 00070 * that it solves the original system before equilibration. 00071 * 00072 * Arguments 00073 * ========= 00074 * 00075 * FACT (input) CHARACTER*1 00076 * Specifies whether or not the factored form of the matrix A is 00077 * supplied on entry, and if not, whether the matrix A should be 00078 * equilibrated before it is factored. 00079 * = 'F': On entry, AF and IPIV contain the factored form of A. 00080 * If EQUED is not 'N', the matrix A has been 00081 * equilibrated with scaling factors given by R and C. 00082 * A, AF, and IPIV are not modified. 00083 * = 'N': The matrix A will be copied to AF and factored. 00084 * = 'E': The matrix A will be equilibrated if necessary, then 00085 * copied to AF and factored. 00086 * 00087 * TRANS (input) CHARACTER*1 00088 * Specifies the form of the system of equations: 00089 * = 'N': A * X = B (No transpose) 00090 * = 'T': A**T * X = B (Transpose) 00091 * = 'C': A**H * X = B (Transpose) 00092 * 00093 * N (input) INTEGER 00094 * The number of linear equations, i.e., the order of the 00095 * matrix A. N >= 0. 00096 * 00097 * NRHS (input) INTEGER 00098 * The number of right hand sides, i.e., the number of columns 00099 * of the matrices B and X. NRHS >= 0. 00100 * 00101 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 00102 * On entry, the N-by-N matrix A. If FACT = 'F' and EQUED is 00103 * not 'N', then A must have been equilibrated by the scaling 00104 * factors in R and/or C. A is not modified if FACT = 'F' or 00105 * 'N', or if FACT = 'E' and EQUED = 'N' on exit. 00106 * 00107 * On exit, if EQUED .ne. 'N', A is scaled as follows: 00108 * EQUED = 'R': A := diag(R) * A 00109 * EQUED = 'C': A := A * diag(C) 00110 * EQUED = 'B': A := diag(R) * A * diag(C). 00111 * 00112 * LDA (input) INTEGER 00113 * The leading dimension of the array A. LDA >= max(1,N). 00114 * 00115 * AF (input or output) DOUBLE PRECISION array, dimension (LDAF,N) 00116 * If FACT = 'F', then AF is an input argument and on entry 00117 * contains the factors L and U from the factorization 00118 * A = P*L*U as computed by DGETRF. If EQUED .ne. 'N', then 00119 * AF is the factored form of the equilibrated matrix A. 00120 * 00121 * If FACT = 'N', then AF is an output argument and on exit 00122 * returns the factors L and U from the factorization A = P*L*U 00123 * of the original matrix A. 00124 * 00125 * If FACT = 'E', then AF is an output argument and on exit 00126 * returns the factors L and U from the factorization A = P*L*U 00127 * of the equilibrated matrix A (see the description of A for 00128 * the form of the equilibrated matrix). 00129 * 00130 * LDAF (input) INTEGER 00131 * The leading dimension of the array AF. LDAF >= max(1,N). 00132 * 00133 * IPIV (input or output) INTEGER array, dimension (N) 00134 * If FACT = 'F', then IPIV is an input argument and on entry 00135 * contains the pivot indices from the factorization A = P*L*U 00136 * as computed by DGETRF; row i of the matrix was interchanged 00137 * with row IPIV(i). 00138 * 00139 * If FACT = 'N', then IPIV is an output argument and on exit 00140 * contains the pivot indices from the factorization A = P*L*U 00141 * of the original matrix A. 00142 * 00143 * If FACT = 'E', then IPIV is an output argument and on exit 00144 * contains the pivot indices from the factorization A = P*L*U 00145 * of the equilibrated matrix A. 00146 * 00147 * EQUED (input or output) CHARACTER*1 00148 * Specifies the form of equilibration that was done. 00149 * = 'N': No equilibration (always true if FACT = 'N'). 00150 * = 'R': Row equilibration, i.e., A has been premultiplied by 00151 * diag(R). 00152 * = 'C': Column equilibration, i.e., A has been postmultiplied 00153 * by diag(C). 00154 * = 'B': Both row and column equilibration, i.e., A has been 00155 * replaced by diag(R) * A * diag(C). 00156 * EQUED is an input argument if FACT = 'F'; otherwise, it is an 00157 * output argument. 00158 * 00159 * R (input or output) DOUBLE PRECISION array, dimension (N) 00160 * The row scale factors for A. If EQUED = 'R' or 'B', A is 00161 * multiplied on the left by diag(R); if EQUED = 'N' or 'C', R 00162 * is not accessed. R is an input argument if FACT = 'F'; 00163 * otherwise, R is an output argument. If FACT = 'F' and 00164 * EQUED = 'R' or 'B', each element of R must be positive. 00165 * 00166 * C (input or output) DOUBLE PRECISION array, dimension (N) 00167 * The column scale factors for A. If EQUED = 'C' or 'B', A is 00168 * multiplied on the right by diag(C); if EQUED = 'N' or 'R', C 00169 * is not accessed. C is an input argument if FACT = 'F'; 00170 * otherwise, C is an output argument. If FACT = 'F' and 00171 * EQUED = 'C' or 'B', each element of C must be positive. 00172 * 00173 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) 00174 * On entry, the N-by-NRHS right hand side matrix B. 00175 * On exit, 00176 * if EQUED = 'N', B is not modified; 00177 * if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by 00178 * diag(R)*B; 00179 * if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is 00180 * overwritten by diag(C)*B. 00181 * 00182 * LDB (input) INTEGER 00183 * The leading dimension of the array B. LDB >= max(1,N). 00184 * 00185 * X (output) DOUBLE PRECISION array, dimension (LDX,NRHS) 00186 * If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X 00187 * to the original system of equations. Note that A and B are 00188 * modified on exit if EQUED .ne. 'N', and the solution to the 00189 * equilibrated system is inv(diag(C))*X if TRANS = 'N' and 00190 * EQUED = 'C' or 'B', or inv(diag(R))*X if TRANS = 'T' or 'C' 00191 * and EQUED = 'R' or 'B'. 00192 * 00193 * LDX (input) INTEGER 00194 * The leading dimension of the array X. LDX >= max(1,N). 00195 * 00196 * RCOND (output) DOUBLE PRECISION 00197 * The estimate of the reciprocal condition number of the matrix 00198 * A after equilibration (if done). If RCOND is less than the 00199 * machine precision (in particular, if RCOND = 0), the matrix 00200 * is singular to working precision. This condition is 00201 * indicated by a return code of INFO > 0. 00202 * 00203 * FERR (output) DOUBLE PRECISION array, dimension (NRHS) 00204 * The estimated forward error bound for each solution vector 00205 * X(j) (the j-th column of the solution matrix X). 00206 * If XTRUE is the true solution corresponding to X(j), FERR(j) 00207 * is an estimated upper bound for the magnitude of the largest 00208 * element in (X(j) - XTRUE) divided by the magnitude of the 00209 * largest element in X(j). The estimate is as reliable as 00210 * the estimate for RCOND, and is almost always a slight 00211 * overestimate of the true error. 00212 * 00213 * BERR (output) DOUBLE PRECISION array, dimension (NRHS) 00214 * The componentwise relative backward error of each solution 00215 * vector X(j) (i.e., the smallest relative change in 00216 * any element of A or B that makes X(j) an exact solution). 00217 * 00218 * WORK (workspace/output) DOUBLE PRECISION array, dimension (4*N) 00219 * On exit, WORK(1) contains the reciprocal pivot growth 00220 * factor norm(A)/norm(U). The "max absolute element" norm is 00221 * used. If WORK(1) is much less than 1, then the stability 00222 * of the LU factorization of the (equilibrated) matrix A 00223 * could be poor. This also means that the solution X, condition 00224 * estimator RCOND, and forward error bound FERR could be 00225 * unreliable. If factorization fails with 0<INFO<=N, then 00226 * WORK(1) contains the reciprocal pivot growth factor for the 00227 * leading INFO columns of A. 00228 * 00229 * IWORK (workspace) INTEGER array, dimension (N) 00230 * 00231 * INFO (output) INTEGER 00232 * = 0: successful exit 00233 * < 0: if INFO = -i, the i-th argument had an illegal value 00234 * > 0: if INFO = i, and i is 00235 * <= N: U(i,i) is exactly zero. The factorization has 00236 * been completed, but the factor U is exactly 00237 * singular, so the solution and error bounds 00238 * could not be computed. RCOND = 0 is returned. 00239 * = N+1: U is nonsingular, but RCOND is less than machine 00240 * precision, meaning that the matrix is singular 00241 * to working precision. Nevertheless, the 00242 * solution and error bounds are computed because 00243 * there are a number of situations where the 00244 * computed solution can be more accurate than the 00245 * value of RCOND would suggest. 00246 * 00247 * ===================================================================== 00248 * 00249 * .. Parameters .. 00250 DOUBLE PRECISION ZERO, ONE 00251 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00252 * .. 00253 * .. Local Scalars .. 00254 LOGICAL COLEQU, EQUIL, NOFACT, NOTRAN, ROWEQU 00255 CHARACTER NORM 00256 INTEGER I, INFEQU, J 00257 DOUBLE PRECISION AMAX, ANORM, BIGNUM, COLCND, RCMAX, RCMIN, 00258 $ ROWCND, RPVGRW, SMLNUM 00259 * .. 00260 * .. External Functions .. 00261 LOGICAL LSAME 00262 DOUBLE PRECISION DLAMCH, DLANGE, DLANTR 00263 EXTERNAL LSAME, DLAMCH, DLANGE, DLANTR 00264 * .. 00265 * .. External Subroutines .. 00266 EXTERNAL DGECON, DGEEQU, DGERFS, DGETRF, DGETRS, DLACPY, 00267 $ DLAQGE, XERBLA 00268 * .. 00269 * .. Intrinsic Functions .. 00270 INTRINSIC MAX, MIN 00271 * .. 00272 * .. Executable Statements .. 00273 * 00274 INFO = 0 00275 NOFACT = LSAME( FACT, 'N' ) 00276 EQUIL = LSAME( FACT, 'E' ) 00277 NOTRAN = LSAME( TRANS, 'N' ) 00278 IF( NOFACT .OR. EQUIL ) THEN 00279 EQUED = 'N' 00280 ROWEQU = .FALSE. 00281 COLEQU = .FALSE. 00282 ELSE 00283 ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' ) 00284 COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' ) 00285 SMLNUM = DLAMCH( 'Safe minimum' ) 00286 BIGNUM = ONE / SMLNUM 00287 END IF 00288 * 00289 * Test the input parameters. 00290 * 00291 IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) ) 00292 $ THEN 00293 INFO = -1 00294 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. 00295 $ LSAME( TRANS, 'C' ) ) THEN 00296 INFO = -2 00297 ELSE IF( N.LT.0 ) THEN 00298 INFO = -3 00299 ELSE IF( NRHS.LT.0 ) THEN 00300 INFO = -4 00301 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00302 INFO = -6 00303 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN 00304 INFO = -8 00305 ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT. 00306 $ ( ROWEQU .OR. COLEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN 00307 INFO = -10 00308 ELSE 00309 IF( ROWEQU ) THEN 00310 RCMIN = BIGNUM 00311 RCMAX = ZERO 00312 DO 10 J = 1, N 00313 RCMIN = MIN( RCMIN, R( J ) ) 00314 RCMAX = MAX( RCMAX, R( J ) ) 00315 10 CONTINUE 00316 IF( RCMIN.LE.ZERO ) THEN 00317 INFO = -11 00318 ELSE IF( N.GT.0 ) THEN 00319 ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM ) 00320 ELSE 00321 ROWCND = ONE 00322 END IF 00323 END IF 00324 IF( COLEQU .AND. INFO.EQ.0 ) THEN 00325 RCMIN = BIGNUM 00326 RCMAX = ZERO 00327 DO 20 J = 1, N 00328 RCMIN = MIN( RCMIN, C( J ) ) 00329 RCMAX = MAX( RCMAX, C( J ) ) 00330 20 CONTINUE 00331 IF( RCMIN.LE.ZERO ) THEN 00332 INFO = -12 00333 ELSE IF( N.GT.0 ) THEN 00334 COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM ) 00335 ELSE 00336 COLCND = ONE 00337 END IF 00338 END IF 00339 IF( INFO.EQ.0 ) THEN 00340 IF( LDB.LT.MAX( 1, N ) ) THEN 00341 INFO = -14 00342 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 00343 INFO = -16 00344 END IF 00345 END IF 00346 END IF 00347 * 00348 IF( INFO.NE.0 ) THEN 00349 CALL XERBLA( 'DGESVX', -INFO ) 00350 RETURN 00351 END IF 00352 * 00353 IF( EQUIL ) THEN 00354 * 00355 * Compute row and column scalings to equilibrate the matrix A. 00356 * 00357 CALL DGEEQU( N, N, A, LDA, R, C, ROWCND, COLCND, AMAX, INFEQU ) 00358 IF( INFEQU.EQ.0 ) THEN 00359 * 00360 * Equilibrate the matrix. 00361 * 00362 CALL DLAQGE( N, N, A, LDA, R, C, ROWCND, COLCND, AMAX, 00363 $ EQUED ) 00364 ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' ) 00365 COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' ) 00366 END IF 00367 END IF 00368 * 00369 * Scale the right hand side. 00370 * 00371 IF( NOTRAN ) THEN 00372 IF( ROWEQU ) THEN 00373 DO 40 J = 1, NRHS 00374 DO 30 I = 1, N 00375 B( I, J ) = R( I )*B( I, J ) 00376 30 CONTINUE 00377 40 CONTINUE 00378 END IF 00379 ELSE IF( COLEQU ) THEN 00380 DO 60 J = 1, NRHS 00381 DO 50 I = 1, N 00382 B( I, J ) = C( I )*B( I, J ) 00383 50 CONTINUE 00384 60 CONTINUE 00385 END IF 00386 * 00387 IF( NOFACT .OR. EQUIL ) THEN 00388 * 00389 * Compute the LU factorization of A. 00390 * 00391 CALL DLACPY( 'Full', N, N, A, LDA, AF, LDAF ) 00392 CALL DGETRF( N, N, AF, LDAF, IPIV, INFO ) 00393 * 00394 * Return if INFO is non-zero. 00395 * 00396 IF( INFO.GT.0 ) THEN 00397 * 00398 * Compute the reciprocal pivot growth factor of the 00399 * leading rank-deficient INFO columns of A. 00400 * 00401 RPVGRW = DLANTR( 'M', 'U', 'N', INFO, INFO, AF, LDAF, 00402 $ WORK ) 00403 IF( RPVGRW.EQ.ZERO ) THEN 00404 RPVGRW = ONE 00405 ELSE 00406 RPVGRW = DLANGE( 'M', N, INFO, A, LDA, WORK ) / RPVGRW 00407 END IF 00408 WORK( 1 ) = RPVGRW 00409 RCOND = ZERO 00410 RETURN 00411 END IF 00412 END IF 00413 * 00414 * Compute the norm of the matrix A and the 00415 * reciprocal pivot growth factor RPVGRW. 00416 * 00417 IF( NOTRAN ) THEN 00418 NORM = '1' 00419 ELSE 00420 NORM = 'I' 00421 END IF 00422 ANORM = DLANGE( NORM, N, N, A, LDA, WORK ) 00423 RPVGRW = DLANTR( 'M', 'U', 'N', N, N, AF, LDAF, WORK ) 00424 IF( RPVGRW.EQ.ZERO ) THEN 00425 RPVGRW = ONE 00426 ELSE 00427 RPVGRW = DLANGE( 'M', N, N, A, LDA, WORK ) / RPVGRW 00428 END IF 00429 * 00430 * Compute the reciprocal of the condition number of A. 00431 * 00432 CALL DGECON( NORM, N, AF, LDAF, ANORM, RCOND, WORK, IWORK, INFO ) 00433 * 00434 * Compute the solution matrix X. 00435 * 00436 CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 00437 CALL DGETRS( TRANS, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO ) 00438 * 00439 * Use iterative refinement to improve the computed solution and 00440 * compute error bounds and backward error estimates for it. 00441 * 00442 CALL DGERFS( TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, 00443 $ LDX, FERR, BERR, WORK, IWORK, INFO ) 00444 * 00445 * Transform the solution matrix X to a solution of the original 00446 * system. 00447 * 00448 IF( NOTRAN ) THEN 00449 IF( COLEQU ) THEN 00450 DO 80 J = 1, NRHS 00451 DO 70 I = 1, N 00452 X( I, J ) = C( I )*X( I, J ) 00453 70 CONTINUE 00454 80 CONTINUE 00455 DO 90 J = 1, NRHS 00456 FERR( J ) = FERR( J ) / COLCND 00457 90 CONTINUE 00458 END IF 00459 ELSE IF( ROWEQU ) THEN 00460 DO 110 J = 1, NRHS 00461 DO 100 I = 1, N 00462 X( I, J ) = R( I )*X( I, J ) 00463 100 CONTINUE 00464 110 CONTINUE 00465 DO 120 J = 1, NRHS 00466 FERR( J ) = FERR( J ) / ROWCND 00467 120 CONTINUE 00468 END IF 00469 * 00470 WORK( 1 ) = RPVGRW 00471 * 00472 * Set INFO = N+1 if the matrix is singular to working precision. 00473 * 00474 IF( RCOND.LT.DLAMCH( 'Epsilon' ) ) 00475 $ INFO = N + 1 00476 RETURN 00477 * 00478 * End of DGESVX 00479 * 00480 END
1.7.3 
