Loading...
Searching...
No Matches
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
defaultmatrixhelper.hh
Dune
Definition affinegeometry.hh:22