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    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