使用非线程安全随机数生成器更正C中pi蒙特卡罗的OpenMP杂注

Question

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

第一：这个SO线程对我没有帮助。

我自己的尝试是以下#pragmaomp语句。我认为I、x和y变量应该由每个线程初始化，并且应该是私有的。z是圆中所有命中数的总和，因此它应该在for循环后的隐含barriere之后求和。

认为主要问题是随机数生成器的静态var。我做了一个关键的部分，在那里调用函数，这样每次只有一个线程可以执行它。但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: %lfn", pi);;
    return 0;
}

这个RNG实际上是一个包含的文件，但由于我不确定我是否创建了正确的头文件，我将它集成到了另一个程序文件中，所以我只有一个.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);
}

我也尝试过使静态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龙卷风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，它仍然很糟糕。它只是比代码中使用的要好一点。

Answer 2

尝试下面的代码。它为每个线程创建一个私有状态。我在rand_r函数中做了类似的操作。为什么OpenMP的计算比单线程多花费100倍的时间？

编辑：我使用Hristo的一些建议更新了我的代码我（第一次）使用了threadprivate。我还使用了一个更好的rand函数，它可以更好地估计pi，但它仍然不够好。

一件奇怪的事情是，我必须在threadprivate之后定义函数rng_int，否则我会得到一个错误"error:'state'在第一次使用后声明为threadprivate'"。我可能应该问一个关于这个的问题。

//gcc -O3 -Wall -pedantic -fopenmp main.c
#include <omp.h>
#include <stdio.h>
#define RNG_MOD 0x80000000
int state;
int rng_int(void);
double rng_doub(double range);
int main() {
    int i, numIn, n;
    double x, y, pi;
    n = 1<<30;
    numIn = 0;
    #pragma omp threadprivate(state)
    #pragma omp parallel private(x, y) reduction(+:numIn) 
    {
        state = 25234 + 17 * omp_get_thread_num();
        #pragma omp for
        for (i = 0; i <= n; i++) {
            x = (double)rng_doub(1.0);
            y = (double)rng_doub(1.0);
            if (x*x + y*y <= 1) numIn++;
        }
    }
    pi = 4.*numIn / n;
    printf("asdf pi %fn", pi);
    return 0;
}
int rng_int(void) {
   // & 0x7fffffff is equivalent to modulo with RNG_MOD = 2^31
   return (state = (state * 1103515245 + 12345) & 0x7fffffff);
}
double rng_doub(double range) {
    return ((double)rng_int()) / (((double)RNG_MOD)/range);
}

您可以在上查看结果（以及编辑和运行代码）http://coliru.stacked-crooked.com/a/23c1753a1b7d1b0d