如何并行生成随机数？

Question

我想使用 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.：我还没有接受任何答案，因为我需要确定这个想法是好的。

Answer 1

我将在这里发布我发布到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 只是为了确保没有任何东西被优化掉，并展示一个小点；你会得到不同数量的线程，因为每个线程都有自己的线程数作为种子；如果您使用相同数量的线程一次又一次地运行相同的代码，出于相同的原因，您将获得相同的总和。无论如何，时间安排（在没有其他用户的 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 倍以上的加速。

Answer 2

您不能在多个线程中使用 C rand() 函数；这会导致未定义的行为。某些实现可能会给您锁定（这会使其变慢）；其他人可能允许线程破坏彼此的状态，可能会导致程序崩溃或只是给出“坏”随机数。

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

Answer 3

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

Answer 4

似乎rand 在 Linux 上的所有线程之间具有全局共享状态，而在 Windows 上则具有线程本地存储状态。由于必要的同步，Linux 上的共享状态导致您的速度变慢。

我认为 C 库中没有一种可移植的方式来在多个线程上并行使用 RNG，因此您需要另一种方式。您可以使用Mersenne Twister。正如marcog所说，您需要以不同的方式初始化每个线程的种子。

Answer 5

在 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 个不同的伪随机数，没有任何互斥锁或竞争条件。

Answer 6

随机数的生成速度非常快，因此通常内存会成为瓶颈。通过在多个线程之间划分此任务，您会产生额外的通信和同步开销（并且不同内核的缓存的同步并不便宜）。

最好使用具有更好random()函数的单线程。