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()の

外側は常に1を返す並列領域のうちomp_get_num_threadsを使用

。代わりにomp_get_max_threads()を使用すると、手作業で上書きされない限り、今後のparallel領域のスレッド数が返されます。特にthreadvecにはエントリが1つしかありません。

は正常に動作しません並列領域で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); 

Allthoughをドキュメントからそれがスレッドセーフであるならば、それは明らかに100％ではないので、単にシリアルループバージョンを使用します。

適切にデフォルトでは、異なる

あなたのRNGをシード、すべてのRNGは同じ番号を播種しているので、明らかに彼らは正確に同一の配列を返します。スレッド番号で適切にシードします。 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; 
}