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