algorithms.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
#ifndef DUNE_GEOMETRY_UTILITY_ALGORITHMS_HH
#define DUNE_GEOMETRY_UTILITY_ALGORITHMS_HH

#include <algorithm>
#include <cmath>
#include <limits>
#include <optional>
#include <type_traits>

#include <dune/common/debugstream.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/ftraits.hh>
#include <dune/common/fvector.hh>
#include <dune/geometry/utility/defaultmatrixhelper.hh>

namespace Dune {
namespace Impl {

template <class R = double>
struct GaussNewtonOptions
{
  int maxIt = 100;

  R absTol = []{ using std::sqrt; return sqrt(std::numeric_limits<R>::epsilon()); }();

  int maxInnerIt = 10;

  R theta = 0.5;
};


enum class GaussNewtonErrorCode
{
  OK = 0,                   //< A solution is found
  JACOBIAN_NOT_INVERTIBLE,  //< The Jacobian is not invertible at the current point
  STAGNATION,               //< No reduction of the residul norm possible
  TOLERANCE_NOT_REACHED     //< The break tolerance for the resodual norm is not reached
};


template <class F, class DF, class Domain,
          class Range = std::invoke_result_t<F, Domain>,
          class R = typename Dune::FieldTraits<Domain>::real_type>
GaussNewtonErrorCode gaussNewton (const F& f, const DF& df, Range y, Domain& x0,
                                  GaussNewtonOptions<R> opts = {})
{
  Domain x = x0;
  Domain dx{};
  Range dy = f(x0) - y;
  R resNorm0 = dy.two_norm();
  R resNorm = 0;

  for (int i = 0; i < opts.maxIt; ++i)
  {
    // Get descent direction dx: (J^T*J)dx = J^T*dy
    const bool invertible = FieldMatrixHelper<R>::xTRightInvA(df(x), dy, dx);

    // break if jacobian is not invertible
    if (!invertible)
      return GaussNewtonErrorCode::JACOBIAN_NOT_INVERTIBLE;

    // line-search procedure to update x with correction dx
    R alpha = 1;
    for (int j = 0; j < opts.maxInnerIt; ++j) {
      x = x0 - alpha * dx;
      dy = f(x) - y;
      resNorm = dy.two_norm();

      if (resNorm < resNorm0)
        break;

      alpha *= opts.theta;
    }

    // cannot reduce the residual
    if (!(resNorm < resNorm0))
      return GaussNewtonErrorCode::STAGNATION;

    x0 = x;
    resNorm0 = resNorm;

    // break if tolerance is reached.
    if (resNorm < opts.absTol)
      return GaussNewtonErrorCode::OK;
  }

  // tolerance could not be reached
  if (!(resNorm < opts.absTol))
    return GaussNewtonErrorCode::TOLERANCE_NOT_REACHED;

  return GaussNewtonErrorCode::OK;
}

} // end namespace Impl
} // end namespace Dune

#endif // DUNE_GEOMETRY_UTILITY_ALGORITHMS_HH