Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <gsl/gsl_integration.h>
7 : #include <limits>
8 : #include <memory>
9 : #include <vector>
10 :
11 : /// \cond
12 : namespace PUP {
13 : class er;
14 : } // namespace PUP
15 : /// \endcond
16 :
17 : /// \ingroup NumericalAlgorithmsGroup
18 : /// Numerical integration algorithms
19 1 : namespace integration {
20 :
21 : namespace detail {
22 :
23 : // The GSL functions require the integrand to have this particular function
24 : // signature. In particular, any extra parameters to the functions must be
25 : // passed as a void* and re-interpreted appropriately.
26 : template <typename IntegrandType>
27 : double integrand(const double x, void* const params) {
28 : // NOLINTNEXTLINE(cppcoreguidelines-pro-type-reinterpret-cast)
29 : auto* const function = reinterpret_cast<IntegrandType*>(params);
30 : return (*function)(x);
31 : }
32 :
33 : class GslQuadAdaptiveImpl {
34 : public:
35 : explicit GslQuadAdaptiveImpl(size_t max_intervals);
36 :
37 : GslQuadAdaptiveImpl() = default;
38 : GslQuadAdaptiveImpl(const GslQuadAdaptiveImpl&) = delete;
39 : GslQuadAdaptiveImpl& operator=(const GslQuadAdaptiveImpl&) = delete;
40 : GslQuadAdaptiveImpl(GslQuadAdaptiveImpl&&) = default;
41 : GslQuadAdaptiveImpl& operator=(GslQuadAdaptiveImpl&& rhs) = default;
42 : ~GslQuadAdaptiveImpl() = default;
43 :
44 : // NOLINTNEXTLINE(google-runtime-references)
45 : void pup(PUP::er& p);
46 :
47 : protected:
48 : template <typename IntegrandType>
49 : gsl_function make_gsl_integrand(IntegrandType& integrand) const {
50 : gsl_function gsl_integrand;
51 : gsl_integrand.function = &detail::integrand<IntegrandType>;
52 : gsl_integrand.params = &integrand;
53 : return gsl_integrand;
54 : }
55 :
56 : struct GslIntegrationWorkspaceDeleter {
57 : void operator()(gsl_integration_workspace* workspace) const;
58 : };
59 :
60 : size_t max_intervals_ = 0;
61 : std::unique_ptr<gsl_integration_workspace, GslIntegrationWorkspaceDeleter>
62 : workspace_{};
63 :
64 : private:
65 : void initialize();
66 : };
67 :
68 : void check_status_code(int status_code);
69 :
70 : void disable_gsl_error_handling();
71 :
72 : } // namespace detail
73 :
74 : /// Each type specifies which algorithm from GSL should be used. It should be
75 : /// chosen according to the problem.
76 1 : enum class GslIntegralType {
77 : /// gsl_integration_qag()
78 : StandardGaussKronrod,
79 : /// gsl_integration_qags()
80 : IntegrableSingularitiesPresent,
81 : /// gsl_integration_qagp()
82 : IntegrableSingularitiesKnown,
83 : /// gsl_integration_qagi()
84 : InfiniteInterval,
85 : /// gsl_integration_qagiu()
86 : UpperBoundaryInfinite,
87 : /// gsl_integration_qagil()
88 : LowerBoundaryInfinite
89 : };
90 :
91 : /*!
92 : * \brief A wrapper around the GSL adaptive integration procedures
93 : *
94 : * All templates take upper bounds to the absolute and relative error;
95 : * `tolerance_abs` and `tolerance_rel(default = 0.0)`, respectively. To compute
96 : * to a specified absolute error, set `tolerance_rel` to zero. To compute to a
97 : * specified relative error, set `tolerance_abs` to zero. For details on the
98 : * algorithm see the GSL documentation on `gsl_integration`.
99 : *
100 : * Here is an example how to use this class. For the function:
101 : *
102 : * \snippet Test_GslQuadAdaptive.cpp integrated_function
103 : *
104 : * the integration should look like:
105 : *
106 : * \snippet Test_GslQuadAdaptive.cpp integration_example
107 : */
108 : template <GslIntegralType TheIntegralType>
109 1 : class GslQuadAdaptive;
110 :
111 : /*!
112 : * \brief Use Gauss-Kronrod rule to integrate a 1D-function
113 : *
114 : * The algorithm for "StandardGaussKronrod" uses the QAG algorithm to employ an
115 : * adaptive Gauss-Kronrod n-points integration rule. Its function takes a
116 : * `lower_boundary` and `upper_boundary` to the integration region, an upper
117 : * bound for the absolute error `tolerance_abs` and a `key`. The latter is an
118 : * index to the array [15, 21, 31, 41, 51, 61], where each element denotes how
119 : * many function evaluations take place in each subinterval. A higher-order rule
120 : * serves better for smooth functions, whereas a lower-order rule saves time for
121 : * functions with local difficulties, such as discontinuities.
122 : */
123 : template <>
124 1 : class GslQuadAdaptive<GslIntegralType::StandardGaussKronrod>
125 : : public detail::GslQuadAdaptiveImpl {
126 : public:
127 : using detail::GslQuadAdaptiveImpl::GslQuadAdaptiveImpl;
128 : template <typename IntegrandType>
129 0 : double operator()(IntegrandType&& integrand, const double lower_boundary,
130 : const double upper_boundary, const double tolerance_abs,
131 : const int key, const double tolerance_rel = 0.) const {
132 : double result = std::numeric_limits<double>::signaling_NaN();
133 : double error = std::numeric_limits<double>::signaling_NaN();
134 : auto gsl_integrand = make_gsl_integrand(integrand);
135 : detail::disable_gsl_error_handling();
136 : const int status_code =
137 : gsl_integration_qag(&gsl_integrand, lower_boundary, upper_boundary,
138 : tolerance_abs, tolerance_rel, this->max_intervals_,
139 : key, this->workspace_.get(), &result, &error);
140 : detail::check_status_code(status_code);
141 : return result;
142 : }
143 : };
144 :
145 : /*!
146 : * \brief Integrates a 1D-function with singularities
147 : *
148 : * The algorithm for "IntegrableSingularitiesPresent" concentrates new,
149 : * increasingly smaller subintervals around an unknown singularity and makes
150 : * successive approximations to the integral which should converge towards a
151 : * limit. The integration region is defined by `lower_boundary` and
152 : * `upper_boundary`.
153 : */
154 : template <>
155 1 : class GslQuadAdaptive<GslIntegralType::IntegrableSingularitiesPresent>
156 : : public detail::GslQuadAdaptiveImpl {
157 : public:
158 : using detail::GslQuadAdaptiveImpl::GslQuadAdaptiveImpl;
159 : template <typename IntegrandType>
160 0 : double operator()(IntegrandType&& integrand, const double lower_boundary,
161 : const double upper_boundary, const double tolerance_abs,
162 : const double tolerance_rel = 0.) const {
163 : double result = std::numeric_limits<double>::signaling_NaN();
164 : double error = std::numeric_limits<double>::signaling_NaN();
165 : auto gsl_integrand = make_gsl_integrand(integrand);
166 : detail::disable_gsl_error_handling();
167 : const int status_code =
168 : gsl_integration_qags(&gsl_integrand, lower_boundary, upper_boundary,
169 : tolerance_abs, tolerance_rel, this->max_intervals_,
170 : this->workspace_.get(), &result, &error);
171 : detail::check_status_code(status_code);
172 : return result;
173 : }
174 : };
175 :
176 : /*!
177 : * \brief Integrates a 1D-function where singularities are known
178 : *
179 : * The algorithm for "IntegrableSingularitiesKnown" uses user-defined
180 : * subintervals given by a vector of doubles `points`, where each element
181 : * denotes an interval boundary.
182 : */
183 : template <>
184 1 : class GslQuadAdaptive<GslIntegralType::IntegrableSingularitiesKnown>
185 : : public detail::GslQuadAdaptiveImpl {
186 : public:
187 : using detail::GslQuadAdaptiveImpl::GslQuadAdaptiveImpl;
188 : template <typename IntegrandType>
189 0 : double operator()(IntegrandType&& integrand,
190 : const std::vector<double>& points,
191 : const double tolerance_abs,
192 : const double tolerance_rel = 0.) const {
193 : double result = std::numeric_limits<double>::signaling_NaN();
194 : double error = std::numeric_limits<double>::signaling_NaN();
195 : auto gsl_integrand = make_gsl_integrand(integrand);
196 : detail::disable_gsl_error_handling();
197 : // The const_cast is necessary because GSL has a weird interface where
198 : // the number of singularities does not change, but it doesn't take
199 : // the argument as a `const double*`. However, the first thing
200 : // `gsl_integration_qagp` does internally is forward to another function
201 : // that does take `points` by `const double*`. If `gsl_integration_qagp`
202 : // were to change the size of `points` this code would be severely broken
203 : // because `std::vector` allocates with `new`, while GSL would likely use
204 : // `malloc` (or some other C allocator). Mixing (de)allocators in such a way
205 : // is undefined behavior.
206 : const int status_code = gsl_integration_qagp(
207 : &gsl_integrand, const_cast<double*>(points.data()), // NOLINT
208 : points.size(), tolerance_abs, tolerance_rel, this->max_intervals_,
209 : this->workspace_.get(), &result, &error);
210 : detail::check_status_code(status_code);
211 : return result;
212 : }
213 : };
214 :
215 : /*!
216 : * \brief Integrates a 1D-function over the interval \f$ (-\infty, +\infty) \f$
217 : *
218 : * The algorithm for "InfiniteInterval" uses the semi-open interval
219 : * \f$ (0, 1] \f$ to map to an infinite interval \f$ (-\infty, +\infty) \f$. Its
220 : * function takes no parameters other than limits `tolerance_abs` and optional
221 : * `tolerance_rel` for the absolute error and the relative error.
222 : */
223 : template <>
224 1 : class GslQuadAdaptive<GslIntegralType::InfiniteInterval>
225 : : public detail::GslQuadAdaptiveImpl {
226 : public:
227 : using detail::GslQuadAdaptiveImpl::GslQuadAdaptiveImpl;
228 : template <typename IntegrandType>
229 0 : double operator()(IntegrandType&& integrand, const double tolerance_abs,
230 : const double tolerance_rel = 0.) const {
231 : double result = std::numeric_limits<double>::signaling_NaN();
232 : double error = std::numeric_limits<double>::signaling_NaN();
233 : auto gsl_integrand = make_gsl_integrand(integrand);
234 : detail::disable_gsl_error_handling();
235 : const int status_code = gsl_integration_qagi(
236 : &gsl_integrand, tolerance_abs, tolerance_rel, this->max_intervals_,
237 : this->workspace_.get(), &result, &error);
238 : detail::check_status_code(status_code);
239 : return result;
240 : }
241 : };
242 :
243 : /*!
244 : * \brief Integrates a 1D-function over the interval \f$ [a, +\infty) \f$
245 : *
246 : * The algorithm for "UpperBoundaryInfinite" maps the semi-open interval
247 : * \f$ (0, 1] \f$ to a semi-infinite interval \f$ [a, +\infty) \f$, where
248 : * \f$ a \f$ is given by `lower_boundary`.
249 : */
250 : template <>
251 1 : class GslQuadAdaptive<GslIntegralType::UpperBoundaryInfinite>
252 : : public detail::GslQuadAdaptiveImpl {
253 : public:
254 : using detail::GslQuadAdaptiveImpl::GslQuadAdaptiveImpl;
255 : template <typename IntegrandType>
256 0 : double operator()(IntegrandType&& integrand, const double lower_boundary,
257 : const double tolerance_abs,
258 : const double tolerance_rel = 0.) const {
259 : double result = std::numeric_limits<double>::signaling_NaN();
260 : double error = std::numeric_limits<double>::signaling_NaN();
261 : auto gsl_integrand = make_gsl_integrand(integrand);
262 : detail::disable_gsl_error_handling();
263 : const int status_code = gsl_integration_qagiu(
264 : &gsl_integrand, lower_boundary, tolerance_abs, tolerance_rel,
265 : this->max_intervals_, this->workspace_.get(), &result, &error);
266 : detail::check_status_code(status_code);
267 : return result;
268 : }
269 : };
270 :
271 : /*!
272 : * \brief Integrates a 1D-function over the interval \f$ (-\infty, b] \f$
273 : *
274 : * The algorithm for "LowerBoundaryInfinite" maps the semi-open interval
275 : * \f$ (0, 1] \f$ to a semi-infinite interval \f$ (-\infty, b] \f$, where
276 : * \f$ b \f$ is given by `upper_boundary`.
277 : */
278 : template <>
279 1 : class GslQuadAdaptive<GslIntegralType::LowerBoundaryInfinite>
280 : : public detail::GslQuadAdaptiveImpl {
281 : public:
282 : using detail::GslQuadAdaptiveImpl::GslQuadAdaptiveImpl;
283 : template <typename IntegrandType>
284 0 : double operator()(IntegrandType&& integrand, const double upper_boundary,
285 : const double tolerance_abs,
286 : const double tolerance_rel = 0.) const {
287 : double result = std::numeric_limits<double>::signaling_NaN();
288 : double error = std::numeric_limits<double>::signaling_NaN();
289 : auto gsl_integrand = make_gsl_integrand(integrand);
290 : detail::disable_gsl_error_handling();
291 : const int status_code = gsl_integration_qagil(
292 : &gsl_integrand, upper_boundary, tolerance_abs, tolerance_rel,
293 : this->max_intervals_, this->workspace_.get(), &result, &error);
294 : detail::check_status_code(status_code);
295 : return result;
296 : }
297 : };
298 : } // namespace integration
|
