cpbt05.c
Go to the documentation of this file.
00001 /* cpbt05.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Table of constant values */
00017 
00018 static integer c__1 = 1;
00019 
00020 /* Subroutine */ int cpbt05_(char *uplo, integer *n, integer *kd, integer *
00021         nrhs, complex *ab, integer *ldab, complex *b, integer *ldb, complex *
00022         x, integer *ldx, complex *xact, integer *ldxact, real *ferr, real *
00023         berr, real *reslts)
00024 {
00025     /* System generated locals */
00026     integer ab_dim1, ab_offset, b_dim1, b_offset, x_dim1, x_offset, xact_dim1,
00027              xact_offset, i__1, i__2, i__3, i__4, i__5;
00028     real r__1, r__2, r__3, r__4;
00029     complex q__1, q__2;
00030 
00031     /* Builtin functions */
00032     double r_imag(complex *);
00033 
00034     /* Local variables */
00035     integer i__, j, k, nz;
00036     real eps, tmp, diff, axbi;
00037     integer imax;
00038     real unfl, ovfl;
00039     extern logical lsame_(char *, char *);
00040     logical upper;
00041     real xnorm;
00042     extern integer icamax_(integer *, complex *, integer *);
00043     extern doublereal slamch_(char *);
00044     real errbnd;
00045 
00046 
00047 /*  -- LAPACK test routine (version 3.1) -- */
00048 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00049 /*     November 2006 */
00050 
00051 /*     .. Scalar Arguments .. */
00052 /*     .. */
00053 /*     .. Array Arguments .. */
00054 /*     .. */
00055 
00056 /*  Purpose */
00057 /*  ======= */
00058 
00059 /*  CPBT05 tests the error bounds from iterative refinement for the */
00060 /*  computed solution to a system of equations A*X = B, where A is a */
00061 /*  Hermitian band matrix. */
00062 
00063 /*  RESLTS(1) = test of the error bound */
00064 /*            = norm(X - XACT) / ( norm(X) * FERR ) */
00065 
00066 /*  A large value is returned if this ratio is not less than one. */
00067 
00068 /*  RESLTS(2) = residual from the iterative refinement routine */
00069 /*            = the maximum of BERR / ( NZ*EPS + (*) ), where */
00070 /*              (*) = NZ*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i ) */
00071 /*              and NZ = max. number of nonzeros in any row of A, plus 1 */
00072 
00073 /*  Arguments */
00074 /*  ========= */
00075 
00076 /*  UPLO    (input) CHARACTER*1 */
00077 /*          Specifies whether the upper or lower triangular part of the */
00078 /*          Hermitian matrix A is stored. */
00079 /*          = 'U':  Upper triangular */
00080 /*          = 'L':  Lower triangular */
00081 
00082 /*  N       (input) INTEGER */
00083 /*          The number of rows of the matrices X, B, and XACT, and the */
00084 /*          order of the matrix A.  N >= 0. */
00085 
00086 /*  KD      (input) INTEGER */
00087 /*          The number of super-diagonals of the matrix A if UPLO = 'U', */
00088 /*          or the number of sub-diagonals if UPLO = 'L'.  KD >= 0. */
00089 
00090 /*  NRHS    (input) INTEGER */
00091 /*          The number of columns of the matrices X, B, and XACT. */
00092 /*          NRHS >= 0. */
00093 
00094 /*  AB      (input) COMPLEX array, dimension (LDAB,N) */
00095 /*          The upper or lower triangle of the Hermitian band matrix A, */
00096 /*          stored in the first KD+1 rows of the array.  The j-th column */
00097 /*          of A is stored in the j-th column of the array AB as follows: */
00098 /*          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; */
00099 /*          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd). */
00100 
00101 /*  LDAB    (input) INTEGER */
00102 /*          The leading dimension of the array AB.  LDAB >= KD+1. */
00103 
00104 /*  B       (input) COMPLEX array, dimension (LDB,NRHS) */
00105 /*          The right hand side vectors for the system of linear */
00106 /*          equations. */
00107 
00108 /*  LDB     (input) INTEGER */
00109 /*          The leading dimension of the array B.  LDB >= max(1,N). */
00110 
00111 /*  X       (input) COMPLEX array, dimension (LDX,NRHS) */
00112 /*          The computed solution vectors.  Each vector is stored as a */
00113 /*          column of the matrix X. */
00114 
00115 /*  LDX     (input) INTEGER */
00116 /*          The leading dimension of the array X.  LDX >= max(1,N). */
00117 
00118 /*  XACT    (input) COMPLEX array, dimension (LDX,NRHS) */
00119 /*          The exact solution vectors.  Each vector is stored as a */
00120 /*          column of the matrix XACT. */
00121 
00122 /*  LDXACT  (input) INTEGER */
00123 /*          The leading dimension of the array XACT.  LDXACT >= max(1,N). */
00124 
00125 /*  FERR    (input) REAL array, dimension (NRHS) */
00126 /*          The estimated forward error bounds for each solution vector */
00127 /*          X.  If XTRUE is the true solution, FERR bounds the magnitude */
00128 /*          of the largest entry in (X - XTRUE) divided by the magnitude */
00129 /*          of the largest entry in X. */
00130 
00131 /*  BERR    (input) REAL array, dimension (NRHS) */
00132 /*          The componentwise relative backward error of each solution */
00133 /*          vector (i.e., the smallest relative change in any entry of A */
00134 /*          or B that makes X an exact solution). */
00135 
00136 /*  RESLTS  (output) REAL array, dimension (2) */
00137 /*          The maximum over the NRHS solution vectors of the ratios: */
00138 /*          RESLTS(1) = norm(X - XACT) / ( norm(X) * FERR ) */
00139 /*          RESLTS(2) = BERR / ( NZ*EPS + (*) ) */
00140 
00141 /*  ===================================================================== */
00142 
00143 /*     .. Parameters .. */
00144 /*     .. */
00145 /*     .. Local Scalars .. */
00146 /*     .. */
00147 /*     .. External Functions .. */
00148 /*     .. */
00149 /*     .. Intrinsic Functions .. */
00150 /*     .. */
00151 /*     .. Statement Functions .. */
00152 /*     .. */
00153 /*     .. Statement Function definitions .. */
00154 /*     .. */
00155 /*     .. Executable Statements .. */
00156 
00157 /*     Quick exit if N = 0 or NRHS = 0. */
00158 
00159     /* Parameter adjustments */
00160     ab_dim1 = *ldab;
00161     ab_offset = 1 + ab_dim1;
00162     ab -= ab_offset;
00163     b_dim1 = *ldb;
00164     b_offset = 1 + b_dim1;
00165     b -= b_offset;
00166     x_dim1 = *ldx;
00167     x_offset = 1 + x_dim1;
00168     x -= x_offset;
00169     xact_dim1 = *ldxact;
00170     xact_offset = 1 + xact_dim1;
00171     xact -= xact_offset;
00172     --ferr;
00173     --berr;
00174     --reslts;
00175 
00176     /* Function Body */
00177     if (*n <= 0 || *nrhs <= 0) {
00178         reslts[1] = 0.f;
00179         reslts[2] = 0.f;
00180         return 0;
00181     }
00182 
00183     eps = slamch_("Epsilon");
00184     unfl = slamch_("Safe minimum");
00185     ovfl = 1.f / unfl;
00186     upper = lsame_(uplo, "U");
00187 /* Computing MAX */
00188     i__1 = *kd, i__2 = *n - 1;
00189     nz = (max(i__1,i__2) << 1) + 1;
00190 
00191 /*     Test 1:  Compute the maximum of */
00192 /*        norm(X - XACT) / ( norm(X) * FERR ) */
00193 /*     over all the vectors X and XACT using the infinity-norm. */
00194 
00195     errbnd = 0.f;
00196     i__1 = *nrhs;
00197     for (j = 1; j <= i__1; ++j) {
00198         imax = icamax_(n, &x[j * x_dim1 + 1], &c__1);
00199 /* Computing MAX */
00200         i__2 = imax + j * x_dim1;
00201         r__3 = (r__1 = x[i__2].r, dabs(r__1)) + (r__2 = r_imag(&x[imax + j * 
00202                 x_dim1]), dabs(r__2));
00203         xnorm = dmax(r__3,unfl);
00204         diff = 0.f;
00205         i__2 = *n;
00206         for (i__ = 1; i__ <= i__2; ++i__) {
00207             i__3 = i__ + j * x_dim1;
00208             i__4 = i__ + j * xact_dim1;
00209             q__2.r = x[i__3].r - xact[i__4].r, q__2.i = x[i__3].i - xact[i__4]
00210                     .i;
00211             q__1.r = q__2.r, q__1.i = q__2.i;
00212 /* Computing MAX */
00213             r__3 = diff, r__4 = (r__1 = q__1.r, dabs(r__1)) + (r__2 = r_imag(&
00214                     q__1), dabs(r__2));
00215             diff = dmax(r__3,r__4);
00216 /* L10: */
00217         }
00218 
00219         if (xnorm > 1.f) {
00220             goto L20;
00221         } else if (diff <= ovfl * xnorm) {
00222             goto L20;
00223         } else {
00224             errbnd = 1.f / eps;
00225             goto L30;
00226         }
00227 
00228 L20:
00229         if (diff / xnorm <= ferr[j]) {
00230 /* Computing MAX */
00231             r__1 = errbnd, r__2 = diff / xnorm / ferr[j];
00232             errbnd = dmax(r__1,r__2);
00233         } else {
00234             errbnd = 1.f / eps;
00235         }
00236 L30:
00237         ;
00238     }
00239     reslts[1] = errbnd;
00240 
00241 /*     Test 2:  Compute the maximum of BERR / ( NZ*EPS + (*) ), where */
00242 /*     (*) = NZ*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i ) */
00243 
00244     i__1 = *nrhs;
00245     for (k = 1; k <= i__1; ++k) {
00246         i__2 = *n;
00247         for (i__ = 1; i__ <= i__2; ++i__) {
00248             i__3 = i__ + k * b_dim1;
00249             tmp = (r__1 = b[i__3].r, dabs(r__1)) + (r__2 = r_imag(&b[i__ + k *
00250                      b_dim1]), dabs(r__2));
00251             if (upper) {
00252 /* Computing MAX */
00253                 i__3 = i__ - *kd;
00254                 i__4 = i__ - 1;
00255                 for (j = max(i__3,1); j <= i__4; ++j) {
00256                     i__3 = *kd + 1 - i__ + j + i__ * ab_dim1;
00257                     i__5 = j + k * x_dim1;
00258                     tmp += ((r__1 = ab[i__3].r, dabs(r__1)) + (r__2 = r_imag(&
00259                             ab[*kd + 1 - i__ + j + i__ * ab_dim1]), dabs(r__2)
00260                             )) * ((r__3 = x[i__5].r, dabs(r__3)) + (r__4 = 
00261                             r_imag(&x[j + k * x_dim1]), dabs(r__4)));
00262 /* L40: */
00263                 }
00264                 i__4 = *kd + 1 + i__ * ab_dim1;
00265                 i__3 = i__ + k * x_dim1;
00266                 tmp += (r__1 = ab[i__4].r, dabs(r__1)) * ((r__2 = x[i__3].r, 
00267                         dabs(r__2)) + (r__3 = r_imag(&x[i__ + k * x_dim1]), 
00268                         dabs(r__3)));
00269 /* Computing MIN */
00270                 i__3 = i__ + *kd;
00271                 i__4 = min(i__3,*n);
00272                 for (j = i__ + 1; j <= i__4; ++j) {
00273                     i__3 = *kd + 1 + i__ - j + j * ab_dim1;
00274                     i__5 = j + k * x_dim1;
00275                     tmp += ((r__1 = ab[i__3].r, dabs(r__1)) + (r__2 = r_imag(&
00276                             ab[*kd + 1 + i__ - j + j * ab_dim1]), dabs(r__2)))
00277                              * ((r__3 = x[i__5].r, dabs(r__3)) + (r__4 = 
00278                             r_imag(&x[j + k * x_dim1]), dabs(r__4)));
00279 /* L50: */
00280                 }
00281             } else {
00282 /* Computing MAX */
00283                 i__4 = i__ - *kd;
00284                 i__3 = i__ - 1;
00285                 for (j = max(i__4,1); j <= i__3; ++j) {
00286                     i__4 = i__ + 1 - j + j * ab_dim1;
00287                     i__5 = j + k * x_dim1;
00288                     tmp += ((r__1 = ab[i__4].r, dabs(r__1)) + (r__2 = r_imag(&
00289                             ab[i__ + 1 - j + j * ab_dim1]), dabs(r__2))) * ((
00290                             r__3 = x[i__5].r, dabs(r__3)) + (r__4 = r_imag(&x[
00291                             j + k * x_dim1]), dabs(r__4)));
00292 /* L60: */
00293                 }
00294                 i__3 = i__ * ab_dim1 + 1;
00295                 i__4 = i__ + k * x_dim1;
00296                 tmp += (r__1 = ab[i__3].r, dabs(r__1)) * ((r__2 = x[i__4].r, 
00297                         dabs(r__2)) + (r__3 = r_imag(&x[i__ + k * x_dim1]), 
00298                         dabs(r__3)));
00299 /* Computing MIN */
00300                 i__4 = i__ + *kd;
00301                 i__3 = min(i__4,*n);
00302                 for (j = i__ + 1; j <= i__3; ++j) {
00303                     i__4 = j + 1 - i__ + i__ * ab_dim1;
00304                     i__5 = j + k * x_dim1;
00305                     tmp += ((r__1 = ab[i__4].r, dabs(r__1)) + (r__2 = r_imag(&
00306                             ab[j + 1 - i__ + i__ * ab_dim1]), dabs(r__2))) * (
00307                             (r__3 = x[i__5].r, dabs(r__3)) + (r__4 = r_imag(&
00308                             x[j + k * x_dim1]), dabs(r__4)));
00309 /* L70: */
00310                 }
00311             }
00312             if (i__ == 1) {
00313                 axbi = tmp;
00314             } else {
00315                 axbi = dmin(axbi,tmp);
00316             }
00317 /* L80: */
00318         }
00319 /* Computing MAX */
00320         r__1 = axbi, r__2 = nz * unfl;
00321         tmp = berr[k] / (nz * eps + nz * unfl / dmax(r__1,r__2));
00322         if (k == 1) {
00323             reslts[2] = tmp;
00324         } else {
00325             reslts[2] = dmax(reslts[2],tmp);
00326         }
00327 /* L90: */
00328     }
00329 
00330     return 0;
00331 
00332 /*     End of CPBT05 */
00333 
00334 } /* cpbt05_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:32