src/NonLinearOptimization/chkder.h

0001 #define chkder_log10e 0.43429448190325182765
0002 #define chkder_factor 100.
0003
0004 namespace Eigen {
0005
0006 namespace internal {
0007
0008 template<typename Scalar>
0009 void chkder(
0010         const Matrix< Scalar, Dynamic, 1 >  &x,
0011         const Matrix< Scalar, Dynamic, 1 >  &fvec,
0012         const Matrix< Scalar, Dynamic, Dynamic > &fjac,
0013         Matrix< Scalar, Dynamic, 1 >  &xp,
0014         const Matrix< Scalar, Dynamic, 1 >  &fvecp,
0015         int mode,
0016         Matrix< Scalar, Dynamic, 1 >  &err
0017         )
0018 {
0019     using std::sqrt;
0020     using std::abs;
0021     using std::log;
0022
0023     typedef DenseIndex Index;
0024
0025     const Scalar eps = sqrt(NumTraits<Scalar>::epsilon());
0026     const Scalar epsf = chkder_factor * NumTraits<Scalar>::epsilon();
0027     const Scalar epslog = chkder_log10e * log(eps);
0028     Scalar temp;
0029
0030     const Index m = fvec.size(), n = x.size();
0031
0032     if (mode != 2) {
0033         /* mode = 1. */
0034         xp.resize(n);
0035         for (Index j = 0; j < n; ++j) {
0036             temp = eps * abs(x[j]);
0037             if (temp == 0.)
0038                 temp = eps;
0039             xp[j] = x[j] + temp;
0040         }
0041     }
0042     else {
0043         /* mode = 2. */
0044         err.setZero(m);
0045         for (Index j = 0; j < n; ++j) {
0046             temp = abs(x[j]);
0047             if (temp == 0.)
0048                 temp = 1.;
0049             err += temp * fjac.col(j);
0050         }
0051         for (Index i = 0; i < m; ++i) {
0052             temp = 1.;
0053             if (fvec[i] != 0. && fvecp[i] != 0. && abs(fvecp[i] - fvec[i]) >= epsf * abs(fvec[i]))
0054                 temp = eps * abs((fvecp[i] - fvec[i]) / eps - err[i]) / (abs(fvec[i]) + abs(fvecp[i]));
0055             err[i] = 1.;
0056             if (temp > NumTraits<Scalar>::epsilon() && temp < eps)
0057                 err[i] = (chkder_log10e * log(temp) - epslog) / epslog;
0058             if (temp >= eps)
0059                 err[i] = 0.;
0060         }
0061     }
0062 }
0063
0064 } // end namespace internal
0065
0066 } // end namespace Eigen