Parallelise numerical integration with OpenMP in C++

Question

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.

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 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{
private:
    double B;
    double x0;
    double xt_pos;
    double y0;
    double yt_pos;
    double invLtot;
    double invtau;
    int n_lim;
public:
    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 performs 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;
}

Answer 1

1. Variable t has a type double in Composite_at_t function, but int in PDFfunction. This may result in loss of significant digits, and there are unnecessary type conversions. A related comment is that it is a good practice not to use integer literals(/constants) in expressions to calculate a double type variable. To give an example const double one_third=1/3; is 0, because 1 and 3 are both integers, so the result of this integer division is zero, but 1.0/3.0 is 0.3333333...

2. Use const keyword wherever possible. It helps the compiler to better optimize the code. For more details please read e.g. What kind of optimization does const offer in C/C++?.

3. To parallelize with OpenMP the easiest solution is to put a #pragma omp parallel for reduction(+:res) line in PDFfunction before the for loop. This is, however, not the best alternative, because you call PDFfunction twice in a row.

4. A better solution is to merge Composite_at_t and PDFfunction. In this case only one loop is needed. If you change the order of multiplications in line res +=... the n * M_PI * invL and exp(-1.0 * (n * n) * c) expressions are enough to calculate only once. I have also simplified some of your expressions (I hope I have not made any mistakes), so your Composite_at_t function should be something like:

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

const double c =  0.5 * M_PI * M_PI * t * invtau;
const double b2 = 4.0 * invLtot * invLtot * t * exp(-2.0 * t * B);

double resx = 0;
double resy = 0;

#pragma omp parallel for reduction(+:resx, resy)
for(int n = 1; n <= n_lim; ++n){
  const double tmp1 = exp(-1.0 * (n * n) * c);
  const double tmp2 = M_PI * n * invLtot;  
  resx += tmp1 * cos(tmp2* xt_pos) * cos(tmp2 * x0);
  resy += tmp1 * cos(tmp2* yt_pos) * cos(tmp2 * y0);
}

return  b2 * (1.0 + 2.0 * resx) * (1.0 + 2.0 * resy);
}