使用 OpenMP 并行化 c++ Monte Carlo 集成的最佳方法是什么？答案

【问题标题】：What is the best way to parellelize a c++ MonteCarlo Integration using OpenMP?使用 OpenMP 并行化 c++ Monte Carlo 集成的最佳方法是什么？
【发布时间】：2018-05-23 03:34:45
【问题描述】：

我正在尝试编写一个 c++ 代码来计算积分 b/n 两个数字 (a,b)。

我编写了以下顺序代码，但我需要知道的是并行化此代码的最佳方法是什么，以便以更快的方式生成随机数。这个c++线程中的随机数生成器安全吗？

我使用的积分公式是：

I = sum(f(xi))*dx 和 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;
}

【问题讨论】：

您应该将前三个语句移出循环。
如果您使用局部变量在单个函数中执行所有操作，那么是的，这是线程安全的，尽管通常您会问自己在多个线程使用时究竟需要什么是“安全的” ，在这里你完全没有分享任何东西。您必须显示您打算在多个线程之间分发/共享的确切内容，才能提出该问题。
@SeanF std::random_device rd;我认为每个线程都有一个种子，但是所有线程只有一个随机设备。这对线程安全有什么问题吗？
@habl 有一个问题，表明它是线程安全的，但第二次调用可能会阻塞，直到第一次调用完成，所以以这种方式使用它并没有多大意义，尽管它可能取决于实现 stackoverflow.com/questions/42157381/…

标签： c++ multithreading parallel-processing openmp

【解决方案1】：

这里我重写了你的代码以使用 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;
}

这样编译

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

【讨论】：

即使代码正确。将生成器保留在循环内部是一个非常糟糕的主意。
现在您在 RNG 上有一个竞争条件。
@Zulan 和 Jorge 如果我把 std::mt19937 gen(rd()); std::uniform_real_distribution dis(0.0, 1.0); for 循环内的这两行？我应该如何防止 std::random_device rd; 上的竞争条件？谢谢！