线程安全的随机数生成器用于蒙特卡罗积分

3

我正在尝试编写一个可以快速计算随机数并可应用于多个线程的程序。我的当前代码如下:

/* Approximating PI using a Monte-Carlo method. */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <omp.h>
#define N 1000000000  /* As lareg as possible for increased accuracy */

double random_function(void);

int main(void)
{
   int i = 0;
    double X, Y;
   double count_inside_temp = 0.0, count_inside = 0.0;
   unsigned int th_id = omp_get_thread_num();
   #pragma omp parallel private(i, X, Y) firstprivate(count_inside_temp)
   {
      srand(th_id);
      #pragma omp for schedule(static)
      for (i = 0; i <= N; i++) {
         X = 2.0 * random_function() - 1.0;
         Y = 2.0 * random_function() - 1.0;
         if ((X * X) + (Y * Y) < 1.0) {
        count_inside_temp += 1.0;
     }
  }
  #pragma omp atomic
      count_inside += count_inside_temp;
   }
   printf("Approximation to PI is = %.10lf\n", (count_inside * 4.0)/ N);
   return 0;
}

double random_function(void)
{
   return ((double) rand() / (double) RAND_MAX);
}

这段代码可以运行,但是从资源管理器的观察来看,我知道它没有使用所有的线程。rand()函数在多线程代码中能正常工作吗?如果不能,有没有好的替代方案?非常感谢。Jack

1个回答

7

rand()函数是否线程安全?也许是,也许不是:

rand()函数不一定是可重入的。不需要可重入性的函数也不需要线程安全。

一个测试和良好的学习练习是将对rand()的调用替换为一个固定的整数,看看会发生什么。

我认为伪随机数生成器就像一个黑盒子,它以一个整数作为输入并返回一个整数作为输出。对于任何给定的输入,输出总是相同的,但在数字序列中没有模式,并且该序列在可能的输出范围内均匀分布。(这个模型并不完全准确,但它能用。)您使用这个黑盒子的方法是选择一个起始数(种子),在应用程序中使用输出值,并将其用作下一次调用随机数生成器的输入。设计API有两种常见方法:

  1. 有两个函数,一个用于设置初始种子(例如:srand(seed)),另一个用于检索序列中的下一个值(例如:rand())。PRNG的状态存储在一种全局变量中。生成新的随机数要么不是线程安全的(很难确定,但输出流将不能被再现),要么在多线程代码中速度较慢(您最终将在状态值周围进行序列化)。
  2. 一个接口,其中PRNG状态暴露给应用程序员。这里通常有三个函数:init_prng(seed),它返回一些不透明表示PRNG状态的东西,get_prng(state),它返回一个随机数并更改状态变量,并且destroy_peng(state),它只是清理已分配的内存等。具有此类型API的PRNG应该都是线程安全的,并且可以在没有锁定的情况下并行运行(因为您负责管理(现在是线程本地的)状态变量)。
我通常使用Fortran编写代码,并使用Ladd's实现的Mersenne Twister PRNG(那个链接值得阅读)。在C语言中有许多适合控制状态的PRNG。PRNG看起来不错,使用它(并在并行区域和私有状态变量内进行初始化和销毁调用)应该能够使您获得相当不错的加速效果。
最后,如果您一次请求整个随机数序列(例如编译器可以矢量化PRNG内部),通常PRNG可以更好地执行。因此,库经常具有类似于get_prng_array(state)函数,它们会像将get_prng放在循环中填充数组元素一样更快地返回一个充满随机数的数组。这将是第二种优化方法(并且需要在并行for循环内添加一个for循环。显然,您不想在此过程中耗尽每个线程的堆栈空间!

Ladd的实现链接已失效,有没有其他Fortran实现的建议? - lalmei
一个Fortran语言的Mersenne Twister多流实现:http://theo.phys.sci.hiroshima-u.ac.jp/~ishikawa/PRNG/mt_stream_en.html - M. S. B.

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接