SUBROUTINE DSYSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, $ LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK, $ IWORK, INFO ) * * -- LAPACK driver routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. CHARACTER FACT, UPLO INTEGER INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS DOUBLE PRECISION RCOND * .. * .. Array Arguments .. INTEGER IPIV( * ), IWORK( * ) DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ), $ BERR( * ), FERR( * ), WORK( * ), X( LDX, * ) * .. * * Purpose * ======= * * DSYSVX uses the diagonal pivoting factorization to compute the * solution to a real system of linear equations A * X = B, * where A is an N-by-N symmetric matrix and X and B are N-by-NRHS * matrices. * * Error bounds on the solution and a condition estimate are also * provided. * * Description * =========== * * The following steps are performed: * * 1. If FACT = 'N', the diagonal pivoting method is used to factor A. * The form of the factorization is * A = U * D * U**T, if UPLO = 'U', or * A = L * D * L**T, if UPLO = 'L', * where U (or L) is a product of permutation and unit upper (lower) * triangular matrices, and D is symmetric and block diagonal with * 1-by-1 and 2-by-2 diagonal blocks. * * 2. If some D(i,i)=0, so that D is exactly singular, then the routine * returns with INFO = i. Otherwise, the factored form of A is used * to estimate the condition number of the matrix A. If the * reciprocal of the condition number is less than machine precision, * INFO = N+1 is returned as a warning, but the routine still goes on * to solve for X and compute error bounds as described below. * * 3. The system of equations is solved for X using the factored form * of A. * * 4. Iterative refinement is applied to improve the computed solution * matrix and calculate error bounds and backward error estimates * for it. * * Arguments * ========= * * FACT (input) CHARACTER*1 * Specifies whether or not the factored form of A has been * supplied on entry. * = 'F': On entry, AF and IPIV contain the factored form of * A. AF and IPIV will not be modified. * = 'N': The matrix A will be copied to AF and factored. * * UPLO (input) CHARACTER*1 * = 'U': Upper triangle of A is stored; * = 'L': Lower triangle of A is stored. * * N (input) INTEGER * The number of linear equations, i.e., the order of the * matrix A. N >= 0. * * NRHS (input) INTEGER * The number of right hand sides, i.e., the number of columns * of the matrices B and X. NRHS >= 0. * * A (input) DOUBLE PRECISION array, dimension (LDA,N) * The symmetric matrix A. If UPLO = 'U', the leading N-by-N * upper triangular part of A contains the upper triangular part * of the matrix A, and the strictly lower triangular part of A * is not referenced. If UPLO = 'L', the leading N-by-N lower * triangular part of A contains the lower triangular part of * the matrix A, and the strictly upper triangular part of A is * not referenced. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * AF (input or output) DOUBLE PRECISION array, dimension (LDAF,N) * If FACT = 'F', then AF is an input argument and on entry * contains the block diagonal matrix D and the multipliers used * to obtain the factor U or L from the factorization * A = U*D*U**T or A = L*D*L**T as computed by DSYTRF. * * If FACT = 'N', then AF is an output argument and on exit * returns the block diagonal matrix D and the multipliers used * to obtain the factor U or L from the factorization * A = U*D*U**T or A = L*D*L**T. * * LDAF (input) INTEGER * The leading dimension of the array AF. LDAF >= max(1,N). * * IPIV (input or output) INTEGER array, dimension (N) * If FACT = 'F', then IPIV is an input argument and on entry * contains details of the interchanges and the block structure * of D, as determined by DSYTRF. * If IPIV(k) > 0, then rows and columns k and IPIV(k) were * interchanged and D(k,k) is a 1-by-1 diagonal block. * If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and * columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) * is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) = * IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were * interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block. * * If FACT = 'N', then IPIV is an output argument and on exit * contains details of the interchanges and the block structure * of D, as determined by DSYTRF. * * B (input) DOUBLE PRECISION array, dimension (LDB,NRHS) * The N-by-NRHS right hand side matrix B. * * LDB (input) INTEGER * The leading dimension of the array B. LDB >= max(1,N). * * X (output) DOUBLE PRECISION array, dimension (LDX,NRHS) * If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X. * * LDX (input) INTEGER * The leading dimension of the array X. LDX >= max(1,N). * * RCOND (output) DOUBLE PRECISION * The estimate of the reciprocal condition number of the matrix * A. If RCOND is less than the machine precision (in * particular, if RCOND = 0), the matrix is singular to working * precision. This condition is indicated by a return code of * INFO > 0. * * FERR (output) DOUBLE PRECISION array, dimension (NRHS) * The estimated forward error bound for each solution vector * X(j) (the j-th column of the solution matrix X). * If XTRUE is the true solution corresponding to X(j), FERR(j) * is an estimated upper bound for the magnitude of the largest * element in (X(j) - XTRUE) divided by the magnitude of the * largest element in X(j). The estimate is as reliable as * the estimate for RCOND, and is almost always a slight * overestimate of the true error. * * BERR (output) DOUBLE PRECISION array, dimension (NRHS) * The componentwise relative backward error of each solution * vector X(j) (i.e., the smallest relative change in * any element of A or B that makes X(j) an exact solution). * * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. * * LWORK (input) INTEGER * The length of WORK. LWORK >= max(1,3*N), and for best * performance, when FACT = 'N', LWORK >= max(1,3*N,N*NB), where * NB is the optimal blocksize for DSYTRF. * * If LWORK = -1, then a workspace query is assumed; the routine * only calculates the optimal size of the WORK array, returns * this value as the first entry of the WORK array, and no error * message related to LWORK is issued by XERBLA. * * IWORK (workspace) INTEGER array, dimension (N) * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * > 0: if INFO = i, and i is * <= N: D(i,i) is exactly zero. The factorization * has been completed but the factor D is exactly * singular, so the solution and error bounds could * not be computed. RCOND = 0 is returned. * = N+1: D is nonsingular, but RCOND is less than machine * precision, meaning that the matrix is singular * to working precision. Nevertheless, the * solution and error bounds are computed because * there are a number of situations where the * computed solution can be more accurate than the * value of RCOND would suggest. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. * .. Local Scalars .. LOGICAL LQUERY, NOFACT INTEGER LWKOPT, NB DOUBLE PRECISION ANORM * .. * .. External Functions .. LOGICAL LSAME INTEGER ILAENV DOUBLE PRECISION DLAMCH, DLANSY EXTERNAL LSAME, ILAENV, DLAMCH, DLANSY * .. * .. External Subroutines .. EXTERNAL DLACPY, DSYCON, DSYRFS, DSYTRF, DSYTRS, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 NOFACT = LSAME( FACT, 'N' ) LQUERY = ( LWORK.EQ.-1 ) IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN INFO = -1 ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) $ THEN INFO = -2 ELSE IF( N.LT.0 ) THEN INFO = -3 ELSE IF( NRHS.LT.0 ) THEN INFO = -4 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -6 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN INFO = -8 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN INFO = -11 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN INFO = -13 ELSE IF( LWORK.LT.MAX( 1, 3*N ) .AND. .NOT.LQUERY ) THEN INFO = -18 END IF * IF( INFO.EQ.0 ) THEN LWKOPT = MAX( 1, 3*N ) IF( NOFACT ) THEN NB = ILAENV( 1, 'DSYTRF', UPLO, N, -1, -1, -1 ) LWKOPT = MAX( LWKOPT, N*NB ) END IF WORK( 1 ) = LWKOPT END IF * IF( INFO.NE.0 ) THEN CALL XERBLA( 'DSYSVX', -INFO ) RETURN ELSE IF( LQUERY ) THEN RETURN END IF * IF( NOFACT ) THEN * * Compute the factorization A = U*D*U' or A = L*D*L'. * CALL DLACPY( UPLO, N, N, A, LDA, AF, LDAF ) CALL DSYTRF( UPLO, N, AF, LDAF, IPIV, WORK, LWORK, INFO ) * * Return if INFO is non-zero. * IF( INFO.GT.0 )THEN RCOND = ZERO RETURN END IF END IF * * Compute the norm of the matrix A. * ANORM = DLANSY( 'I', UPLO, N, A, LDA, WORK ) * * Compute the reciprocal of the condition number of A. * CALL DSYCON( UPLO, N, AF, LDAF, IPIV, ANORM, RCOND, WORK, IWORK, $ INFO ) * * Compute the solution vectors X. * CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) CALL DSYTRS( UPLO, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO ) * * Use iterative refinement to improve the computed solutions and * compute error bounds and backward error estimates for them. * CALL DSYRFS( UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, $ LDX, FERR, BERR, WORK, IWORK, INFO ) * * Set INFO = N+1 if the matrix is singular to working precision. * IF( RCOND.LT.DLAMCH( 'Epsilon' ) ) $ INFO = N + 1 * WORK( 1 ) = LWKOPT * RETURN * * End of DSYSVX * END