chkder.h
1 #define chkder_log10e 0.43429448190325182765
2 #define chkder_factor 100.
3 
4 namespace Eigen {
5 
6 namespace internal {
7 
8 template<typename Scalar>
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