001:SUBROUTINESPOSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED, 002: $ S, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, 003: $ IWORK, INFO ) 004:*005:* -- LAPACK driver routine (version 3.2) --006:* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..007:* November 2006008:*009:* .. Scalar Arguments ..010: CHARACTER EQUED, FACT, UPLO 011: INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS 012: REAL RCOND 013:* ..014:* .. Array Arguments ..015: INTEGERIWORK( * ) 016: REALA( LDA, * ),AF( LDAF, * ),B( LDB, * ), 017: $BERR( * ),FERR( * ),S( * ),WORK( * ), 018: $X( LDX, * ) 019:* ..020:*021:* Purpose022:* =======023:*024:* SPOSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to025:* compute the solution to a real system of linear equations026:* A * X = B,027:* where A is an N-by-N symmetric positive definite matrix and X and B028:* 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:* diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B041:* Whether or not the system will be equilibrated depends on the042:* scaling of the matrix A, but if equilibration is used, A is043:* overwritten by diag(S)*A*diag(S) and B by diag(S)*B.044:*045:* 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to046:* factor the matrix A (after equilibration if FACT = 'E') as047:* A = U**T* U, if UPLO = 'U', or048:* A = L * L**T, if UPLO = 'L',049:* where U is an upper triangular matrix and L is a lower triangular050:* matrix.051:*052:* 3. If the leading i-by-i principal minor is not positive definite,053:* then the routine returns with INFO = i. Otherwise, the factored054:* form of A is used to estimate the condition number of the matrix055:* A. If the reciprocal of the condition number is less than machine056:* precision, INFO = N+1 is returned as a warning, but the routine057:* still goes on to solve for X and compute error bounds as058:* described below.059:*060:* 4. The system of equations is solved for X using the factored form061:* of A.062:*063:* 5. Iterative refinement is applied to improve the computed solution064:* matrix and calculate error bounds and backward error estimates065:* for it.066:*067:* 6. If equilibration was used, the matrix X is premultiplied by068:* diag(S) so that it solves the original system before069:* equilibration.070:*071:* Arguments072:* =========073:*074:* FACT (input) CHARACTER*1075:* Specifies whether or not the factored form of the matrix A is076:* supplied on entry, and if not, whether the matrix A should be077:* equilibrated before it is factored.078:* = 'F': On entry, AF contains the factored form of A.079:* If EQUED = 'Y', the matrix A has been equilibrated080:* with scaling factors given by S. A and AF will not081:* be modified.082:* = 'N': The matrix A will be copied to AF and factored.083:* = 'E': The matrix A will be equilibrated if necessary, then084:* copied to AF and factored.085:*086:* UPLO (input) CHARACTER*1087:* = 'U': Upper triangle of A is stored;088:* = 'L': Lower triangle of A is stored.089:*090:* N (input) INTEGER091:* The number of linear equations, i.e., the order of the092:* matrix A. N >= 0.093:*094:* NRHS (input) INTEGER095:* The number of right hand sides, i.e., the number of columns096:* of the matrices B and X. NRHS >= 0.097:*098:* A (input/output) REAL array, dimension (LDA,N)099:* On entry, the symmetric matrix A, except if FACT = 'F' and100:* EQUED = 'Y', then A must contain the equilibrated matrix101:* diag(S)*A*diag(S). If UPLO = 'U', the leading102:* N-by-N upper triangular part of A contains the upper103:* triangular part of the matrix A, and the strictly lower104:* triangular part of A is not referenced. If UPLO = 'L', the105:* leading N-by-N lower triangular part of A contains the lower106:* triangular part of the matrix A, and the strictly upper107:* triangular part of A is not referenced. A is not modified if108:* FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit.109:*110:* On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by111:* diag(S)*A*diag(S).112:*113:* LDA (input) INTEGER114:* The leading dimension of the array A. LDA >= max(1,N).115:*116:* AF (input or output) REAL array, dimension (LDAF,N)117:* If FACT = 'F', then AF is an input argument and on entry118:* contains the triangular factor U or L from the Cholesky119:* factorization A = U**T*U or A = L*L**T, in the same storage120:* format as A. If EQUED .ne. 'N', then AF is the factored form121:* of the equilibrated matrix diag(S)*A*diag(S).122:*123:* If FACT = 'N', then AF is an output argument and on exit124:* returns the triangular factor U or L from the Cholesky125:* factorization A = U**T*U or A = L*L**T of the original126:* matrix A.127:*128:* If FACT = 'E', then AF is an output argument and on exit129:* returns the triangular factor U or L from the Cholesky130:* factorization A = U**T*U or A = L*L**T of the equilibrated131:* matrix A (see the description of A for the form of the132:* equilibrated matrix).133:*134:* LDAF (input) INTEGER135:* The leading dimension of the array AF. LDAF >= max(1,N).136:*137:* EQUED (input or output) CHARACTER*1138:* Specifies the form of equilibration that was done.139:* = 'N': No equilibration (always true if FACT = 'N').140:* = 'Y': Equilibration was done, i.e., A has been replaced by141:* diag(S) * A * diag(S).142:* EQUED is an input argument if FACT = 'F'; otherwise, it is an143:* output argument.144:*145:* S (input or output) REAL array, dimension (N)146:* The scale factors for A; not accessed if EQUED = 'N'. S is147:* an input argument if FACT = 'F'; otherwise, S is an output148:* argument. If FACT = 'F' and EQUED = 'Y', each element of S149:* must be positive.150:*151:* B (input/output) REAL array, dimension (LDB,NRHS)152:* On entry, the N-by-NRHS right hand side matrix B.153:* On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y',154:* B is overwritten by diag(S) * B.155:*156:* LDB (input) INTEGER157:* The leading dimension of the array B. LDB >= max(1,N).158:*159:* X (output) REAL array, dimension (LDX,NRHS)160:* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to161:* the original system of equations. Note that if EQUED = 'Y',162:* A and B are modified on exit, and the solution to the163:* equilibrated system is inv(diag(S))*X.164:*165:* LDX (input) INTEGER166:* The leading dimension of the array X. LDX >= max(1,N).167:*168:* RCOND (output) REAL169:* The estimate of the reciprocal condition number of the matrix170:* A after equilibration (if done). If RCOND is less than the171:* machine precision (in particular, if RCOND = 0), the matrix172:* is singular to working precision. This condition is173:* indicated by a return code of INFO > 0.174:*175:* FERR (output) REAL array, dimension (NRHS)176:* The estimated forward error bound for each solution vector177:* X(j) (the j-th column of the solution matrix X).178:* If XTRUE is the true solution corresponding to X(j), FERR(j)179:* is an estimated upper bound for the magnitude of the largest180:* element in (X(j) - XTRUE) divided by the magnitude of the181:* largest element in X(j). The estimate is as reliable as182:* the estimate for RCOND, and is almost always a slight183:* overestimate of the true error.184:*185:* BERR (output) REAL array, dimension (NRHS)186:* The componentwise relative backward error of each solution187:* vector X(j) (i.e., the smallest relative change in188:* any element of A or B that makes X(j) an exact solution).189:*190:* WORK (workspace) REAL array, dimension (3*N)191:*192:* IWORK (workspace) INTEGER array, dimension (N)193:*194:* INFO (output) INTEGER195:* = 0: successful exit196:* < 0: if INFO = -i, the i-th argument had an illegal value197:* > 0: if INFO = i, and i is198:* <= N: the leading minor of order i of A is199:* not positive definite, so the factorization200:* could not be completed, and the solution has not201:* been computed. RCOND = 0 is returned.202:* = N+1: U is nonsingular, but RCOND is less than machine203:* precision, meaning that the matrix is singular204:* to working precision. Nevertheless, the205:* solution and error bounds are computed because206:* there are a number of situations where the207:* computed solution can be more accurate than the208:* value of RCOND would suggest.209:*210:* =====================================================================211:*212:* .. Parameters ..213: REAL ZERO, ONE 214:PARAMETER( ZERO = 0.0E+0, ONE = 1.0E+0 ) 215:* ..216:* .. Local Scalars ..217:LOGICALEQUIL, NOFACT, RCEQU 218: INTEGER I, INFEQU, J 219: REAL AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM 220:* ..221:* .. External Functions ..222:LOGICALLSAME 223: REAL SLAMCH, SLANSY 224:EXTERNALLSAME, SLAMCH, SLANSY 225:* ..226:* .. External Subroutines ..227:EXTERNALSLACPY, SLAQSY, SPOCON, SPOEQU, SPORFS, SPOTRF, 228: $ SPOTRS, XERBLA 229:* ..230:* .. Intrinsic Functions ..231:INTRINSICMAX, MIN 232:* ..233:* .. Executable Statements ..234:*235: INFO = 0 236: NOFACT =LSAME( FACT, 'N' ) 237: EQUIL =LSAME( FACT, 'E' ) 238:IF( NOFACT .OR. EQUIL )THEN239: EQUED = 'N' 240: RCEQU = .FALSE. 241:ELSE242: RCEQU =LSAME( EQUED, 'Y' ) 243: SMLNUM =SLAMCH( 'Safe minimum' ) 244: BIGNUM = ONE / SMLNUM 245:ENDIF246:*247:* Test the input parameters.248:*249:IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) ) 250: $THEN251: INFO = -1 252:ELSEIF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) 253: $THEN254: INFO = -2 255:ELSEIF( N.LT.0 )THEN256: INFO = -3 257:ELSEIF( NRHS.LT.0 )THEN258: INFO = -4 259:ELSEIF( LDA.LT.MAX( 1, N ) )THEN260: INFO = -6 261:ELSEIF( LDAF.LT.MAX( 1, N ) )THEN262: INFO = -8 263:ELSEIF(LSAME( FACT, 'F' ) .AND. .NOT. 264: $ ( RCEQU .OR.LSAME( EQUED, 'N' ) ) )THEN265: INFO = -9 266:ELSE267:IF( RCEQU )THEN268: SMIN = BIGNUM 269: SMAX = ZERO 270:DO10 J = 1, N 271: SMIN =MIN( SMIN,S( J ) ) 272: SMAX =MAX( SMAX,S( J ) ) 273: 10CONTINUE274:IF( SMIN.LE.ZERO )THEN275: INFO = -10 276:ELSEIF( N.GT.0 )THEN277: SCOND =MAX( SMIN, SMLNUM ) /MIN( SMAX, BIGNUM ) 278:ELSE279: SCOND = ONE 280:ENDIF281:ENDIF282:IF( INFO.EQ.0 )THEN283:IF( LDB.LT.MAX( 1, N ) )THEN284: INFO = -12 285:ELSEIF( LDX.LT.MAX( 1, N ) )THEN286: INFO = -14 287:ENDIF288:ENDIF289:ENDIF290:*291:IF( INFO.NE.0 )THEN292:CALLXERBLA( 'SPOSVX', -INFO ) 293:RETURN294:ENDIF295:*296:IF( EQUIL )THEN297:*298:* Compute row and column scalings to equilibrate the matrix A.299:*300:CALLSPOEQU( N, A, LDA, S, SCOND, AMAX, INFEQU ) 301:IF( INFEQU.EQ.0 )THEN302:*303:* Equilibrate the matrix.304:*305:CALLSLAQSY( UPLO, N, A, LDA, S, SCOND, AMAX, EQUED ) 306: RCEQU =LSAME( EQUED, 'Y' ) 307:ENDIF308:ENDIF309:*310:* Scale the right hand side.311:*312:IF( RCEQU )THEN313:DO30 J = 1, NRHS 314:DO20 I = 1, N 315:B( I, J ) =S( I )*B( I, J ) 316: 20CONTINUE317: 30CONTINUE318:ENDIF319:*320:IF( NOFACT .OR. EQUIL )THEN321:*322:* Compute the Cholesky factorization A = U'*U or A = L*L'.323:*324:CALLSLACPY( UPLO, N, N, A, LDA, AF, LDAF ) 325:CALLSPOTRF( UPLO, N, AF, LDAF, INFO ) 326:*327:* Return if INFO is non-zero.328:*329:IF( INFO.GT.0 )THEN330: RCOND = ZERO 331:RETURN332:ENDIF333:ENDIF334:*335:* Compute the norm of the matrix A.336:*337: ANORM =SLANSY( '1', UPLO, N, A, LDA, WORK ) 338:*339:* Compute the reciprocal of the condition number of A.340:*341:CALLSPOCON( UPLO, N, AF, LDAF, ANORM, RCOND, WORK, IWORK, INFO ) 342:*343:* Compute the solution matrix X.344:*345:CALLSLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 346:CALLSPOTRS( UPLO, N, NRHS, AF, LDAF, X, LDX, INFO ) 347:*348:* Use iterative refinement to improve the computed solution and349:* compute error bounds and backward error estimates for it.350:*351:CALLSPORFS( UPLO, N, NRHS, A, LDA, AF, LDAF, B, LDB, X, LDX, 352: $ FERR, BERR, WORK, IWORK, INFO ) 353:*354:* Transform the solution matrix X to a solution of the original355:* system.356:*357:IF( RCEQU )THEN358:DO50 J = 1, NRHS 359:DO40 I = 1, N 360:X( I, J ) =S( I )*X( I, J ) 361: 40CONTINUE362: 50CONTINUE363:DO60 J = 1, NRHS 364:FERR( J ) =FERR( J ) / SCOND 365: 60CONTINUE366:ENDIF367:*368:* Set INFO = N+1 if the matrix is singular to working precision.369:*370:IF( RCOND.LT.SLAMCH( 'Epsilon' ) ) 371: $ INFO = N + 1 372:*373:RETURN374:*375:* End of SPOSVX376:*377:END378: