001:SUBROUTINEDGESVX( FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, 002: $ EQUED, R, C, B, LDB, X, LDX, RCOND, FERR, BERR, 003: $ WORK, IWORK, INFO ) 004:*005:* -- LAPACK driver routine (version 3.2) --006:* -- LAPACK is a software package provided by Univ. of Tennessee, --007:* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--008:* November 2006009:*010:* .. Scalar Arguments ..011: CHARACTER EQUED, FACT, TRANS 012: INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS 013: DOUBLE PRECISION RCOND 014:* ..015:* .. Array Arguments ..016: INTEGERIPIV( * ),IWORK( * ) 017: DOUBLE PRECISIONA( LDA, * ),AF( LDAF, * ),B( LDB, * ), 018: $BERR( * ),C( * ),FERR( * ),R( * ), 019: $WORK( * ),X( LDX, * ) 020:* ..021:*022:* Purpose023:* =======024:*025:* DGESVX uses the LU factorization to compute the solution to a real026:* system of linear equations027:* A * X = B,028:* where A is an N-by-N matrix and X and B are N-by-NRHS matrices.029:*030:* Error bounds on the solution and a condition estimate are also031:* provided.032:*033:* Description034:* ===========035:*036:* The following steps are performed:037:*038:* 1. If FACT = 'E', real scaling factors are computed to equilibrate039:* the system:040:* TRANS = 'N': diag(R)*A*diag(C) *inv(diag(C))*X = diag(R)*B041:* TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B042:* TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B043:* Whether or not the system will be equilibrated depends on the044:* scaling of the matrix A, but if equilibration is used, A is045:* overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N')046:* or diag(C)*B (if TRANS = 'T' or 'C').047:*048:* 2. If FACT = 'N' or 'E', the LU decomposition is used to factor the049:* matrix A (after equilibration if FACT = 'E') as050:* A = P * L * U,051:* where P is a permutation matrix, L is a unit lower triangular052:* matrix, and U is upper triangular.053:*054:* 3. If some U(i,i)=0, so that U is exactly singular, then the routine055:* returns with INFO = i. Otherwise, the factored form of A is used056:* to estimate the condition number of the matrix A. If the057:* reciprocal of the condition number is less than machine precision,058:* INFO = N+1 is returned as a warning, but the routine still goes on059:* to solve for X and compute error bounds as described below.060:*061:* 4. The system of equations is solved for X using the factored form062:* of A.063:*064:* 5. Iterative refinement is applied to improve the computed solution065:* matrix and calculate error bounds and backward error estimates066:* for it.067:*068:* 6. If equilibration was used, the matrix X is premultiplied by069:* diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so070:* that it solves the original system before equilibration.071:*072:* Arguments073:* =========074:*075:* FACT (input) CHARACTER*1076:* Specifies whether or not the factored form of the matrix A is077:* supplied on entry, and if not, whether the matrix A should be078:* equilibrated before it is factored.079:* = 'F': On entry, AF and IPIV contain the factored form of A.080:* If EQUED is not 'N', the matrix A has been081:* equilibrated with scaling factors given by R and C.082:* A, AF, and IPIV are not modified.083:* = 'N': The matrix A will be copied to AF and factored.084:* = 'E': The matrix A will be equilibrated if necessary, then085:* copied to AF and factored.086:*087:* TRANS (input) CHARACTER*1088:* Specifies the form of the system of equations:089:* = 'N': A * X = B (No transpose)090:* = 'T': A**T * X = B (Transpose)091:* = 'C': A**H * X = B (Transpose)092:*093:* N (input) INTEGER094:* The number of linear equations, i.e., the order of the095:* matrix A. N >= 0.096:*097:* NRHS (input) INTEGER098:* The number of right hand sides, i.e., the number of columns099:* of the matrices B and X. NRHS >= 0.100:*101:* A (input/output) DOUBLE PRECISION array, dimension (LDA,N)102:* On entry, the N-by-N matrix A. If FACT = 'F' and EQUED is103:* not 'N', then A must have been equilibrated by the scaling104:* factors in R and/or C. A is not modified if FACT = 'F' or105:* 'N', or if FACT = 'E' and EQUED = 'N' on exit.106:*107:* On exit, if EQUED .ne. 'N', A is scaled as follows:108:* EQUED = 'R': A := diag(R) * A109:* EQUED = 'C': A := A * diag(C)110:* EQUED = 'B': A := diag(R) * A * diag(C).111:*112:* LDA (input) INTEGER113:* The leading dimension of the array A. LDA >= max(1,N).114:*115:* AF (input or output) DOUBLE PRECISION array, dimension (LDAF,N)116:* If FACT = 'F', then AF is an input argument and on entry117:* contains the factors L and U from the factorization118:* A = P*L*U as computed by DGETRF. If EQUED .ne. 'N', then119:* AF is the factored form of the equilibrated matrix A.120:*121:* If FACT = 'N', then AF is an output argument and on exit122:* returns the factors L and U from the factorization A = P*L*U123:* of the original matrix A.124:*125:* If FACT = 'E', then AF is an output argument and on exit126:* returns the factors L and U from the factorization A = P*L*U127:* of the equilibrated matrix A (see the description of A for128:* the form of the equilibrated matrix).129:*130:* LDAF (input) INTEGER131:* The leading dimension of the array AF. LDAF >= max(1,N).132:*133:* IPIV (input or output) INTEGER array, dimension (N)134:* If FACT = 'F', then IPIV is an input argument and on entry135:* contains the pivot indices from the factorization A = P*L*U136:* as computed by DGETRF; row i of the matrix was interchanged137:* with row IPIV(i).138:*139:* If FACT = 'N', then IPIV is an output argument and on exit140:* contains the pivot indices from the factorization A = P*L*U141:* of the original matrix A.142:*143:* If FACT = 'E', then IPIV is an output argument and on exit144:* contains the pivot indices from the factorization A = P*L*U145:* of the equilibrated matrix A.146:*147:* EQUED (input or output) CHARACTER*1148:* Specifies the form of equilibration that was done.149:* = 'N': No equilibration (always true if FACT = 'N').150:* = 'R': Row equilibration, i.e., A has been premultiplied by151:* diag(R).152:* = 'C': Column equilibration, i.e., A has been postmultiplied153:* by diag(C).154:* = 'B': Both row and column equilibration, i.e., A has been155:* replaced by diag(R) * A * diag(C).156:* EQUED is an input argument if FACT = 'F'; otherwise, it is an157:* output argument.158:*159:* R (input or output) DOUBLE PRECISION array, dimension (N)160:* The row scale factors for A. If EQUED = 'R' or 'B', A is161:* multiplied on the left by diag(R); if EQUED = 'N' or 'C', R162:* is not accessed. R is an input argument if FACT = 'F';163:* otherwise, R is an output argument. If FACT = 'F' and164:* EQUED = 'R' or 'B', each element of R must be positive.165:*166:* C (input or output) DOUBLE PRECISION array, dimension (N)167:* The column scale factors for A. If EQUED = 'C' or 'B', A is168:* multiplied on the right by diag(C); if EQUED = 'N' or 'R', C169:* is not accessed. C is an input argument if FACT = 'F';170:* otherwise, C is an output argument. If FACT = 'F' and171:* EQUED = 'C' or 'B', each element of C must be positive.172:*173:* B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)174:* On entry, the N-by-NRHS right hand side matrix B.175:* On exit,176:* if EQUED = 'N', B is not modified;177:* if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by178:* diag(R)*B;179:* if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is180:* overwritten by diag(C)*B.181:*182:* LDB (input) INTEGER183:* The leading dimension of the array B. LDB >= max(1,N).184:*185:* X (output) DOUBLE PRECISION array, dimension (LDX,NRHS)186:* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X187:* to the original system of equations. Note that A and B are188:* modified on exit if EQUED .ne. 'N', and the solution to the189:* equilibrated system is inv(diag(C))*X if TRANS = 'N' and190:* EQUED = 'C' or 'B', or inv(diag(R))*X if TRANS = 'T' or 'C'191:* and EQUED = 'R' or 'B'.192:*193:* LDX (input) INTEGER194:* The leading dimension of the array X. LDX >= max(1,N).195:*196:* RCOND (output) DOUBLE PRECISION197:* The estimate of the reciprocal condition number of the matrix198:* A after equilibration (if done). If RCOND is less than the199:* machine precision (in particular, if RCOND = 0), the matrix200:* is singular to working precision. This condition is201:* indicated by a return code of INFO > 0.202:*203:* FERR (output) DOUBLE PRECISION array, dimension (NRHS)204:* The estimated forward error bound for each solution vector205:* X(j) (the j-th column of the solution matrix X).206:* If XTRUE is the true solution corresponding to X(j), FERR(j)207:* is an estimated upper bound for the magnitude of the largest208:* element in (X(j) - XTRUE) divided by the magnitude of the209:* largest element in X(j). The estimate is as reliable as210:* the estimate for RCOND, and is almost always a slight211:* overestimate of the true error.212:*213:* BERR (output) DOUBLE PRECISION array, dimension (NRHS)214:* The componentwise relative backward error of each solution215:* vector X(j) (i.e., the smallest relative change in216:* any element of A or B that makes X(j) an exact solution).217:*218:* WORK (workspace/output) DOUBLE PRECISION array, dimension (4*N)219:* On exit, WORK(1) contains the reciprocal pivot growth220:* factor norm(A)/norm(U). The "max absolute element" norm is221:* used. If WORK(1) is much less than 1, then the stability222:* of the LU factorization of the (equilibrated) matrix A223:* could be poor. This also means that the solution X, condition224:* estimator RCOND, and forward error bound FERR could be225:* unreliable. If factorization fails with 0<INFO<=N, then226:* WORK(1) contains the reciprocal pivot growth factor for the227:* leading INFO columns of A.228:*229:* IWORK (workspace) INTEGER array, dimension (N)230:*231:* INFO (output) INTEGER232:* = 0: successful exit233:* < 0: if INFO = -i, the i-th argument had an illegal value234:* > 0: if INFO = i, and i is235:* <= N: U(i,i) is exactly zero. The factorization has236:* been completed, but the factor U is exactly237:* singular, so the solution and error bounds238:* could not be computed. RCOND = 0 is returned.239:* = N+1: U is nonsingular, but RCOND is less than machine240:* precision, meaning that the matrix is singular241:* to working precision. Nevertheless, the242:* solution and error bounds are computed because243:* there are a number of situations where the244:* computed solution can be more accurate than the245:* value of RCOND would suggest.246:*247:* =====================================================================248:*249:* .. Parameters ..250: DOUBLE PRECISION ZERO, ONE 251:PARAMETER( ZERO = 0.0D+0, ONE = 1.0D+0 ) 252:* ..253:* .. Local Scalars ..254:LOGICALCOLEQU, EQUIL, NOFACT, NOTRAN, ROWEQU 255: CHARACTER NORM 256: INTEGER I, INFEQU, J 257: DOUBLE PRECISION AMAX, ANORM, BIGNUM, COLCND, RCMAX, RCMIN, 258: $ ROWCND, RPVGRW, SMLNUM 259:* ..260:* .. External Functions ..261:LOGICALLSAME 262: DOUBLE PRECISION DLAMCH, DLANGE, DLANTR 263:EXTERNALLSAME, DLAMCH, DLANGE, DLANTR 264:* ..265:* .. External Subroutines ..266:EXTERNALDGECON, DGEEQU, DGERFS, DGETRF, DGETRS, DLACPY, 267: $ DLAQGE, XERBLA 268:* ..269:* .. Intrinsic Functions ..270:INTRINSICMAX, MIN 271:* ..272:* .. Executable Statements ..273:*274: INFO = 0 275: NOFACT =LSAME( FACT, 'N' ) 276: EQUIL =LSAME( FACT, 'E' ) 277: NOTRAN =LSAME( TRANS, 'N' ) 278:IF( NOFACT .OR. EQUIL )THEN279: EQUED = 'N' 280: ROWEQU = .FALSE. 281: COLEQU = .FALSE. 282:ELSE283: ROWEQU =LSAME( EQUED, 'R' ) .OR.LSAME( EQUED, 'B' ) 284: COLEQU =LSAME( EQUED, 'C' ) .OR.LSAME( EQUED, 'B' ) 285: SMLNUM =DLAMCH( 'Safe minimum' ) 286: BIGNUM = ONE / SMLNUM 287:ENDIF288:*289:* Test the input parameters.290:*291:IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) ) 292: $THEN293: INFO = -1 294:ELSEIF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. 295: $LSAME( TRANS, 'C' ) )THEN296: INFO = -2 297:ELSEIF( N.LT.0 )THEN298: INFO = -3 299:ELSEIF( NRHS.LT.0 )THEN300: INFO = -4 301:ELSEIF( LDA.LT.MAX( 1, N ) )THEN302: INFO = -6 303:ELSEIF( LDAF.LT.MAX( 1, N ) )THEN304: INFO = -8 305:ELSEIF(LSAME( FACT, 'F' ) .AND. .NOT. 306: $ ( ROWEQU .OR. COLEQU .OR.LSAME( EQUED, 'N' ) ) )THEN307: INFO = -10 308:ELSE309:IF( ROWEQU )THEN310: RCMIN = BIGNUM 311: RCMAX = ZERO 312:DO10 J = 1, N 313: RCMIN =MIN( RCMIN,R( J ) ) 314: RCMAX =MAX( RCMAX,R( J ) ) 315: 10CONTINUE316:IF( RCMIN.LE.ZERO )THEN317: INFO = -11 318:ELSEIF( N.GT.0 )THEN319: ROWCND =MAX( RCMIN, SMLNUM ) /MIN( RCMAX, BIGNUM ) 320:ELSE321: ROWCND = ONE 322:ENDIF323:ENDIF324:IF( COLEQU .AND. INFO.EQ.0 )THEN325: RCMIN = BIGNUM 326: RCMAX = ZERO 327:DO20 J = 1, N 328: RCMIN =MIN( RCMIN,C( J ) ) 329: RCMAX =MAX( RCMAX,C( J ) ) 330: 20CONTINUE331:IF( RCMIN.LE.ZERO )THEN332: INFO = -12 333:ELSEIF( N.GT.0 )THEN334: COLCND =MAX( RCMIN, SMLNUM ) /MIN( RCMAX, BIGNUM ) 335:ELSE336: COLCND = ONE 337:ENDIF338:ENDIF339:IF( INFO.EQ.0 )THEN340:IF( LDB.LT.MAX( 1, N ) )THEN341: INFO = -14 342:ELSEIF( LDX.LT.MAX( 1, N ) )THEN343: INFO = -16 344:ENDIF345:ENDIF346:ENDIF347:*348:IF( INFO.NE.0 )THEN349:CALLXERBLA( 'DGESVX', -INFO ) 350:RETURN351:ENDIF352:*353:IF( EQUIL )THEN354:*355:* Compute row and column scalings to equilibrate the matrix A.356:*357:CALLDGEEQU( N, N, A, LDA, R, C, ROWCND, COLCND, AMAX, INFEQU ) 358:IF( INFEQU.EQ.0 )THEN359:*360:* Equilibrate the matrix.361:*362:CALLDLAQGE( N, N, A, LDA, R, C, ROWCND, COLCND, AMAX, 363: $ EQUED ) 364: ROWEQU =LSAME( EQUED, 'R' ) .OR.LSAME( EQUED, 'B' ) 365: COLEQU =LSAME( EQUED, 'C' ) .OR.LSAME( EQUED, 'B' ) 366:ENDIF367:ENDIF368:*369:* Scale the right hand side.370:*371:IF( NOTRAN )THEN372:IF( ROWEQU )THEN373:DO40 J = 1, NRHS 374:DO30 I = 1, N 375:B( I, J ) =R( I )*B( I, J ) 376: 30CONTINUE377: 40CONTINUE378:ENDIF379:ELSEIF( COLEQU )THEN380:DO60 J = 1, NRHS 381:DO50 I = 1, N 382:B( I, J ) =C( I )*B( I, J ) 383: 50CONTINUE384: 60CONTINUE385:ENDIF386:*387:IF( NOFACT .OR. EQUIL )THEN388:*389:* Compute the LU factorization of A.390:*391:CALLDLACPY( 'Full', N, N, A, LDA, AF, LDAF ) 392:CALLDGETRF( N, N, AF, LDAF, IPIV, INFO ) 393:*394:* Return if INFO is non-zero.395:*396:IF( INFO.GT.0 )THEN397:*398:* Compute the reciprocal pivot growth factor of the399:* leading rank-deficient INFO columns of A.400:*401: RPVGRW =DLANTR( 'M', 'U', 'N', INFO, INFO, AF, LDAF, 402: $ WORK ) 403:IF( RPVGRW.EQ.ZERO )THEN404: RPVGRW = ONE 405:ELSE406: RPVGRW =DLANGE( 'M', N, INFO, A, LDA, WORK ) / RPVGRW 407:ENDIF408:WORK( 1 ) = RPVGRW 409: RCOND = ZERO 410:RETURN411:ENDIF412:ENDIF413:*414:* Compute the norm of the matrix A and the415:* reciprocal pivot growth factor RPVGRW.416:*417:IF( NOTRAN )THEN418: NORM = '1' 419:ELSE420: NORM = 'I' 421:ENDIF422: ANORM =DLANGE( NORM, N, N, A, LDA, WORK ) 423: RPVGRW =DLANTR( 'M', 'U', 'N', N, N, AF, LDAF, WORK ) 424:IF( RPVGRW.EQ.ZERO )THEN425: RPVGRW = ONE 426:ELSE427: RPVGRW =DLANGE( 'M', N, N, A, LDA, WORK ) / RPVGRW 428:ENDIF429:*430:* Compute the reciprocal of the condition number of A.431:*432:CALLDGECON( NORM, N, AF, LDAF, ANORM, RCOND, WORK, IWORK, INFO ) 433:*434:* Compute the solution matrix X.435:*436:CALLDLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 437:CALLDGETRS( TRANS, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO ) 438:*439:* Use iterative refinement to improve the computed solution and440:* compute error bounds and backward error estimates for it.441:*442:CALLDGERFS( TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, 443: $ LDX, FERR, BERR, WORK, IWORK, INFO ) 444:*445:* Transform the solution matrix X to a solution of the original446:* system.447:*448:IF( NOTRAN )THEN449:IF( COLEQU )THEN450:DO80 J = 1, NRHS 451:DO70 I = 1, N 452:X( I, J ) =C( I )*X( I, J ) 453: 70CONTINUE454: 80CONTINUE455:DO90 J = 1, NRHS 456:FERR( J ) =FERR( J ) / COLCND 457: 90CONTINUE458:ENDIF459:ELSEIF( ROWEQU )THEN460:DO110 J = 1, NRHS 461:DO100 I = 1, N 462:X( I, J ) =R( I )*X( I, J ) 463: 100CONTINUE464: 110CONTINUE465:DO120 J = 1, NRHS 466:FERR( J ) =FERR( J ) / ROWCND 467: 120CONTINUE468:ENDIF469:*470:WORK( 1 ) = RPVGRW 471:*472:* Set INFO = N+1 if the matrix is singular to working precision.473:*474:IF( RCOND.LT.DLAMCH( 'Epsilon' ) ) 475: $ INFO = N + 1 476:RETURN477:*478:* End of DGESVX479:*480:END481: