1 /* ./src_f77/dspsvx.f -- translated by f2c (version 20030320).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include <punc/vf2c.h>
7
8 /* Table of constant values */
9
10 static integer c__1 = 1;
11
dspsvx_(char * fact,char * uplo,integer * n,integer * nrhs,doublereal * ap,doublereal * afp,integer * ipiv,doublereal * b,integer * ldb,doublereal * x,integer * ldx,doublereal * rcond,doublereal * ferr,doublereal * berr,doublereal * work,integer * iwork,integer * info,ftnlen fact_len,ftnlen uplo_len)12 /* Subroutine */ int dspsvx_(char *fact, char *uplo, integer *n, integer *
13 nrhs, doublereal *ap, doublereal *afp, integer *ipiv, doublereal *b,
14 integer *ldb, doublereal *x, integer *ldx, doublereal *rcond,
15 doublereal *ferr, doublereal *berr, doublereal *work, integer *iwork,
16 integer *info, ftnlen fact_len, ftnlen uplo_len)
17 {
18 /* System generated locals */
19 integer b_dim1, b_offset, x_dim1, x_offset, i__1;
20
21 /* Local variables */
22 extern logical lsame_(char *, char *, ftnlen, ftnlen);
23 static doublereal anorm;
24 extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
25 doublereal *, integer *);
26 extern doublereal dlamch_(char *, ftnlen);
27 static logical nofact;
28 extern /* Subroutine */ int dlacpy_(char *, integer *, integer *,
29 doublereal *, integer *, doublereal *, integer *, ftnlen),
30 xerbla_(char *, integer *, ftnlen);
31 extern doublereal dlansp_(char *, char *, integer *, doublereal *,
32 doublereal *, ftnlen, ftnlen);
33 extern /* Subroutine */ int dspcon_(char *, integer *, doublereal *,
34 integer *, doublereal *, doublereal *, doublereal *, integer *,
35 integer *, ftnlen), dsprfs_(char *, integer *, integer *,
36 doublereal *, doublereal *, integer *, doublereal *, integer *,
37 doublereal *, integer *, doublereal *, doublereal *, doublereal *,
38 integer *, integer *, ftnlen), dsptrf_(char *, integer *,
39 doublereal *, integer *, integer *, ftnlen), dsptrs_(char *,
40 integer *, integer *, doublereal *, integer *, doublereal *,
41 integer *, integer *, ftnlen);
42
43
44 /* -- LAPACK driver routine (version 3.0) -- */
45 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
46 /* Courant Institute, Argonne National Lab, and Rice University */
47 /* June 30, 1999 */
48
49 /* .. Scalar Arguments .. */
50 /* .. */
51 /* .. Array Arguments .. */
52 /* .. */
53
54 /* Purpose */
55 /* ======= */
56
57 /* DSPSVX uses the diagonal pivoting factorization A = U*D*U**T or */
58 /* A = L*D*L**T to compute the solution to a real system of linear */
59 /* equations A * X = B, where A is an N-by-N symmetric matrix stored */
60 /* in packed format and X and B are N-by-NRHS matrices. */
61
62 /* Error bounds on the solution and a condition estimate are also */
63 /* provided. */
64
65 /* Description */
66 /* =========== */
67
68 /* The following steps are performed: */
69
70 /* 1. If FACT = 'N', the diagonal pivoting method is used to factor A as */
71 /* A = U * D * U**T, if UPLO = 'U', or */
72 /* A = L * D * L**T, if UPLO = 'L', */
73 /* where U (or L) is a product of permutation and unit upper (lower) */
74 /* triangular matrices and D is symmetric and block diagonal with */
75 /* 1-by-1 and 2-by-2 diagonal blocks. */
76
77 /* 2. If some D(i,i)=0, so that D is exactly singular, then the routine */
78 /* returns with INFO = i. Otherwise, the factored form of A is used */
79 /* to estimate the condition number of the matrix A. If the */
80 /* reciprocal of the condition number is less than machine precision, */
81 /* INFO = N+1 is returned as a warning, but the routine still goes on */
82 /* to solve for X and compute error bounds as described below. */
83
84 /* 3. The system of equations is solved for X using the factored form */
85 /* of A. */
86
87 /* 4. Iterative refinement is applied to improve the computed solution */
88 /* matrix and calculate error bounds and backward error estimates */
89 /* for it. */
90
91 /* Arguments */
92 /* ========= */
93
94 /* FACT (input) CHARACTER*1 */
95 /* Specifies whether or not the factored form of A has been */
96 /* supplied on entry. */
97 /* = 'F': On entry, AFP and IPIV contain the factored form of */
98 /* A. AP, AFP and IPIV will not be modified. */
99 /* = 'N': The matrix A will be copied to AFP and factored. */
100
101 /* UPLO (input) CHARACTER*1 */
102 /* = 'U': Upper triangle of A is stored; */
103 /* = 'L': Lower triangle of A is stored. */
104
105 /* N (input) INTEGER */
106 /* The number of linear equations, i.e., the order of the */
107 /* matrix A. N >= 0. */
108
109 /* NRHS (input) INTEGER */
110 /* The number of right hand sides, i.e., the number of columns */
111 /* of the matrices B and X. NRHS >= 0. */
112
113 /* AP (input) DOUBLE PRECISION array, dimension (N*(N+1)/2) */
114 /* The upper or lower triangle of the symmetric matrix A, packed */
115 /* columnwise in a linear array. The j-th column of A is stored */
116 /* in the array AP as follows: */
117 /* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; */
118 /* if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n. */
119 /* See below for further details. */
120
121 /* AFP (input or output) DOUBLE PRECISION array, dimension */
122 /* (N*(N+1)/2) */
123 /* If FACT = 'F', then AFP is an input argument and on entry */
124 /* contains the block diagonal matrix D and the multipliers used */
125 /* to obtain the factor U or L from the factorization */
126 /* A = U*D*U**T or A = L*D*L**T as computed by DSPTRF, stored as */
127 /* a packed triangular matrix in the same storage format as A. */
128
129 /* If FACT = 'N', then AFP is an output argument and on exit */
130 /* contains the block diagonal matrix D and the multipliers used */
131 /* to obtain the factor U or L from the factorization */
132 /* A = U*D*U**T or A = L*D*L**T as computed by DSPTRF, stored as */
133 /* a packed triangular matrix in the same storage format as A. */
134
135 /* IPIV (input or output) INTEGER array, dimension (N) */
136 /* If FACT = 'F', then IPIV is an input argument and on entry */
137 /* contains details of the interchanges and the block structure */
138 /* of D, as determined by DSPTRF. */
139 /* If IPIV(k) > 0, then rows and columns k and IPIV(k) were */
140 /* interchanged and D(k,k) is a 1-by-1 diagonal block. */
141 /* If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and */
142 /* columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) */
143 /* is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) = */
144 /* IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were */
145 /* interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block. */
146
147 /* If FACT = 'N', then IPIV is an output argument and on exit */
148 /* contains details of the interchanges and the block structure */
149 /* of D, as determined by DSPTRF. */
150
151 /* B (input) DOUBLE PRECISION array, dimension (LDB,NRHS) */
152 /* The N-by-NRHS right hand side matrix B. */
153
154 /* LDB (input) INTEGER */
155 /* The leading dimension of the array B. LDB >= max(1,N). */
156
157 /* X (output) DOUBLE PRECISION array, dimension (LDX,NRHS) */
158 /* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X. */
159
160 /* LDX (input) INTEGER */
161 /* The leading dimension of the array X. LDX >= max(1,N). */
162
163 /* RCOND (output) DOUBLE PRECISION */
164 /* The estimate of the reciprocal condition number of the matrix */
165 /* A. If RCOND is less than the machine precision (in */
166 /* particular, if RCOND = 0), the matrix is singular to working */
167 /* precision. This condition is indicated by a return code of */
168 /* INFO > 0. */
169
170 /* FERR (output) DOUBLE PRECISION array, dimension (NRHS) */
171 /* The estimated forward error bound for each solution vector */
172 /* X(j) (the j-th column of the solution matrix X). */
173 /* If XTRUE is the true solution corresponding to X(j), FERR(j) */
174 /* is an estimated upper bound for the magnitude of the largest */
175 /* element in (X(j) - XTRUE) divided by the magnitude of the */
176 /* largest element in X(j). The estimate is as reliable as */
177 /* the estimate for RCOND, and is almost always a slight */
178 /* overestimate of the true error. */
179
180 /* BERR (output) DOUBLE PRECISION array, dimension (NRHS) */
181 /* The componentwise relative backward error of each solution */
182 /* vector X(j) (i.e., the smallest relative change in */
183 /* any element of A or B that makes X(j) an exact solution). */
184
185 /* WORK (workspace) DOUBLE PRECISION array, dimension (3*N) */
186
187 /* IWORK (workspace) INTEGER array, dimension (N) */
188
189 /* INFO (output) INTEGER */
190 /* = 0: successful exit */
191 /* < 0: if INFO = -i, the i-th argument had an illegal value */
192 /* > 0: if INFO = i, and i is */
193 /* <= N: D(i,i) is exactly zero. The factorization */
194 /* has been completed but the factor D is exactly */
195 /* singular, so the solution and error bounds could */
196 /* not be computed. RCOND = 0 is returned. */
197 /* = N+1: D is nonsingular, but RCOND is less than machine */
198 /* precision, meaning that the matrix is singular */
199 /* to working precision. Nevertheless, the */
200 /* solution and error bounds are computed because */
201 /* there are a number of situations where the */
202 /* computed solution can be more accurate than the */
203 /* value of RCOND would suggest. */
204
205 /* Further Details */
206 /* =============== */
207
208 /* The packed storage scheme is illustrated by the following example */
209 /* when N = 4, UPLO = 'U': */
210
211 /* Two-dimensional storage of the symmetric matrix A: */
212
213 /* a11 a12 a13 a14 */
214 /* a22 a23 a24 */
215 /* a33 a34 (aij = aji) */
216 /* a44 */
217
218 /* Packed storage of the upper triangle of A: */
219
220 /* AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ] */
221
222 /* ===================================================================== */
223
224 /* .. Parameters .. */
225 /* .. */
226 /* .. Local Scalars .. */
227 /* .. */
228 /* .. External Functions .. */
229 /* .. */
230 /* .. External Subroutines .. */
231 /* .. */
232 /* .. Intrinsic Functions .. */
233 /* .. */
234 /* .. Executable Statements .. */
235
236 /* Test the input parameters. */
237
238 /* Parameter adjustments */
239 --ap;
240 --afp;
241 --ipiv;
242 b_dim1 = *ldb;
243 b_offset = 1 + b_dim1;
244 b -= b_offset;
245 x_dim1 = *ldx;
246 x_offset = 1 + x_dim1;
247 x -= x_offset;
248 --ferr;
249 --berr;
250 --work;
251 --iwork;
252
253 /* Function Body */
254 *info = 0;
255 nofact = lsame_(fact, "N", (ftnlen)1, (ftnlen)1);
256 if (! nofact && ! lsame_(fact, "F", (ftnlen)1, (ftnlen)1)) {
257 *info = -1;
258 } else if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo,
259 "L", (ftnlen)1, (ftnlen)1)) {
260 *info = -2;
261 } else if (*n < 0) {
262 *info = -3;
263 } else if (*nrhs < 0) {
264 *info = -4;
265 } else if (*ldb < max(1,*n)) {
266 *info = -9;
267 } else if (*ldx < max(1,*n)) {
268 *info = -11;
269 }
270 if (*info != 0) {
271 i__1 = -(*info);
272 xerbla_("DSPSVX", &i__1, (ftnlen)6);
273 return 0;
274 }
275
276 if (nofact) {
277
278 /* Compute the factorization A = U*D*U' or A = L*D*L'. */
279
280 i__1 = *n * (*n + 1) / 2;
281 dcopy_(&i__1, &ap[1], &c__1, &afp[1], &c__1);
282 dsptrf_(uplo, n, &afp[1], &ipiv[1], info, (ftnlen)1);
283
284 /* Return if INFO is non-zero. */
285
286 if (*info != 0) {
287 if (*info > 0) {
288 *rcond = 0.;
289 }
290 return 0;
291 }
292 }
293
294 /* Compute the norm of the matrix A. */
295
296 anorm = dlansp_("I", uplo, n, &ap[1], &work[1], (ftnlen)1, (ftnlen)1);
297
298 /* Compute the reciprocal of the condition number of A. */
299
300 dspcon_(uplo, n, &afp[1], &ipiv[1], &anorm, rcond, &work[1], &iwork[1],
301 info, (ftnlen)1);
302
303 /* Set INFO = N+1 if the matrix is singular to working precision. */
304
305 if (*rcond < dlamch_("Epsilon", (ftnlen)7)) {
306 *info = *n + 1;
307 }
308
309 /* Compute the solution vectors X. */
310
311 dlacpy_("Full", n, nrhs, &b[b_offset], ldb, &x[x_offset], ldx, (ftnlen)4);
312 dsptrs_(uplo, n, nrhs, &afp[1], &ipiv[1], &x[x_offset], ldx, info, (
313 ftnlen)1);
314
315 /* Use iterative refinement to improve the computed solutions and */
316 /* compute error bounds and backward error estimates for them. */
317
318 dsprfs_(uplo, n, nrhs, &ap[1], &afp[1], &ipiv[1], &b[b_offset], ldb, &x[
319 x_offset], ldx, &ferr[1], &berr[1], &work[1], &iwork[1], info, (
320 ftnlen)1);
321
322 return 0;
323
324 /* End of DSPSVX */
325
326 } /* dspsvx_ */
327
328