Dune Core Modules (unstable)

algorithms.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
4// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
5#ifndef DUNE_GEOMETRY_UTILITY_ALGORITHMS_HH
6#define DUNE_GEOMETRY_UTILITY_ALGORITHMS_HH
7
8#include <algorithm>
9#include <cmath>
10#include <limits>
11#include <optional>
12#include <type_traits>
13
18#include <dune/geometry/affinegeometry.hh> // for FieldMatrixHelper
19
20namespace Dune {
21namespace Impl {
22
23template <class R = double>
24struct GaussNewtonOptions
25{
27 int maxIt = 100;
28
30 R absTol = []{ using std::sqrt; return sqrt(std::numeric_limits<R>::epsilon()); }();
31
33 int maxInnerIt = 10;
34
36 R theta = 0.5;
37};
38
39
41enum class GaussNewtonErrorCode
42{
43 OK = 0, //< A solution is found
44 JACOBIAN_NOT_INVERTIBLE, //< The Jacobian is not invertible at the current point
45 STAGNATION, //< No reduction of the residul norm possible
46 TOLERANCE_NOT_REACHED //< The break tolerance for the resodual norm is not reached
47};
48
49
62template <class F, class DF, class Domain,
63 class Range = std::invoke_result_t<F, Domain>,
64 class R = typename Dune::FieldTraits<Domain>::real_type>
65GaussNewtonErrorCode gaussNewton (const F& f, const DF& df, Range y, Domain& x0,
66 GaussNewtonOptions<R> opts = {})
67{
68 Domain x = x0;
69 Domain dx{};
70 Range dy = f(x0) - y;
71 R resNorm0 = dy.two_norm();
72 R resNorm = 0;
73
74 for (int i = 0; i < opts.maxIt; ++i)
75 {
76 // Get descent direction dx: (J^T*J)dx = J^T*dy
77 const bool invertible = FieldMatrixHelper<R>::xTRightInvA(df(x), dy, dx);
78
79 // break if jacobian is not invertible
80 if (!invertible)
81 return GaussNewtonErrorCode::JACOBIAN_NOT_INVERTIBLE;
82
83 // line-search procedure to update x with correction dx
84 R alpha = 1;
85 for (int j = 0; j < opts.maxInnerIt; ++j) {
86 x = x0 - alpha * dx;
87 dy = f(x) - y;
88 resNorm = dy.two_norm();
89
90 if (resNorm < resNorm0)
91 break;
92
93 alpha *= opts.theta;
94 }
95
96 // cannot reduce the residual
97 if (!(resNorm < resNorm0))
98 return GaussNewtonErrorCode::STAGNATION;
99
100 x0 = x;
101 resNorm0 = resNorm;
102
103 // break if tolerance is reached.
104 if (resNorm < opts.absTol)
105 return GaussNewtonErrorCode::OK;
106 }
107
108 // tolerance could not be reached
109 if (!(resNorm < opts.absTol))
110 return GaussNewtonErrorCode::TOLERANCE_NOT_REACHED;
111
112 return GaussNewtonErrorCode::OK;
113}
114
115} // end namespace Impl
116} // end namespace Dune
117
118#endif // DUNE_GEOMETRY_UTILITY_ALGORITHMS_HH
An implementation of the Geometry interface for affine geometries.
Defines several output streams for messages of different importance.
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Type traits to determine the type of reals (when working with complex numbers)
Implements a vector constructed from a given type representing a field and a compile-time given size.
Dune namespace.
Definition: alignedallocator.hh:13
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Dec 22, 23:30, 2024)