00001 /* dpot05.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 dpot05_(char *uplo, integer *n, integer *nrhs, 00021 doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal * 00022 x, integer *ldx, doublereal *xact, integer *ldxact, doublereal *ferr, 00023 doublereal *berr, doublereal *reslts) 00024 { 00025 /* System generated locals */ 00026 integer a_dim1, a_offset, b_dim1, b_offset, x_dim1, x_offset, xact_dim1, 00027 xact_offset, i__1, i__2, i__3; 00028 doublereal d__1, d__2, d__3; 00029 00030 /* Local variables */ 00031 integer i__, j, k; 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 /* DPOT05 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 n by n 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 / ( (n+1)*EPS + (*) ), where */ 00066 /* (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i ) */ 00067 00068 /* Arguments */ 00069 /* ========= */ 00070 00071 /* UPLO (input) CHARACTER*1 */ 00072 /* Specifies whether the upper or lower triangular part of the */ 00073 /* symmetric matrix A is stored. */ 00074 /* = 'U': Upper triangular */ 00075 /* = 'L': Lower triangular */ 00076 00077 /* N (input) INTEGER */ 00078 /* The number of rows of the matrices X, B, and XACT, and the */ 00079 /* order of the matrix A. N >= 0. */ 00080 00081 /* NRHS (input) INTEGER */ 00082 /* The number of columns of the matrices X, B, and XACT. */ 00083 /* NRHS >= 0. */ 00084 00085 /* A (input) DOUBLE PRECISION array, dimension (LDA,N) */ 00086 /* The symmetric matrix A. If UPLO = 'U', the leading n by n */ 00087 /* upper triangular part of A contains the upper triangular part */ 00088 /* of the matrix A, and the strictly lower triangular part of A */ 00089 /* is not referenced. If UPLO = 'L', the leading n by n lower */ 00090 /* triangular part of A contains the lower triangular part of */ 00091 /* the matrix A, and the strictly upper triangular part of A is */ 00092 /* not referenced. */ 00093 00094 /* LDA (input) INTEGER */ 00095 /* The leading dimension of the array A. LDA >= max(1,N). */ 00096 00097 /* B (input) DOUBLE PRECISION array, dimension (LDB,NRHS) */ 00098 /* The right hand side vectors for the system of linear */ 00099 /* equations. */ 00100 00101 /* LDB (input) INTEGER */ 00102 /* The leading dimension of the array B. LDB >= max(1,N). */ 00103 00104 /* X (input) DOUBLE PRECISION array, dimension (LDX,NRHS) */ 00105 /* The computed solution vectors. Each vector is stored as a */ 00106 /* column of the matrix X. */ 00107 00108 /* LDX (input) INTEGER */ 00109 /* The leading dimension of the array X. LDX >= max(1,N). */ 00110 00111 /* XACT (input) DOUBLE PRECISION array, dimension (LDX,NRHS) */ 00112 /* The exact solution vectors. Each vector is stored as a */ 00113 /* column of the matrix XACT. */ 00114 00115 /* LDXACT (input) INTEGER */ 00116 /* The leading dimension of the array XACT. LDXACT >= max(1,N). */ 00117 00118 /* FERR (input) DOUBLE PRECISION array, dimension (NRHS) */ 00119 /* The estimated forward error bounds for each solution vector */ 00120 /* X. If XTRUE is the true solution, FERR bounds the magnitude */ 00121 /* of the largest entry in (X - XTRUE) divided by the magnitude */ 00122 /* of the largest entry in X. */ 00123 00124 /* BERR (input) DOUBLE PRECISION array, dimension (NRHS) */ 00125 /* The componentwise relative backward error of each solution */ 00126 /* vector (i.e., the smallest relative change in any entry of A */ 00127 /* or B that makes X an exact solution). */ 00128 00129 /* RESLTS (output) DOUBLE PRECISION array, dimension (2) */ 00130 /* The maximum over the NRHS solution vectors of the ratios: */ 00131 /* RESLTS(1) = norm(X - XACT) / ( norm(X) * FERR ) */ 00132 /* RESLTS(2) = BERR / ( (n+1)*EPS + (*) ) */ 00133 00134 /* ===================================================================== */ 00135 00136 /* .. Parameters .. */ 00137 /* .. */ 00138 /* .. Local Scalars .. */ 00139 /* .. */ 00140 /* .. External Functions .. */ 00141 /* .. */ 00142 /* .. Intrinsic Functions .. */ 00143 /* .. */ 00144 /* .. Executable Statements .. */ 00145 00146 /* Quick exit if N = 0 or NRHS = 0. */ 00147 00148 /* Parameter adjustments */ 00149 a_dim1 = *lda; 00150 a_offset = 1 + a_dim1; 00151 a -= a_offset; 00152 b_dim1 = *ldb; 00153 b_offset = 1 + b_dim1; 00154 b -= b_offset; 00155 x_dim1 = *ldx; 00156 x_offset = 1 + x_dim1; 00157 x -= x_offset; 00158 xact_dim1 = *ldxact; 00159 xact_offset = 1 + xact_dim1; 00160 xact -= xact_offset; 00161 --ferr; 00162 --berr; 00163 --reslts; 00164 00165 /* Function Body */ 00166 if (*n <= 0 || *nrhs <= 0) { 00167 reslts[1] = 0.; 00168 reslts[2] = 0.; 00169 return 0; 00170 } 00171 00172 eps = dlamch_("Epsilon"); 00173 unfl = dlamch_("Safe minimum"); 00174 ovfl = 1. / unfl; 00175 upper = lsame_(uplo, "U"); 00176 00177 /* Test 1: Compute the maximum of */ 00178 /* norm(X - XACT) / ( norm(X) * FERR ) */ 00179 /* over all the vectors X and XACT using the infinity-norm. */ 00180 00181 errbnd = 0.; 00182 i__1 = *nrhs; 00183 for (j = 1; j <= i__1; ++j) { 00184 imax = idamax_(n, &x[j * x_dim1 + 1], &c__1); 00185 /* Computing MAX */ 00186 d__2 = (d__1 = x[imax + j * x_dim1], abs(d__1)); 00187 xnorm = max(d__2,unfl); 00188 diff = 0.; 00189 i__2 = *n; 00190 for (i__ = 1; i__ <= i__2; ++i__) { 00191 /* Computing MAX */ 00192 d__2 = diff, d__3 = (d__1 = x[i__ + j * x_dim1] - xact[i__ + j * 00193 xact_dim1], abs(d__1)); 00194 diff = max(d__2,d__3); 00195 /* L10: */ 00196 } 00197 00198 if (xnorm > 1.) { 00199 goto L20; 00200 } else if (diff <= ovfl * xnorm) { 00201 goto L20; 00202 } else { 00203 errbnd = 1. / eps; 00204 goto L30; 00205 } 00206 00207 L20: 00208 if (diff / xnorm <= ferr[j]) { 00209 /* Computing MAX */ 00210 d__1 = errbnd, d__2 = diff / xnorm / ferr[j]; 00211 errbnd = max(d__1,d__2); 00212 } else { 00213 errbnd = 1. / eps; 00214 } 00215 L30: 00216 ; 00217 } 00218 reslts[1] = errbnd; 00219 00220 /* Test 2: Compute the maximum of BERR / ( (n+1)*EPS + (*) ), where */ 00221 /* (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i ) */ 00222 00223 i__1 = *nrhs; 00224 for (k = 1; k <= i__1; ++k) { 00225 i__2 = *n; 00226 for (i__ = 1; i__ <= i__2; ++i__) { 00227 tmp = (d__1 = b[i__ + k * b_dim1], abs(d__1)); 00228 if (upper) { 00229 i__3 = i__; 00230 for (j = 1; j <= i__3; ++j) { 00231 tmp += (d__1 = a[j + i__ * a_dim1], abs(d__1)) * (d__2 = 00232 x[j + k * x_dim1], abs(d__2)); 00233 /* L40: */ 00234 } 00235 i__3 = *n; 00236 for (j = i__ + 1; j <= i__3; ++j) { 00237 tmp += (d__1 = a[i__ + j * a_dim1], abs(d__1)) * (d__2 = 00238 x[j + k * x_dim1], abs(d__2)); 00239 /* L50: */ 00240 } 00241 } else { 00242 i__3 = i__ - 1; 00243 for (j = 1; j <= i__3; ++j) { 00244 tmp += (d__1 = a[i__ + j * a_dim1], abs(d__1)) * (d__2 = 00245 x[j + k * x_dim1], abs(d__2)); 00246 /* L60: */ 00247 } 00248 i__3 = *n; 00249 for (j = i__; j <= i__3; ++j) { 00250 tmp += (d__1 = a[j + i__ * a_dim1], abs(d__1)) * (d__2 = 00251 x[j + k * x_dim1], abs(d__2)); 00252 /* L70: */ 00253 } 00254 } 00255 if (i__ == 1) { 00256 axbi = tmp; 00257 } else { 00258 axbi = min(axbi,tmp); 00259 } 00260 /* L80: */ 00261 } 00262 /* Computing MAX */ 00263 d__1 = axbi, d__2 = (*n + 1) * unfl; 00264 tmp = berr[k] / ((*n + 1) * eps + (*n + 1) * unfl / max(d__1,d__2)); 00265 if (k == 1) { 00266 reslts[2] = tmp; 00267 } else { 00268 reslts[2] = max(reslts[2],tmp); 00269 } 00270 /* L90: */ 00271 } 00272 00273 return 0; 00274 00275 /* End of DPOT05 */ 00276 00277 } /* dpot05_ */