Thread: 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?