GSL+OMP：C++ 中的线程安全随机数生成器

Question

我有一个代码，我试图在其中并行执行。

#include<iostream>
#include<omp.h>
#include<math.h>
#include<cstdlib>
#include<iterator>
#include<string.h>
#include<vector>
#include<map>
#include<time.h>
#include<gsl/gsl_rng.h>
#include<gsl/gsl_randist.h>

gsl_rng ** threadvec = new gsl_rng*[omp_get_num_threads()];
using namespace std;

int main(){
   clock_t begin = omp_get_wtime();
   vector<double> PopVals;
   map<int, vector<double> > BigMap;
   int Num1 = 100; 
   double randval;
   int Num2 = 10; 
   #pragma omp parallel
   {
       gsl_rng_env_setup();     
       for (int b = 0; b < omp_get_num_threads(); b++)
           threadvec[b] = gsl_rng_alloc(gsl_rng_taus);  
   }
   for( int i = 0; i < Num1; i++){
       PopVals.resize(Num2);
       #pragma omp parallel for
          for( int j = 0; j < Num2; j++){   
              randval = gsl_rng_uniform(threadvec[omp_get_thread_num()]);   
              PopVals[j] = randval; 
          }
       BigMap.insert(make_pair(i,PopVals));
       PopVals.clear();
   }

map<int,vector<double> >::iterator it = BigMap.find(Num1-1);
vector<double> OutVals = it->second; 

for (int i = 0; i < Num2; i++)
    cout << endl << OutVals[i] << endl; 

for (int b = 0; b < omp_get_num_threads(); b++)
        gsl_rng_free(threadvec[b]);

clock_t end = omp_get_wtime(); 
double elapsed_time = double(end - begin);
cout << endl << "Time taken to run: " << elapsed_time <<  " secs" << endl;
}

当我 运行 这样做时，有 8 个线程并行执行嵌套循环，但我一直看到每个线程的随机数相同。我将这种行为归因于每次迭代都没有设置种子。如果有人能指出，我如何以线程安全的方式在循环的每次迭代中生成唯一的随机数，那就太好了。

以上代码输出0.793816，10次。然而，我希望内部循环中的每个值都有唯一的数字。

谢谢。

Answer 1

这里有多个问题。

使用 omp_get_num_threads 个平行区域

在平行区域之外，omp_get_num_threads() 总是 returns 1。请改用 omp_get_max_threads()，它将 return 任何即将到来的 parallel 区域的线程数，除非手动覆盖。特别是 threadvec 只有一个条目。

不要在并行区域初始化环境

在并行区域中调用 gsl_rng_env_setup 将无法正常工作。您还试图通过所有线程分配整个 rngs 向量...只需删除并行区域并正确使用 omp_get_max_threads() 即可。或者你也可以这样做：

gsl_rng_env_setup(); // serial
#pragma omp parallel
threadvec[omp_get_thread_num()] = gsl_rng_alloc(gsl_rng_taus);

尽管从文档中还不能 100% 清楚这是否是线程安全的，所以只需使用串行循环版本。

以不同的方式正确播种你的 rngs

默认情况下，所有 rng 都使用相同的数字作为种子，因此很明显它们将 return 完全相同的序列。用线程号正确地播种它们，例如gsl_rng_set(threadvec[b], b * 101);。请注意，Tausworthe 生成器很奇怪。当用 0 或 1.

播种时，那些特定的会生成相同的数字序列

隐式共享变量

您的变量 randval 是在并行区域之外定义的，因此它是隐式共享的。您可以强制它是私有的，但最好尽可能在本地声明变量。这使得对 OpenMP 代码的推理变得更加容易。

最后看起来像这样：

#include <cstdlib>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <iostream>
#include <iterator>
#include <map>
#include <math.h>
#include <omp.h>
#include <string.h>
#include <time.h>
#include <vector>

// DO NOT using namespace std;

int main() {
  clock_t begin = omp_get_wtime();
  std::vector<double> PopVals;
  std::map<int, std::vector<double>> BigMap;
  constexpr int Num1 = 100;
  constexpr int Num2 = 10;
  gsl_rng_env_setup();
  gsl_rng **threadvec = new gsl_rng *[omp_get_max_threads()];
  for (int b = 0; b < omp_get_max_threads(); b++) {
    threadvec[b] = gsl_rng_alloc(gsl_rng_taus);
    gsl_rng_set(threadvec[b], b * 101);
  }
  for (int i = 0; i < Num1; i++) {
    PopVals.resize(Num2);
    #pragma omp parallel for
    for (int j = 0; j < Num2; j++) {
      double randval = gsl_rng_uniform(threadvec[omp_get_thread_num()]);
      PopVals[j] = randval;
    }
    BigMap.insert(std::make_pair(i, PopVals));
    PopVals.clear();
  }

  std::map<int, std::vector<double>>::iterator it = BigMap.find(Num1 - 1);
  std::vector<double> OutVals = it->second;

  for (int i = 0; i < Num2; i++)
    std::cout << std::endl << OutVals[i] << std::endl;

  for (int b = 0; b < omp_get_max_threads(); b++)
    gsl_rng_free(threadvec[b]);

  clock_t end = omp_get_wtime();
  double elapsed_time = double(end - begin);
  std::cout << std::endl << "Time taken to run: " << elapsed_time << " secs" << std::endl;
  delete[] threadvec;
}