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