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