001:SUBROUTINESPPSVX( FACT, UPLO, N, NRHS, AP, AFP, EQUED, S, B, LDB, 002: $ X, LDX, RCOND, FERR, BERR, WORK, IWORK, INFO ) 003:*004:* -- LAPACK driver routine (version 3.2) --005:* -- LAPACK is a software package provided by Univ. of Tennessee, --006:* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--007:* November 2006008:*009:* .. Scalar Arguments ..010: CHARACTER EQUED, FACT, UPLO 011: INTEGER INFO, LDB, LDX, N, NRHS 012: REAL RCOND 013:* ..014:* .. Array Arguments ..015: INTEGERIWORK( * ) 016: REALAFP( * ),AP( * ),B( LDB, * ),BERR( * ), 017: $FERR( * ),S( * ),WORK( * ),X( LDX, * ) 018:* ..019:*020:* Purpose021:* =======022:*023:* SPPSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to024:* compute the solution to a real system of linear equations025:* A * X = B,026:* where A is an N-by-N symmetric positive definite matrix stored in027:* packed format and X and B are N-by-NRHS matrices.028:*029:* Error bounds on the solution and a condition estimate are also030:* provided.031:*032:* Description033:* ===========034:*035:* The following steps are performed:036:*037:* 1. If FACT = 'E', real scaling factors are computed to equilibrate038:* the system:039:* diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B040:* Whether or not the system will be equilibrated depends on the041:* scaling of the matrix A, but if equilibration is used, A is042:* overwritten by diag(S)*A*diag(S) and B by diag(S)*B.043:*044:* 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to045:* factor the matrix A (after equilibration if FACT = 'E') as046:* A = U**T* U, if UPLO = 'U', or047:* A = L * L**T, if UPLO = 'L',048:* where U is an upper triangular matrix and L is a lower triangular049:* matrix.050:*051:* 3. If the leading i-by-i principal minor is not positive definite,052:* then the routine returns with INFO = i. Otherwise, the factored053:* form of A is used to estimate the condition number of the matrix054:* A. If the reciprocal of the condition number is less than machine055:* precision, INFO = N+1 is returned as a warning, but the routine056:* still goes on to solve for X and compute error bounds as057:* described below.058:*059:* 4. The system of equations is solved for X using the factored form060:* of A.061:*062:* 5. Iterative refinement is applied to improve the computed solution063:* matrix and calculate error bounds and backward error estimates064:* for it.065:*066:* 6. If equilibration was used, the matrix X is premultiplied by067:* diag(S) so that it solves the original system before068:* equilibration.069:*070:* Arguments071:* =========072:*073:* FACT (input) CHARACTER*1074:* Specifies whether or not the factored form of the matrix A is075:* supplied on entry, and if not, whether the matrix A should be076:* equilibrated before it is factored.077:* = 'F': On entry, AFP contains the factored form of A.078:* If EQUED = 'Y', the matrix A has been equilibrated079:* with scaling factors given by S. AP and AFP will not080:* be modified.081:* = 'N': The matrix A will be copied to AFP and factored.082:* = 'E': The matrix A will be equilibrated if necessary, then083:* copied to AFP and factored.084:*085:* UPLO (input) CHARACTER*1086:* = 'U': Upper triangle of A is stored;087:* = 'L': Lower triangle of A is stored.088:*089:* N (input) INTEGER090:* The number of linear equations, i.e., the order of the091:* matrix A. N >= 0.092:*093:* NRHS (input) INTEGER094:* The number of right hand sides, i.e., the number of columns095:* of the matrices B and X. NRHS >= 0.096:*097:* AP (input/output) REAL array, dimension (N*(N+1)/2)098:* On entry, the upper or lower triangle of the symmetric matrix099:* A, packed columnwise in a linear array, except if FACT = 'F'100:* and EQUED = 'Y', then A must contain the equilibrated matrix101:* diag(S)*A*diag(S). The j-th column of A is stored in the102:* array AP as follows:103:* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;104:* if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.105:* See below for further details. A is not modified if106:* FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit.107:*108:* On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by109:* diag(S)*A*diag(S).110:*111:* AFP (input or output) REAL array, dimension112:* (N*(N+1)/2)113:* If FACT = 'F', then AFP is an input argument and on entry114:* contains the triangular factor U or L from the Cholesky115:* factorization A = U'*U or A = L*L', in the same storage116:* format as A. If EQUED .ne. 'N', then AFP is the factored117:* form of the equilibrated matrix A.118:*119:* If FACT = 'N', then AFP is an output argument and on exit120:* returns the triangular factor U or L from the Cholesky121:* factorization A = U'*U or A = L*L' of the original matrix A.122:*123:* If FACT = 'E', then AFP is an output argument and on exit124:* returns the triangular factor U or L from the Cholesky125:* factorization A = U'*U or A = L*L' of the equilibrated126:* matrix A (see the description of AP for the form of the127:* equilibrated matrix).128:*129:* EQUED (input or output) CHARACTER*1130:* Specifies the form of equilibration that was done.131:* = 'N': No equilibration (always true if FACT = 'N').132:* = 'Y': Equilibration was done, i.e., A has been replaced by133:* diag(S) * A * diag(S).134:* EQUED is an input argument if FACT = 'F'; otherwise, it is an135:* output argument.136:*137:* S (input or output) REAL array, dimension (N)138:* The scale factors for A; not accessed if EQUED = 'N'. S is139:* an input argument if FACT = 'F'; otherwise, S is an output140:* argument. If FACT = 'F' and EQUED = 'Y', each element of S141:* must be positive.142:*143:* B (input/output) REAL array, dimension (LDB,NRHS)144:* On entry, the N-by-NRHS right hand side matrix B.145:* On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y',146:* B is overwritten by diag(S) * B.147:*148:* LDB (input) INTEGER149:* The leading dimension of the array B. LDB >= max(1,N).150:*151:* X (output) REAL array, dimension (LDX,NRHS)152:* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to153:* the original system of equations. Note that if EQUED = 'Y',154:* A and B are modified on exit, and the solution to the155:* equilibrated system is inv(diag(S))*X.156:*157:* LDX (input) INTEGER158:* The leading dimension of the array X. LDX >= max(1,N).159:*160:* RCOND (output) REAL161:* The estimate of the reciprocal condition number of the matrix162:* A after equilibration (if done). If RCOND is less than the163:* machine precision (in particular, if RCOND = 0), the matrix164:* is singular to working precision. This condition is165:* indicated by a return code of INFO > 0.166:*167:* FERR (output) REAL array, dimension (NRHS)168:* The estimated forward error bound for each solution vector169:* X(j) (the j-th column of the solution matrix X).170:* If XTRUE is the true solution corresponding to X(j), FERR(j)171:* is an estimated upper bound for the magnitude of the largest172:* element in (X(j) - XTRUE) divided by the magnitude of the173:* largest element in X(j). The estimate is as reliable as174:* the estimate for RCOND, and is almost always a slight175:* overestimate of the true error.176:*177:* BERR (output) REAL array, dimension (NRHS)178:* The componentwise relative backward error of each solution179:* vector X(j) (i.e., the smallest relative change in180:* any element of A or B that makes X(j) an exact solution).181:*182:* WORK (workspace) REAL array, dimension (3*N)183:*184:* IWORK (workspace) INTEGER array, dimension (N)185:*186:* INFO (output) INTEGER187:* = 0: successful exit188:* < 0: if INFO = -i, the i-th argument had an illegal value189:* > 0: if INFO = i, and i is190:* <= N: the leading minor of order i of A is191:* not positive definite, so the factorization192:* could not be completed, and the solution has not193:* been computed. RCOND = 0 is returned.194:* = N+1: U is nonsingular, but RCOND is less than machine195:* precision, meaning that the matrix is singular196:* to working precision. Nevertheless, the197:* solution and error bounds are computed because198:* there are a number of situations where the199:* computed solution can be more accurate than the200:* value of RCOND would suggest.201:*202:* Further Details203:* ===============204:*205:* The packed storage scheme is illustrated by the following example206:* when N = 4, UPLO = 'U':207:*208:* Two-dimensional storage of the symmetric matrix A:209:*210:* a11 a12 a13 a14211:* a22 a23 a24212:* a33 a34 (aij = conjg(aji))213:* a44214:*215:* Packed storage of the upper triangle of A:216:*217:* AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]218:*219:* =====================================================================220:*221:* .. Parameters ..222: REAL ZERO, ONE 223:PARAMETER( ZERO = 0.0E+0, ONE = 1.0E+0 ) 224:* ..225:* .. Local Scalars ..226:LOGICALEQUIL, NOFACT, RCEQU 227: INTEGER I, INFEQU, J 228: REAL AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM 229:* ..230:* .. External Functions ..231:LOGICALLSAME 232: REAL SLAMCH, SLANSP 233:EXTERNALLSAME, SLAMCH, SLANSP 234:* ..235:* .. External Subroutines ..236:EXTERNALSCOPY, SLACPY, SLAQSP, SPPCON, SPPEQU, SPPRFS, 237: $ SPPTRF, SPPTRS, XERBLA 238:* ..239:* .. Intrinsic Functions ..240:INTRINSICMAX, MIN 241:* ..242:* .. Executable Statements ..243:*244: INFO = 0 245: NOFACT =LSAME( FACT, 'N' ) 246: EQUIL =LSAME( FACT, 'E' ) 247:IF( NOFACT .OR. EQUIL )THEN248: EQUED = 'N' 249: RCEQU = .FALSE. 250:ELSE251: RCEQU =LSAME( EQUED, 'Y' ) 252: SMLNUM =SLAMCH( 'Safe minimum' ) 253: BIGNUM = ONE / SMLNUM 254:ENDIF255:*256:* Test the input parameters.257:*258:IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) ) 259: $THEN260: INFO = -1 261:ELSEIF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) 262: $THEN263: INFO = -2 264:ELSEIF( N.LT.0 )THEN265: INFO = -3 266:ELSEIF( NRHS.LT.0 )THEN267: INFO = -4 268:ELSEIF(LSAME( FACT, 'F' ) .AND. .NOT. 269: $ ( RCEQU .OR.LSAME( EQUED, 'N' ) ) )THEN270: INFO = -7 271:ELSE272:IF( RCEQU )THEN273: SMIN = BIGNUM 274: SMAX = ZERO 275:DO10 J = 1, N 276: SMIN =MIN( SMIN,S( J ) ) 277: SMAX =MAX( SMAX,S( J ) ) 278: 10CONTINUE279:IF( SMIN.LE.ZERO )THEN280: INFO = -8 281:ELSEIF( N.GT.0 )THEN282: SCOND =MAX( SMIN, SMLNUM ) /MIN( SMAX, BIGNUM ) 283:ELSE284: SCOND = ONE 285:ENDIF286:ENDIF287:IF( INFO.EQ.0 )THEN288:IF( LDB.LT.MAX( 1, N ) )THEN289: INFO = -10 290:ELSEIF( LDX.LT.MAX( 1, N ) )THEN291: INFO = -12 292:ENDIF293:ENDIF294:ENDIF295:*296:IF( INFO.NE.0 )THEN297:CALLXERBLA( 'SPPSVX', -INFO ) 298:RETURN299:ENDIF300:*301:IF( EQUIL )THEN302:*303:* Compute row and column scalings to equilibrate the matrix A.304:*305:CALLSPPEQU( UPLO, N, AP, S, SCOND, AMAX, INFEQU ) 306:IF( INFEQU.EQ.0 )THEN307:*308:* Equilibrate the matrix.309:*310:CALLSLAQSP( UPLO, N, AP, S, SCOND, AMAX, EQUED ) 311: RCEQU =LSAME( EQUED, 'Y' ) 312:ENDIF313:ENDIF314:*315:* Scale the right-hand side.316:*317:IF( RCEQU )THEN318:DO30 J = 1, NRHS 319:DO20 I = 1, N 320:B( I, J ) =S( I )*B( I, J ) 321: 20CONTINUE322: 30CONTINUE323:ENDIF324:*325:IF( NOFACT .OR. EQUIL )THEN326:*327:* Compute the Cholesky factorization A = U'*U or A = L*L'.328:*329:CALLSCOPY( N*( N+1 ) / 2, AP, 1, AFP, 1 ) 330:CALLSPPTRF( UPLO, N, AFP, INFO ) 331:*332:* Return if INFO is non-zero.333:*334:IF( INFO.GT.0 )THEN335: RCOND = ZERO 336:RETURN337:ENDIF338:ENDIF339:*340:* Compute the norm of the matrix A.341:*342: ANORM =SLANSP( 'I', UPLO, N, AP, WORK ) 343:*344:* Compute the reciprocal of the condition number of A.345:*346:CALLSPPCON( UPLO, N, AFP, ANORM, RCOND, WORK, IWORK, INFO ) 347:*348:* Compute the solution matrix X.349:*350:CALLSLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 351:CALLSPPTRS( UPLO, N, NRHS, AFP, X, LDX, INFO ) 352:*353:* Use iterative refinement to improve the computed solution and354:* compute error bounds and backward error estimates for it.355:*356:CALLSPPRFS( UPLO, N, NRHS, AP, AFP, B, LDB, X, LDX, FERR, BERR, 357: $ WORK, IWORK, INFO ) 358:*359:* Transform the solution matrix X to a solution of the original360:* system.361:*362:IF( RCEQU )THEN363:DO50 J = 1, NRHS 364:DO40 I = 1, N 365:X( I, J ) =S( I )*X( I, J ) 366: 40CONTINUE367: 50CONTINUE368:DO60 J = 1, NRHS 369:FERR( J ) =FERR( J ) / SCOND 370: 60CONTINUE371:ENDIF372:*373:* Set INFO = N+1 if the matrix is singular to working precision.374:*375:IF( RCOND.LT.SLAMCH( 'Epsilon' ) ) 376: $ INFO = N + 1 377:*378:RETURN379:*380:* End of SPPSVX381:*382:END383: