1 /* ./src_f77/ctprfs.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
ctprfs_(char * uplo,char * trans,char * diag,integer * n,integer * nrhs,complex * ap,complex * b,integer * ldb,complex * x,integer * ldx,real * ferr,real * berr,complex * work,real * rwork,integer * info,ftnlen uplo_len,ftnlen trans_len,ftnlen diag_len)12 /* Subroutine */ int ctprfs_(char *uplo, char *trans, char *diag, integer *n,
13 integer *nrhs, complex *ap, complex *b, integer *ldb, complex *x,
14 integer *ldx, real *ferr, real *berr, complex *work, real *rwork,
15 integer *info, ftnlen uplo_len, ftnlen trans_len, ftnlen diag_len)
16 {
17 /* System generated locals */
18 integer b_dim1, b_offset, x_dim1, x_offset, i__1, i__2, i__3, i__4, i__5;
19 real r__1, r__2, r__3, r__4;
20 complex q__1;
21
22 /* Builtin functions */
23 double r_imag(complex *);
24
25 /* Local variables */
26 static integer i__, j, k;
27 static real s;
28 static integer kc;
29 static real 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 *), ctpmv_(char *, char *, char *,
38 integer *, complex *, complex *, integer *, ftnlen, ftnlen,
39 ftnlen);
40 static logical upper;
41 extern /* Subroutine */ int ctpsv_(char *, char *, char *, integer *,
42 complex *, complex *, integer *, ftnlen, ftnlen, ftnlen), clacon_(
43 integer *, complex *, complex *, real *, integer *);
44 extern doublereal slamch_(char *, ftnlen);
45 static real safmin;
46 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
47 static logical notran;
48 static char transn[1], transt[1];
49 static logical nounit;
50 static real lstres;
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 /* CTPRFS provides error bounds and backward error estimates for the */
67 /* solution to a system of linear equations with a triangular packed */
68 /* coefficient matrix. */
69
70 /* The solution matrix X must be computed by CTPTRS or some other */
71 /* means before entering this routine. CTPRFS does not do iterative */
72 /* refinement because doing so cannot improve the backward error. */
73
74 /* Arguments */
75 /* ========= */
76
77 /* UPLO (input) CHARACTER*1 */
78 /* = 'U': A is upper triangular; */
79 /* = 'L': A is lower triangular. */
80
81 /* TRANS (input) CHARACTER*1 */
82 /* Specifies the form of the system of equations: */
83 /* = 'N': A * X = B (No transpose) */
84 /* = 'T': A**T * X = B (Transpose) */
85 /* = 'C': A**H * X = B (Conjugate transpose) */
86
87 /* DIAG (input) CHARACTER*1 */
88 /* = 'N': A is non-unit triangular; */
89 /* = 'U': A is unit triangular. */
90
91 /* N (input) INTEGER */
92 /* The order of the matrix A. N >= 0. */
93
94 /* NRHS (input) INTEGER */
95 /* The number of right hand sides, i.e., the number of columns */
96 /* of the matrices B and X. NRHS >= 0. */
97
98 /* AP (input) COMPLEX array, dimension (N*(N+1)/2) */
99 /* The upper or lower triangular matrix A, packed columnwise in */
100 /* a linear array. The j-th column of A is stored in the array */
101 /* AP as follows: */
102 /* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; */
103 /* if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. */
104 /* If DIAG = 'U', the diagonal elements of A are not referenced */
105 /* and are assumed to be 1. */
106
107 /* B (input) COMPLEX array, dimension (LDB,NRHS) */
108 /* The right hand side matrix B. */
109
110 /* LDB (input) INTEGER */
111 /* The leading dimension of the array B. LDB >= max(1,N). */
112
113 /* X (input) COMPLEX array, dimension (LDX,NRHS) */
114 /* The solution matrix X. */
115
116 /* LDX (input) INTEGER */
117 /* The leading dimension of the array X. LDX >= max(1,N). */
118
119 /* FERR (output) REAL array, dimension (NRHS) */
120 /* The estimated forward error bound for each solution vector */
121 /* X(j) (the j-th column of the solution matrix X). */
122 /* If XTRUE is the true solution corresponding to X(j), FERR(j) */
123 /* is an estimated upper bound for the magnitude of the largest */
124 /* element in (X(j) - XTRUE) divided by the magnitude of the */
125 /* largest element in X(j). The estimate is as reliable as */
126 /* the estimate for RCOND, and is almost always a slight */
127 /* overestimate of the true error. */
128
129 /* BERR (output) REAL array, dimension (NRHS) */
130 /* The componentwise relative backward error of each solution */
131 /* vector X(j) (i.e., the smallest relative change in */
132 /* any element of A or B that makes X(j) an exact solution). */
133
134 /* WORK (workspace) COMPLEX array, dimension (2*N) */
135
136 /* RWORK (workspace) REAL array, dimension (N) */
137
138 /* INFO (output) INTEGER */
139 /* = 0: successful exit */
140 /* < 0: if INFO = -i, the i-th argument had an illegal value */
141
142 /* ===================================================================== */
143
144 /* .. Parameters .. */
145 /* .. */
146 /* .. Local Scalars .. */
147 /* .. */
148 /* .. External Subroutines .. */
149 /* .. */
150 /* .. Intrinsic Functions .. */
151 /* .. */
152 /* .. External Functions .. */
153 /* .. */
154 /* .. Statement Functions .. */
155 /* .. */
156 /* .. Statement Function definitions .. */
157 /* .. */
158 /* .. Executable Statements .. */
159
160 /* Test the input parameters. */
161
162 /* Parameter adjustments */
163 --ap;
164 b_dim1 = *ldb;
165 b_offset = 1 + b_dim1;
166 b -= b_offset;
167 x_dim1 = *ldx;
168 x_offset = 1 + x_dim1;
169 x -= x_offset;
170 --ferr;
171 --berr;
172 --work;
173 --rwork;
174
175 /* Function Body */
176 *info = 0;
177 upper = lsame_(uplo, "U", (ftnlen)1, (ftnlen)1);
178 notran = lsame_(trans, "N", (ftnlen)1, (ftnlen)1);
179 nounit = lsame_(diag, "N", (ftnlen)1, (ftnlen)1);
180
181 if (! upper && ! lsame_(uplo, "L", (ftnlen)1, (ftnlen)1)) {
182 *info = -1;
183 } else if (! notran && ! lsame_(trans, "T", (ftnlen)1, (ftnlen)1) && !
184 lsame_(trans, "C", (ftnlen)1, (ftnlen)1)) {
185 *info = -2;
186 } else if (! nounit && ! lsame_(diag, "U", (ftnlen)1, (ftnlen)1)) {
187 *info = -3;
188 } else if (*n < 0) {
189 *info = -4;
190 } else if (*nrhs < 0) {
191 *info = -5;
192 } else if (*ldb < max(1,*n)) {
193 *info = -8;
194 } else if (*ldx < max(1,*n)) {
195 *info = -10;
196 }
197 if (*info != 0) {
198 i__1 = -(*info);
199 xerbla_("CTPRFS", &i__1, (ftnlen)6);
200 return 0;
201 }
202
203 /* Quick return if possible */
204
205 if (*n == 0 || *nrhs == 0) {
206 i__1 = *nrhs;
207 for (j = 1; j <= i__1; ++j) {
208 ferr[j] = 0.f;
209 berr[j] = 0.f;
210 /* L10: */
211 }
212 return 0;
213 }
214
215 if (notran) {
216 *(unsigned char *)transn = 'N';
217 *(unsigned char *)transt = 'C';
218 } else {
219 *(unsigned char *)transn = 'C';
220 *(unsigned char *)transt = 'N';
221 }
222
223 /* NZ = maximum number of nonzero elements in each row of A, plus 1 */
224
225 nz = *n + 1;
226 eps = slamch_("Epsilon", (ftnlen)7);
227 safmin = slamch_("Safe minimum", (ftnlen)12);
228 safe1 = nz * safmin;
229 safe2 = safe1 / eps;
230
231 /* Do for each right hand side */
232
233 i__1 = *nrhs;
234 for (j = 1; j <= i__1; ++j) {
235
236 /* Compute residual R = B - op(A) * X, */
237 /* where op(A) = A, A**T, or A**H, depending on TRANS. */
238
239 ccopy_(n, &x[j * x_dim1 + 1], &c__1, &work[1], &c__1);
240 ctpmv_(uplo, trans, diag, n, &ap[1], &work[1], &c__1, (ftnlen)1, (
241 ftnlen)1, (ftnlen)1);
242 q__1.r = -1.f, q__1.i = -0.f;
243 caxpy_(n, &q__1, &b[j * b_dim1 + 1], &c__1, &work[1], &c__1);
244
245 /* Compute componentwise relative backward error from formula */
246
247 /* max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) ) */
248
249 /* where abs(Z) is the componentwise absolute value of the matrix */
250 /* or vector Z. If the i-th component of the denominator is less */
251 /* than SAFE2, then SAFE1 is added to the i-th components of the */
252 /* numerator and denominator before dividing. */
253
254 i__2 = *n;
255 for (i__ = 1; i__ <= i__2; ++i__) {
256 i__3 = i__ + j * b_dim1;
257 rwork[i__] = (r__1 = b[i__3].r, dabs(r__1)) + (r__2 = r_imag(&b[
258 i__ + j * b_dim1]), dabs(r__2));
259 /* L20: */
260 }
261
262 if (notran) {
263
264 /* Compute abs(A)*abs(X) + abs(B). */
265
266 if (upper) {
267 kc = 1;
268 if (nounit) {
269 i__2 = *n;
270 for (k = 1; k <= i__2; ++k) {
271 i__3 = k + j * x_dim1;
272 xk = (r__1 = x[i__3].r, dabs(r__1)) + (r__2 = r_imag(&
273 x[k + j * x_dim1]), dabs(r__2));
274 i__3 = k;
275 for (i__ = 1; i__ <= i__3; ++i__) {
276 i__4 = kc + i__ - 1;
277 rwork[i__] += ((r__1 = ap[i__4].r, dabs(r__1)) + (
278 r__2 = r_imag(&ap[kc + i__ - 1]), dabs(
279 r__2))) * xk;
280 /* L30: */
281 }
282 kc += k;
283 /* L40: */
284 }
285 } else {
286 i__2 = *n;
287 for (k = 1; k <= i__2; ++k) {
288 i__3 = k + j * x_dim1;
289 xk = (r__1 = x[i__3].r, dabs(r__1)) + (r__2 = r_imag(&
290 x[k + j * x_dim1]), dabs(r__2));
291 i__3 = k - 1;
292 for (i__ = 1; i__ <= i__3; ++i__) {
293 i__4 = kc + i__ - 1;
294 rwork[i__] += ((r__1 = ap[i__4].r, dabs(r__1)) + (
295 r__2 = r_imag(&ap[kc + i__ - 1]), dabs(
296 r__2))) * xk;
297 /* L50: */
298 }
299 rwork[k] += xk;
300 kc += k;
301 /* L60: */
302 }
303 }
304 } else {
305 kc = 1;
306 if (nounit) {
307 i__2 = *n;
308 for (k = 1; k <= i__2; ++k) {
309 i__3 = k + j * x_dim1;
310 xk = (r__1 = x[i__3].r, dabs(r__1)) + (r__2 = r_imag(&
311 x[k + j * x_dim1]), dabs(r__2));
312 i__3 = *n;
313 for (i__ = k; i__ <= i__3; ++i__) {
314 i__4 = kc + i__ - k;
315 rwork[i__] += ((r__1 = ap[i__4].r, dabs(r__1)) + (
316 r__2 = r_imag(&ap[kc + i__ - k]), dabs(
317 r__2))) * xk;
318 /* L70: */
319 }
320 kc = kc + *n - k + 1;
321 /* L80: */
322 }
323 } else {
324 i__2 = *n;
325 for (k = 1; k <= i__2; ++k) {
326 i__3 = k + j * x_dim1;
327 xk = (r__1 = x[i__3].r, dabs(r__1)) + (r__2 = r_imag(&
328 x[k + j * x_dim1]), dabs(r__2));
329 i__3 = *n;
330 for (i__ = k + 1; i__ <= i__3; ++i__) {
331 i__4 = kc + i__ - k;
332 rwork[i__] += ((r__1 = ap[i__4].r, dabs(r__1)) + (
333 r__2 = r_imag(&ap[kc + i__ - k]), dabs(
334 r__2))) * xk;
335 /* L90: */
336 }
337 rwork[k] += xk;
338 kc = kc + *n - k + 1;
339 /* L100: */
340 }
341 }
342 }
343 } else {
344
345 /* Compute abs(A**H)*abs(X) + abs(B). */
346
347 if (upper) {
348 kc = 1;
349 if (nounit) {
350 i__2 = *n;
351 for (k = 1; k <= i__2; ++k) {
352 s = 0.f;
353 i__3 = k;
354 for (i__ = 1; i__ <= i__3; ++i__) {
355 i__4 = kc + i__ - 1;
356 i__5 = i__ + j * x_dim1;
357 s += ((r__1 = ap[i__4].r, dabs(r__1)) + (r__2 =
358 r_imag(&ap[kc + i__ - 1]), dabs(r__2))) *
359 ((r__3 = x[i__5].r, dabs(r__3)) + (r__4 =
360 r_imag(&x[i__ + j * x_dim1]), dabs(r__4)))
361 ;
362 /* L110: */
363 }
364 rwork[k] += s;
365 kc += k;
366 /* L120: */
367 }
368 } else {
369 i__2 = *n;
370 for (k = 1; k <= i__2; ++k) {
371 i__3 = k + j * x_dim1;
372 s = (r__1 = x[i__3].r, dabs(r__1)) + (r__2 = r_imag(&
373 x[k + j * x_dim1]), dabs(r__2));
374 i__3 = k - 1;
375 for (i__ = 1; i__ <= i__3; ++i__) {
376 i__4 = kc + i__ - 1;
377 i__5 = i__ + j * x_dim1;
378 s += ((r__1 = ap[i__4].r, dabs(r__1)) + (r__2 =
379 r_imag(&ap[kc + i__ - 1]), dabs(r__2))) *
380 ((r__3 = x[i__5].r, dabs(r__3)) + (r__4 =
381 r_imag(&x[i__ + j * x_dim1]), dabs(r__4)))
382 ;
383 /* L130: */
384 }
385 rwork[k] += s;
386 kc += k;
387 /* L140: */
388 }
389 }
390 } else {
391 kc = 1;
392 if (nounit) {
393 i__2 = *n;
394 for (k = 1; k <= i__2; ++k) {
395 s = 0.f;
396 i__3 = *n;
397 for (i__ = k; i__ <= i__3; ++i__) {
398 i__4 = kc + i__ - k;
399 i__5 = i__ + j * x_dim1;
400 s += ((r__1 = ap[i__4].r, dabs(r__1)) + (r__2 =
401 r_imag(&ap[kc + i__ - k]), dabs(r__2))) *
402 ((r__3 = x[i__5].r, dabs(r__3)) + (r__4 =
403 r_imag(&x[i__ + j * x_dim1]), dabs(r__4)))
404 ;
405 /* L150: */
406 }
407 rwork[k] += s;
408 kc = kc + *n - k + 1;
409 /* L160: */
410 }
411 } else {
412 i__2 = *n;
413 for (k = 1; k <= i__2; ++k) {
414 i__3 = k + j * x_dim1;
415 s = (r__1 = x[i__3].r, dabs(r__1)) + (r__2 = r_imag(&
416 x[k + j * x_dim1]), dabs(r__2));
417 i__3 = *n;
418 for (i__ = k + 1; i__ <= i__3; ++i__) {
419 i__4 = kc + i__ - k;
420 i__5 = i__ + j * x_dim1;
421 s += ((r__1 = ap[i__4].r, dabs(r__1)) + (r__2 =
422 r_imag(&ap[kc + i__ - k]), dabs(r__2))) *
423 ((r__3 = x[i__5].r, dabs(r__3)) + (r__4 =
424 r_imag(&x[i__ + j * x_dim1]), dabs(r__4)))
425 ;
426 /* L170: */
427 }
428 rwork[k] += s;
429 kc = kc + *n - k + 1;
430 /* L180: */
431 }
432 }
433 }
434 }
435 s = 0.f;
436 i__2 = *n;
437 for (i__ = 1; i__ <= i__2; ++i__) {
438 if (rwork[i__] > safe2) {
439 /* Computing MAX */
440 i__3 = i__;
441 r__3 = s, r__4 = ((r__1 = work[i__3].r, dabs(r__1)) + (r__2 =
442 r_imag(&work[i__]), dabs(r__2))) / rwork[i__];
443 s = dmax(r__3,r__4);
444 } else {
445 /* Computing MAX */
446 i__3 = i__;
447 r__3 = s, r__4 = ((r__1 = work[i__3].r, dabs(r__1)) + (r__2 =
448 r_imag(&work[i__]), dabs(r__2)) + safe1) / (rwork[i__]
449 + safe1);
450 s = dmax(r__3,r__4);
451 }
452 /* L190: */
453 }
454 berr[j] = s;
455
456 /* Bound error from formula */
457
458 /* norm(X - XTRUE) / norm(X) .le. FERR = */
459 /* norm( abs(inv(op(A)))* */
460 /* ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X) */
461
462 /* where */
463 /* norm(Z) is the magnitude of the largest component of Z */
464 /* inv(op(A)) is the inverse of op(A) */
465 /* abs(Z) is the componentwise absolute value of the matrix or */
466 /* vector Z */
467 /* NZ is the maximum number of nonzeros in any row of A, plus 1 */
468 /* EPS is machine epsilon */
469
470 /* The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B)) */
471 /* is incremented by SAFE1 if the i-th component of */
472 /* abs(op(A))*abs(X) + abs(B) is less than SAFE2. */
473
474 /* Use CLACON to estimate the infinity-norm of the matrix */
475 /* inv(op(A)) * diag(W), */
476 /* where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) */
477
478 i__2 = *n;
479 for (i__ = 1; i__ <= i__2; ++i__) {
480 if (rwork[i__] > safe2) {
481 i__3 = i__;
482 rwork[i__] = (r__1 = work[i__3].r, dabs(r__1)) + (r__2 =
483 r_imag(&work[i__]), dabs(r__2)) + nz * eps * rwork[
484 i__];
485 } else {
486 i__3 = i__;
487 rwork[i__] = (r__1 = work[i__3].r, dabs(r__1)) + (r__2 =
488 r_imag(&work[i__]), dabs(r__2)) + nz * eps * rwork[
489 i__] + safe1;
490 }
491 /* L200: */
492 }
493
494 kase = 0;
495 L210:
496 clacon_(n, &work[*n + 1], &work[1], &ferr[j], &kase);
497 if (kase != 0) {
498 if (kase == 1) {
499
500 /* Multiply by diag(W)*inv(op(A)**H). */
501
502 ctpsv_(uplo, transt, diag, n, &ap[1], &work[1], &c__1, (
503 ftnlen)1, (ftnlen)1, (ftnlen)1);
504 i__2 = *n;
505 for (i__ = 1; i__ <= i__2; ++i__) {
506 i__3 = i__;
507 i__4 = i__;
508 i__5 = i__;
509 q__1.r = rwork[i__4] * work[i__5].r, q__1.i = rwork[i__4]
510 * work[i__5].i;
511 work[i__3].r = q__1.r, work[i__3].i = q__1.i;
512 /* L220: */
513 }
514 } else {
515
516 /* Multiply by inv(op(A))*diag(W). */
517
518 i__2 = *n;
519 for (i__ = 1; i__ <= i__2; ++i__) {
520 i__3 = i__;
521 i__4 = i__;
522 i__5 = i__;
523 q__1.r = rwork[i__4] * work[i__5].r, q__1.i = rwork[i__4]
524 * work[i__5].i;
525 work[i__3].r = q__1.r, work[i__3].i = q__1.i;
526 /* L230: */
527 }
528 ctpsv_(uplo, transn, diag, n, &ap[1], &work[1], &c__1, (
529 ftnlen)1, (ftnlen)1, (ftnlen)1);
530 }
531 goto L210;
532 }
533
534 /* Normalize error. */
535
536 lstres = 0.f;
537 i__2 = *n;
538 for (i__ = 1; i__ <= i__2; ++i__) {
539 /* Computing MAX */
540 i__3 = i__ + j * x_dim1;
541 r__3 = lstres, r__4 = (r__1 = x[i__3].r, dabs(r__1)) + (r__2 =
542 r_imag(&x[i__ + j * x_dim1]), dabs(r__2));
543 lstres = dmax(r__3,r__4);
544 /* L240: */
545 }
546 if (lstres != 0.f) {
547 ferr[j] /= lstres;
548 }
549
550 /* L250: */
551 }
552
553 return 0;
554
555 /* End of CTPRFS */
556
557 } /* ctprfs_ */
558
559