GaussLegendreQuadratureQ.hh

Go to the documentation of this file.

#ifndef NPSTAT_GAUSSLEGENDREQUADRATUREQ_HH_
#define NPSTAT_GAUSSLEGENDREQUADRATUREQ_HH_
 
/*!
// \file GaussLegendreQuadratureQ.hh
//
// \brief Gauss-Legendre quadratures in quad precision
//
// Author: I. Volobouev
//
// July 2020
*/
 
#include <vector>
#include <utility>
 
#include "npstat/nm/lapack_double.h"
#include "npstat/nm/SimpleFunctors.hh"
 
namespace npstat {
    /** 
    // Gauss-Legendre quadrature. Internally, operations are performed
    // in lapack_double precision.
    */
    class GaussLegendreQuadratureQ
    {
    public:
        enum {interval_quadrature = 1};
 
        /**
        // At the moment, the following numbers of points are supported:
        // 2, 4, 6, 8, 10, 12, 16, 32, 64, 100, 128, 256, 512, 1024.
        //
        // If an unsupported number of points is given in the
        // constructor, std::invalid_argument exception will be thrown.
        */
        explicit GaussLegendreQuadratureQ(unsigned npoints);
 
        /** Return the number of quadrature points */
        inline unsigned npoints() const {return npoints_;}
 
        /**
        // The abscissae are returned for positive points only,
        // so the buffer length should be at least npoints/2.
        */
        void getAbscissae(lapack_double* abscissae, unsigned len) const;
 
        /** Return abscissae for all points */
        void getAllAbscissae(lapack_double* abscissae, unsigned len) const;
 
        /**
        // The weights are returned for positive points only,
        // so the buffer length should be at least npoints/2.
        */
        void getWeights(lapack_double* weights, unsigned len) const;
 
        /** Return weights for all points */
        void getAllWeights(lapack_double* weights, unsigned len) const;
 
        /** Perform the quadrature on [a, b] interval */
        template <typename FcnResult, typename FcnArg>
        lapack_double integrate(const Functor1<FcnResult,FcnArg>& fcn,
                                lapack_double a, lapack_double b) const;
 
        /**
        // This method splits the interval [a, b] into "nsplit"
        // subintervals of equal length, applies Gauss-Legendre
        // quadrature to each subinterval, and sums the results.
        */
        template <typename FcnResult, typename FcnArg>
        lapack_double integrate(const Functor1<FcnResult,FcnArg>& fcn,
                                lapack_double a, lapack_double b,
                                unsigned nsplit) const;
 
        /**
        // Weighted integration points on the given interval, suitable
        // for constructing orthogonal polynomials w.r.t. the given
        // weight function (in particular, by ContOrthoPoly1D class).
        // Naturally, rule with "npoints" points must be able to calculate
        // polynomial normalization integrals exactly.
        */
        template <class Functor>
        std::vector<std::pair<lapack_double,lapack_double> >
        weightedIntegrationPoints(
            const Functor& weight, lapack_double a, lapack_double b) const;
 
        /** Check if the rule with the given number of points is supported */
        static bool isAllowed(unsigned npoints);
 
        /** The complete set of allowed rules, in the increasing order */
        static std::vector<unsigned> allowedNPonts();
 
        /**
        // Minimum number of points, among the supported rules, which
        // integrates a polynomial with the given degree exactly.
        // Returns 0 if the degree is out of range.
        */
        static unsigned minimalExactRule(unsigned polyDegree);
 
    private:
        GaussLegendreQuadratureQ();
 
        const lapack_double* a_;
        const lapack_double* w_;
        mutable std::vector<std::pair<lapack_double, lapack_double> > buf_;
        unsigned npoints_;
    };
}
 
#include "npstat/nm/GaussLegendreQuadratureQ.icc"
 
#endif // NPSTAT_GAUSSLEGENDREQUADRATUREQ_HH_