What is the best way to parellelize a c++ MonteCarlo Integration using OpenMP?

Question

I am trying to write a c++ code which calculates an Integration b/n two numbers (a,b).

I wrote the following sequential code, but what I need to know is what is the best way to parellelize this code in order to generate the random numbers in a faster way. And is this random number generator in c++ thread safe?

The integration formula I am using is:

I = sum(f(xi))*dx and dx =(b-a)/n

double fun(double x) //f(x) = x;
{
    return x;
}

double MonteCarloIntegration (double a, double b, int n)
{
    if(a > b){
        return MonteCarloIntegration(b, a, n);
    }    
    double sum = 0.0;
    double r= 0.0;    

    for (int i = 1; i <= n; i++) 
    {
        std::random_device rd; 
        std::mt19937 gen(rd());  
        std::uniform_real_distribution<double> dis(0.0, 1.0);

        r = dis(gen);
        sum = sum + fun(a+((b-a)*r));
    }
    sum = ((b-a)/n)*sum; 
    return sum;
}

int main(int argc,char * argv[]) { 

    if (argc < 2) {
        std::cerr << "use: " << argv[0] 
                                     << " Numer_of_Random_samples (n) \n";
        std::cerr << " Example:\n  " << argv[0] << " 1000000 \n\n";
        return -1;
    }

    double b = 4.0; //lower bound
    double a = 7.0; //upper bound

    int n = atoi(argv[1]);

    std::cout <<MonteCarloIntegration(a,b,n);   

    return 0;
}

If you're doing everything inside a single function with local variables, then yes, that is thread safe, although typically you would be asking yourself what exactly needs to be "safe" when used by multiple threads, and here you share absolutely nothing. You would have to show what exactly you intend to distribute/share amongst multiple threads for that question to even be asked. — Sean F, May 23 '18 at 03:56
@SeanF std::random_device rd; I think the threads will have one seed per each but, only one random device for all them. Does this have any problem for the thread safety? — habl, May 23 '18 at 04:15
@habl there is a question on that, suggesting it is thread safe but a second call likely blocks until first call completes so not much point using it that way, although it probably depends on implementation https://stackoverflow.com/questions/42157381/thread-safety-of-stdrandom-device — Sean F, May 23 '18 at 12:06

Jorge · Answer 1 · 2018-05-24T07:49:08.970

Here I have rewritten your code to use OpenMP

#include <random>
#include <iostream>

double fun(double x) //f(x) = x;
{
    return x;
}

double MonteCarloIntegration (double a, double b, int n)
{
  std::random_device rd;
  std::mt19937 gen(rd());
  std::uniform_real_distribution<double> dis(0.0, 1.0);

  if(a > b)
    {
    return MonteCarloIntegration(b, a, n);
    }
  double sum = 0.0;

#pragma omp parallel for reduction(+:sum)
  for (int i = 1; i <= n; i++)
    {

    double r = dis(gen);
    sum = sum + fun(a+((b-a)*r));
    }
  sum = ((b-a)/n)*sum;
  return sum;
}

int main(int argc,char * argv[]) {

  if (argc < 2)
    {
    std::cerr << "use: " << argv[0]
              << " Numer_of_Random_samples (n) \n";
    std::cerr << " Example:\n  " << argv[0] << " 1000000 \n\n";
    return -1;
    }

  double b = 4.0; //lower bound
  double a = 7.0; //upper bound

  int n = atoi(argv[1]);

  std::cout << MonteCarloIntegration(a,b,n) << std::endl;

  return 0;
}

Compile this way

g++ -O3 -fopenmp integrate-mc.cxx -std=c++11 -o integrate-mc

Even if the code is correct. It is a really bad idea to keep the generator inside of the loop. — Zulan, May 23 '18 at 23:34
@Zulan and Jorge what if I put **std::mt19937 gen(rd()); std::uniform_real_distribution dis(0.0, 1.0);** this two lines inside the for loop? How should I prevent the race condition on **std::random_device rd;** ? Thank you! — habl, May 26 '18 at 01:22

What is the best way to parellelize a c++ MonteCarlo Integration using OpenMP?

1 Answers1