SpECTRE Documentation Coverage Report
Current view: top level - NumericalAlgorithms/Integration - GslQuadAdaptive.hpp Hit Total Coverage
Commit: 2db722c93a8e9b106e406b439b79c8e05c2057fb Lines: 8 15 53.3 %
Date: 2021-03-03 22:01:00
Legend: Lines: hit not hit

          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             : 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) noexcept;
      36             : 
      37             :   GslQuadAdaptiveImpl() = default;
      38             :   GslQuadAdaptiveImpl(const GslQuadAdaptiveImpl&) = delete;
      39             :   GslQuadAdaptiveImpl& operator=(const GslQuadAdaptiveImpl&) = delete;
      40             :   GslQuadAdaptiveImpl(GslQuadAdaptiveImpl&&) noexcept = default;
      41             :   GslQuadAdaptiveImpl& operator=(GslQuadAdaptiveImpl&& rhs) noexcept = default;
      42             :   ~GslQuadAdaptiveImpl() noexcept = default;
      43             : 
      44             :   void pup(PUP::er& p) noexcept;  // NOLINT(google-runtime-references)
      45             : 
      46             :  protected:
      47             :   template <typename IntegrandType>
      48             :   gsl_function* gsl_integrand(IntegrandType&& integrand) const noexcept {
      49             :     gsl_integrand_.function = &detail::integrand<IntegrandType>;
      50             :     gsl_integrand_.params = &integrand;
      51             :     return &gsl_integrand_;
      52             :   }
      53             : 
      54             :   struct GslIntegrationWorkspaceDeleter {
      55             :     void operator()(gsl_integration_workspace* workspace) const noexcept;
      56             :   };
      57             : 
      58             :   size_t max_intervals_ = 0;
      59             :   std::unique_ptr<gsl_integration_workspace, GslIntegrationWorkspaceDeleter>
      60             :       workspace_{};
      61             : 
      62             :  private:
      63             :   void initialize() noexcept;
      64             : 
      65             :   mutable gsl_function gsl_integrand_{};
      66             : };
      67             : 
      68             : void check_status_code(int status_code);
      69             : 
      70             : void disable_gsl_error_handling() noexcept;
      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             : class GslQuadAdaptive<GslIntegralType::StandardGaussKronrod>
     125           1 :     : 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             :     detail::disable_gsl_error_handling();
     135             :     const int status_code = gsl_integration_qag(
     136             :         gsl_integrand(std::forward<IntegrandType>(integrand)), lower_boundary,
     137             :         upper_boundary, tolerance_abs, tolerance_rel, this->max_intervals_, key,
     138             :         this->workspace_.get(), &result, &error);
     139             :     detail::check_status_code(status_code);
     140             :     return result;
     141             :   }
     142             : };
     143             : 
     144             : /*!
     145             :  * \brief Integrates a 1D-function with singularities
     146             :  *
     147             :  * The algorithm for "IntegrableSingularitiesPresent" concentrates new,
     148             :  * increasingly smaller subintervals around an unknown singularity and makes
     149             :  * successive approximations to the integral which should converge towards a
     150             :  * limit. The integration region is defined by `lower_boundary` and
     151             :  * `upper_boundary`.
     152             :  */
     153             : template <>
     154             : class GslQuadAdaptive<GslIntegralType::IntegrableSingularitiesPresent>
     155           1 :     : public detail::GslQuadAdaptiveImpl {
     156             :  public:
     157             :   using detail::GslQuadAdaptiveImpl::GslQuadAdaptiveImpl;
     158             :   template <typename IntegrandType>
     159           0 :   double operator()(IntegrandType&& integrand, const double lower_boundary,
     160             :                     const double upper_boundary, const double tolerance_abs,
     161             :                     const double tolerance_rel = 0.) const {
     162             :     double result = std::numeric_limits<double>::signaling_NaN();
     163             :     double error = std::numeric_limits<double>::signaling_NaN();
     164             :     detail::disable_gsl_error_handling();
     165             :     const int status_code = gsl_integration_qags(
     166             :         gsl_integrand(std::forward<IntegrandType>(integrand)), lower_boundary,
     167             :         upper_boundary, tolerance_abs, tolerance_rel, this->max_intervals_,
     168             :         this->workspace_.get(), &result, &error);
     169             :     detail::check_status_code(status_code);
     170             :     return result;
     171             :   }
     172             : };
     173             : 
     174             : /*!
     175             :  * \brief Integrates a 1D-function where singularities are known
     176             :  *
     177             :  * The algorithm for "IntegrableSingularitiesKnown" uses user-defined
     178             :  * subintervals given by a vector of doubles `points`, where each element
     179             :  * denotes an interval boundary.
     180             :  */
     181             : template <>
     182             : class GslQuadAdaptive<GslIntegralType::IntegrableSingularitiesKnown>
     183           1 :     : public detail::GslQuadAdaptiveImpl {
     184             :  public:
     185             :   using detail::GslQuadAdaptiveImpl::GslQuadAdaptiveImpl;
     186             :   template <typename IntegrandType>
     187           0 :   double operator()(IntegrandType&& integrand,
     188             :                     const std::vector<double>& points,
     189             :                     const double tolerance_abs,
     190             :                     const double tolerance_rel = 0.) const {
     191             :     double result = std::numeric_limits<double>::signaling_NaN();
     192             :     double error = std::numeric_limits<double>::signaling_NaN();
     193             :     detail::disable_gsl_error_handling();
     194             :     // The const_cast is necessary because GSL has a weird interface where
     195             :     // the number of singularities does not change, but it doesn't take
     196             :     // the argument as a `const double*`. However, the first thing
     197             :     // `gsl_integration_qagp` does internally is forward to another function
     198             :     // that does take `points` by `const double*`. If `gsl_integration_qagp`
     199             :     // were to change the size of `points` this code would be severely broken
     200             :     // because `std::vector` allocates with `new`, while GSL would likely use
     201             :     // `malloc` (or some other C allocator). Mixing (de)allocators in such a way
     202             :     // is undefined behavior.
     203             :     const int status_code = gsl_integration_qagp(
     204             :         gsl_integrand(std::forward<IntegrandType>(integrand)),
     205             :         const_cast<double*>(points.data()),  // NOLINT
     206             :         points.size(), tolerance_abs, tolerance_rel, this->max_intervals_,
     207             :         this->workspace_.get(), &result, &error);
     208             :     detail::check_status_code(status_code);
     209             :     return result;
     210             :   }
     211             : };
     212             : 
     213             : /*!
     214             :  * \brief Integrates a 1D-function over the interval \f$ (-\infty, +\infty) \f$
     215             :  *
     216             :  * The algorithm for "InfiniteInterval" uses the semi-open interval
     217             :  * \f$ (0, 1] \f$ to map to an infinite interval \f$ (-\infty, +\infty) \f$. Its
     218             :  * function takes no parameters other than limits `tolerance_abs` and optional
     219             :  * `tolerance_rel` for the absolute error and the relative error.
     220             :  */
     221             : template <>
     222             : class GslQuadAdaptive<GslIntegralType::InfiniteInterval>
     223           1 :     : public detail::GslQuadAdaptiveImpl {
     224             :  public:
     225             :   using detail::GslQuadAdaptiveImpl::GslQuadAdaptiveImpl;
     226             :   template <typename IntegrandType>
     227           0 :   double operator()(IntegrandType&& integrand, const double tolerance_abs,
     228             :                     const double tolerance_rel = 0.) const {
     229             :     double result = std::numeric_limits<double>::signaling_NaN();
     230             :     double error = std::numeric_limits<double>::signaling_NaN();
     231             :     detail::disable_gsl_error_handling();
     232             :     const int status_code = gsl_integration_qagi(
     233             :         gsl_integrand(std::forward<IntegrandType>(integrand)), tolerance_abs,
     234             :         tolerance_rel, this->max_intervals_, this->workspace_.get(), &result,
     235             :         &error);
     236             :     detail::check_status_code(status_code);
     237             :     return result;
     238             :   }
     239             : };
     240             : 
     241             : /*!
     242             :  * \brief Integrates a 1D-function over the interval \f$ [a, +\infty) \f$
     243             :  *
     244             :  * The algorithm for "UpperBoundaryInfinite" maps the semi-open interval
     245             :  * \f$ (0, 1] \f$ to a semi-infinite interval \f$ [a, +\infty) \f$, where
     246             :  * \f$ a \f$ is given by `lower_boundary`.
     247             :  */
     248             : template <>
     249             : class GslQuadAdaptive<GslIntegralType::UpperBoundaryInfinite>
     250           1 :     : public detail::GslQuadAdaptiveImpl {
     251             :  public:
     252             :   using detail::GslQuadAdaptiveImpl::GslQuadAdaptiveImpl;
     253             :   template <typename IntegrandType>
     254           0 :   double operator()(IntegrandType&& integrand, const double lower_boundary,
     255             :                     const double tolerance_abs,
     256             :                     const double tolerance_rel = 0.) const {
     257             :     double result = std::numeric_limits<double>::signaling_NaN();
     258             :     double error = std::numeric_limits<double>::signaling_NaN();
     259             :     detail::disable_gsl_error_handling();
     260             :     const int status_code = gsl_integration_qagiu(
     261             :         gsl_integrand(std::forward<IntegrandType>(integrand)), lower_boundary,
     262             :         tolerance_abs, tolerance_rel, this->max_intervals_,
     263             :         this->workspace_.get(), &result, &error);
     264             :     detail::check_status_code(status_code);
     265             :     return result;
     266             :   }
     267             : };
     268             : 
     269             : /*!
     270             :  * \brief Integrates a 1D-function over the interval \f$ (-\infty, b] \f$
     271             :  *
     272             :  * The algorithm for "LowerBoundaryInfinite" maps the semi-open interval
     273             :  * \f$ (0, 1] \f$ to a semi-infinite interval \f$ (-\infty, b] \f$, where
     274             :  * \f$ b \f$ is given by `upper_boundary`.
     275             :  */
     276             : template <>
     277             : class GslQuadAdaptive<GslIntegralType::LowerBoundaryInfinite>
     278           1 :     : public detail::GslQuadAdaptiveImpl {
     279             :  public:
     280             :   using detail::GslQuadAdaptiveImpl::GslQuadAdaptiveImpl;
     281             :   template <typename IntegrandType>
     282           0 :   double operator()(IntegrandType&& integrand, const double upper_boundary,
     283             :                     const double tolerance_abs,
     284             :                     const double tolerance_rel = 0.) const {
     285             :     double result = std::numeric_limits<double>::signaling_NaN();
     286             :     double error = std::numeric_limits<double>::signaling_NaN();
     287             :     detail::disable_gsl_error_handling();
     288             :     const int status_code = gsl_integration_qagil(
     289             :         gsl_integrand(std::forward<IntegrandType>(integrand)), upper_boundary,
     290             :         tolerance_abs, tolerance_rel, this->max_intervals_,
     291             :         this->workspace_.get(), &result, &error);
     292             :     detail::check_status_code(status_code);
     293             :     return result;
     294             :   }
     295             : };
     296             : }  // namespace integration

Generated by: LCOV version 1.14