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


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