I want to generate a reproducible sequence of random numbers in parallel with OpenMP. I wrote a small code that creates for each thread its own RNG with different seed. But when I start generating random numbers in parallel, something strange happens. With uniform distribution I get reproducible results, but with normal distribution i get different results!
More specifically, every time I run the program, the array seems to be shuffled.
My question: Why this might be happening? And how do I resolve this issue?
#include <iostream>
#include <random>
#include <omp.h>
int main() {
std::normal_distribution<> normal(0,1);
std::uniform_real_distribution<> uniform(0,1);
std::mt19937 *rng = new std::mt19937[omp_get_max_threads()];
for(int i = 0; i < omp_get_max_threads(); ++i) {
rng[i].seed(1+i);
}
int N = 100;
double* arr = new double[N];
#pragma omp parallel for
for(int i = 0; i < N; ++i) {
arr[i] = normal(rng[omp_get_thread_num()]);
}
for(int i = 0; i < N; ++i) {
std::cout << arr[i] << std::endl;
}
delete[] arr;
delete[] rng;
return 0;
}
Some implementations of std::normal_distribution generate variates two-at-a-time. They store the second variate in the distribution object and return it the next time it's called. Having multiple threads accessing the same normal distribution is the likely the reason for seeing non-deterministic results.
std::uniform_real_distribution tends not to have any hidden state and hence call the underlying bit generator every time so doesn't exhibit this behaviour.
The existence of RandomNumberDistribution.reset() method hints that it's not safe to share objects between threads.
As suggested in the comments, I'd suggest editing your code so that it doesn't call rand() if you want deterministic results. The way the code is written and what you're asking has different intents and caused me (and other posters) to get confused.