thread-safety - 使用非线程安全随机数生成器更正 C 中 pi monte carlo 的 OpenMP pragma

Question

我需要一些帮助来通过给定的随机数生成器使用 openmp 的蒙特卡洛方法并行化 pi 计算，这不是线程安全的。

第一个:This所以线程没有帮助我。

我自己的尝试是以下#pragma omp 语句。我认为 i、x 和 y 变量应该由每个线程初始化，而不应该是私有(private)的。 z 是圆圈中所有命中的总和，因此它应该在 for 循环之后的隐含屏障之后求和。

认为主要问题是随机数生成器的静态变量。我做了一个调用函数的关键部分，这样每次只有一个线程可以执行它。但是 Pi 解决方案不会随着更高的值而扩展。

注意:我不应该使用另一个 RNG，但可以对其进行一些小改动。

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

    int i, z = 0, threads = 8, iters = 100000;
    double x,y, pi;

    #pragma omp parallel firstprivate(i,x,y) reduction(+:z) num_threads(threads)
        for (i=0; i<iters; ++i) {
            #pragma omp critical
            {
                x = rng_doub(1.0);
                y = rng_doub(1.0);
            }
            if ((x*x+y*y) <= 1.0)
                z++;
        }

    pi = ((double) z / (double) (iters*threads))*4.0;
    printf("Pi: %lf\n", pi);;
    return 0;
}

这个RNG实际上是一个include文件，但是我不确定我创建的头文件是否正确，所以我将它集成到另一个程序文件中，所以我只有一个.c文件。

#define RNG_MOD 741025

int rng_int(void) {
    static int state = 0;

    return (state = (1366 * state + 150889) % RNG_MOD);
}

double rng_doub(double range) {
    return ((double) rng_int()) / (double) ((RNG_MOD - 1)/range);
}

我也试过将static int状态设为全局，但它并没有改变我的结果，也许我做错了。所以请你能帮我做出正确的改变吗？非常感谢!

Answer 1

您的原始线性全等 PRNG 的循环长度为 49400，因此您只能获得 29700 个唯一测试点。这是用于任何类型的蒙特卡罗模拟的糟糕生成器。即使你做了 100000000 次试验，你也不会更接近 Pi 的真实值，因为你只是一遍又一遍地重复相同的点，结果 z 的最终值和iters 简单地乘以相同的常量，最后在除法期间取消。

Z 玻色子引入的每线程种子稍微改善了这种情况，唯一点的数量随着 OpenMP 线程总数的增加而增加。增加不是线性的，因为如果一个 PRNG 的种子落在另一个 PRNG 的序列中，则两个 PRNG 产生相同的序列，移位不超过 49400 个元素。给定周期长度，每个 PRNG 覆盖总输出范围的 49400/RNG_MOD = 6.7%，这就是两个 PRNG 同步的概率。总共有 RNG_MOD/49400 = 15 个可能的唯一序列。这基本上意味着在最好的播种情况下，您将无法超过 30 个线程，因为任何其他线程只会重复其他一些线程的结果。乘数 2 是因为每个点使用序列中的两个元素，因此如果将序列移动一个元素，可能会得到一组不同的点。

最终的解决方案是完全放弃您的 PRNG 并坚持使用类似 Mersenne twister MT19937 的东西，其循环长度为 2¹⁹⁹³⁷ − 1 和非常强大的播种算法。如果您不能像您在问题中所述那样使用另一个 PRNG，请至少修改 LCG 的常量以匹配 rand() 中使用的常量:

int rng_int(void) {
   static int state = 1;

   // & 0x7fffffff is equivalent to modulo with RNG_MOD = 2^31
   return (state = (state * 1103515245 + 12345) & 0x7fffffff);
}

请注意 rand() 不是一个好的 PRNG - 它仍然很糟糕。它只比您的代码中使用的好一点。