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