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