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