I would like to parallelise with OpenMP a one-dimensional integral using the integrate function implemented in the boost library. I found a rather obscure example (for my c++ skills) in the boost web site: https://www.boost.org/doc/libs/1_56_0/libs/numeric/odeint/doc/html/boost_numeric_odeint/tutorial/parallel_computation_with_openmp_and_mpi.html

I was wondering whether someone has experience with that, and may point out how I could implement this approach on my study case.

Also, do you think that the for loop within PDFfunction could create problems in the parallelisation of the integral?

PDFfunction represents the probability density function that is used within the integrand function:

double PDFfunction(double invL, int t, double invtau, double x0, double x, int n_lim) {

const double c =  M_PI * (M_PI/4) * ((2 * t) * invtau);

double res = 0;

for(int n = 1; n <= n_lim; ++n){

  res += exp(-1 * (n * n) * c) * cos((n * M_PI * x) * invL) * cos((n * M_PI * x0) * invL);


return invL + ((2 * invL) * res);

Composite_at_t is a function that makes use of the PDFfunction to compute pbA and pbB:

double Composite_at_t(double t, double B, double x0, double xt_pos, double y0, double yt_pos, double invLtot, double invtau, int n_lim) {

double pbA = PDFfunction(invLtot, t, invtau, x0, xt_pos, n_lim);
double pbB = PDFfunction(invLtot, t, invtau, y0, yt_pos, n_lim);
double b1 = 2 * (2 * t) * exp(-2 * t * B);
return pbA * pbB * b1;

Composite_at_tCLASS is a Func class which computes a first integral over variable t.

class Composite_at_tCLASS: public Func{
    double B;
    double x0;
    double xt_pos;
    double y0;
    double yt_pos;
    double invLtot;
    double invtau;
    int n_lim;
    Composite_at_tCLASS(double B_, double x0_, double xt_pos_, double y0_, double yt_pos_, double invLtot_, double invtau_, int n_lim_) : B(B_), x0(x0_), xt_pos(xt_pos_), y0(y0_), yt_pos(yt_pos_), invLtot(invLtot_), invtau(invtau_), n_lim(n_lim_) {}
    double operator()(const double& t) const{
        return Composite_at_t(t, B, x0, xt_pos, y0, yt_pos, invLtot, invtau, n_lim);

integrate_CompositeCLASS is the actual function that uses the class Composite_at_tCLASS and perform the integral over t, between 0 and time_lim.

double integrate_CompositeCLASS(double B, double x0, double xt_pos, double y0, double yt_pos, double invLtot, double invtau, int n_lim, double time_lim){

    Composite_at_tCLASS f(B, x0, xt_pos, y0, yt_pos, invLtot, invtau, n_lim);
    double err_est;
    int err_code;
    double res = integrate(f, 0, time_lim, err_est, err_code);
    return res;


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.