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