为什么这个随机数生成器不是线程安全的？

Question

我使用 rand() 函数生成 0,1 之间的伪随机数用于模拟目的，但是当我决定让我的 C++ 代码并行运行（通过 OpenMP）时，我注意到 rand() 不是线程- 安全，也不是很统一。

所以我转而使用（所谓的）更统一的生成器，这些生成器在其他问题的许多答案中都出现了。看起来像这样

double rnd(const double & min, const double & max) {
    static thread_local mt19937* generator = nullptr;
    if (!generator) generator = new mt19937(clock() + omp_get_thread_num());
    uniform_real_distribution<double> distribution(min, max);
    return fabs(distribution(*generator));
}

但我在模拟的原始问题中发现了许多科学错误。既违反rand()的结果又违反常识的问题。

所以我写了一个代码，用这个函数生成 500k 个随机数，计算它们的平均值，然后重复 200 次并绘制结果。

double SUM=0;
for(r=0; r<=10; r+=0.05){   
    #pragma omp parallel for ordered schedule(static)
    for(w=1; w<=500000; w++){   
        double a;
        a=rnd(0,1);
        SUM=SUM+a;
    } 
    SUM=SUM/w_max;
    ft<<r<<'\t'<<SUM<<'\n';
    SUM=0;
}   

我们知道，如果不是 500k，我可以无限次执行它，它应该是一条值为 0.5 的简单线。但是 500k 我们会有 0.5 左右的波动。

单线程运行代码时，结果可以接受：

但这是 2 个线程的结果：

3 个线程：

4 个线程：

我现在没有我的 8 线程 CPU，但结果还是值得的。

如您所见，它们都不统一，并且在平均值附近波动很大。

那么这个伪随机生成器也是线程不安全的吗？

还是我在某个地方犯了错误？

Answer 1

我将对您的测试输出提出三个意见：

• 它的方差比一个好的随机源的平均值要强得多。您通过与单线程结果进行比较自己观察到这一点。

• 计算出的平均值随着线程数的增加而减少，并且永远不会达到原来的 0.5（即不仅方差更高，而且平均值也发生了变化）。

• 数据中存在时间关系，在 4 线程情况下尤其明显。

所有这些都可以通过代码中存在的竞争条件来解释：您从多个线程分配给SUM。递增双精度不是原子操作（即使在 x86 上，您可能会在寄存器上进行原子读取和写入）。两个线程可以读取当前值（例如 10），增加它（例如都加 0.5），然后将值写回内存。现在您有两个线程写入 10.5 而不是正确的 11。

尝试并发（不同步）写入SUM 的线程越多，丢失的更改就越多。这解释了所有观察结果：

• 线程在每次单独运行中相互竞争的程度决定了丢失了多少结果，这可能因运行而异。

• 比赛越多（例如线程越多），平均值越低，因为丢失的结果越多。您永远无法超过统计平均值 0.5，因为您只会丢失写入。

• 随着线程和调度程序的“稳定”，差异会减少。这与您应该在基准测试时“预热”测试的原因类似。

不用说，这是未定义的行为。它只是在 x86 CPU 上显示良性行为，但这不是 C++ 标准向您保证的。如您所知，double 的各个字节可能会被不同的线程同时写入，从而导致完全垃圾。

正确的解决方案是在本地添加双打线程，然后（通过同步）最后将它们添加在一起。 OMP 有针对此特定目的的缩减条款。

对于整数类型，您可以使用std::atomic<IntegralType>::fetch_add()。 std::atomic<double> 存在，但（在 C++20 之前）提到的函数（和其他）仅可用于整数类型。

Answer 2

问题不在于你的 RNG，而在于你的总结。 SUM 上只有一个竞争条件。要解决此问题，请使用减少，例如将编译指示更改为：

#pragma omp parallel for ordered schedule(static) reduction(+:SUM)

请注意，将thread_local 与 OpenMP 一起使用在技术上不是定义的行为。它可能会在实践中起作用，但 OpenMP 和 C++11 线程概念之间的交互并没有得到很好的定义（另请参见 this question）。因此，对您来说安全的 OpenMP 替代方案是：

static mt19937* generator = nullptr;
#pragma omp threadprivate(generator)