如何使用 openMP 使这段代码线程安全？ Monte Carlo二维积分

Question

我正在尝试使用 openmp 学习并行化。
此代码执行二维 Monte Carlo 积分。对于 N = 9 亿

，程序执行（使用 omp）时间约为 25 秒

#include <cmath>
#include <cstdint>
#include <iostream>
#include <omp.h>
#include <random>

double func(double x, double y) {
    return sqrt(x+y);
}

using std::cout;
using std::endl;


const double a = 0.0, b = 1.0, c = 3.0, d = 5.0;
const uint64_t N = 900000000; 

int main() {
    double sum = 0.0;
    std::random_device rd; 

#pragma omp parallel
    {
        uint32_t seed;

#pragma omp critical
        {
            seed = rd();
        }

        std::mt19937 gen(seed);
        std::uniform_real_distribution<double> dist(0.0, 1.0);
#pragma omp for reduction(+ : sum)
        for (int64_t i = 0; i < N; i++) { 
         
            double x = a + (b - a) * dist(gen);
            double y = c + (d - c) * dist(gen);
            
            
            sum += func(x, y);
            // sum - local for each thread
            // initialized to 0 in each thread
            // and will be added to the global sum when exiting the loop
            
        }
    }

    // The result is calculated after the parallel section

    double ans = sum * (b - a) * (d - c)/ N;
    cout << ans<<endl;
    return 0;
}

如何转换成线程安全的代码（教授说不是线程安全的）？

Answer 1

据我所知，您的教授是错误的。代码是线程安全的。减少条款写对了。

我遇到的唯一问题是您初始化随机数生成器的方式。

每个线程调用一次 random_device 效率低下
这增加了生日问题的可能性。两个线程有一个（非常小的）机会以相同的随机状态开始，因为它们接收到相同的随机值。

创建多个随机序列的更有效方法是将每个线程的种子递增 1。一个好的随机数生成器会将其转换为截然不同的随机值。参见此处的示例： https://www.johndcook.com/blog/2016/01/29/random-number-generator-seed-mistakes/

另一个小问题是您可以将缩减移动到并行部分的末尾。然后您可以将 for 循环声明为 nowait 以避免在循环结束时出现一个隐式 openmp 障碍。这种优化是可能的，因为您不需要并行部分内的最终总和。但是，此更改只是一种优化。它不会改变您代码的正确性。

代码看起来像这样：

#include <cmath>
// using std::sqrt
#include <random>
// using std::random_device, std::mt19937, std::uniform_real_distribution
#include <iostream>
// using std::cout
#include <omp.h>
// using omp_get_thread_num


int main() {
    const double a = 0.0, b = 1.0, c = 3.0, d = 5.0;
    const int64_t N = 900000000;
    double sum = 0.0;
    std::random_device rd;
    // initial seed
    const std::random_device::result_type seed = rd();
#   pragma omp parallel reduction(+ : sum)
    {
        std::uniform_real_distribution<double> dist(0.0, 1.0);
        /*
        * We add the thread number to the seed. This is better than using multiple
        * random numbers because it avoids the birthday problem. See for example
        * https://www.johndcook.com/blog/2016/01/29/random-number-generator-seed-mistakes/
        */
        std::mt19937 gen(seed + static_cast<unsigned>(omp_get_thread_num()));
#       pragma omp for nowait
        for(int64_t i = 0; i < N; ++i) { 
            const double x = a + (b - a) * dist(gen);
            const double y = c + (d - c) * dist(gen);
            sum += std::sqrt(x + y);
        }
    }
    const double ans = sum * (b - a) * (d - c)/ N;
    std::cout << ans << '\n';
    return 0;
}

另请注意如何缩进编译指示以提高可读性。只需将 # 保留在行首即可。

与类似问题的关系

正如 Victor 提出的那样，这里有一个类似的问题：

问题是这段代码是否可以线程安全地使用：

int f() {
    std::random_device seeder;
    std::mt19937 engine(seeder());
    std::uniform_int_distribution<int> dist(1, 6);
    return dist(engine);
}

答案是肯定的，正如那里的答案所指出的：

No C++ std type uses global data in a non-thread-safe way.

但是，创建一个新的 mt19937 只是为了计算一个整数是没有意义的。所以那里的答案提出了多次重用 mt19937 实例的方法。但是重用需要线程安全因为

By default, one instance of a type cannot be accessed from two threads without synchronization.

最简单的答案是使用线程本地实例。我们在这里使用类似的方法，在并行部分中每个线程使用一个实例。

我的解决方案与那里提出的解决方案之间存在一些主要差异：

另一个解决方案仍然存在生日问题。但是，除非您创建大量线程，否则这个问题是小问题，所以这是可以原谅的（恕我直言）。此外，您可以通过使用不同的初始化策略来避免这种情况。例如，这会起作用：

static std::atomic<std::default_random_engine::result_type> seed(
      std::default_random_engine{}());
thread_local std::mt19937 rnd(
      seed.fetch_add(1, std::memory_order_relaxed));

那里的 mt19937 可能会活得更久，如果你重用线程，它会被更多地重用。这可能是一个好处。但是在这里的问题中，线程没有被重用。此外，我们在丢弃 RNG 之前使用了很长时间，因此使用时间更长的好处不如其他问题重要。
可能存在与线程本地数据相关的隐藏成本。 mt19937 具有与之关联的大量内存。如果不需要，我建议不要过度使用线程本地存储。另外，与任何隐藏的全局状态一样，它有点代码味道。对初始化的直接控制可能是优先的，例如用于测试

正在初始化梅森扭曲函数

我还应该提到，仅从单个整数种子值初始化具有巨大内部状态的梅森扭曲器可能不太理想。 C++'s solution of seed_seq may also be flawed 但这是我们拥有的最好的。所以这是我编写的代码大纲：

    std::array<std::mt19937::result_type, std::mt19937::state_size> seed;
    std::random_device rdev;
    std::generate(seed.begin(), seed.end(), std::ref(rdev));
#   pragma omp parallel firstprivate(seed)
    {
      seed.front() += static_cast<unsigned>(omp_get_thread_num());
      std::seed_seq sseq(seed.begin(), seed.end());
      std::mt19937 rnd(sseq);
    }

使用 MT19937 Monte Carlo

我还需要在 Wikipedia

上提到这个小笔记

Multiple instances that differ only in seed value (but not other parameters) are not generally appropriate for Monte-Carlo simulations that require independent random number generators, though there exists a method for choosing multiple sets of parameter values.[44][45]

[44] Makoto Matsumoto; Takuji Nishimura. "Dynamic Creation of Pseudorandom Number Generators" (PDF). Retrieved 19 July 2015.

[45] Hiroshi Haramoto; Makoto Matsumoto; Takuji Nishimura; François Panneton; Pierre L'Ecuyer. "Efficient Jump Ahead for F2-Linear Random Number Generators" (PDF). Retrieved 12 Nov 2015

如何使用 openMP 使这段代码线程安全？ Monte Carlo二维积分

How to make this code thread safe with openMP? Monte Carlo two-dimensional integration

c++

multithreading

parallel-processing

openmp

montecarlo

与类似问题的关系

正在初始化梅森扭曲函数

使用 MT19937 Monte Carlo