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


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