1 #define chkder_log10e 0.43429448190325182765 2 #define chkder_factor 100. 3 4 namespace Eigen { 5 6 namespace internal { 7 8 template<typename Scalar> 9 void chkder( 10 const Matrix< Scalar, Dynamic, 1 > &x, 11 const Matrix< Scalar, Dynamic, 1 > &fvec, 12 const Matrix< Scalar, Dynamic, Dynamic > &fjac, 13 Matrix< Scalar, Dynamic, 1 > &xp, 14 const Matrix< Scalar, Dynamic, 1 > &fvecp, 15 int mode, 16 Matrix< Scalar, Dynamic, 1 > &err 17 ) 18 { 19 typedef DenseIndex Index; 20 21 const Scalar eps = sqrt(NumTraits<Scalar>::epsilon()); 22 const Scalar epsf = chkder_factor * NumTraits<Scalar>::epsilon(); 23 const Scalar epslog = chkder_log10e * log(eps); 24 Scalar temp; 25 26 const Index m = fvec.size(), n = x.size(); 27 28 if (mode != 2) { 29 /* mode = 1. */ 30 xp.resize(n); 31 for (Index j = 0; j < n; ++j) { 32 temp = eps * abs(x[j]); 33 if (temp == 0.) 34 temp = eps; 35 xp[j] = x[j] + temp; 36 } 37 } 38 else { 39 /* mode = 2. */ 40 err.setZero(m); 41 for (Index j = 0; j < n; ++j) { 42 temp = abs(x[j]); 43 if (temp == 0.) 44 temp = 1.; 45 err += temp * fjac.col(j); 46 } 47 for (Index i = 0; i < m; ++i) { 48 temp = 1.; 49 if (fvec[i] != 0. && fvecp[i] != 0. && abs(fvecp[i] - fvec[i]) >= epsf * abs(fvec[i])) 50 temp = eps * abs((fvecp[i] - fvec[i]) / eps - err[i]) / (abs(fvec[i]) + abs(fvecp[i])); 51 err[i] = 1.; 52 if (temp > NumTraits<Scalar>::epsilon() && temp < eps) 53 err[i] = (chkder_log10e * log(temp) - epslog) / epslog; 54 if (temp >= eps) 55 err[i] = 0.; 56 } 57 } 58 } 59 60 } // end namespace internal 61 62 } // end namespace Eigen 63
