001:SUBROUTINESSYSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, 002: $ LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK, 003: $ 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 FACT, UPLO 012: INTEGER INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS 013: REAL RCOND 014:* ..015:* .. Array Arguments ..016: INTEGERIPIV( * ),IWORK( * ) 017: REALA( LDA, * ),AF( LDAF, * ),B( LDB, * ), 018: $BERR( * ),FERR( * ),WORK( * ),X( LDX, * ) 019:* ..020:*021:* Purpose022:* =======023:*024:* SSYSVX uses the diagonal pivoting factorization to compute the025:* solution to a real system of linear equations A * X = B,026:* where A is an N-by-N symmetric matrix and X and B are N-by-NRHS027:* 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 = 'N', the diagonal pivoting method is used to factor A.038:* The form of the factorization is039:* A = U * D * U**T, if UPLO = 'U', or040:* A = L * D * L**T, if UPLO = 'L',041:* where U (or L) is a product of permutation and unit upper (lower)042:* triangular matrices, and D is symmetric and block diagonal with043:* 1-by-1 and 2-by-2 diagonal blocks.044:*045:* 2. If some D(i,i)=0, so that D is exactly singular, then the routine046:* returns with INFO = i. Otherwise, the factored form of A is used047:* to estimate the condition number of the matrix A. If the048:* reciprocal of the condition number is less than machine precision,049:* INFO = N+1 is returned as a warning, but the routine still goes on050:* to solve for X and compute error bounds as described below.051:*052:* 3. The system of equations is solved for X using the factored form053:* of A.054:*055:* 4. Iterative refinement is applied to improve the computed solution056:* matrix and calculate error bounds and backward error estimates057:* for it.058:*059:* Arguments060:* =========061:*062:* FACT (input) CHARACTER*1063:* Specifies whether or not the factored form of A has been064:* supplied on entry.065:* = 'F': On entry, AF and IPIV contain the factored form of066:* A. AF and IPIV will not be modified.067:* = 'N': The matrix A will be copied to AF and factored.068:*069:* UPLO (input) CHARACTER*1070:* = 'U': Upper triangle of A is stored;071:* = 'L': Lower triangle of A is stored.072:*073:* N (input) INTEGER074:* The number of linear equations, i.e., the order of the075:* matrix A. N >= 0.076:*077:* NRHS (input) INTEGER078:* The number of right hand sides, i.e., the number of columns079:* of the matrices B and X. NRHS >= 0.080:*081:* A (input) REAL array, dimension (LDA,N)082:* The symmetric matrix A. If UPLO = 'U', the leading N-by-N083:* upper triangular part of A contains the upper triangular part084:* of the matrix A, and the strictly lower triangular part of A085:* is not referenced. If UPLO = 'L', the leading N-by-N lower086:* triangular part of A contains the lower triangular part of087:* the matrix A, and the strictly upper triangular part of A is088:* not referenced.089:*090:* LDA (input) INTEGER091:* The leading dimension of the array A. LDA >= max(1,N).092:*093:* AF (input or output) REAL array, dimension (LDAF,N)094:* If FACT = 'F', then AF is an input argument and on entry095:* contains the block diagonal matrix D and the multipliers used096:* to obtain the factor U or L from the factorization097:* A = U*D*U**T or A = L*D*L**T as computed by SSYTRF.098:*099:* If FACT = 'N', then AF is an output argument and on exit100:* returns the block diagonal matrix D and the multipliers used101:* to obtain the factor U or L from the factorization102:* A = U*D*U**T or A = L*D*L**T.103:*104:* LDAF (input) INTEGER105:* The leading dimension of the array AF. LDAF >= max(1,N).106:*107:* IPIV (input or output) INTEGER array, dimension (N)108:* If FACT = 'F', then IPIV is an input argument and on entry109:* contains details of the interchanges and the block structure110:* of D, as determined by SSYTRF.111:* If IPIV(k) > 0, then rows and columns k and IPIV(k) were112:* interchanged and D(k,k) is a 1-by-1 diagonal block.113:* If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and114:* columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)115:* is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =116:* IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were117:* interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.118:*119:* If FACT = 'N', then IPIV is an output argument and on exit120:* contains details of the interchanges and the block structure121:* of D, as determined by SSYTRF.122:*123:* B (input) REAL array, dimension (LDB,NRHS)124:* The N-by-NRHS right hand side matrix B.125:*126:* LDB (input) INTEGER127:* The leading dimension of the array B. LDB >= max(1,N).128:*129:* X (output) REAL array, dimension (LDX,NRHS)130:* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.131:*132:* LDX (input) INTEGER133:* The leading dimension of the array X. LDX >= max(1,N).134:*135:* RCOND (output) REAL136:* The estimate of the reciprocal condition number of the matrix137:* A. If RCOND is less than the machine precision (in138:* particular, if RCOND = 0), the matrix is singular to working139:* precision. This condition is indicated by a return code of140:* INFO > 0.141:*142:* FERR (output) REAL array, dimension (NRHS)143:* The estimated forward error bound for each solution vector144:* X(j) (the j-th column of the solution matrix X).145:* If XTRUE is the true solution corresponding to X(j), FERR(j)146:* is an estimated upper bound for the magnitude of the largest147:* element in (X(j) - XTRUE) divided by the magnitude of the148:* largest element in X(j). The estimate is as reliable as149:* the estimate for RCOND, and is almost always a slight150:* overestimate of the true error.151:*152:* BERR (output) REAL array, dimension (NRHS)153:* The componentwise relative backward error of each solution154:* vector X(j) (i.e., the smallest relative change in155:* any element of A or B that makes X(j) an exact solution).156:*157:* WORK (workspace/output) REAL array, dimension (MAX(1,LWORK))158:* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.159:*160:* LWORK (input) INTEGER161:* The length of WORK. LWORK >= max(1,3*N), and for best162:* performance, when FACT = 'N', LWORK >= max(1,3*N,N*NB), where163:* NB is the optimal blocksize for SSYTRF.164:*165:* If LWORK = -1, then a workspace query is assumed; the routine166:* only calculates the optimal size of the WORK array, returns167:* this value as the first entry of the WORK array, and no error168:* message related to LWORK is issued by XERBLA.169:*170:* IWORK (workspace) INTEGER array, dimension (N)171:*172:* INFO (output) INTEGER173:* = 0: successful exit174:* < 0: if INFO = -i, the i-th argument had an illegal value175:* > 0: if INFO = i, and i is176:* <= N: D(i,i) is exactly zero. The factorization177:* has been completed but the factor D is exactly178:* singular, so the solution and error bounds could179:* not be computed. RCOND = 0 is returned.180:* = N+1: D is nonsingular, but RCOND is less than machine181:* precision, meaning that the matrix is singular182:* to working precision. Nevertheless, the183:* solution and error bounds are computed because184:* there are a number of situations where the185:* computed solution can be more accurate than the186:* value of RCOND would suggest.187:*188:* =====================================================================189:*190:* .. Parameters ..191: REAL ZERO 192:PARAMETER( ZERO = 0.0E+0 ) 193:* ..194:* .. Local Scalars ..195:LOGICALLQUERY, NOFACT 196: INTEGER LWKOPT, NB 197: REAL ANORM 198:* ..199:* .. External Functions ..200:LOGICALLSAME 201: INTEGER ILAENV 202: REAL SLAMCH, SLANSY 203:EXTERNALILAENV, LSAME, SLAMCH, SLANSY 204:* ..205:* .. External Subroutines ..206:EXTERNALSLACPY, SSYCON, SSYRFS, SSYTRF, SSYTRS, XERBLA 207:* ..208:* .. Intrinsic Functions ..209:INTRINSICMAX 210:* ..211:* .. Executable Statements ..212:*213:* Test the input parameters.214:*215: INFO = 0 216: NOFACT =LSAME( FACT, 'N' ) 217: LQUERY = ( LWORK.EQ.-1 ) 218:IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) )THEN219: INFO = -1 220:ELSEIF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) 221: $THEN222: INFO = -2 223:ELSEIF( N.LT.0 )THEN224: INFO = -3 225:ELSEIF( NRHS.LT.0 )THEN226: INFO = -4 227:ELSEIF( LDA.LT.MAX( 1, N ) )THEN228: INFO = -6 229:ELSEIF( LDAF.LT.MAX( 1, N ) )THEN230: INFO = -8 231:ELSEIF( LDB.LT.MAX( 1, N ) )THEN232: INFO = -11 233:ELSEIF( LDX.LT.MAX( 1, N ) )THEN234: INFO = -13 235:ELSEIF( LWORK.LT.MAX( 1, 3*N ) .AND. .NOT.LQUERY )THEN236: INFO = -18 237:ENDIF238:*239:IF( INFO.EQ.0 )THEN240: LWKOPT =MAX( 1, 3*N ) 241:IF( NOFACT )THEN242: NB =ILAENV( 1, 'SSYTRF', UPLO, N, -1, -1, -1 ) 243: LWKOPT =MAX( LWKOPT, N*NB ) 244:ENDIF245:WORK( 1 ) = LWKOPT 246:ENDIF247:*248:IF( INFO.NE.0 )THEN249:CALLXERBLA( 'SSYSVX', -INFO ) 250:RETURN251:ELSEIF( LQUERY )THEN252:RETURN253:ENDIF254:*255:IF( NOFACT )THEN256:*257:* Compute the factorization A = U*D*U' or A = L*D*L'.258:*259:CALLSLACPY( UPLO, N, N, A, LDA, AF, LDAF ) 260:CALLSSYTRF( UPLO, N, AF, LDAF, IPIV, WORK, LWORK, INFO ) 261:*262:* Return if INFO is non-zero.263:*264:IF( INFO.GT.0 )THEN265: RCOND = ZERO 266:RETURN267:ENDIF268:ENDIF269:*270:* Compute the norm of the matrix A.271:*272: ANORM =SLANSY( 'I', UPLO, N, A, LDA, WORK ) 273:*274:* Compute the reciprocal of the condition number of A.275:*276:CALLSSYCON( UPLO, N, AF, LDAF, IPIV, ANORM, RCOND, WORK, IWORK, 277: $ INFO ) 278:*279:* Compute the solution vectors X.280:*281:CALLSLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 282:CALLSSYTRS( UPLO, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO ) 283:*284:* Use iterative refinement to improve the computed solutions and285:* compute error bounds and backward error estimates for them.286:*287:CALLSSYRFS( UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, 288: $ LDX, FERR, BERR, WORK, IWORK, INFO ) 289:*290:* Set INFO = N+1 if the matrix is singular to working precision.291:*292:IF( RCOND.LT.SLAMCH( 'Epsilon' ) ) 293: $ INFO = N + 1 294:*295:WORK( 1 ) = LWKOPT 296:*297:RETURN298:*299:* End of SSYSVX300:*301:END302: