GSL + OMP: générateurs de nombres aléatoires de thread sécurisé en C++

Question

J'ai un code dans lequel j'essaie d'exécuter en parallèle.

#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; 
} 

Quand je lance cela, il y a 8 fils d'exécution de la boucle imbriquée en parallèle, mais je continue à voir le même nombre aléatoire pour chaque thread. J'ai attribué ce comportement à l'absence de mise en graine, à chaque itération. Ce serait génial si quelqu'un pouvait le faire remarquer, comment puis-je générer des nombres aléatoires uniques dans chaque itération de la boucle d'une manière thread sûre.

La sortie du code ci-dessus est de 0,793816, 10 fois. Alors que je veux des nombres uniques pour chacune des valeurs de la boucle interne.

Merci.

Answer 1

Il y a plusieurs problèmes ici.

L'utilisation omp_get_num_threads de régions parallèles

En dehors d'une région parallèle, omp_get_num_threads() retourne toujours 1. Utilisez omp_get_max_threads() à la place, il retournera le nombre de threads pour toute région à venir parallel à moins d'être surchargé manuellement. Surtout threadvec a une seule entrée.

Ne pas initialiser l'environnement dans une région parallèle

Appel gsl_rng_env_setup dans une région parallèle ne fonctionnera pas correctement. Aussi vous essayez d'allouer le vecteur entier de rngs par tous les threads ... Supprimez simplement la région parallèle et utilisez omp_get_max_threads() correctement. Ou vous pouvez aussi le faire:

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

Allthough de la documentation, il est pas 100% clair si cela est threadsafe, donc il suffit d'utiliser la version en boucle série.

correctement vos ensemencer RNG différemment

Par défaut, tous les RNG sont ensemencés avec le même nombre, donc évidemment ils retourneront exactement la même séquence. Les graver correctement avec le numéro de fil, par ex. gsl_rng_set(threadvec[b], b * 101);. Notez que les générateurs Tausworthe sont étranges. Ceux-ci génèrent la même séquence de nombres lorsqu'ils sont ensemencés avec 0 ou 1.

des variables partagées Implicitement

votre variable randval est définie en dehors de la région parallèle, d'où il est implicitement partagé. Vous pouvez le forcer à être privé, mais il est préférable de déclarer les variables aussi localement que possible. Cela rend le raisonnement sur le code OpenMP beaucoup plus facile.

À la fin, il ressemble à ceci:

#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; 
}