如何并行生成随机数？

Question

如何并行生成随机数？

17

我想使用OpenMP并行生成伪随机数，类似于这样：

int i;
#pragma omp parallel for
for (i=0;i<100;i++)
{
    printf("%d %d %d\n",i,omp_get_thread_num(),rand());
} 
return 0;

我已在Windows上进行了测试，获得了巨大的加速，但每个线程生成的数字完全相同。我还在Linux上进行了测试，发现并行版本在8核处理器上约比顺序版本慢10倍，但每个线程生成的数字不同。

有没有办法既能提升速度又能生成不同的数字？

编辑 27.11.2010
我认为通过借鉴Jonathan Dursi的帖子中的想法，我已经解决了这个问题。似乎以下代码在Linux和Windows上都能快速运行。数字也是伪随机的。您对此有何看法？

int seed[10];

int main(int argc, char **argv) 
{
int i,s;
for (i=0;i<10;i++)
    seed[i] = rand();

#pragma omp parallel private(s)
{
    s = seed[omp_get_thread_num()];
    #pragma omp for
    for (i=0;i<1000;i++)
    {
        printf("%d %d %d\n",i,omp_get_thread_num(),s);
        s=(s*17931+7391); // those numbers should be choosen more carefully
    }
    seed[omp_get_thread_num()] = s;
}
return 0; 
}

PS.：我还没有接受任何答案，因为我需要确保这个想法是好的。

- Tomek Tarczynski

rand 是一个非常低质量的伪随机数生成器，只有在需要兼容性时才应该使用它（例如为了复制使用此相同糟糕的 PRNG 的模拟运行）。大多数操作系统/库提供更好的 PRNG（例如 FreeBSD 有 random，lrand48，arc4random 等）。 - Dave C

此外，考虑使用基于计数器的伪随机数生成器，例如在论文《并行随机数：像1、2、3一样简单》中描述的那些。链接 - Dave C

6个回答

9

在多个线程中不能使用C语言的rand()函数，这会导致未定义的行为。一些实现可能会提供锁定（这会使它变慢）；其他实现可能允许线程覆盖彼此的状态，可能会导致程序崩溃或产生“错误”的随机数。

要解决这个问题，可以编写自己的PRNG实现，或者使用一个现有的PRNG实现，允许调用者存储并传递状态给PRNG迭代器函数。

- R.. GitHub STOP HELPING ICE

非常好。一般的想法是rand()（和其他函数）会锁定访问，有效地强制线程等待（使实现接近单线程甚至更慢），但从未指出可能会“破坏彼此状态，可能导致程序崩溃”的可能性，我曾亲眼目睹过这种情况！ - SChepurin

4

让每个线程根据其线程ID设置不同的种子，例如srand(omp_get_thread_num() * 1000)；

- moinudin

2

如果没有一些逻辑检查种子是否在所有线程上初始化，那么几乎可以肯定这不会消除Linux上的减速。 - Axel Gneiting

说明：http://software.intel.com/en-us/blogs/2009/11/05/use-of-rand-in-openmp-parallel-sections/ - chrisaycock

@Axel 这可能是因为rand()有一个原子操作，它会锁定。你需要寻找一个非锁定的随机数生成器。 - moinudin

我尝试使用rand_r()函数来测试可重入版本是否更快（无需锁定），但在我的系统上花费的时间与之前相同。 - chrisaycock

rand 不一定会进行锁定，也不应该这样做。从多个线程调用它会导致未定义的行为。 - R.. GitHub STOP HELPING ICE

Axel是正确的，它没有解决Linux上的减速问题。 - Tomek Tarczynski

4

似乎在Linux上，rand具有全局共享状态，并且在Windows上具有线程本地存储状态。由于必须进行同步，Linux上的共享状态会导致速度变慢。我认为，在C库中没有一种可移植的方法可以在多个线程上并行使用RNG，因此您需要另一个。您可以使用Mersenne Twister。正如marcog所说，您需要为每个线程不同地初始化种子。

- Axel Gneiting

实际上，除了使用自己的互斥锁来包装对 rand() 的调用之外，没有其他可移植的方法......但这样做会很不划算。 - R.. GitHub STOP HELPING ICE

1

rand_r()是可移植的（在POSIX 1.c中）且可重入的。 - Jonathan Dursi

我需要更仔细地研究Mersenne Twister，因为这种方法并不像大多数伪随机数生成器那样显而易见。 - Tomek Tarczynski

Jonathan，可重入的rand函数无法帮助进行并行数值生成，因为它需要同步。Tomek，你是在使用纯C还是C++？ - Axel Gneiting

我阅读了rand_r的规范，发现我错了。种子不是全局状态，而是作为函数参数给出的。所以它可以工作，但无论如何都不具备可移植性。 - Axel Gneiting

3

在Linux/Unix系统中，您可以使用以下命令：

long jrand48(unsigned short xsubi[3]);

xsubi[3]是随机数生成器的状态，编码方式如下：

#include<stdio.h>
#include<stdlib.h>
#include <algorithm> 
int main() {
  unsigned short *xsub;
#pragma omp parallel private(xsub)
  {  
    xsub = new unsigned short[3];
    xsub[0]=xsub[1]=xsub[2]= 3+omp_get_thread_num();
    int j;
#pragma omp for
    for(j=0;j<10;j++) 
      printf("%d [%d] %ld\n", j, omp_get_thread_num(), jrand48(xsub));
  }
}

编译使用

g++-mp-4.4 -Wall -Wextra -O2 -march=native -fopenmp -D_GLIBCXX_PARALLEL jrand.cc -o jrand

(将g++-mp-4.4替换为您需要调用的g++版本4.4或4.3) 然后您会得到

$ ./jrand 
0 [0] 1344229389
1 [0] 1845350537
2 [0] 229759373
3 [0] 1219688060
4 [0] -553792943
5 [1] 360650087
6 [1] -404254894
7 [1] 1678400333
8 [1] 1373359290
9 [1] 171280263

即在没有任何互斥锁或竞态条件的情况下生成10个不同的伪随机数。

- Riko Jacob

你能详细解释一下你的答案吗？我从未听说过jrand48，我认为这个函数不在任何标准库中。 - Tomek Tarczynski

jrand48属于drand48()和lrand48()的“家族”。它是“标准C库(libc，-lc)”的一部分，因此需要包含头文件<stdlib.h>。 - Riko Jacob

0

随机数可以非常快速地生成，因此通常内存会成为瓶颈。通过在多个线程之间分配此任务，您会创建额外的通信和同步开销（不同核心缓存的同步也不便宜）。

最好使用单个线程和更好的random()函数。

- ruslik

1

这对我来说可能不是一个好的解决方案，因为我的程序会生成很多随机数，并且需要并发处理。 - Tomek Tarczynski

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Jonathan Dursi · Accepted Answer

我将在此处发布我在Concurrent random number generation中发布的内容：

我认为你正在寻找rand_r()，它明确将当前RNG状态作为参数。然后每个线程应该有自己的种子数据副本（无论您是否希望每个线程以相同的种子开始或不同的种子取决于您要做什么，在这里您需要它们是不同的，否则您会一遍又一遍地得到相同的行）。这里有一些关于rand_r()和线程安全性的讨论：whether rand_r is real thread safe?。

所以说，假设您想让每个线程的种子从其线程编号开始（这可能不是您想要的，因为它会在每次使用相同数量的线程运行时给出相同的结果，但只是一个例子）：

#pragma omp parallel default(none)
{
    int i;
    unsigned int myseed = omp_get_thread_num();
    #pragma omp for
    for(i=0; i<100; i++)
            printf("%d %d %d\n",i,omp_get_thread_num(),rand_r(&myseed));
}

编辑: 只是为了好玩，检查了一下上面的代码是否可以获得任何加速。完整代码如下：

#define NRANDS 1000000
int main(int argc, char **argv) {

    struct timeval t;
    int a[NRANDS];

    tick(&t);
    #pragma omp parallel default(none) shared(a)
    {
        int i;
        unsigned int myseed = omp_get_thread_num();
        #pragma omp for
        for(i=0; i<NRANDS; i++)
                a[i] = rand_r(&myseed);
    }
    double sum = 0.;
    double time=tock(&t);
    for (long int i=0; i<NRANDS; i++) {
        sum += a[i];
    }
    printf("Time = %lf, sum = %lf\n", time, sum);

    return 0;
}

其中tick和tock只是gettimeofday()的包装器，而tock()返回秒数差。为了确保没有任何优化被忽略，并演示一个小点，会打印出sum；因为每个线程都有自己的threadnum作为种子，所以使用不同数量的线程将得到不同的数字；如果您使用相同数量的线程多次运行相同的代码，则出于同样的原因，您将获得相同的总和。无论如何，在8核nehalem框架上运行（没有其他用户），时间如下：

$ export OMP_NUM_THREADS=1
$ ./rand
Time = 0.008639, sum = 1074808568711883.000000

$ export OMP_NUM_THREADS=2
$ ./rand
Time = 0.006274, sum = 1074093295878604.000000

$ export OMP_NUM_THREADS=4
$ ./rand
Time = 0.005335, sum = 1073422298606608.000000

$ export OMP_NUM_THREADS=8
$ ./rand
Time = 0.004163, sum = 1073971133482410.000000

因此，加速效果不是很明显；正如@ruslik所指出的那样，这并不是一个真正的计算密集型过程，其他问题，如内存带宽开始发挥作用。因此，在8个核心上只有略微超过2倍的加速效果。