1 #define chkder_log10e 0.43429448190325182765
2 #define chkder_factor 100.
3
4 namespace Eigen {
5
6 namespace internal {
7
8 template<typename Scalar>
chkder(const Matrix<Scalar,Dynamic,1> & x,const Matrix<Scalar,Dynamic,1> & fvec,const Matrix<Scalar,Dynamic,Dynamic> & fjac,Matrix<Scalar,Dynamic,1> & xp,const Matrix<Scalar,Dynamic,1> & fvecp,int mode,Matrix<Scalar,Dynamic,1> & err)9 void 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