chkder.h
00001 
00002 #define chkder_log10e 0.43429448190325182765
00003 #define chkder_factor 100.
00004 
00005 namespace Eigen { 
00006 
00007 namespace internal {
00008 
00009 template<typename Scalar>
00010 void chkder(
00011         const Matrix< Scalar, Dynamic, 1 >  &x,
00012         const Matrix< Scalar, Dynamic, 1 >  &fvec,
00013         const Matrix< Scalar, Dynamic, Dynamic > &fjac,
00014         Matrix< Scalar, Dynamic, 1 >  &xp,
00015         const Matrix< Scalar, Dynamic, 1 >  &fvecp,
00016         int mode,
00017         Matrix< Scalar, Dynamic, 1 >  &err
00018         )
00019 {
00020     typedef DenseIndex Index;
00021 
00022     const Scalar eps = sqrt(NumTraits<Scalar>::epsilon());
00023     const Scalar epsf = chkder_factor * NumTraits<Scalar>::epsilon();
00024     const Scalar epslog = chkder_log10e * log(eps);
00025     Scalar temp;
00026 
00027     const Index m = fvec.size(), n = x.size();
00028 
00029     if (mode != 2) {
00030         /* mode = 1. */
00031         xp.resize(n);
00032         for (Index j = 0; j < n; ++j) {
00033             temp = eps * abs(x[j]);
00034             if (temp == 0.)
00035                 temp = eps;
00036             xp[j] = x[j] + temp;
00037         }
00038     }
00039     else {
00040         /* mode = 2. */
00041         err.setZero(m); 
00042         for (Index j = 0; j < n; ++j) {
00043             temp = abs(x[j]);
00044             if (temp == 0.)
00045                 temp = 1.;
00046             err += temp * fjac.col(j);
00047         }
00048         for (Index i = 0; i < m; ++i) {
00049             temp = 1.;
00050             if (fvec[i] != 0. && fvecp[i] != 0. && abs(fvecp[i] - fvec[i]) >= epsf * abs(fvec[i]))
00051                 temp = eps * abs((fvecp[i] - fvec[i]) / eps - err[i]) / (abs(fvec[i]) + abs(fvecp[i]));
00052             err[i] = 1.;
00053             if (temp > NumTraits<Scalar>::epsilon() && temp < eps)
00054                 err[i] = (chkder_log10e * log(temp) - epslog) / epslog;
00055             if (temp >= eps)
00056                 err[i] = 0.;
00057         }
00058     }
00059 }
00060 
00061 } // end namespace internal
00062 
00063 } // end namespace Eigen