contrib/Eigen/src/NonLinearOptimization/chkder.h

// SPDX-FileCopyrightText: The Eigen Authors
// SPDX-License-Identifier: MPL-2.0

#ifndef EIGEN_NONLINEAROPTIMIZATION_CHKDER_H
#define EIGEN_NONLINEAROPTIMIZATION_CHKDER_H

#define chkder_log10e 0.43429448190325182765
#define chkder_factor 100.

// IWYU pragma: private
#include "./InternalHeaderCheck.h"

namespace Eigen {

namespace internal {

template <typename Scalar>
void 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) {
  using std::abs;
  using std::log;
  using std::sqrt;

  typedef DenseIndex Index;

  const Scalar eps = sqrt(NumTraits<Scalar>::epsilon());
  const Scalar epsf = chkder_factor * NumTraits<Scalar>::epsilon();
  const Scalar epslog = chkder_log10e * log(eps);
  Scalar temp;

  const Index m = fvec.size(), n = x.size();

  if (mode != 2) {
    /* mode = 1. */
    xp.resize(n);
    for (Index j = 0; j < n; ++j) {
      temp = eps * abs(x[j]);
      if (temp == 0.) temp = eps;
      xp[j] = x[j] + temp;
    }
  } else {
    /* mode = 2. */
    err.setZero(m);
    for (Index j = 0; j < n; ++j) {
      temp = abs(x[j]);
      if (temp == 0.) temp = 1.;
      err += temp * fjac.col(j);
    }
    for (Index i = 0; i < m; ++i) {
      temp = 1.;
      if (fvec[i] != 0. && fvecp[i] != 0. && abs(fvecp[i] - fvec[i]) >= epsf * abs(fvec[i]))
        temp = eps * abs((fvecp[i] - fvec[i]) / eps - err[i]) / (abs(fvec[i]) + abs(fvecp[i]));
      err[i] = 1.;
      if (temp > NumTraits<Scalar>::epsilon() && temp < eps) err[i] = (chkder_log10e * log(temp) - epslog) / epslog;
      if (temp >= eps) err[i] = 0.;
    }
  }
}

}  // end namespace internal

}  // end namespace Eigen

#undef chkder_log10e
#undef chkder_factor

#endif  // EIGEN_NONLINEAROPTIMIZATION_CHKDER_H