1 /* ./src_f77/sporfs.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 static real c_b12 = -1.f;
12 static real c_b14 = 1.f;
13
sporfs_(char * uplo,integer * n,integer * nrhs,real * a,integer * lda,real * af,integer * ldaf,real * b,integer * ldb,real * x,integer * ldx,real * ferr,real * berr,real * work,integer * iwork,integer * info,ftnlen uplo_len)14 /* Subroutine */ int sporfs_(char *uplo, integer *n, integer *nrhs, real *a,
15 integer *lda, real *af, integer *ldaf, real *b, integer *ldb, real *x,
16 integer *ldx, real *ferr, real *berr, real *work, integer *iwork,
17 integer *info, ftnlen uplo_len)
18 {
19 /* System generated locals */
20 integer a_dim1, a_offset, af_dim1, af_offset, b_dim1, b_offset, x_dim1,
21 x_offset, i__1, i__2, i__3;
22 real r__1, r__2, r__3;
23
24 /* Local variables */
25 static integer i__, j, k;
26 static real s, xk;
27 static integer nz;
28 static real eps;
29 static integer kase;
30 static real safe1, safe2;
31 extern logical lsame_(char *, char *, ftnlen, ftnlen);
32 static integer count;
33 static logical upper;
34 extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *,
35 integer *), saxpy_(integer *, real *, real *, integer *, real *,
36 integer *), ssymv_(char *, integer *, real *, real *, integer *,
37 real *, integer *, real *, real *, integer *, ftnlen);
38 extern doublereal slamch_(char *, ftnlen);
39 static real safmin;
40 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen), slacon_(
41 integer *, real *, real *, integer *, real *, integer *);
42 static real lstres;
43 extern /* Subroutine */ int spotrs_(char *, integer *, integer *, real *,
44 integer *, real *, integer *, integer *, ftnlen);
45
46
47 /* -- LAPACK routine (version 3.0) -- */
48 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
49 /* Courant Institute, Argonne National Lab, and Rice University */
50 /* September 30, 1994 */
51
52 /* .. Scalar Arguments .. */
53 /* .. */
54 /* .. Array Arguments .. */
55 /* .. */
56
57 /* Purpose */
58 /* ======= */
59
60 /* SPORFS improves the computed solution to a system of linear */
61 /* equations when the coefficient matrix is symmetric positive definite, */
62 /* and provides error bounds and backward error estimates for the */
63 /* solution. */
64
65 /* Arguments */
66 /* ========= */
67
68 /* UPLO (input) CHARACTER*1 */
69 /* = 'U': Upper triangle of A is stored; */
70 /* = 'L': Lower triangle of A is stored. */
71
72 /* N (input) INTEGER */
73 /* The order of the matrix A. N >= 0. */
74
75 /* NRHS (input) INTEGER */
76 /* The number of right hand sides, i.e., the number of columns */
77 /* of the matrices B and X. NRHS >= 0. */
78
79 /* A (input) REAL array, dimension (LDA,N) */
80 /* The symmetric matrix A. If UPLO = 'U', the leading N-by-N */
81 /* upper triangular part of A contains the upper triangular part */
82 /* of the matrix A, and the strictly lower triangular part of A */
83 /* is not referenced. If UPLO = 'L', the leading N-by-N lower */
84 /* triangular part of A contains the lower triangular part of */
85 /* the matrix A, and the strictly upper triangular part of A is */
86 /* not referenced. */
87
88 /* LDA (input) INTEGER */
89 /* The leading dimension of the array A. LDA >= max(1,N). */
90
91 /* AF (input) REAL array, dimension (LDAF,N) */
92 /* The triangular factor U or L from the Cholesky factorization */
93 /* A = U**T*U or A = L*L**T, as computed by SPOTRF. */
94
95 /* LDAF (input) INTEGER */
96 /* The leading dimension of the array AF. LDAF >= max(1,N). */
97
98 /* B (input) REAL array, dimension (LDB,NRHS) */
99 /* The right hand side matrix B. */
100
101 /* LDB (input) INTEGER */
102 /* The leading dimension of the array B. LDB >= max(1,N). */
103
104 /* X (input/output) REAL array, dimension (LDX,NRHS) */
105 /* On entry, the solution matrix X, as computed by SPOTRS. */
106 /* On exit, the improved solution matrix X. */
107
108 /* LDX (input) INTEGER */
109 /* The leading dimension of the array X. LDX >= max(1,N). */
110
111 /* FERR (output) REAL array, dimension (NRHS) */
112 /* The estimated forward error bound for each solution vector */
113 /* X(j) (the j-th column of the solution matrix X). */
114 /* If XTRUE is the true solution corresponding to X(j), FERR(j) */
115 /* is an estimated upper bound for the magnitude of the largest */
116 /* element in (X(j) - XTRUE) divided by the magnitude of the */
117 /* largest element in X(j). The estimate is as reliable as */
118 /* the estimate for RCOND, and is almost always a slight */
119 /* overestimate of the true error. */
120
121 /* BERR (output) REAL array, dimension (NRHS) */
122 /* The componentwise relative backward error of each solution */
123 /* vector X(j) (i.e., the smallest relative change in */
124 /* any element of A or B that makes X(j) an exact solution). */
125
126 /* WORK (workspace) REAL array, dimension (3*N) */
127
128 /* IWORK (workspace) INTEGER array, dimension (N) */
129
130 /* INFO (output) INTEGER */
131 /* = 0: successful exit */
132 /* < 0: if INFO = -i, the i-th argument had an illegal value */
133
134 /* Internal Parameters */
135 /* =================== */
136
137 /* ITMAX is the maximum number of steps of iterative refinement. */
138
139 /* ===================================================================== */
140
141 /* .. Parameters .. */
142 /* .. */
143 /* .. Local Scalars .. */
144 /* .. */
145 /* .. External Subroutines .. */
146 /* .. */
147 /* .. Intrinsic Functions .. */
148 /* .. */
149 /* .. External Functions .. */
150 /* .. */
151 /* .. Executable Statements .. */
152
153 /* Test the input parameters. */
154
155 /* Parameter adjustments */
156 a_dim1 = *lda;
157 a_offset = 1 + a_dim1;
158 a -= a_offset;
159 af_dim1 = *ldaf;
160 af_offset = 1 + af_dim1;
161 af -= af_offset;
162 b_dim1 = *ldb;
163 b_offset = 1 + b_dim1;
164 b -= b_offset;
165 x_dim1 = *ldx;
166 x_offset = 1 + x_dim1;
167 x -= x_offset;
168 --ferr;
169 --berr;
170 --work;
171 --iwork;
172
173 /* Function Body */
174 *info = 0;
175 upper = lsame_(uplo, "U", (ftnlen)1, (ftnlen)1);
176 if (! upper && ! lsame_(uplo, "L", (ftnlen)1, (ftnlen)1)) {
177 *info = -1;
178 } else if (*n < 0) {
179 *info = -2;
180 } else if (*nrhs < 0) {
181 *info = -3;
182 } else if (*lda < max(1,*n)) {
183 *info = -5;
184 } else if (*ldaf < max(1,*n)) {
185 *info = -7;
186 } else if (*ldb < max(1,*n)) {
187 *info = -9;
188 } else if (*ldx < max(1,*n)) {
189 *info = -11;
190 }
191 if (*info != 0) {
192 i__1 = -(*info);
193 xerbla_("SPORFS", &i__1, (ftnlen)6);
194 return 0;
195 }
196
197 /* Quick return if possible */
198
199 if (*n == 0 || *nrhs == 0) {
200 i__1 = *nrhs;
201 for (j = 1; j <= i__1; ++j) {
202 ferr[j] = 0.f;
203 berr[j] = 0.f;
204 /* L10: */
205 }
206 return 0;
207 }
208
209 /* NZ = maximum number of nonzero elements in each row of A, plus 1 */
210
211 nz = *n + 1;
212 eps = slamch_("Epsilon", (ftnlen)7);
213 safmin = slamch_("Safe minimum", (ftnlen)12);
214 safe1 = nz * safmin;
215 safe2 = safe1 / eps;
216
217 /* Do for each right hand side */
218
219 i__1 = *nrhs;
220 for (j = 1; j <= i__1; ++j) {
221
222 count = 1;
223 lstres = 3.f;
224 L20:
225
226 /* Loop until stopping criterion is satisfied. */
227
228 /* Compute residual R = B - A * X */
229
230 scopy_(n, &b[j * b_dim1 + 1], &c__1, &work[*n + 1], &c__1);
231 ssymv_(uplo, n, &c_b12, &a[a_offset], lda, &x[j * x_dim1 + 1], &c__1,
232 &c_b14, &work[*n + 1], &c__1, (ftnlen)1);
233
234 /* Compute componentwise relative backward error from formula */
235
236 /* max(i) ( abs(R(i)) / ( abs(A)*abs(X) + abs(B) )(i) ) */
237
238 /* where abs(Z) is the componentwise absolute value of the matrix */
239 /* or vector Z. If the i-th component of the denominator is less */
240 /* than SAFE2, then SAFE1 is added to the i-th components of the */
241 /* numerator and denominator before dividing. */
242
243 i__2 = *n;
244 for (i__ = 1; i__ <= i__2; ++i__) {
245 work[i__] = (r__1 = b[i__ + j * b_dim1], dabs(r__1));
246 /* L30: */
247 }
248
249 /* Compute abs(A)*abs(X) + abs(B). */
250
251 if (upper) {
252 i__2 = *n;
253 for (k = 1; k <= i__2; ++k) {
254 s = 0.f;
255 xk = (r__1 = x[k + j * x_dim1], dabs(r__1));
256 i__3 = k - 1;
257 for (i__ = 1; i__ <= i__3; ++i__) {
258 work[i__] += (r__1 = a[i__ + k * a_dim1], dabs(r__1)) *
259 xk;
260 s += (r__1 = a[i__ + k * a_dim1], dabs(r__1)) * (r__2 = x[
261 i__ + j * x_dim1], dabs(r__2));
262 /* L40: */
263 }
264 work[k] = work[k] + (r__1 = a[k + k * a_dim1], dabs(r__1)) *
265 xk + s;
266 /* L50: */
267 }
268 } else {
269 i__2 = *n;
270 for (k = 1; k <= i__2; ++k) {
271 s = 0.f;
272 xk = (r__1 = x[k + j * x_dim1], dabs(r__1));
273 work[k] += (r__1 = a[k + k * a_dim1], dabs(r__1)) * xk;
274 i__3 = *n;
275 for (i__ = k + 1; i__ <= i__3; ++i__) {
276 work[i__] += (r__1 = a[i__ + k * a_dim1], dabs(r__1)) *
277 xk;
278 s += (r__1 = a[i__ + k * a_dim1], dabs(r__1)) * (r__2 = x[
279 i__ + j * x_dim1], dabs(r__2));
280 /* L60: */
281 }
282 work[k] += s;
283 /* L70: */
284 }
285 }
286 s = 0.f;
287 i__2 = *n;
288 for (i__ = 1; i__ <= i__2; ++i__) {
289 if (work[i__] > safe2) {
290 /* Computing MAX */
291 r__2 = s, r__3 = (r__1 = work[*n + i__], dabs(r__1)) / work[
292 i__];
293 s = dmax(r__2,r__3);
294 } else {
295 /* Computing MAX */
296 r__2 = s, r__3 = ((r__1 = work[*n + i__], dabs(r__1)) + safe1)
297 / (work[i__] + safe1);
298 s = dmax(r__2,r__3);
299 }
300 /* L80: */
301 }
302 berr[j] = s;
303
304 /* Test stopping criterion. Continue iterating if */
305 /* 1) The residual BERR(J) is larger than machine epsilon, and */
306 /* 2) BERR(J) decreased by at least a factor of 2 during the */
307 /* last iteration, and */
308 /* 3) At most ITMAX iterations tried. */
309
310 if (berr[j] > eps && berr[j] * 2.f <= lstres && count <= 5) {
311
312 /* Update solution and try again. */
313
314 spotrs_(uplo, n, &c__1, &af[af_offset], ldaf, &work[*n + 1], n,
315 info, (ftnlen)1);
316 saxpy_(n, &c_b14, &work[*n + 1], &c__1, &x[j * x_dim1 + 1], &c__1)
317 ;
318 lstres = berr[j];
319 ++count;
320 goto L20;
321 }
322
323 /* Bound error from formula */
324
325 /* norm(X - XTRUE) / norm(X) .le. FERR = */
326 /* norm( abs(inv(A))* */
327 /* ( abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) ))) / norm(X) */
328
329 /* where */
330 /* norm(Z) is the magnitude of the largest component of Z */
331 /* inv(A) is the inverse of A */
332 /* abs(Z) is the componentwise absolute value of the matrix or */
333 /* vector Z */
334 /* NZ is the maximum number of nonzeros in any row of A, plus 1 */
335 /* EPS is machine epsilon */
336
337 /* The i-th component of abs(R)+NZ*EPS*(abs(A)*abs(X)+abs(B)) */
338 /* is incremented by SAFE1 if the i-th component of */
339 /* abs(A)*abs(X) + abs(B) is less than SAFE2. */
340
341 /* Use SLACON to estimate the infinity-norm of the matrix */
342 /* inv(A) * diag(W), */
343 /* where W = abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) ))) */
344
345 i__2 = *n;
346 for (i__ = 1; i__ <= i__2; ++i__) {
347 if (work[i__] > safe2) {
348 work[i__] = (r__1 = work[*n + i__], dabs(r__1)) + nz * eps *
349 work[i__];
350 } else {
351 work[i__] = (r__1 = work[*n + i__], dabs(r__1)) + nz * eps *
352 work[i__] + safe1;
353 }
354 /* L90: */
355 }
356
357 kase = 0;
358 L100:
359 slacon_(n, &work[(*n << 1) + 1], &work[*n + 1], &iwork[1], &ferr[j], &
360 kase);
361 if (kase != 0) {
362 if (kase == 1) {
363
364 /* Multiply by diag(W)*inv(A'). */
365
366 spotrs_(uplo, n, &c__1, &af[af_offset], ldaf, &work[*n + 1],
367 n, info, (ftnlen)1);
368 i__2 = *n;
369 for (i__ = 1; i__ <= i__2; ++i__) {
370 work[*n + i__] = work[i__] * work[*n + i__];
371 /* L110: */
372 }
373 } else if (kase == 2) {
374
375 /* Multiply by inv(A)*diag(W). */
376
377 i__2 = *n;
378 for (i__ = 1; i__ <= i__2; ++i__) {
379 work[*n + i__] = work[i__] * work[*n + i__];
380 /* L120: */
381 }
382 spotrs_(uplo, n, &c__1, &af[af_offset], ldaf, &work[*n + 1],
383 n, info, (ftnlen)1);
384 }
385 goto L100;
386 }
387
388 /* Normalize error. */
389
390 lstres = 0.f;
391 i__2 = *n;
392 for (i__ = 1; i__ <= i__2; ++i__) {
393 /* Computing MAX */
394 r__2 = lstres, r__3 = (r__1 = x[i__ + j * x_dim1], dabs(r__1));
395 lstres = dmax(r__2,r__3);
396 /* L130: */
397 }
398 if (lstres != 0.f) {
399 ferr[j] /= lstres;
400 }
401
402 /* L140: */
403 }
404
405 return 0;
406
407 /* End of SPORFS */
408
409 } /* sporfs_ */
410
411