在 OMP 中使用线程生成相同的随机数

Question

我正在尝试使用 OMP 对一些代码进行多线程处理。目前我的顺序版本使用 rand() 生成一组具有一致种子的随机数，以便它们在每次运行时返回相同的结果。我想并行化我的代码，但 rand() 不是线程安全的。有人可以告诉我如何使用在线程上工作的随机数生成器，这样我就可以在每次测试时生成相同的数据集，类似于使用带有 rand() 的种子。我的并行化代码如下：

    long linkCnt;
    long j;
    long t;
    srand(randInit);

    linkCnt=0; //Keep track of # of outgoing edges
 #pragma omp parallel for schedule(runtime) private(t)
    for(j = 0; j < G->N; j++)
    {

        if(i == j){
            t = NO_CONN;
        } else {
            t = (rand() % ((MAX_EDGE_WEIGTH-1) * 2)+1); //50% of having no connection
            if(t > MAX_EDGE_WEIGTH){
                //t = INF; //Like no connection
                t = NO_CONN; //Like no connection
            } else {
                linkCnt++;
                G->visited[j] = VISITED; //Do this to find isolated nods that have no incomming edge
            }
        }

        G->node[i][j] = t;
    }

Answer 1

似乎有几个问题在这里混为一谈。

首先，rand() 函数的非线程安全特性意味着从不同线程同时调用 rand() 可能会产生与顺序调用不同的值。用一个简单的例子来解释这一点可能是最容易的，所以让我们看一下 PCG 的 32 位版本，因为它既好又简单。它有一个 32 位的状态，并生成这样的 32 位数字：

static uint32_t state = ...;

static uint32_t generate(void) {
  uint32_t s = state;
  uint32_t v = ((s >> ((s >> 28) + 4)) ^ s) * (277803737U);
  v ^= v >> 22;
  state = state * 747796405U + 1729U;
  return v;
}

现在想想如果两个线程几乎同时调用generate() 会发生什么。也许它们都为state 得到相同的值，因此两次生成相同的随机数。也许一个人在另一个人阅读之前更新了state，所以他们得到了不同的值。

我们可以通过使用互斥体保护generate() 函数来消除该问题，或者在32 位PGC 的情况下（这就是为什么I use it 用于可再现数字）使用原子。如果我们这样做，那么我们将始终以相同的顺序获得相同的数字。

问题的第二部分是当调用者的顺序在您的代码中混淆时会发生什么。假设您有两个线程（称为 A 和 B），它们每个都必须运行循环的两次迭代。即使您从线程安全源获取随机数，调用的顺序也可能是 AABB、ABAB、ABBA、BBAA、BABA 或 BAAB，每个调用顺序都会导致您的代码生成不同的结果。

有几种直接的方法可以解决这个问题。首先，您可以使用同步原语来确保每次迭代都以您想要的顺序调用generate。最简单的方法可能是使用队列，但是您会浪费大量时间进行同步，并且会失去一些并行处理的机会（并且您必须大量重组代码）。

如果您的迭代次数相对较少，您可能会考虑提前生成一个数组。想想：

int i;
int nums[LEN];
for (i = 0 ; i < LEN ; i++)
  nums[i] = generate();
#pragma omp parallel for ...
for (i = 0 ; i < LEN ; i++) {
  do_stuff(nums[i]);
}

不过，更好的解决方案可能是完全放弃生成随机数的想法，转而使用散列。 https://burtleburtle.net/bob/hash/integer.html 有一些选择。例如

for (int i = 0 ; i < LEN ; i++) {
  do_stuff(hash(i));
}

当然你可以加入一些盐，甚至可以使用rand() 来生成你的盐。

Answer 2

这是一种基于块的方法，它将问题空间划分为 N/BLOCK_SIZE 个块，并使用您的 randInit + 每个块的块号重新播种 RNG。无论您拥有多少线程，这都会提供可重现的输出。它还为 N + x 序列生成相同的初始 N 个数字。这只要你保持相同的 BLOCK_SIZE。

一个好的块大小可能类似于典型的 N / (max_num_procs * 2)。但仍有试验空间。

#include <omp.h>
#include <stdio.h>
#include <stdlib.h>

#define N_DEFAULT  48   //Default number of nodes 
#define BLOCK_SIZE 12   //BLOCK SIZE number of nodes per block.
                        //Changes this changes reseed frequencey, 
                        //.. and hence the generated sequence
#define randInit   42   //Had to be something.

int main(int argc , char* argv[])
{
    int N=N_DEFAULT;
    if (argc >1)
        N=atoi(argv[1]); 

    int rands[N];// keep our random numbers for sequential debug output
    int n=BLOCK_SIZE;
    int num_blocks=(N+BLOCK_SIZE-1)/ BLOCK_SIZE;  // ceil(N/BLOCK_SIZE)
    int nt=omp_get_max_threads();   


    printf("          N: %d\n",N);
    printf("     Blocks: %d, (size: %d)\n",num_blocks,n);
    printf("    threads: %d\n",nt);

    //Parallel random generation
    #pragma omp parallel for schedule(runtime) 
    for (int J=0;J<num_blocks;J++)
    {
        int block_seed=randInit+J;      // unique block seed
        int start = J * n;
        int end= start+n > N?N:start+n;

        int tid = omp_get_thread_num(); // Just for debug
        printf("thread %d: works on block %d (%d - %d )\n",tid,J,start,end);

        for (int j=start; j < end;j++)
        {           
            int t=rand_r(&block_seed); //rand_r provides thread safe (re-entrant rand)
            rands[j]=t;
        }
    }

    //Output for debug single thread
    for (int j=0; j < N;j++)
    {
        printf("%d : %d \n",j,rands[j]);
    }   
    return 0;
}

输出不同的N，线程数如下所示。

N: 24                                   N: 27
Blocks: 3, (size: 8)                    Blocks: 4, (size: 8)
threads: 4                              threads: 1
-------------------------------------|-------------------------------
thread 1: works on block 1 (8 - 16 )    thread 0: works on block 0 (0 - 8 )
thread 2: works on block 2 (16 - 24 )   thread 0: works on block 1 (8 - 16 )
thread 0: works on block 0 (0 - 8 )     thread 0: works on block 2 (16 - 24 )
                                        thread 0: works on block 3 (24 - 27 )
-------------------------------------|-------------------------------
0 : 681191333                           0 : 681191333
1 : 928546885                           1 : 928546885
2 : 1457394273                          2 : 1457394273
3 : 941445650                           3 : 941445650
4 : 2129613237                          4 : 2129613237
5 : 1661015563                          5 : 1661015563
6 : 2071432601                          6 : 2071432601
7 : 222443696                           7 : 222443696
8 : 1156886562                          8 : 1156886562
9 : 398918689                           9 : 398918689
10 : 170756699                         10 : 170756699
11 : 703115845                         11 : 703115845
12 : 1424182583                        12 : 1424182583
13 : 1516198481                        13 : 1516198481
14 : 1740837599                        14 : 1740837599
15 : 1148851528                        15 : 1148851528
16 : 1633630368                        16 : 1633630368
17 : 2015727614                        17 : 2015727614
18 : 1031602773                        18 : 1031602773
19 : 463737465                         19 : 463737465
20 : 720848057                         20 : 720848057
21 : 1369285272                        21 : 1369285272
22 : 1411290150                        22 : 1411290150
23 : 2074210785                        23 : 2074210785
-------------------------------------|-------------------------------
                                       24 : 2109326622
                                       25 : 1486099418
                                       26 : 1892448847