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()始终返回1。请改用omp_get_max_threads()，它将返回任何即将到来的parallel区域的线程数，除非手动覆盖。特别是threadvec只有一个条目。

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

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

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

尽管从文档中看，这是否是线程安全的，但并不是 100% 清楚，因此只需使用串行循环版本。

以不同的方式正确播种您的 rng

默认情况下，所有 rng 都使用相同的数字播种，因此显然它们将返回完全相同的序列。用线号正确播种，例如gsl_rng_set(threadvec[b], b * 101);.请注意，陶斯沃思发电机很奇怪。当播种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;
}