1 SUBROUTINE DPOSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED,
2 $ S, B, LDB, X, LDX, RCOND, FERR, BERR, WORK,
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 EQUED, FACT, UPLO
12 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS
13 DOUBLE PRECISION RCOND
14 * ..
15 * .. Array Arguments ..
16 INTEGER IWORK( * )
17 DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
18 $ BERR( * ), FERR( * ), S( * ), WORK( * ),
19 $ X( LDX, * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * DPOSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to
26 * compute the solution to a real system of linear equations
27 * A * X = B,
28 * where A is an N-by-N symmetric positive definite matrix and X and B
29 * are N-by-NRHS matrices.
30 *
31 * Error bounds on the solution and a condition estimate are also
32 * provided.
33 *
34 * Description
35 * ===========
36 *
37 * The following steps are performed:
38 *
39 * 1. If FACT = 'E', real scaling factors are computed to equilibrate
40 * the system:
41 * diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B
42 * Whether or not the system will be equilibrated depends on the
43 * scaling of the matrix A, but if equilibration is used, A is
44 * overwritten by diag(S)*A*diag(S) and B by diag(S)*B.
45 *
46 * 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
47 * factor the matrix A (after equilibration if FACT = 'E') as
48 * A = U**T* U, if UPLO = 'U', or
49 * A = L * L**T, if UPLO = 'L',
50 * where U is an upper triangular matrix and L is a lower triangular
51 * matrix.
52 *
53 * 3. If the leading i-by-i principal minor is not positive definite,
54 * then the routine returns with INFO = i. Otherwise, the factored
55 * form of A is used to estimate the condition number of the matrix
56 * A. If the reciprocal of the condition number is less than machine
57 * precision, INFO = N+1 is returned as a warning, but the routine
58 * still goes on to solve for X and compute error bounds as
59 * described below.
60 *
61 * 4. The system of equations is solved for X using the factored form
62 * of A.
63 *
64 * 5. Iterative refinement is applied to improve the computed solution
65 * matrix and calculate error bounds and backward error estimates
66 * for it.
67 *
68 * 6. If equilibration was used, the matrix X is premultiplied by
69 * diag(S) so that it solves the original system before
70 * equilibration.
71 *
72 * Arguments
73 * =========
74 *
75 * FACT (input) CHARACTER*1
76 * Specifies whether or not the factored form of the matrix A is
77 * supplied on entry, and if not, whether the matrix A should be
78 * equilibrated before it is factored.
79 * = 'F': On entry, AF contains the factored form of A.
80 * If EQUED = 'Y', the matrix A has been equilibrated
81 * with scaling factors given by S. A and AF will not
82 * be modified.
83 * = 'N': The matrix A will be copied to AF and factored.
84 * = 'E': The matrix A will be equilibrated if necessary, then
85 * copied to AF and factored.
86 *
87 * UPLO (input) CHARACTER*1
88 * = 'U': Upper triangle of A is stored;
89 * = 'L': Lower triangle of A is stored.
90 *
91 * N (input) INTEGER
92 * The number of linear equations, i.e., the order of the
93 * matrix A. N >= 0.
94 *
95 * NRHS (input) INTEGER
96 * The number of right hand sides, i.e., the number of columns
97 * of the matrices B and X. NRHS >= 0.
98 *
99 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
100 * On entry, the symmetric matrix A, except if FACT = 'F' and
101 * EQUED = 'Y', then A must contain the equilibrated matrix
102 * diag(S)*A*diag(S). If UPLO = 'U', the leading
103 * N-by-N upper triangular part of A contains the upper
104 * triangular part of the matrix A, and the strictly lower
105 * triangular part of A is not referenced. If UPLO = 'L', the
106 * leading N-by-N lower triangular part of A contains the lower
107 * triangular part of the matrix A, and the strictly upper
108 * triangular part of A is not referenced. A is not modified if
109 * FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit.
110 *
111 * On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
112 * diag(S)*A*diag(S).
113 *
114 * LDA (input) INTEGER
115 * The leading dimension of the array A. LDA >= max(1,N).
116 *
117 * AF (input or output) DOUBLE PRECISION array, dimension (LDAF,N)
118 * If FACT = 'F', then AF is an input argument and on entry
119 * contains the triangular factor U or L from the Cholesky
120 * factorization A = U**T*U or A = L*L**T, in the same storage
121 * format as A. If EQUED .ne. 'N', then AF is the factored form
122 * of the equilibrated matrix diag(S)*A*diag(S).
123 *
124 * If FACT = 'N', then AF is an output argument and on exit
125 * returns the triangular factor U or L from the Cholesky
126 * factorization A = U**T*U or A = L*L**T of the original
127 * matrix A.
128 *
129 * If FACT = 'E', then AF is an output argument and on exit
130 * returns the triangular factor U or L from the Cholesky
131 * factorization A = U**T*U or A = L*L**T of the equilibrated
132 * matrix A (see the description of A for the form of the
133 * equilibrated matrix).
134 *
135 * LDAF (input) INTEGER
136 * The leading dimension of the array AF. LDAF >= max(1,N).
137 *
138 * EQUED (input or output) CHARACTER*1
139 * Specifies the form of equilibration that was done.
140 * = 'N': No equilibration (always true if FACT = 'N').
141 * = 'Y': Equilibration was done, i.e., A has been replaced by
142 * diag(S) * A * diag(S).
143 * EQUED is an input argument if FACT = 'F'; otherwise, it is an
144 * output argument.
145 *
146 * S (input or output) DOUBLE PRECISION array, dimension (N)
147 * The scale factors for A; not accessed if EQUED = 'N'. S is
148 * an input argument if FACT = 'F'; otherwise, S is an output
149 * argument. If FACT = 'F' and EQUED = 'Y', each element of S
150 * must be positive.
151 *
152 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
153 * On entry, the N-by-NRHS right hand side matrix B.
154 * On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y',
155 * B is overwritten by diag(S) * B.
156 *
157 * LDB (input) INTEGER
158 * The leading dimension of the array B. LDB >= max(1,N).
159 *
160 * X (output) DOUBLE PRECISION array, dimension (LDX,NRHS)
161 * If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to
162 * the original system of equations. Note that if EQUED = 'Y',
163 * A and B are modified on exit, and the solution to the
164 * equilibrated system is inv(diag(S))*X.
165 *
166 * LDX (input) INTEGER
167 * The leading dimension of the array X. LDX >= max(1,N).
168 *
169 * RCOND (output) DOUBLE PRECISION
170 * The estimate of the reciprocal condition number of the matrix
171 * A after equilibration (if done). If RCOND is less than the
172 * machine precision (in particular, if RCOND = 0), the matrix
173 * is singular to working precision. This condition is
174 * indicated by a return code of INFO > 0.
175 *
176 * FERR (output) DOUBLE PRECISION array, dimension (NRHS)
177 * The estimated forward error bound for each solution vector
178 * X(j) (the j-th column of the solution matrix X).
179 * If XTRUE is the true solution corresponding to X(j), FERR(j)
180 * is an estimated upper bound for the magnitude of the largest
181 * element in (X(j) - XTRUE) divided by the magnitude of the
182 * largest element in X(j). The estimate is as reliable as
183 * the estimate for RCOND, and is almost always a slight
184 * overestimate of the true error.
185 *
186 * BERR (output) DOUBLE PRECISION array, dimension (NRHS)
187 * The componentwise relative backward error of each solution
188 * vector X(j) (i.e., the smallest relative change in
189 * any element of A or B that makes X(j) an exact solution).
190 *
191 * WORK (workspace) DOUBLE PRECISION array, dimension (3*N)
192 *
193 * IWORK (workspace) INTEGER array, dimension (N)
194 *
195 * INFO (output) INTEGER
196 * = 0: successful exit
197 * < 0: if INFO = -i, the i-th argument had an illegal value
198 * > 0: if INFO = i, and i is
199 * <= N: the leading minor of order i of A is
200 * not positive definite, so the factorization
201 * could not be completed, and the solution has not
202 * been computed. RCOND = 0 is returned.
203 * = N+1: U is nonsingular, but RCOND is less than machine
204 * precision, meaning that the matrix is singular
205 * to working precision. Nevertheless, the
206 * solution and error bounds are computed because
207 * there are a number of situations where the
208 * computed solution can be more accurate than the
209 * value of RCOND would suggest.
210 *
211 * =====================================================================
212 *
213 * .. Parameters ..
214 DOUBLE PRECISION ZERO, ONE
215 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
216 * ..
217 * .. Local Scalars ..
218 LOGICAL EQUIL, NOFACT, RCEQU
219 INTEGER I, INFEQU, J
220 DOUBLE PRECISION AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM
221 * ..
222 * .. External Functions ..
223 LOGICAL LSAME
224 DOUBLE PRECISION DLAMCH, DLANSY
225 EXTERNAL LSAME, DLAMCH, DLANSY
226 * ..
227 * .. External Subroutines ..
228 EXTERNAL DLACPY, DLAQSY, DPOCON, DPOEQU, DPORFS, DPOTRF,
229 $ DPOTRS, XERBLA
230 * ..
231 * .. Intrinsic Functions ..
232 INTRINSIC MAX, MIN
233 * ..
234 * .. Executable Statements ..
235 *
236 INFO = 0
237 NOFACT = LSAME( FACT, 'N' )
238 EQUIL = LSAME( FACT, 'E' )
239 IF( NOFACT .OR. EQUIL ) THEN
240 EQUED = 'N'
241 RCEQU = .FALSE.
242 ELSE
243 RCEQU = LSAME( EQUED, 'Y' )
244 SMLNUM = DLAMCH( 'Safe minimum' )
245 BIGNUM = ONE / SMLNUM
246 END IF
247 *
248 * Test the input parameters.
249 *
250 IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) )
251 $ THEN
252 INFO = -1
253 ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) )
254 $ THEN
255 INFO = -2
256 ELSE IF( N.LT.0 ) THEN
257 INFO = -3
258 ELSE IF( NRHS.LT.0 ) THEN
259 INFO = -4
260 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
261 INFO = -6
262 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
263 INFO = -8
264 ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT.
265 $ ( RCEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN
266 INFO = -9
267 ELSE
268 IF( RCEQU ) THEN
269 SMIN = BIGNUM
270 SMAX = ZERO
271 DO 10 J = 1, N
272 SMIN = MIN( SMIN, S( J ) )
273 SMAX = MAX( SMAX, S( J ) )
274 10 CONTINUE
275 IF( SMIN.LE.ZERO ) THEN
276 INFO = -10
277 ELSE IF( N.GT.0 ) THEN
278 SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM )
279 ELSE
280 SCOND = ONE
281 END IF
282 END IF
283 IF( INFO.EQ.0 ) THEN
284 IF( LDB.LT.MAX( 1, N ) ) THEN
285 INFO = -12
286 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
287 INFO = -14
288 END IF
289 END IF
290 END IF
291 *
292 IF( INFO.NE.0 ) THEN
293 CALL XERBLA( 'DPOSVX', -INFO )
294 RETURN
295 END IF
296 *
297 IF( EQUIL ) THEN
298 *
299 * Compute row and column scalings to equilibrate the matrix A.
300 *
301 CALL DPOEQU( N, A, LDA, S, SCOND, AMAX, INFEQU )
302 IF( INFEQU.EQ.0 ) THEN
303 *
304 * Equilibrate the matrix.
305 *
306 CALL DLAQSY( UPLO, N, A, LDA, S, SCOND, AMAX, EQUED )
307 RCEQU = LSAME( EQUED, 'Y' )
308 END IF
309 END IF
310 *
311 * Scale the right hand side.
312 *
313 IF( RCEQU ) THEN
314 DO 30 J = 1, NRHS
315 DO 20 I = 1, N
316 B( I, J ) = S( I )*B( I, J )
317 20 CONTINUE
318 30 CONTINUE
319 END IF
320 *
321 IF( NOFACT .OR. EQUIL ) THEN
322 *
323 * Compute the Cholesky factorization A = U**T *U or A = L*L**T.
324 *
325 CALL DLACPY( UPLO, N, N, A, LDA, AF, LDAF )
326 CALL DPOTRF( UPLO, N, AF, LDAF, INFO )
327 *
328 * Return if INFO is non-zero.
329 *
330 IF( INFO.GT.0 )THEN
331 RCOND = ZERO
332 RETURN
333 END IF
334 END IF
335 *
336 * Compute the norm of the matrix A.
337 *
338 ANORM = DLANSY( '1', UPLO, N, A, LDA, WORK )
339 *
340 * Compute the reciprocal of the condition number of A.
341 *
342 CALL DPOCON( UPLO, N, AF, LDAF, ANORM, RCOND, WORK, IWORK, INFO )
343 *
344 * Compute the solution matrix X.
345 *
346 CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
347 CALL DPOTRS( UPLO, N, NRHS, AF, LDAF, X, LDX, INFO )
348 *
349 * Use iterative refinement to improve the computed solution and
350 * compute error bounds and backward error estimates for it.
351 *
352 CALL DPORFS( UPLO, N, NRHS, A, LDA, AF, LDAF, B, LDB, X, LDX,
353 $ FERR, BERR, WORK, IWORK, INFO )
354 *
355 * Transform the solution matrix X to a solution of the original
356 * system.
357 *
358 IF( RCEQU ) THEN
359 DO 50 J = 1, NRHS
360 DO 40 I = 1, N
361 X( I, J ) = S( I )*X( I, J )
362 40 CONTINUE
363 50 CONTINUE
364 DO 60 J = 1, NRHS
365 FERR( J ) = FERR( J ) / SCOND
366 60 CONTINUE
367 END IF
368 *
369 * Set INFO = N+1 if the matrix is singular to working precision.
370 *
371 IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
372 $ INFO = N + 1
373 *
374 RETURN
375 *
376 * End of DPOSVX
377 *
378 END
2 $ S, B, LDB, X, LDX, RCOND, FERR, BERR, WORK,
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 EQUED, FACT, UPLO
12 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS
13 DOUBLE PRECISION RCOND
14 * ..
15 * .. Array Arguments ..
16 INTEGER IWORK( * )
17 DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
18 $ BERR( * ), FERR( * ), S( * ), WORK( * ),
19 $ X( LDX, * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * DPOSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to
26 * compute the solution to a real system of linear equations
27 * A * X = B,
28 * where A is an N-by-N symmetric positive definite matrix and X and B
29 * are N-by-NRHS matrices.
30 *
31 * Error bounds on the solution and a condition estimate are also
32 * provided.
33 *
34 * Description
35 * ===========
36 *
37 * The following steps are performed:
38 *
39 * 1. If FACT = 'E', real scaling factors are computed to equilibrate
40 * the system:
41 * diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B
42 * Whether or not the system will be equilibrated depends on the
43 * scaling of the matrix A, but if equilibration is used, A is
44 * overwritten by diag(S)*A*diag(S) and B by diag(S)*B.
45 *
46 * 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
47 * factor the matrix A (after equilibration if FACT = 'E') as
48 * A = U**T* U, if UPLO = 'U', or
49 * A = L * L**T, if UPLO = 'L',
50 * where U is an upper triangular matrix and L is a lower triangular
51 * matrix.
52 *
53 * 3. If the leading i-by-i principal minor is not positive definite,
54 * then the routine returns with INFO = i. Otherwise, the factored
55 * form of A is used to estimate the condition number of the matrix
56 * A. If the reciprocal of the condition number is less than machine
57 * precision, INFO = N+1 is returned as a warning, but the routine
58 * still goes on to solve for X and compute error bounds as
59 * described below.
60 *
61 * 4. The system of equations is solved for X using the factored form
62 * of A.
63 *
64 * 5. Iterative refinement is applied to improve the computed solution
65 * matrix and calculate error bounds and backward error estimates
66 * for it.
67 *
68 * 6. If equilibration was used, the matrix X is premultiplied by
69 * diag(S) so that it solves the original system before
70 * equilibration.
71 *
72 * Arguments
73 * =========
74 *
75 * FACT (input) CHARACTER*1
76 * Specifies whether or not the factored form of the matrix A is
77 * supplied on entry, and if not, whether the matrix A should be
78 * equilibrated before it is factored.
79 * = 'F': On entry, AF contains the factored form of A.
80 * If EQUED = 'Y', the matrix A has been equilibrated
81 * with scaling factors given by S. A and AF will not
82 * be modified.
83 * = 'N': The matrix A will be copied to AF and factored.
84 * = 'E': The matrix A will be equilibrated if necessary, then
85 * copied to AF and factored.
86 *
87 * UPLO (input) CHARACTER*1
88 * = 'U': Upper triangle of A is stored;
89 * = 'L': Lower triangle of A is stored.
90 *
91 * N (input) INTEGER
92 * The number of linear equations, i.e., the order of the
93 * matrix A. N >= 0.
94 *
95 * NRHS (input) INTEGER
96 * The number of right hand sides, i.e., the number of columns
97 * of the matrices B and X. NRHS >= 0.
98 *
99 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
100 * On entry, the symmetric matrix A, except if FACT = 'F' and
101 * EQUED = 'Y', then A must contain the equilibrated matrix
102 * diag(S)*A*diag(S). If UPLO = 'U', the leading
103 * N-by-N upper triangular part of A contains the upper
104 * triangular part of the matrix A, and the strictly lower
105 * triangular part of A is not referenced. If UPLO = 'L', the
106 * leading N-by-N lower triangular part of A contains the lower
107 * triangular part of the matrix A, and the strictly upper
108 * triangular part of A is not referenced. A is not modified if
109 * FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit.
110 *
111 * On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
112 * diag(S)*A*diag(S).
113 *
114 * LDA (input) INTEGER
115 * The leading dimension of the array A. LDA >= max(1,N).
116 *
117 * AF (input or output) DOUBLE PRECISION array, dimension (LDAF,N)
118 * If FACT = 'F', then AF is an input argument and on entry
119 * contains the triangular factor U or L from the Cholesky
120 * factorization A = U**T*U or A = L*L**T, in the same storage
121 * format as A. If EQUED .ne. 'N', then AF is the factored form
122 * of the equilibrated matrix diag(S)*A*diag(S).
123 *
124 * If FACT = 'N', then AF is an output argument and on exit
125 * returns the triangular factor U or L from the Cholesky
126 * factorization A = U**T*U or A = L*L**T of the original
127 * matrix A.
128 *
129 * If FACT = 'E', then AF is an output argument and on exit
130 * returns the triangular factor U or L from the Cholesky
131 * factorization A = U**T*U or A = L*L**T of the equilibrated
132 * matrix A (see the description of A for the form of the
133 * equilibrated matrix).
134 *
135 * LDAF (input) INTEGER
136 * The leading dimension of the array AF. LDAF >= max(1,N).
137 *
138 * EQUED (input or output) CHARACTER*1
139 * Specifies the form of equilibration that was done.
140 * = 'N': No equilibration (always true if FACT = 'N').
141 * = 'Y': Equilibration was done, i.e., A has been replaced by
142 * diag(S) * A * diag(S).
143 * EQUED is an input argument if FACT = 'F'; otherwise, it is an
144 * output argument.
145 *
146 * S (input or output) DOUBLE PRECISION array, dimension (N)
147 * The scale factors for A; not accessed if EQUED = 'N'. S is
148 * an input argument if FACT = 'F'; otherwise, S is an output
149 * argument. If FACT = 'F' and EQUED = 'Y', each element of S
150 * must be positive.
151 *
152 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
153 * On entry, the N-by-NRHS right hand side matrix B.
154 * On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y',
155 * B is overwritten by diag(S) * B.
156 *
157 * LDB (input) INTEGER
158 * The leading dimension of the array B. LDB >= max(1,N).
159 *
160 * X (output) DOUBLE PRECISION array, dimension (LDX,NRHS)
161 * If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to
162 * the original system of equations. Note that if EQUED = 'Y',
163 * A and B are modified on exit, and the solution to the
164 * equilibrated system is inv(diag(S))*X.
165 *
166 * LDX (input) INTEGER
167 * The leading dimension of the array X. LDX >= max(1,N).
168 *
169 * RCOND (output) DOUBLE PRECISION
170 * The estimate of the reciprocal condition number of the matrix
171 * A after equilibration (if done). If RCOND is less than the
172 * machine precision (in particular, if RCOND = 0), the matrix
173 * is singular to working precision. This condition is
174 * indicated by a return code of INFO > 0.
175 *
176 * FERR (output) DOUBLE PRECISION array, dimension (NRHS)
177 * The estimated forward error bound for each solution vector
178 * X(j) (the j-th column of the solution matrix X).
179 * If XTRUE is the true solution corresponding to X(j), FERR(j)
180 * is an estimated upper bound for the magnitude of the largest
181 * element in (X(j) - XTRUE) divided by the magnitude of the
182 * largest element in X(j). The estimate is as reliable as
183 * the estimate for RCOND, and is almost always a slight
184 * overestimate of the true error.
185 *
186 * BERR (output) DOUBLE PRECISION array, dimension (NRHS)
187 * The componentwise relative backward error of each solution
188 * vector X(j) (i.e., the smallest relative change in
189 * any element of A or B that makes X(j) an exact solution).
190 *
191 * WORK (workspace) DOUBLE PRECISION array, dimension (3*N)
192 *
193 * IWORK (workspace) INTEGER array, dimension (N)
194 *
195 * INFO (output) INTEGER
196 * = 0: successful exit
197 * < 0: if INFO = -i, the i-th argument had an illegal value
198 * > 0: if INFO = i, and i is
199 * <= N: the leading minor of order i of A is
200 * not positive definite, so the factorization
201 * could not be completed, and the solution has not
202 * been computed. RCOND = 0 is returned.
203 * = N+1: U is nonsingular, but RCOND is less than machine
204 * precision, meaning that the matrix is singular
205 * to working precision. Nevertheless, the
206 * solution and error bounds are computed because
207 * there are a number of situations where the
208 * computed solution can be more accurate than the
209 * value of RCOND would suggest.
210 *
211 * =====================================================================
212 *
213 * .. Parameters ..
214 DOUBLE PRECISION ZERO, ONE
215 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
216 * ..
217 * .. Local Scalars ..
218 LOGICAL EQUIL, NOFACT, RCEQU
219 INTEGER I, INFEQU, J
220 DOUBLE PRECISION AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM
221 * ..
222 * .. External Functions ..
223 LOGICAL LSAME
224 DOUBLE PRECISION DLAMCH, DLANSY
225 EXTERNAL LSAME, DLAMCH, DLANSY
226 * ..
227 * .. External Subroutines ..
228 EXTERNAL DLACPY, DLAQSY, DPOCON, DPOEQU, DPORFS, DPOTRF,
229 $ DPOTRS, XERBLA
230 * ..
231 * .. Intrinsic Functions ..
232 INTRINSIC MAX, MIN
233 * ..
234 * .. Executable Statements ..
235 *
236 INFO = 0
237 NOFACT = LSAME( FACT, 'N' )
238 EQUIL = LSAME( FACT, 'E' )
239 IF( NOFACT .OR. EQUIL ) THEN
240 EQUED = 'N'
241 RCEQU = .FALSE.
242 ELSE
243 RCEQU = LSAME( EQUED, 'Y' )
244 SMLNUM = DLAMCH( 'Safe minimum' )
245 BIGNUM = ONE / SMLNUM
246 END IF
247 *
248 * Test the input parameters.
249 *
250 IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) )
251 $ THEN
252 INFO = -1
253 ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) )
254 $ THEN
255 INFO = -2
256 ELSE IF( N.LT.0 ) THEN
257 INFO = -3
258 ELSE IF( NRHS.LT.0 ) THEN
259 INFO = -4
260 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
261 INFO = -6
262 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
263 INFO = -8
264 ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT.
265 $ ( RCEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN
266 INFO = -9
267 ELSE
268 IF( RCEQU ) THEN
269 SMIN = BIGNUM
270 SMAX = ZERO
271 DO 10 J = 1, N
272 SMIN = MIN( SMIN, S( J ) )
273 SMAX = MAX( SMAX, S( J ) )
274 10 CONTINUE
275 IF( SMIN.LE.ZERO ) THEN
276 INFO = -10
277 ELSE IF( N.GT.0 ) THEN
278 SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM )
279 ELSE
280 SCOND = ONE
281 END IF
282 END IF
283 IF( INFO.EQ.0 ) THEN
284 IF( LDB.LT.MAX( 1, N ) ) THEN
285 INFO = -12
286 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
287 INFO = -14
288 END IF
289 END IF
290 END IF
291 *
292 IF( INFO.NE.0 ) THEN
293 CALL XERBLA( 'DPOSVX', -INFO )
294 RETURN
295 END IF
296 *
297 IF( EQUIL ) THEN
298 *
299 * Compute row and column scalings to equilibrate the matrix A.
300 *
301 CALL DPOEQU( N, A, LDA, S, SCOND, AMAX, INFEQU )
302 IF( INFEQU.EQ.0 ) THEN
303 *
304 * Equilibrate the matrix.
305 *
306 CALL DLAQSY( UPLO, N, A, LDA, S, SCOND, AMAX, EQUED )
307 RCEQU = LSAME( EQUED, 'Y' )
308 END IF
309 END IF
310 *
311 * Scale the right hand side.
312 *
313 IF( RCEQU ) THEN
314 DO 30 J = 1, NRHS
315 DO 20 I = 1, N
316 B( I, J ) = S( I )*B( I, J )
317 20 CONTINUE
318 30 CONTINUE
319 END IF
320 *
321 IF( NOFACT .OR. EQUIL ) THEN
322 *
323 * Compute the Cholesky factorization A = U**T *U or A = L*L**T.
324 *
325 CALL DLACPY( UPLO, N, N, A, LDA, AF, LDAF )
326 CALL DPOTRF( UPLO, N, AF, LDAF, INFO )
327 *
328 * Return if INFO is non-zero.
329 *
330 IF( INFO.GT.0 )THEN
331 RCOND = ZERO
332 RETURN
333 END IF
334 END IF
335 *
336 * Compute the norm of the matrix A.
337 *
338 ANORM = DLANSY( '1', UPLO, N, A, LDA, WORK )
339 *
340 * Compute the reciprocal of the condition number of A.
341 *
342 CALL DPOCON( UPLO, N, AF, LDAF, ANORM, RCOND, WORK, IWORK, INFO )
343 *
344 * Compute the solution matrix X.
345 *
346 CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
347 CALL DPOTRS( UPLO, N, NRHS, AF, LDAF, X, LDX, INFO )
348 *
349 * Use iterative refinement to improve the computed solution and
350 * compute error bounds and backward error estimates for it.
351 *
352 CALL DPORFS( UPLO, N, NRHS, A, LDA, AF, LDAF, B, LDB, X, LDX,
353 $ FERR, BERR, WORK, IWORK, INFO )
354 *
355 * Transform the solution matrix X to a solution of the original
356 * system.
357 *
358 IF( RCEQU ) THEN
359 DO 50 J = 1, NRHS
360 DO 40 I = 1, N
361 X( I, J ) = S( I )*X( I, J )
362 40 CONTINUE
363 50 CONTINUE
364 DO 60 J = 1, NRHS
365 FERR( J ) = FERR( J ) / SCOND
366 60 CONTINUE
367 END IF
368 *
369 * Set INFO = N+1 if the matrix is singular to working precision.
370 *
371 IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
372 $ INFO = N + 1
373 *
374 RETURN
375 *
376 * End of DPOSVX
377 *
378 END