1 SUBROUTINE DSYSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B,
2 $ LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK,
3 $ IWORK, INFO )
4 *
5 * -- LAPACK driver routine (version 3.3.1) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * -- April 2011 --
9 *
10 * .. Scalar Arguments ..
11 CHARACTER FACT, UPLO
12 INTEGER INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS
13 DOUBLE PRECISION RCOND
14 * ..
15 * .. Array Arguments ..
16 INTEGER IPIV( * ), IWORK( * )
17 DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
18 $ BERR( * ), FERR( * ), WORK( * ), X( LDX, * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * DSYSVX uses the diagonal pivoting factorization to compute the
25 * solution to a real system of linear equations A * X = B,
26 * where A is an N-by-N symmetric matrix and X and B are N-by-NRHS
27 * matrices.
28 *
29 * Error bounds on the solution and a condition estimate are also
30 * provided.
31 *
32 * Description
33 * ===========
34 *
35 * The following steps are performed:
36 *
37 * 1. If FACT = 'N', the diagonal pivoting method is used to factor A.
38 * The form of the factorization is
39 * A = U * D * U**T, if UPLO = 'U', or
40 * A = L * D * L**T, if UPLO = 'L',
41 * where U (or L) is a product of permutation and unit upper (lower)
42 * triangular matrices, and D is symmetric and block diagonal with
43 * 1-by-1 and 2-by-2 diagonal blocks.
44 *
45 * 2. If some D(i,i)=0, so that D is exactly singular, then the routine
46 * returns with INFO = i. Otherwise, the factored form of A is used
47 * to estimate the condition number of the matrix A. If the
48 * reciprocal of the condition number is less than machine precision,
49 * INFO = N+1 is returned as a warning, but the routine still goes on
50 * to solve for X and compute error bounds as described below.
51 *
52 * 3. The system of equations is solved for X using the factored form
53 * of A.
54 *
55 * 4. Iterative refinement is applied to improve the computed solution
56 * matrix and calculate error bounds and backward error estimates
57 * for it.
58 *
59 * Arguments
60 * =========
61 *
62 * FACT (input) CHARACTER*1
63 * Specifies whether or not the factored form of A has been
64 * supplied on entry.
65 * = 'F': On entry, AF and IPIV contain the factored form of
66 * A. AF and IPIV will not be modified.
67 * = 'N': The matrix A will be copied to AF and factored.
68 *
69 * UPLO (input) CHARACTER*1
70 * = 'U': Upper triangle of A is stored;
71 * = 'L': Lower triangle of A is stored.
72 *
73 * N (input) INTEGER
74 * The number of linear equations, i.e., the order of the
75 * matrix A. N >= 0.
76 *
77 * NRHS (input) INTEGER
78 * The number of right hand sides, i.e., the number of columns
79 * of the matrices B and X. NRHS >= 0.
80 *
81 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
82 * The symmetric matrix A. If UPLO = 'U', the leading N-by-N
83 * upper triangular part of A contains the upper triangular part
84 * of the matrix A, and the strictly lower triangular part of A
85 * is not referenced. If UPLO = 'L', the leading N-by-N lower
86 * triangular part of A contains the lower triangular part of
87 * the matrix A, and the strictly upper triangular part of A is
88 * not referenced.
89 *
90 * LDA (input) INTEGER
91 * The leading dimension of the array A. LDA >= max(1,N).
92 *
93 * AF (input or output) DOUBLE PRECISION array, dimension (LDAF,N)
94 * If FACT = 'F', then AF is an input argument and on entry
95 * contains the block diagonal matrix D and the multipliers used
96 * to obtain the factor U or L from the factorization
97 * A = U*D*U**T or A = L*D*L**T as computed by DSYTRF.
98 *
99 * If FACT = 'N', then AF is an output argument and on exit
100 * returns the block diagonal matrix D and the multipliers used
101 * to obtain the factor U or L from the factorization
102 * A = U*D*U**T or A = L*D*L**T.
103 *
104 * LDAF (input) INTEGER
105 * 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 entry
109 * contains details of the interchanges and the block structure
110 * of D, as determined by DSYTRF.
111 * If IPIV(k) > 0, then rows and columns k and IPIV(k) were
112 * 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 and
114 * 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) were
117 * 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 exit
120 * contains details of the interchanges and the block structure
121 * of D, as determined by DSYTRF.
122 *
123 * B (input) DOUBLE PRECISION array, dimension (LDB,NRHS)
124 * The N-by-NRHS right hand side matrix B.
125 *
126 * LDB (input) INTEGER
127 * The leading dimension of the array B. LDB >= max(1,N).
128 *
129 * X (output) DOUBLE PRECISION array, dimension (LDX,NRHS)
130 * If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
131 *
132 * LDX (input) INTEGER
133 * The leading dimension of the array X. LDX >= max(1,N).
134 *
135 * RCOND (output) DOUBLE PRECISION
136 * The estimate of the reciprocal condition number of the matrix
137 * A. If RCOND is less than the machine precision (in
138 * particular, if RCOND = 0), the matrix is singular to working
139 * precision. This condition is indicated by a return code of
140 * INFO > 0.
141 *
142 * FERR (output) DOUBLE PRECISION array, dimension (NRHS)
143 * The estimated forward error bound for each solution vector
144 * 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 largest
147 * element in (X(j) - XTRUE) divided by the magnitude of the
148 * largest element in X(j). The estimate is as reliable as
149 * the estimate for RCOND, and is almost always a slight
150 * overestimate of the true error.
151 *
152 * BERR (output) DOUBLE PRECISION array, dimension (NRHS)
153 * The componentwise relative backward error of each solution
154 * vector X(j) (i.e., the smallest relative change in
155 * any element of A or B that makes X(j) an exact solution).
156 *
157 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
158 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
159 *
160 * LWORK (input) INTEGER
161 * The length of WORK. LWORK >= max(1,3*N), and for best
162 * performance, when FACT = 'N', LWORK >= max(1,3*N,N*NB), where
163 * NB is the optimal blocksize for DSYTRF.
164 *
165 * If LWORK = -1, then a workspace query is assumed; the routine
166 * only calculates the optimal size of the WORK array, returns
167 * this value as the first entry of the WORK array, and no error
168 * message related to LWORK is issued by XERBLA.
169 *
170 * IWORK (workspace) INTEGER array, dimension (N)
171 *
172 * INFO (output) INTEGER
173 * = 0: successful exit
174 * < 0: if INFO = -i, the i-th argument had an illegal value
175 * > 0: if INFO = i, and i is
176 * <= N: D(i,i) is exactly zero. The factorization
177 * has been completed but the factor D is exactly
178 * singular, so the solution and error bounds could
179 * not be computed. RCOND = 0 is returned.
180 * = N+1: D is nonsingular, but RCOND is less than machine
181 * precision, meaning that the matrix is singular
182 * to working precision. Nevertheless, the
183 * solution and error bounds are computed because
184 * there are a number of situations where the
185 * computed solution can be more accurate than the
186 * value of RCOND would suggest.
187 *
188 * =====================================================================
189 *
190 * .. Parameters ..
191 DOUBLE PRECISION ZERO
192 PARAMETER ( ZERO = 0.0D+0 )
193 * ..
194 * .. Local Scalars ..
195 LOGICAL LQUERY, NOFACT
196 INTEGER LWKOPT, NB
197 DOUBLE PRECISION ANORM
198 * ..
199 * .. External Functions ..
200 LOGICAL LSAME
201 INTEGER ILAENV
202 DOUBLE PRECISION DLAMCH, DLANSY
203 EXTERNAL LSAME, ILAENV, DLAMCH, DLANSY
204 * ..
205 * .. External Subroutines ..
206 EXTERNAL DLACPY, DSYCON, DSYRFS, DSYTRF, DSYTRS, XERBLA
207 * ..
208 * .. Intrinsic Functions ..
209 INTRINSIC MAX
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' ) ) THEN
219 INFO = -1
220 ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) )
221 $ THEN
222 INFO = -2
223 ELSE IF( N.LT.0 ) THEN
224 INFO = -3
225 ELSE IF( NRHS.LT.0 ) THEN
226 INFO = -4
227 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
228 INFO = -6
229 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
230 INFO = -8
231 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
232 INFO = -11
233 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
234 INFO = -13
235 ELSE IF( LWORK.LT.MAX( 1, 3*N ) .AND. .NOT.LQUERY ) THEN
236 INFO = -18
237 END IF
238 *
239 IF( INFO.EQ.0 ) THEN
240 LWKOPT = MAX( 1, 3*N )
241 IF( NOFACT ) THEN
242 NB = ILAENV( 1, 'DSYTRF', UPLO, N, -1, -1, -1 )
243 LWKOPT = MAX( LWKOPT, N*NB )
244 END IF
245 WORK( 1 ) = LWKOPT
246 END IF
247 *
248 IF( INFO.NE.0 ) THEN
249 CALL XERBLA( 'DSYSVX', -INFO )
250 RETURN
251 ELSE IF( LQUERY ) THEN
252 RETURN
253 END IF
254 *
255 IF( NOFACT ) THEN
256 *
257 * Compute the factorization A = U*D*U**T or A = L*D*L**T.
258 *
259 CALL DLACPY( UPLO, N, N, A, LDA, AF, LDAF )
260 CALL DSYTRF( UPLO, N, AF, LDAF, IPIV, WORK, LWORK, INFO )
261 *
262 * Return if INFO is non-zero.
263 *
264 IF( INFO.GT.0 )THEN
265 RCOND = ZERO
266 RETURN
267 END IF
268 END IF
269 *
270 * Compute the norm of the matrix A.
271 *
272 ANORM = DLANSY( 'I', UPLO, N, A, LDA, WORK )
273 *
274 * Compute the reciprocal of the condition number of A.
275 *
276 CALL DSYCON( UPLO, N, AF, LDAF, IPIV, ANORM, RCOND, WORK, IWORK,
277 $ INFO )
278 *
279 * Compute the solution vectors X.
280 *
281 CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
282 CALL DSYTRS( UPLO, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO )
283 *
284 * Use iterative refinement to improve the computed solutions and
285 * compute error bounds and backward error estimates for them.
286 *
287 CALL DSYRFS( 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.DLAMCH( 'Epsilon' ) )
293 $ INFO = N + 1
294 *
295 WORK( 1 ) = LWKOPT
296 *
297 RETURN
298 *
299 * End of DSYSVX
300 *
301 END
2 $ LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK,
3 $ IWORK, INFO )
4 *
5 * -- LAPACK driver routine (version 3.3.1) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * -- April 2011 --
9 *
10 * .. Scalar Arguments ..
11 CHARACTER FACT, UPLO
12 INTEGER INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS
13 DOUBLE PRECISION RCOND
14 * ..
15 * .. Array Arguments ..
16 INTEGER IPIV( * ), IWORK( * )
17 DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
18 $ BERR( * ), FERR( * ), WORK( * ), X( LDX, * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * DSYSVX uses the diagonal pivoting factorization to compute the
25 * solution to a real system of linear equations A * X = B,
26 * where A is an N-by-N symmetric matrix and X and B are N-by-NRHS
27 * matrices.
28 *
29 * Error bounds on the solution and a condition estimate are also
30 * provided.
31 *
32 * Description
33 * ===========
34 *
35 * The following steps are performed:
36 *
37 * 1. If FACT = 'N', the diagonal pivoting method is used to factor A.
38 * The form of the factorization is
39 * A = U * D * U**T, if UPLO = 'U', or
40 * A = L * D * L**T, if UPLO = 'L',
41 * where U (or L) is a product of permutation and unit upper (lower)
42 * triangular matrices, and D is symmetric and block diagonal with
43 * 1-by-1 and 2-by-2 diagonal blocks.
44 *
45 * 2. If some D(i,i)=0, so that D is exactly singular, then the routine
46 * returns with INFO = i. Otherwise, the factored form of A is used
47 * to estimate the condition number of the matrix A. If the
48 * reciprocal of the condition number is less than machine precision,
49 * INFO = N+1 is returned as a warning, but the routine still goes on
50 * to solve for X and compute error bounds as described below.
51 *
52 * 3. The system of equations is solved for X using the factored form
53 * of A.
54 *
55 * 4. Iterative refinement is applied to improve the computed solution
56 * matrix and calculate error bounds and backward error estimates
57 * for it.
58 *
59 * Arguments
60 * =========
61 *
62 * FACT (input) CHARACTER*1
63 * Specifies whether or not the factored form of A has been
64 * supplied on entry.
65 * = 'F': On entry, AF and IPIV contain the factored form of
66 * A. AF and IPIV will not be modified.
67 * = 'N': The matrix A will be copied to AF and factored.
68 *
69 * UPLO (input) CHARACTER*1
70 * = 'U': Upper triangle of A is stored;
71 * = 'L': Lower triangle of A is stored.
72 *
73 * N (input) INTEGER
74 * The number of linear equations, i.e., the order of the
75 * matrix A. N >= 0.
76 *
77 * NRHS (input) INTEGER
78 * The number of right hand sides, i.e., the number of columns
79 * of the matrices B and X. NRHS >= 0.
80 *
81 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
82 * The symmetric matrix A. If UPLO = 'U', the leading N-by-N
83 * upper triangular part of A contains the upper triangular part
84 * of the matrix A, and the strictly lower triangular part of A
85 * is not referenced. If UPLO = 'L', the leading N-by-N lower
86 * triangular part of A contains the lower triangular part of
87 * the matrix A, and the strictly upper triangular part of A is
88 * not referenced.
89 *
90 * LDA (input) INTEGER
91 * The leading dimension of the array A. LDA >= max(1,N).
92 *
93 * AF (input or output) DOUBLE PRECISION array, dimension (LDAF,N)
94 * If FACT = 'F', then AF is an input argument and on entry
95 * contains the block diagonal matrix D and the multipliers used
96 * to obtain the factor U or L from the factorization
97 * A = U*D*U**T or A = L*D*L**T as computed by DSYTRF.
98 *
99 * If FACT = 'N', then AF is an output argument and on exit
100 * returns the block diagonal matrix D and the multipliers used
101 * to obtain the factor U or L from the factorization
102 * A = U*D*U**T or A = L*D*L**T.
103 *
104 * LDAF (input) INTEGER
105 * 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 entry
109 * contains details of the interchanges and the block structure
110 * of D, as determined by DSYTRF.
111 * If IPIV(k) > 0, then rows and columns k and IPIV(k) were
112 * 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 and
114 * 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) were
117 * 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 exit
120 * contains details of the interchanges and the block structure
121 * of D, as determined by DSYTRF.
122 *
123 * B (input) DOUBLE PRECISION array, dimension (LDB,NRHS)
124 * The N-by-NRHS right hand side matrix B.
125 *
126 * LDB (input) INTEGER
127 * The leading dimension of the array B. LDB >= max(1,N).
128 *
129 * X (output) DOUBLE PRECISION array, dimension (LDX,NRHS)
130 * If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
131 *
132 * LDX (input) INTEGER
133 * The leading dimension of the array X. LDX >= max(1,N).
134 *
135 * RCOND (output) DOUBLE PRECISION
136 * The estimate of the reciprocal condition number of the matrix
137 * A. If RCOND is less than the machine precision (in
138 * particular, if RCOND = 0), the matrix is singular to working
139 * precision. This condition is indicated by a return code of
140 * INFO > 0.
141 *
142 * FERR (output) DOUBLE PRECISION array, dimension (NRHS)
143 * The estimated forward error bound for each solution vector
144 * 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 largest
147 * element in (X(j) - XTRUE) divided by the magnitude of the
148 * largest element in X(j). The estimate is as reliable as
149 * the estimate for RCOND, and is almost always a slight
150 * overestimate of the true error.
151 *
152 * BERR (output) DOUBLE PRECISION array, dimension (NRHS)
153 * The componentwise relative backward error of each solution
154 * vector X(j) (i.e., the smallest relative change in
155 * any element of A or B that makes X(j) an exact solution).
156 *
157 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
158 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
159 *
160 * LWORK (input) INTEGER
161 * The length of WORK. LWORK >= max(1,3*N), and for best
162 * performance, when FACT = 'N', LWORK >= max(1,3*N,N*NB), where
163 * NB is the optimal blocksize for DSYTRF.
164 *
165 * If LWORK = -1, then a workspace query is assumed; the routine
166 * only calculates the optimal size of the WORK array, returns
167 * this value as the first entry of the WORK array, and no error
168 * message related to LWORK is issued by XERBLA.
169 *
170 * IWORK (workspace) INTEGER array, dimension (N)
171 *
172 * INFO (output) INTEGER
173 * = 0: successful exit
174 * < 0: if INFO = -i, the i-th argument had an illegal value
175 * > 0: if INFO = i, and i is
176 * <= N: D(i,i) is exactly zero. The factorization
177 * has been completed but the factor D is exactly
178 * singular, so the solution and error bounds could
179 * not be computed. RCOND = 0 is returned.
180 * = N+1: D is nonsingular, but RCOND is less than machine
181 * precision, meaning that the matrix is singular
182 * to working precision. Nevertheless, the
183 * solution and error bounds are computed because
184 * there are a number of situations where the
185 * computed solution can be more accurate than the
186 * value of RCOND would suggest.
187 *
188 * =====================================================================
189 *
190 * .. Parameters ..
191 DOUBLE PRECISION ZERO
192 PARAMETER ( ZERO = 0.0D+0 )
193 * ..
194 * .. Local Scalars ..
195 LOGICAL LQUERY, NOFACT
196 INTEGER LWKOPT, NB
197 DOUBLE PRECISION ANORM
198 * ..
199 * .. External Functions ..
200 LOGICAL LSAME
201 INTEGER ILAENV
202 DOUBLE PRECISION DLAMCH, DLANSY
203 EXTERNAL LSAME, ILAENV, DLAMCH, DLANSY
204 * ..
205 * .. External Subroutines ..
206 EXTERNAL DLACPY, DSYCON, DSYRFS, DSYTRF, DSYTRS, XERBLA
207 * ..
208 * .. Intrinsic Functions ..
209 INTRINSIC MAX
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' ) ) THEN
219 INFO = -1
220 ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) )
221 $ THEN
222 INFO = -2
223 ELSE IF( N.LT.0 ) THEN
224 INFO = -3
225 ELSE IF( NRHS.LT.0 ) THEN
226 INFO = -4
227 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
228 INFO = -6
229 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
230 INFO = -8
231 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
232 INFO = -11
233 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
234 INFO = -13
235 ELSE IF( LWORK.LT.MAX( 1, 3*N ) .AND. .NOT.LQUERY ) THEN
236 INFO = -18
237 END IF
238 *
239 IF( INFO.EQ.0 ) THEN
240 LWKOPT = MAX( 1, 3*N )
241 IF( NOFACT ) THEN
242 NB = ILAENV( 1, 'DSYTRF', UPLO, N, -1, -1, -1 )
243 LWKOPT = MAX( LWKOPT, N*NB )
244 END IF
245 WORK( 1 ) = LWKOPT
246 END IF
247 *
248 IF( INFO.NE.0 ) THEN
249 CALL XERBLA( 'DSYSVX', -INFO )
250 RETURN
251 ELSE IF( LQUERY ) THEN
252 RETURN
253 END IF
254 *
255 IF( NOFACT ) THEN
256 *
257 * Compute the factorization A = U*D*U**T or A = L*D*L**T.
258 *
259 CALL DLACPY( UPLO, N, N, A, LDA, AF, LDAF )
260 CALL DSYTRF( UPLO, N, AF, LDAF, IPIV, WORK, LWORK, INFO )
261 *
262 * Return if INFO is non-zero.
263 *
264 IF( INFO.GT.0 )THEN
265 RCOND = ZERO
266 RETURN
267 END IF
268 END IF
269 *
270 * Compute the norm of the matrix A.
271 *
272 ANORM = DLANSY( 'I', UPLO, N, A, LDA, WORK )
273 *
274 * Compute the reciprocal of the condition number of A.
275 *
276 CALL DSYCON( UPLO, N, AF, LDAF, IPIV, ANORM, RCOND, WORK, IWORK,
277 $ INFO )
278 *
279 * Compute the solution vectors X.
280 *
281 CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
282 CALL DSYTRS( UPLO, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO )
283 *
284 * Use iterative refinement to improve the computed solutions and
285 * compute error bounds and backward error estimates for them.
286 *
287 CALL DSYRFS( 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.DLAMCH( 'Epsilon' ) )
293 $ INFO = N + 1
294 *
295 WORK( 1 ) = LWKOPT
296 *
297 RETURN
298 *
299 * End of DSYSVX
300 *
301 END