 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DSYSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, 00002 $ LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK, 00003 $ IWORK, INFO ) 00004 * 00005 * -- LAPACK driver routine (version 3.3.1) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * -- April 2011 -- 00009 * 00010 * .. Scalar Arguments .. 00011 CHARACTER FACT, UPLO 00012 INTEGER INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS 00013 DOUBLE PRECISION RCOND 00014 * .. 00015 * .. Array Arguments .. 00016 INTEGER IPIV( * ), IWORK( * ) 00017 DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 00018 $ BERR( * ), FERR( * ), WORK( * ), X( LDX, * ) 00019 * .. 00020 * 00021 * Purpose 00022 * ======= 00023 * 00024 * DSYSVX uses the diagonal pivoting factorization to compute the 00025 * solution to a real system of linear equations A * X = B, 00026 * where A is an N-by-N symmetric matrix and X and B are N-by-NRHS 00027 * matrices. 00028 * 00029 * Error bounds on the solution and a condition estimate are also 00030 * provided. 00031 * 00032 * Description 00033 * =========== 00034 * 00035 * The following steps are performed: 00036 * 00037 * 1. If FACT = 'N', the diagonal pivoting method is used to factor A. 00038 * The form of the factorization is 00039 * A = U * D * U**T, if UPLO = 'U', or 00040 * A = L * D * L**T, if UPLO = 'L', 00041 * where U (or L) is a product of permutation and unit upper (lower) 00042 * triangular matrices, and D is symmetric and block diagonal with 00043 * 1-by-1 and 2-by-2 diagonal blocks. 00044 * 00045 * 2. If some D(i,i)=0, so that D is exactly singular, then the routine 00046 * returns with INFO = i. Otherwise, the factored form of A is used 00047 * to estimate the condition number of the matrix A. If the 00048 * reciprocal of the condition number is less than machine precision, 00049 * INFO = N+1 is returned as a warning, but the routine still goes on 00050 * to solve for X and compute error bounds as described below. 00051 * 00052 * 3. The system of equations is solved for X using the factored form 00053 * of A. 00054 * 00055 * 4. Iterative refinement is applied to improve the computed solution 00056 * matrix and calculate error bounds and backward error estimates 00057 * for it. 00058 * 00059 * Arguments 00060 * ========= 00061 * 00062 * FACT (input) CHARACTER*1 00063 * Specifies whether or not the factored form of A has been 00064 * supplied on entry. 00065 * = 'F': On entry, AF and IPIV contain the factored form of 00066 * A. AF and IPIV will not be modified. 00067 * = 'N': The matrix A will be copied to AF and factored. 00068 * 00069 * UPLO (input) CHARACTER*1 00070 * = 'U': Upper triangle of A is stored; 00071 * = 'L': Lower triangle of A is stored. 00072 * 00073 * N (input) INTEGER 00074 * The number of linear equations, i.e., the order of the 00075 * matrix A. N >= 0. 00076 * 00077 * NRHS (input) INTEGER 00078 * The number of right hand sides, i.e., the number of columns 00079 * of the matrices B and X. NRHS >= 0. 00080 * 00081 * A (input) DOUBLE PRECISION array, dimension (LDA,N) 00082 * The symmetric matrix A. If UPLO = 'U', the leading N-by-N 00083 * upper triangular part of A contains the upper triangular part 00084 * of the matrix A, and the strictly lower triangular part of A 00085 * is not referenced. If UPLO = 'L', the leading N-by-N lower 00086 * triangular part of A contains the lower triangular part of 00087 * the matrix A, and the strictly upper triangular part of A is 00088 * not referenced. 00089 * 00090 * LDA (input) INTEGER 00091 * The leading dimension of the array A. LDA >= max(1,N). 00092 * 00093 * AF (input or output) DOUBLE PRECISION array, dimension (LDAF,N) 00094 * If FACT = 'F', then AF is an input argument and on entry 00095 * contains the block diagonal matrix D and the multipliers used 00096 * to obtain the factor U or L from the factorization 00097 * A = U*D*U**T or A = L*D*L**T as computed by DSYTRF. 00098 * 00099 * If FACT = 'N', then AF is an output argument and on exit 00100 * returns the block diagonal matrix D and the multipliers used 00101 * to obtain the factor U or L from the factorization 00102 * A = U*D*U**T or A = L*D*L**T. 00103 * 00104 * LDAF (input) INTEGER 00105 * The leading dimension of the array AF. LDAF >= max(1,N). 00106 * 00107 * IPIV (input or output) INTEGER array, dimension (N) 00108 * If FACT = 'F', then IPIV is an input argument and on entry 00109 * contains details of the interchanges and the block structure 00110 * of D, as determined by DSYTRF. 00111 * If IPIV(k) > 0, then rows and columns k and IPIV(k) were 00112 * interchanged and D(k,k) is a 1-by-1 diagonal block. 00113 * If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and 00114 * columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) 00115 * is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) = 00116 * IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were 00117 * interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block. 00118 * 00119 * If FACT = 'N', then IPIV is an output argument and on exit 00120 * contains details of the interchanges and the block structure 00121 * of D, as determined by DSYTRF. 00122 * 00123 * B (input) DOUBLE PRECISION array, dimension (LDB,NRHS) 00124 * The N-by-NRHS right hand side matrix B. 00125 * 00126 * LDB (input) INTEGER 00127 * The leading dimension of the array B. LDB >= max(1,N). 00128 * 00129 * X (output) DOUBLE PRECISION array, dimension (LDX,NRHS) 00130 * If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X. 00131 * 00132 * LDX (input) INTEGER 00133 * The leading dimension of the array X. LDX >= max(1,N). 00134 * 00135 * RCOND (output) DOUBLE PRECISION 00136 * The estimate of the reciprocal condition number of the matrix 00137 * A. If RCOND is less than the machine precision (in 00138 * particular, if RCOND = 0), the matrix is singular to working 00139 * precision. This condition is indicated by a return code of 00140 * INFO > 0. 00141 * 00142 * FERR (output) DOUBLE PRECISION array, dimension (NRHS) 00143 * The estimated forward error bound for each solution vector 00144 * X(j) (the j-th column of the solution matrix X). 00145 * If XTRUE is the true solution corresponding to X(j), FERR(j) 00146 * is an estimated upper bound for the magnitude of the largest 00147 * element in (X(j) - XTRUE) divided by the magnitude of the 00148 * largest element in X(j). The estimate is as reliable as 00149 * the estimate for RCOND, and is almost always a slight 00150 * overestimate of the true error. 00151 * 00152 * BERR (output) DOUBLE PRECISION array, dimension (NRHS) 00153 * The componentwise relative backward error of each solution 00154 * vector X(j) (i.e., the smallest relative change in 00155 * any element of A or B that makes X(j) an exact solution). 00156 * 00157 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 00158 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00159 * 00160 * LWORK (input) INTEGER 00161 * The length of WORK. LWORK >= max(1,3*N), and for best 00162 * performance, when FACT = 'N', LWORK >= max(1,3*N,N*NB), where 00163 * NB is the optimal blocksize for DSYTRF. 00164 * 00165 * If LWORK = -1, then a workspace query is assumed; the routine 00166 * only calculates the optimal size of the WORK array, returns 00167 * this value as the first entry of the WORK array, and no error 00168 * message related to LWORK is issued by XERBLA. 00169 * 00170 * IWORK (workspace) INTEGER array, dimension (N) 00171 * 00172 * INFO (output) INTEGER 00173 * = 0: successful exit 00174 * < 0: if INFO = -i, the i-th argument had an illegal value 00175 * > 0: if INFO = i, and i is 00176 * <= N: D(i,i) is exactly zero. The factorization 00177 * has been completed but the factor D is exactly 00178 * singular, so the solution and error bounds could 00179 * not be computed. RCOND = 0 is returned. 00180 * = N+1: D is nonsingular, but RCOND is less than machine 00181 * precision, meaning that the matrix is singular 00182 * to working precision. Nevertheless, the 00183 * solution and error bounds are computed because 00184 * there are a number of situations where the 00185 * computed solution can be more accurate than the 00186 * value of RCOND would suggest. 00187 * 00188 * ===================================================================== 00189 * 00190 * .. Parameters .. 00191 DOUBLE PRECISION ZERO 00192 PARAMETER ( ZERO = 0.0D+0 ) 00193 * .. 00194 * .. Local Scalars .. 00195 LOGICAL LQUERY, NOFACT 00196 INTEGER LWKOPT, NB 00197 DOUBLE PRECISION ANORM 00198 * .. 00199 * .. External Functions .. 00200 LOGICAL LSAME 00201 INTEGER ILAENV 00202 DOUBLE PRECISION DLAMCH, DLANSY 00203 EXTERNAL LSAME, ILAENV, DLAMCH, DLANSY 00204 * .. 00205 * .. External Subroutines .. 00206 EXTERNAL DLACPY, DSYCON, DSYRFS, DSYTRF, DSYTRS, XERBLA 00207 * .. 00208 * .. Intrinsic Functions .. 00209 INTRINSIC MAX 00210 * .. 00211 * .. Executable Statements .. 00212 * 00213 * Test the input parameters. 00214 * 00215 INFO = 0 00216 NOFACT = LSAME( FACT, 'N' ) 00217 LQUERY = ( LWORK.EQ.-1 ) 00218 IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN 00219 INFO = -1 00220 ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) 00221 $ THEN 00222 INFO = -2 00223 ELSE IF( N.LT.0 ) THEN 00224 INFO = -3 00225 ELSE IF( NRHS.LT.0 ) THEN 00226 INFO = -4 00227 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00228 INFO = -6 00229 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN 00230 INFO = -8 00231 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00232 INFO = -11 00233 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 00234 INFO = -13 00235 ELSE IF( LWORK.LT.MAX( 1, 3*N ) .AND. .NOT.LQUERY ) THEN 00236 INFO = -18 00237 END IF 00238 * 00239 IF( INFO.EQ.0 ) THEN 00240 LWKOPT = MAX( 1, 3*N ) 00241 IF( NOFACT ) THEN 00242 NB = ILAENV( 1, 'DSYTRF', UPLO, N, -1, -1, -1 ) 00243 LWKOPT = MAX( LWKOPT, N*NB ) 00244 END IF 00245 WORK( 1 ) = LWKOPT 00246 END IF 00247 * 00248 IF( INFO.NE.0 ) THEN 00249 CALL XERBLA( 'DSYSVX', -INFO ) 00250 RETURN 00251 ELSE IF( LQUERY ) THEN 00252 RETURN 00253 END IF 00254 * 00255 IF( NOFACT ) THEN 00256 * 00257 * Compute the factorization A = U*D*U**T or A = L*D*L**T. 00258 * 00259 CALL DLACPY( UPLO, N, N, A, LDA, AF, LDAF ) 00260 CALL DSYTRF( UPLO, N, AF, LDAF, IPIV, WORK, LWORK, INFO ) 00261 * 00262 * Return if INFO is non-zero. 00263 * 00264 IF( INFO.GT.0 )THEN 00265 RCOND = ZERO 00266 RETURN 00267 END IF 00268 END IF 00269 * 00270 * Compute the norm of the matrix A. 00271 * 00272 ANORM = DLANSY( 'I', UPLO, N, A, LDA, WORK ) 00273 * 00274 * Compute the reciprocal of the condition number of A. 00275 * 00276 CALL DSYCON( UPLO, N, AF, LDAF, IPIV, ANORM, RCOND, WORK, IWORK, 00277 $ INFO ) 00278 * 00279 * Compute the solution vectors X. 00280 * 00281 CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 00282 CALL DSYTRS( UPLO, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO ) 00283 * 00284 * Use iterative refinement to improve the computed solutions and 00285 * compute error bounds and backward error estimates for them. 00286 * 00287 CALL DSYRFS( UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, 00288 $ LDX, FERR, BERR, WORK, IWORK, INFO ) 00289 * 00290 * Set INFO = N+1 if the matrix is singular to working precision. 00291 * 00292 IF( RCOND.LT.DLAMCH( 'Epsilon' ) ) 00293 $ INFO = N + 1 00294 * 00295 WORK( 1 ) = LWKOPT 00296 * 00297 RETURN 00298 * 00299 * End of DSYSVX 00300 * 00301 END
1.7.3 
