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