001:SUBROUTINEDSYSVX( 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:* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..007:* November 2006008:*009:* .. Scalar Arguments ..010: CHARACTER FACT, UPLO 011: INTEGER INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS 012: DOUBLE PRECISION RCOND 013:* ..014:* .. Array Arguments ..015: INTEGERIPIV( * ),IWORK( * ) 016: DOUBLE PRECISIONA( LDA, * ),AF( LDAF, * ),B( LDB, * ), 017: $BERR( * ),FERR( * ),WORK( * ),X( LDX, * ) 018:* ..019:*020:* Purpose021:* =======022:*023:* DSYSVX uses the diagonal pivoting factorization to compute the024:* solution to a real system of linear equations A * X = B,025:* where A is an N-by-N symmetric matrix and X and B are N-by-NRHS026:* matrices.027:*028:* Error bounds on the solution and a condition estimate are also029:* provided.030:*031:* Description032:* ===========033:*034:* The following steps are performed:035:*036:* 1. If FACT = 'N', the diagonal pivoting method is used to factor A.037:* The form of the factorization is038:* A = U * D * U**T, if UPLO = 'U', or039:* A = L * D * L**T, if UPLO = 'L',040:* where U (or L) is a product of permutation and unit upper (lower)041:* triangular matrices, and D is symmetric and block diagonal with042:* 1-by-1 and 2-by-2 diagonal blocks.043:*044:* 2. If some D(i,i)=0, so that D is exactly singular, then the routine045:* returns with INFO = i. Otherwise, the factored form of A is used046:* to estimate the condition number of the matrix A. If the047:* reciprocal of the condition number is less than machine precision,048:* INFO = N+1 is returned as a warning, but the routine still goes on049:* to solve for X and compute error bounds as described below.050:*051:* 3. The system of equations is solved for X using the factored form052:* of A.053:*054:* 4. Iterative refinement is applied to improve the computed solution055:* matrix and calculate error bounds and backward error estimates056:* for it.057:*058:* Arguments059:* =========060:*061:* FACT (input) CHARACTER*1062:* Specifies whether or not the factored form of A has been063:* supplied on entry.064:* = 'F': On entry, AF and IPIV contain the factored form of065:* A. AF and IPIV will not be modified.066:* = 'N': The matrix A will be copied to AF and factored.067:*068:* UPLO (input) CHARACTER*1069:* = 'U': Upper triangle of A is stored;070:* = 'L': Lower triangle of A is stored.071:*072:* N (input) INTEGER073:* The number of linear equations, i.e., the order of the074:* matrix A. N >= 0.075:*076:* NRHS (input) INTEGER077:* The number of right hand sides, i.e., the number of columns078:* of the matrices B and X. NRHS >= 0.079:*080:* A (input) DOUBLE PRECISION array, dimension (LDA,N)081:* The symmetric matrix A. If UPLO = 'U', the leading N-by-N082:* upper triangular part of A contains the upper triangular part083:* of the matrix A, and the strictly lower triangular part of A084:* is not referenced. If UPLO = 'L', the leading N-by-N lower085:* triangular part of A contains the lower triangular part of086:* the matrix A, and the strictly upper triangular part of A is087:* not referenced.088:*089:* LDA (input) INTEGER090:* The leading dimension of the array A. LDA >= max(1,N).091:*092:* AF (input or output) DOUBLE PRECISION array, dimension (LDAF,N)093:* If FACT = 'F', then AF is an input argument and on entry094:* contains the block diagonal matrix D and the multipliers used095:* to obtain the factor U or L from the factorization096:* A = U*D*U**T or A = L*D*L**T as computed by DSYTRF.097:*098:* If FACT = 'N', then AF is an output argument and on exit099:* returns the block diagonal matrix D and the multipliers used100:* to obtain the factor U or L from the factorization101:* A = U*D*U**T or A = L*D*L**T.102:*103:* LDAF (input) INTEGER104:* The leading dimension of the array AF. LDAF >= max(1,N).105:*106:* IPIV (input or output) INTEGER array, dimension (N)107:* If FACT = 'F', then IPIV is an input argument and on entry108:* contains details of the interchanges and the block structure109:* of D, as determined by DSYTRF.110:* If IPIV(k) > 0, then rows and columns k and IPIV(k) were111:* interchanged and D(k,k) is a 1-by-1 diagonal block.112:* If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and113:* columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)114:* is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =115:* IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were116:* interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.117:*118:* If FACT = 'N', then IPIV is an output argument and on exit119:* contains details of the interchanges and the block structure120:* of D, as determined by DSYTRF.121:*122:* B (input) DOUBLE PRECISION array, dimension (LDB,NRHS)123:* The N-by-NRHS right hand side matrix B.124:*125:* LDB (input) INTEGER126:* The leading dimension of the array B. LDB >= max(1,N).127:*128:* X (output) DOUBLE PRECISION array, dimension (LDX,NRHS)129:* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.130:*131:* LDX (input) INTEGER132:* The leading dimension of the array X. LDX >= max(1,N).133:*134:* RCOND (output) DOUBLE PRECISION135:* The estimate of the reciprocal condition number of the matrix136:* A. If RCOND is less than the machine precision (in137:* particular, if RCOND = 0), the matrix is singular to working138:* precision. This condition is indicated by a return code of139:* INFO > 0.140:*141:* FERR (output) DOUBLE PRECISION array, dimension (NRHS)142:* The estimated forward error bound for each solution vector143:* X(j) (the j-th column of the solution matrix X).144:* If XTRUE is the true solution corresponding to X(j), FERR(j)145:* is an estimated upper bound for the magnitude of the largest146:* element in (X(j) - XTRUE) divided by the magnitude of the147:* largest element in X(j). The estimate is as reliable as148:* the estimate for RCOND, and is almost always a slight149:* overestimate of the true error.150:*151:* BERR (output) DOUBLE PRECISION array, dimension (NRHS)152:* The componentwise relative backward error of each solution153:* vector X(j) (i.e., the smallest relative change in154:* any element of A or B that makes X(j) an exact solution).155:*156:* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))157:* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.158:*159:* LWORK (input) INTEGER160:* The length of WORK. LWORK >= max(1,3*N), and for best161:* performance, when FACT = 'N', LWORK >= max(1,3*N,N*NB), where162:* NB is the optimal blocksize for DSYTRF.163:*164:* If LWORK = -1, then a workspace query is assumed; the routine165:* only calculates the optimal size of the WORK array, returns166:* this value as the first entry of the WORK array, and no error167:* message related to LWORK is issued by XERBLA.168:*169:* IWORK (workspace) INTEGER array, dimension (N)170:*171:* INFO (output) INTEGER172:* = 0: successful exit173:* < 0: if INFO = -i, the i-th argument had an illegal value174:* > 0: if INFO = i, and i is175:* <= N: D(i,i) is exactly zero. The factorization176:* has been completed but the factor D is exactly177:* singular, so the solution and error bounds could178:* not be computed. RCOND = 0 is returned.179:* = N+1: D is nonsingular, but RCOND is less than machine180:* precision, meaning that the matrix is singular181:* to working precision. Nevertheless, the182:* solution and error bounds are computed because183:* there are a number of situations where the184:* computed solution can be more accurate than the185:* value of RCOND would suggest.186:*187:* =====================================================================188:*189:* .. Parameters ..190: DOUBLE PRECISION ZERO 191:PARAMETER( ZERO = 0.0D+0 ) 192:* ..193:* .. Local Scalars ..194:LOGICALLQUERY, NOFACT 195: INTEGER LWKOPT, NB 196: DOUBLE PRECISION ANORM 197:* ..198:* .. External Functions ..199:LOGICALLSAME 200: INTEGER ILAENV 201: DOUBLE PRECISION DLAMCH, DLANSY 202:EXTERNALLSAME, ILAENV, DLAMCH, DLANSY 203:* ..204:* .. External Subroutines ..205:EXTERNALDLACPY, DSYCON, DSYRFS, DSYTRF, DSYTRS, XERBLA 206:* ..207:* .. Intrinsic Functions ..208:INTRINSICMAX 209:* ..210:* .. Executable Statements ..211:*212:* Test the input parameters.213:*214: INFO = 0 215: NOFACT =LSAME( FACT, 'N' ) 216: LQUERY = ( LWORK.EQ.-1 ) 217:IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) )THEN218: INFO = -1 219:ELSEIF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) 220: $THEN221: INFO = -2 222:ELSEIF( N.LT.0 )THEN223: INFO = -3 224:ELSEIF( NRHS.LT.0 )THEN225: INFO = -4 226:ELSEIF( LDA.LT.MAX( 1, N ) )THEN227: INFO = -6 228:ELSEIF( LDAF.LT.MAX( 1, N ) )THEN229: INFO = -8 230:ELSEIF( LDB.LT.MAX( 1, N ) )THEN231: INFO = -11 232:ELSEIF( LDX.LT.MAX( 1, N ) )THEN233: INFO = -13 234:ELSEIF( LWORK.LT.MAX( 1, 3*N ) .AND. .NOT.LQUERY )THEN235: INFO = -18 236:ENDIF237:*238:IF( INFO.EQ.0 )THEN239: LWKOPT =MAX( 1, 3*N ) 240:IF( NOFACT )THEN241: NB =ILAENV( 1, 'DSYTRF', UPLO, N, -1, -1, -1 ) 242: LWKOPT =MAX( LWKOPT, N*NB ) 243:ENDIF244:WORK( 1 ) = LWKOPT 245:ENDIF246:*247:IF( INFO.NE.0 )THEN248:CALLXERBLA( 'DSYSVX', -INFO ) 249:RETURN250:ELSEIF( LQUERY )THEN251:RETURN252:ENDIF253:*254:IF( NOFACT )THEN255:*256:* Compute the factorization A = U*D*U' or A = L*D*L'.257:*258:CALLDLACPY( UPLO, N, N, A, LDA, AF, LDAF ) 259:CALLDSYTRF( UPLO, N, AF, LDAF, IPIV, WORK, LWORK, INFO ) 260:*261:* Return if INFO is non-zero.262:*263:IF( INFO.GT.0 )THEN264: RCOND = ZERO 265:RETURN266:ENDIF267:ENDIF268:*269:* Compute the norm of the matrix A.270:*271: ANORM =DLANSY( 'I', UPLO, N, A, LDA, WORK ) 272:*273:* Compute the reciprocal of the condition number of A.274:*275:CALLDSYCON( UPLO, N, AF, LDAF, IPIV, ANORM, RCOND, WORK, IWORK, 276: $ INFO ) 277:*278:* Compute the solution vectors X.279:*280:CALLDLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 281:CALLDSYTRS( UPLO, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO ) 282:*283:* Use iterative refinement to improve the computed solutions and284:* compute error bounds and backward error estimates for them.285:*286:CALLDSYRFS( UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, 287: $ LDX, FERR, BERR, WORK, IWORK, INFO ) 288:*289:* Set INFO = N+1 if the matrix is singular to working precision.290:*291:IF( RCOND.LT.DLAMCH( 'Epsilon' ) ) 292: $ INFO = N + 1 293:*294:WORK( 1 ) = LWKOPT 295:*296:RETURN297:*298:* End of DSYSVX299:*300:END301: