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