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