1 /* ./src_f77/sposvx.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
sposvx_(char * fact,char * uplo,integer * n,integer * nrhs,real * a,integer * lda,real * af,integer * ldaf,char * equed,real * s,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,ftnlen equed_len)8 /* Subroutine */ int sposvx_(char *fact, char *uplo, integer *n, integer *
9 nrhs, real *a, integer *lda, real *af, integer *ldaf, char *equed,
10 real *s, real *b, integer *ldb, real *x, integer *ldx, real *rcond,
11 real *ferr, real *berr, real *work, integer *iwork, integer *info,
12 ftnlen fact_len, ftnlen uplo_len, ftnlen equed_len)
13 {
14 /* System generated locals */
15 integer a_dim1, a_offset, af_dim1, af_offset, b_dim1, b_offset, x_dim1,
16 x_offset, i__1, i__2;
17 real r__1, r__2;
18
19 /* Local variables */
20 static integer i__, j;
21 static real amax, smin, smax;
22 extern logical lsame_(char *, char *, ftnlen, ftnlen);
23 static real scond, anorm;
24 static logical equil, rcequ;
25 extern doublereal slamch_(char *, ftnlen);
26 static logical nofact;
27 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
28 static real bignum;
29 static integer infequ;
30 extern /* Subroutine */ int slacpy_(char *, integer *, integer *, real *,
31 integer *, real *, integer *, ftnlen), spocon_(char *, integer *,
32 real *, integer *, real *, real *, real *, integer *, integer *,
33 ftnlen);
34 extern doublereal slansy_(char *, char *, integer *, real *, integer *,
35 real *, ftnlen, ftnlen);
36 static real smlnum;
37 extern /* Subroutine */ int slaqsy_(char *, integer *, real *, integer *,
38 real *, real *, real *, char *, ftnlen, ftnlen), spoequ_(integer *
39 , real *, integer *, real *, real *, real *, integer *), sporfs_(
40 char *, integer *, integer *, real *, integer *, real *, integer *
41 , real *, integer *, real *, integer *, real *, real *, real *,
42 integer *, integer *, ftnlen), spotrf_(char *, integer *, real *,
43 integer *, integer *, ftnlen), spotrs_(char *, integer *, integer
44 *, real *, integer *, real *, integer *, integer *, ftnlen);
45
46
47 /* -- LAPACK driver routine (version 3.0) -- */
48 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
49 /* Courant Institute, Argonne National Lab, and Rice University */
50 /* June 30, 1999 */
51
52 /* .. Scalar Arguments .. */
53 /* .. */
54 /* .. Array Arguments .. */
55 /* .. */
56
57 /* Purpose */
58 /* ======= */
59
60 /* SPOSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to */
61 /* compute the solution to a real system of linear equations */
62 /* A * X = B, */
63 /* where A is an N-by-N symmetric positive definite matrix and X and B */
64 /* are N-by-NRHS matrices. */
65
66 /* Error bounds on the solution and a condition estimate are also */
67 /* provided. */
68
69 /* Description */
70 /* =========== */
71
72 /* The following steps are performed: */
73
74 /* 1. If FACT = 'E', real scaling factors are computed to equilibrate */
75 /* the system: */
76 /* diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B */
77 /* Whether or not the system will be equilibrated depends on the */
78 /* scaling of the matrix A, but if equilibration is used, A is */
79 /* overwritten by diag(S)*A*diag(S) and B by diag(S)*B. */
80
81 /* 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to */
82 /* factor the matrix A (after equilibration if FACT = 'E') as */
83 /* A = U**T* U, if UPLO = 'U', or */
84 /* A = L * L**T, if UPLO = 'L', */
85 /* where U is an upper triangular matrix and L is a lower triangular */
86 /* matrix. */
87
88 /* 3. If the leading i-by-i principal minor is not positive definite, */
89 /* then the routine returns with INFO = i. Otherwise, the factored */
90 /* form of A is used to estimate the condition number of the matrix */
91 /* A. If the reciprocal of the condition number is less than machine */
92 /* precision, INFO = N+1 is returned as a warning, but the routine */
93 /* still goes on to solve for X and compute error bounds as */
94 /* described below. */
95
96 /* 4. The system of equations is solved for X using the factored form */
97 /* of A. */
98
99 /* 5. Iterative refinement is applied to improve the computed solution */
100 /* matrix and calculate error bounds and backward error estimates */
101 /* for it. */
102
103 /* 6. If equilibration was used, the matrix X is premultiplied by */
104 /* diag(S) so that it solves the original system before */
105 /* equilibration. */
106
107 /* Arguments */
108 /* ========= */
109
110 /* FACT (input) CHARACTER*1 */
111 /* Specifies whether or not the factored form of the matrix A is */
112 /* supplied on entry, and if not, whether the matrix A should be */
113 /* equilibrated before it is factored. */
114 /* = 'F': On entry, AF contains the factored form of A. */
115 /* If EQUED = 'Y', the matrix A has been equilibrated */
116 /* with scaling factors given by S. A and AF will not */
117 /* be modified. */
118 /* = 'N': The matrix A will be copied to AF and factored. */
119 /* = 'E': The matrix A will be equilibrated if necessary, then */
120 /* copied to AF and factored. */
121
122 /* UPLO (input) CHARACTER*1 */
123 /* = 'U': Upper triangle of A is stored; */
124 /* = 'L': Lower triangle of A is stored. */
125
126 /* N (input) INTEGER */
127 /* The number of linear equations, i.e., the order of the */
128 /* matrix A. N >= 0. */
129
130 /* NRHS (input) INTEGER */
131 /* The number of right hand sides, i.e., the number of columns */
132 /* of the matrices B and X. NRHS >= 0. */
133
134 /* A (input/output) REAL array, dimension (LDA,N) */
135 /* On entry, the symmetric matrix A, except if FACT = 'F' and */
136 /* EQUED = 'Y', then A must contain the equilibrated matrix */
137 /* diag(S)*A*diag(S). If UPLO = 'U', the leading */
138 /* N-by-N upper triangular part of A contains the upper */
139 /* triangular part of the matrix A, and the strictly lower */
140 /* triangular part of A is not referenced. If UPLO = 'L', the */
141 /* leading N-by-N lower triangular part of A contains the lower */
142 /* triangular part of the matrix A, and the strictly upper */
143 /* triangular part of A is not referenced. A is not modified if */
144 /* FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit. */
145
146 /* On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by */
147 /* diag(S)*A*diag(S). */
148
149 /* LDA (input) INTEGER */
150 /* The leading dimension of the array A. LDA >= max(1,N). */
151
152 /* AF (input or output) REAL array, dimension (LDAF,N) */
153 /* If FACT = 'F', then AF is an input argument and on entry */
154 /* contains the triangular factor U or L from the Cholesky */
155 /* factorization A = U**T*U or A = L*L**T, in the same storage */
156 /* format as A. If EQUED .ne. 'N', then AF is the factored form */
157 /* of the equilibrated matrix diag(S)*A*diag(S). */
158
159 /* If FACT = 'N', then AF is an output argument and on exit */
160 /* returns the triangular factor U or L from the Cholesky */
161 /* factorization A = U**T*U or A = L*L**T of the original */
162 /* matrix A. */
163
164 /* If FACT = 'E', then AF is an output argument and on exit */
165 /* returns the triangular factor U or L from the Cholesky */
166 /* factorization A = U**T*U or A = L*L**T of the equilibrated */
167 /* matrix A (see the description of A for the form of the */
168 /* equilibrated matrix). */
169
170 /* LDAF (input) INTEGER */
171 /* The leading dimension of the array AF. LDAF >= max(1,N). */
172
173 /* EQUED (input or output) CHARACTER*1 */
174 /* Specifies the form of equilibration that was done. */
175 /* = 'N': No equilibration (always true if FACT = 'N'). */
176 /* = 'Y': Equilibration was done, i.e., A has been replaced by */
177 /* diag(S) * A * diag(S). */
178 /* EQUED is an input argument if FACT = 'F'; otherwise, it is an */
179 /* output argument. */
180
181 /* S (input or output) REAL array, dimension (N) */
182 /* The scale factors for A; not accessed if EQUED = 'N'. S is */
183 /* an input argument if FACT = 'F'; otherwise, S is an output */
184 /* argument. If FACT = 'F' and EQUED = 'Y', each element of S */
185 /* must be positive. */
186
187 /* B (input/output) REAL array, dimension (LDB,NRHS) */
188 /* On entry, the N-by-NRHS right hand side matrix B. */
189 /* On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y', */
190 /* B is overwritten by diag(S) * B. */
191
192 /* LDB (input) INTEGER */
193 /* The leading dimension of the array B. LDB >= max(1,N). */
194
195 /* X (output) REAL array, dimension (LDX,NRHS) */
196 /* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to */
197 /* the original system of equations. Note that if EQUED = 'Y', */
198 /* A and B are modified on exit, and the solution to the */
199 /* equilibrated system is inv(diag(S))*X. */
200
201 /* LDX (input) INTEGER */
202 /* The leading dimension of the array X. LDX >= max(1,N). */
203
204 /* RCOND (output) REAL */
205 /* The estimate of the reciprocal condition number of the matrix */
206 /* A after equilibration (if done). If RCOND is less than the */
207 /* machine precision (in particular, if RCOND = 0), the matrix */
208 /* is singular to working precision. This condition is */
209 /* indicated by a return code of INFO > 0. */
210
211 /* FERR (output) REAL array, dimension (NRHS) */
212 /* The estimated forward error bound for each solution vector */
213 /* X(j) (the j-th column of the solution matrix X). */
214 /* If XTRUE is the true solution corresponding to X(j), FERR(j) */
215 /* is an estimated upper bound for the magnitude of the largest */
216 /* element in (X(j) - XTRUE) divided by the magnitude of the */
217 /* largest element in X(j). The estimate is as reliable as */
218 /* the estimate for RCOND, and is almost always a slight */
219 /* overestimate of the true error. */
220
221 /* BERR (output) REAL array, dimension (NRHS) */
222 /* The componentwise relative backward error of each solution */
223 /* vector X(j) (i.e., the smallest relative change in */
224 /* any element of A or B that makes X(j) an exact solution). */
225
226 /* WORK (workspace) REAL array, dimension (3*N) */
227
228 /* IWORK (workspace) INTEGER array, dimension (N) */
229
230 /* INFO (output) INTEGER */
231 /* = 0: successful exit */
232 /* < 0: if INFO = -i, the i-th argument had an illegal value */
233 /* > 0: if INFO = i, and i is */
234 /* <= N: the leading minor of order i of A is */
235 /* not positive definite, so the factorization */
236 /* could not be completed, and the solution has not */
237 /* been computed. RCOND = 0 is returned. */
238 /* = N+1: U is nonsingular, but RCOND is less than machine */
239 /* precision, meaning that the matrix is singular */
240 /* to working precision. Nevertheless, the */
241 /* solution and error bounds are computed because */
242 /* there are a number of situations where the */
243 /* computed solution can be more accurate than the */
244 /* value of RCOND would suggest. */
245
246 /* ===================================================================== */
247
248 /* .. Parameters .. */
249 /* .. */
250 /* .. Local Scalars .. */
251 /* .. */
252 /* .. External Functions .. */
253 /* .. */
254 /* .. External Subroutines .. */
255 /* .. */
256 /* .. Intrinsic Functions .. */
257 /* .. */
258 /* .. Executable Statements .. */
259
260 /* Parameter adjustments */
261 a_dim1 = *lda;
262 a_offset = 1 + a_dim1;
263 a -= a_offset;
264 af_dim1 = *ldaf;
265 af_offset = 1 + af_dim1;
266 af -= af_offset;
267 --s;
268 b_dim1 = *ldb;
269 b_offset = 1 + b_dim1;
270 b -= b_offset;
271 x_dim1 = *ldx;
272 x_offset = 1 + x_dim1;
273 x -= x_offset;
274 --ferr;
275 --berr;
276 --work;
277 --iwork;
278
279 /* Function Body */
280 *info = 0;
281 nofact = lsame_(fact, "N", (ftnlen)1, (ftnlen)1);
282 equil = lsame_(fact, "E", (ftnlen)1, (ftnlen)1);
283 if (nofact || equil) {
284 *(unsigned char *)equed = 'N';
285 rcequ = FALSE_;
286 } else {
287 rcequ = lsame_(equed, "Y", (ftnlen)1, (ftnlen)1);
288 smlnum = slamch_("Safe minimum", (ftnlen)12);
289 bignum = 1.f / smlnum;
290 }
291
292 /* Test the input parameters. */
293
294 if (! nofact && ! equil && ! lsame_(fact, "F", (ftnlen)1, (ftnlen)1)) {
295 *info = -1;
296 } else if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo,
297 "L", (ftnlen)1, (ftnlen)1)) {
298 *info = -2;
299 } else if (*n < 0) {
300 *info = -3;
301 } else if (*nrhs < 0) {
302 *info = -4;
303 } else if (*lda < max(1,*n)) {
304 *info = -6;
305 } else if (*ldaf < max(1,*n)) {
306 *info = -8;
307 } else if (lsame_(fact, "F", (ftnlen)1, (ftnlen)1) && ! (rcequ || lsame_(
308 equed, "N", (ftnlen)1, (ftnlen)1))) {
309 *info = -9;
310 } else {
311 if (rcequ) {
312 smin = bignum;
313 smax = 0.f;
314 i__1 = *n;
315 for (j = 1; j <= i__1; ++j) {
316 /* Computing MIN */
317 r__1 = smin, r__2 = s[j];
318 smin = dmin(r__1,r__2);
319 /* Computing MAX */
320 r__1 = smax, r__2 = s[j];
321 smax = dmax(r__1,r__2);
322 /* L10: */
323 }
324 if (smin <= 0.f) {
325 *info = -10;
326 } else if (*n > 0) {
327 scond = dmax(smin,smlnum) / dmin(smax,bignum);
328 } else {
329 scond = 1.f;
330 }
331 }
332 if (*info == 0) {
333 if (*ldb < max(1,*n)) {
334 *info = -12;
335 } else if (*ldx < max(1,*n)) {
336 *info = -14;
337 }
338 }
339 }
340
341 if (*info != 0) {
342 i__1 = -(*info);
343 xerbla_("SPOSVX", &i__1, (ftnlen)6);
344 return 0;
345 }
346
347 if (equil) {
348
349 /* Compute row and column scalings to equilibrate the matrix A. */
350
351 spoequ_(n, &a[a_offset], lda, &s[1], &scond, &amax, &infequ);
352 if (infequ == 0) {
353
354 /* Equilibrate the matrix. */
355
356 slaqsy_(uplo, n, &a[a_offset], lda, &s[1], &scond, &amax, equed, (
357 ftnlen)1, (ftnlen)1);
358 rcequ = lsame_(equed, "Y", (ftnlen)1, (ftnlen)1);
359 }
360 }
361
362 /* Scale the right hand side. */
363
364 if (rcequ) {
365 i__1 = *nrhs;
366 for (j = 1; j <= i__1; ++j) {
367 i__2 = *n;
368 for (i__ = 1; i__ <= i__2; ++i__) {
369 b[i__ + j * b_dim1] = s[i__] * b[i__ + j * b_dim1];
370 /* L20: */
371 }
372 /* L30: */
373 }
374 }
375
376 if (nofact || equil) {
377
378 /* Compute the Cholesky factorization A = U'*U or A = L*L'. */
379
380 slacpy_(uplo, n, n, &a[a_offset], lda, &af[af_offset], ldaf, (ftnlen)
381 1);
382 spotrf_(uplo, n, &af[af_offset], ldaf, info, (ftnlen)1);
383
384 /* Return if INFO is non-zero. */
385
386 if (*info != 0) {
387 if (*info > 0) {
388 *rcond = 0.f;
389 }
390 return 0;
391 }
392 }
393
394 /* Compute the norm of the matrix A. */
395
396 anorm = slansy_("1", uplo, n, &a[a_offset], lda, &work[1], (ftnlen)1, (
397 ftnlen)1);
398
399 /* Compute the reciprocal of the condition number of A. */
400
401 spocon_(uplo, n, &af[af_offset], ldaf, &anorm, rcond, &work[1], &iwork[1],
402 info, (ftnlen)1);
403
404 /* Set INFO = N+1 if the matrix is singular to working precision. */
405
406 if (*rcond < slamch_("Epsilon", (ftnlen)7)) {
407 *info = *n + 1;
408 }
409
410 /* Compute the solution matrix X. */
411
412 slacpy_("Full", n, nrhs, &b[b_offset], ldb, &x[x_offset], ldx, (ftnlen)4);
413 spotrs_(uplo, n, nrhs, &af[af_offset], ldaf, &x[x_offset], ldx, info, (
414 ftnlen)1);
415
416 /* Use iterative refinement to improve the computed solution and */
417 /* compute error bounds and backward error estimates for it. */
418
419 sporfs_(uplo, n, nrhs, &a[a_offset], lda, &af[af_offset], ldaf, &b[
420 b_offset], ldb, &x[x_offset], ldx, &ferr[1], &berr[1], &work[1], &
421 iwork[1], info, (ftnlen)1);
422
423 /* Transform the solution matrix X to a solution of the original */
424 /* system. */
425
426 if (rcequ) {
427 i__1 = *nrhs;
428 for (j = 1; j <= i__1; ++j) {
429 i__2 = *n;
430 for (i__ = 1; i__ <= i__2; ++i__) {
431 x[i__ + j * x_dim1] = s[i__] * x[i__ + j * x_dim1];
432 /* L40: */
433 }
434 /* L50: */
435 }
436 i__1 = *nrhs;
437 for (j = 1; j <= i__1; ++j) {
438 ferr[j] /= scond;
439 /* L60: */
440 }
441 }
442
443 return 0;
444
445 /* End of SPOSVX */
446
447 } /* sposvx_ */
448
449