$treeview $search $mathjax
Eigen-unsupported
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
00001 #define chkder_log10e 0.43429448190325182765 00002 #define chkder_factor 100. 00003 00004 namespace Eigen { 00005 00006 namespace internal { 00007 00008 template<typename Scalar> 00009 void chkder( 00010 const Matrix< Scalar, Dynamic, 1 > &x, 00011 const Matrix< Scalar, Dynamic, 1 > &fvec, 00012 const Matrix< Scalar, Dynamic, Dynamic > &fjac, 00013 Matrix< Scalar, Dynamic, 1 > &xp, 00014 const Matrix< Scalar, Dynamic, 1 > &fvecp, 00015 int mode, 00016 Matrix< Scalar, Dynamic, 1 > &err 00017 ) 00018 { 00019 using std::sqrt; 00020 using std::abs; 00021 using std::log; 00022 00023 typedef DenseIndex Index; 00024 00025 const Scalar eps = sqrt(NumTraits<Scalar>::epsilon()); 00026 const Scalar epsf = chkder_factor * NumTraits<Scalar>::epsilon(); 00027 const Scalar epslog = chkder_log10e * log(eps); 00028 Scalar temp; 00029 00030 const Index m = fvec.size(), n = x.size(); 00031 00032 if (mode != 2) { 00033 /* mode = 1. */ 00034 xp.resize(n); 00035 for (Index j = 0; j < n; ++j) { 00036 temp = eps * abs(x[j]); 00037 if (temp == 0.) 00038 temp = eps; 00039 xp[j] = x[j] + temp; 00040 } 00041 } 00042 else { 00043 /* mode = 2. */ 00044 err.setZero(m); 00045 for (Index j = 0; j < n; ++j) { 00046 temp = abs(x[j]); 00047 if (temp == 0.) 00048 temp = 1.; 00049 err += temp * fjac.col(j); 00050 } 00051 for (Index i = 0; i < m; ++i) { 00052 temp = 1.; 00053 if (fvec[i] != 0. && fvecp[i] != 0. && abs(fvecp[i] - fvec[i]) >= epsf * abs(fvec[i])) 00054 temp = eps * abs((fvecp[i] - fvec[i]) / eps - err[i]) / (abs(fvec[i]) + abs(fvecp[i])); 00055 err[i] = 1.; 00056 if (temp > NumTraits<Scalar>::epsilon() && temp < eps) 00057 err[i] = (chkder_log10e * log(temp) - epslog) / epslog; 00058 if (temp >= eps) 00059 err[i] = 0.; 00060 } 00061 } 00062 } 00063 00064 } // end namespace internal 00065 00066 } // end namespace Eigen