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