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