Generating random numbers in parallel
Hi
I generate a series of random numbers in parallel (using OpenMP), but depending on what number of threads I invoke, I get a different result. From that I conclude that I have made an error somewhere!
Here is the MWE, which generates a number between 0..1 and increments a variable if the generated variable is larger than 0.5:
Code:
#include <random>
typedef std::uniform_real_distribution<double> distr_uni;
#define max_threads 1
using namespace std;
int main(int argc, char* argv[])
{
int reservoir_counter, accepted_tube=0;
double r;
omp_set_num_threads(max_threads);
#pragma omp parallel
{
mt19937 eng(0);
distr_uni uniform(0, 1);
#pragma omp for private(r, reservoir_counter) reduction(+:accepted_tube)
for(reservoir_counter=0; reservoir_counter<100000; reservoir_counter++)
{
r = uniform(eng);
if(r>0.5)
{
accepted_tube++;
}
}
}
cout << accepted_tube << endl;
return 0;
}
When I set max_threads=1 I get 50027, but when max_threads=60 (on a machine that supports it....) I get 50440.
The sensitive RNG and its engine I have declared within the parallelized area, so it's not really clear to me where the error can possibly be.
Can someone spot my error that is apparently there?