使用非线程安全随机数生成器在 C 中为 pi monte carlo 更正 OpenMP 编译指示

Posted 2023-02-16

【中文标题】使用非线程安全随机数生成器在 C 中为 pi monte carlo 更正 OpenMP 编译指示

【英文标题】：Correct OpenMP pragmas for pi monte carlo in C with not thread-safe random number generator

【发布时间】：2013-12-25 11:42:30

【问题描述】：

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

首先：This SO 线程对我没有帮助。

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

认为主要问题是随机数生成器的静态 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: %lf\n", 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 状态设为全局，但它不会改变我的结果，也许我做错了。所以请你能帮我做出正确的改变吗？非常感谢！

【问题讨论】：

能否请您澄清一下这些词的含义：“但是 Pi 解决方案无法随着更高的值扩展。” 除了关键部分将您的序列化这一事实之外线程，你不会得到任何加速，代码对我来说看起来是正确的。

是的，当然。我的意思是，如果我运行更多的迭代，计算出的 pi 应该更接近真实的 pi 值。但是使用这个随机数生成器，我一般看不到这种行为。讲解员说这是因为状态变量的线程不安全。我应该将其设置为全局并使用一个或多个正确的#pragma omp 语句来处理它。但是我尝试了很多组合，但没有任何改变。不知道，如果我将状态变量设为全局变量，而不是静态变量？在这个地方很关键吗？需要共享状态（）？还是更好的线程私有（状态）？已经尝试了很多。

我用你的随机函数更新了我的答案。

您的讲解员可能提到的 OpenMP 结构是 threadprivate。请参阅下面的半答案，了解为什么它不会大大改善解决方案。

@DannyArcher，我使用 Hristo 的一些建议更新了我的答案。

【参考方案1】：

试试下面的代码。它为每个线程创建一个私有状态。我对 rand_r 函数 Why does calculation with OpenMP take 100x more time than with a single thread? 做了类似的事情

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

一个奇怪的事情是我必须在 threadprivate 之后定义函数 rng_int 否则我得到一个错误“错误：'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 %f\n", 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查看结果（并编辑和运行代码）

【讨论】：

非常感谢。到目前为止，这有效。但是如果我增加 n，结果并没有变得更接近 pi。尝试了10000000。而且如果我整合线程数，我也应该增加，结果也不同，但不是正确的方式。认为这取决于 rng 函数本身，所以不是你的错。但现在我想我明白如何处理状态变量了，所以到目前为止我知道的就足够了。再次非常感谢您。

这可能是使用自定义随机数生成器的问题。我对伪随机数生成器知之甚少。我觉得使用rand_r 更舒服。但重点是您希望为每个线程使用不同的私有种子，而不是共享种子。

也可能存在计算状态的溢出。在计算状态时，您可能希望使用 64 位整数。

使用线程私有全局状态而不是将其作为参数传递给每个调用怎么样？

@HristoIliev，我使用 threadprivate 编辑了代码（我第一次使用它）并按照您的建议更改了 rand 函数。我遇到了一个奇怪的错误，我解决了“错误：'state'在第一次使用后声明了'threadprivate'”。也许我会问一个关于它的问题。

【参考方案2】：

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

Z boson 引入的每线程种子稍微改善了这种情况，唯一点的数量随着 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 - 它仍然很糟糕。它比您代码中使用的要好一点。

【讨论】：

这是很好的信息 Hristo (+1)。这就解释了为什么即使使用rand_r 也会产生不好的结果。我会说这是你的错，因为我刚刚复制了你openmp-program-is-slower-than-sequential-one，但是我在另一个答案中看到你为rand_r 写的“请注意，这是一个非常糟糕的蒙特卡罗模拟生成器，erand48 更好，最好的是如果可用，请使用 GNU Scientific Library 中的“Mersenne Twister”类型生成器。我应该阅读更多。

如果你对并行随机数生成感兴趣，你也应该看看在 SC11 上获得最佳论文奖的thesalmons.org/john/random123/papers/random123sc11.pdf。

感谢您的论文。我把它打印出来并阅读它。 Broadwell 处理器将在几个月后问世，并将具有 RDSEED 指令，该指令将生成真正的随机数。我想知道这会产生什么影响。