Каков наилучший способ интеграции C ++ MonteCarlo с использованием OpenMP?

Question

Я пытаюсь написать код на С ++, который вычисляет Интеграцию b / n с двумя числами (a, b).

Я написал следующий последовательный код, но мне нужно знать, каков наилучший способ parellelize этого кода для более быстрой генерации случайных чисел.И является ли этот генератор случайных чисел в потоке c ++ безопасным?

Я использую формулу интегрирования:

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

Answer 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