C++ source
|
C++ template code for Gaussian numerical integration
over n-dimensional cube [-1 .. 1]
developed as C++ template metaprogram
Template parameters:
int ORD — order of the quadrature (must be >1)
int DIM — space dimension (must be >0)
typename FUN — integrand (may be a function or a functional object)
Control parameters:
CYCLE_INSTEAD_OF_STD_ALGO — if defined then cycles are used else std.algorithms
UNROLL_LOOPS_LIMIT — minimal array size to process by the cycle/std.algo (else by unrolled cycle)
NESTED_CYCLE_FOR_DIM_LESS_THAN_4 — if defined then ordinary nested cycles are used for dimensions 1-3
Other parameters:
real_t — real numbers type synonym
Source files:
Header (contains almost all the code)
cpp-file (contains the iterative calculation of the Gaussian quadrature coefficients that is necessary for orders >5)
How to use — a sample:
#include <iostream>
#include "GaussIntegr.h"
// integrand examples:
real_t fun3(real_t x, real_t y, real_t z) {return x*y*z;}
real_t fun3arr(const real_t* xyz) {return xyz[0]*xyz[1]*xyz[2];}
struct ffun3 {real_t operator()(real_t x, real_t y, real_t z)const {return x*y*z;} };
//integrands utilizing a vector data type:
struct vect
{real_t x,y,z;
vect(real_t v=0.): x(v), y(v),z(v) {}
real_t operator[](int idx)const {return idx==0 ? x : idx==1 ? y : z;}
real_t& operator[](int idx) {return idx==0 ? x : idx==1 ? y : z;}
vect& operator+=(const vect& v) {x+=v.x; y+=v.y; z+=v.z; return *this;}
vect& operator*=(real_t v) {x*=v; y*=v; z*=v; return *this;}
vect operator+(const vect& v)const { return (vect(*this))+=v;}
vect operator*(real_t v)const { return (vect(*this))*=v;}
};
real_t vfun3(vect v) {return v[0]*v[1]*v[2];}
vect vfun3v(const vect& v) {vect result(v); result.x+=v.x; result.y+=v.y; result.z+=v.z; return result;}
//----------------------------------
int main()
{
// creation of an integrating object for 3d-space and 4th order of Gauss quadrature:
static const fIntegrate<3,4> integr;
real_t result(0.);
std::cout<<integr.ByPlusAssgn_ArrArg(fun3arr,result)<<'\n'; //this method is absent if defined NESTED_CYCLE_FOR_DIM_LESS_THAN_4
/*arg*/
std::cout<<integr.ByPlusAssgn<real_t>(ffun3(),result)<<'\n';
std::cout<<integr.ByPlusAssgn<vect>(vfun3,result)<<'\n';
/*ret*/
std::cout<<integr.ByPlus_ArrArg<real_t>(fun3arr)<<'\n'; //is absent if defined NESTED_CYCLE_FOR_DIM_LESS_THAN_4
/*arg,ret*/
std::cout<<integr.ByPlus<real_t,real_t>(ffun3())<<'\n';
std::cout<<integr.ByPlus<vect,real_t>(vfun3)<<'\n';
<<
vect v; <<
integr.ByPlusAssgn<vect>(fvunv,v);
v = integr.ByPlus<vect,vect>(fvunv);
return 0;
}
Some test results (for the Gauss order=3):
|