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_ */