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


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:14