phynexis-deps/armadillo-code
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

ensure single-threaded and multi-threaded RNG results are the same

Browse Source
This commit is contained in:
conrad
2022-12-07 13:45:31 +10:00
parent ffd570ef2f
commit de03e259fd
1 changed files with 3 additions and 174 deletions
Download Patch File Download Diff File Expand all files Collapse all files
include/armadillo_bits/arma_rng.hpp
+3 -174
View File
@@ -593,7 +593,7 @@ struct arma_rng::randn
inline
static
void
fill_simple(eT* mem, const uword N)
fill(eT* mem, const uword N)
{
#if defined(ARMA_RNG_ALT)
{
@@ -628,56 +628,6 @@ struct arma_rng::randn
}
inline
static
void
fill(eT* mem, const uword N)
{
#if defined(ARMA_USE_OPENMP)
{
if((N < 1024) || omp_in_parallel()) { arma_rng::randn<eT>::fill_simple(mem, N); return; }
typedef typename std::mt19937_64::result_type local_seed_type;
const uword n_threads = uword( mp_thread_limit::get() );
std::vector< std::mt19937_64 > engine(n_threads);
std::vector< std::normal_distribution<double> > distr(n_threads);
for(uword t=0; t < n_threads; ++t)
{
std::mt19937_64& t_engine = engine[t];
t_engine.seed( local_seed_type(t) + local_seed_type(arma_rng::randi<local_seed_type>()) );
}
const uword chunk_size = N / n_threads;
#pragma omp parallel for schedule(static) num_threads(int(n_threads))
for(uword t=0; t < n_threads; ++t)
{
const uword start = (t+0) * chunk_size;
const uword endp1 = (t+1) * chunk_size;
std::mt19937_64& t_engine = engine[t];
std::normal_distribution<double>& t_distr = distr[t];
for(uword i=start; i < endp1; ++i) { mem[i] = eT( t_distr(t_engine)); }
}
std::mt19937_64& t0_engine = engine[0];
std::normal_distribution<double>& t0_distr = distr[0];
for(uword i=(n_threads*chunk_size); i < N; ++i) { mem[i] = eT( t0_distr(t0_engine)); }
}
#else
{
arma_rng::randn<eT>::fill_simple(mem, N);
}
#endif
}
inline
static
void
@@ -694,7 +644,6 @@ struct arma_rng::randn
mem[i] = (val * sd) + mu;
}
}
};
@@ -745,7 +694,7 @@ struct arma_rng::randn< std::complex<T> >
inline
static
void
fill_simple(std::complex<T>* mem, const uword N)
fill(std::complex<T>* mem, const uword N)
{
#if defined(ARMA_RNG_ALT)
{
@@ -796,68 +745,6 @@ struct arma_rng::randn< std::complex<T> >
}
inline
static
void
fill(std::complex<T>* mem, const uword N)
{
#if defined(ARMA_USE_OPENMP)
{
if((N < 512) || omp_in_parallel()) { arma_rng::randn< std::complex<T> >::fill_simple(mem, N); return; }
typedef typename std::mt19937_64::result_type local_seed_type;
const uword n_threads = uword( mp_thread_limit::get() );
std::vector< std::mt19937_64 > engine(n_threads);
std::vector< std::normal_distribution<double> > distr(n_threads);
for(uword t=0; t < n_threads; ++t)
{
std::mt19937_64& t_engine = engine[t];
t_engine.seed( local_seed_type(t) + local_seed_type(arma_rng::randi<local_seed_type>()) );
}
const uword chunk_size = N / n_threads;
#pragma omp parallel for schedule(static) num_threads(int(n_threads))
for(uword t=0; t < n_threads; ++t)
{
const uword start = (t+0) * chunk_size;
const uword endp1 = (t+1) * chunk_size;
std::mt19937_64& t_engine = engine[t];
std::normal_distribution<double>& t_distr = distr[t];
for(uword i=start; i < endp1; ++i)
{
const T val1 = T( t_distr(t_engine) );
const T val2 = T( t_distr(t_engine) );
mem[i] = std::complex<T>(val1, val2);
}
}
std::mt19937_64& t0_engine = engine[0];
std::normal_distribution<double>& t0_distr = distr[0];
for(uword i=(n_threads*chunk_size); i < N; ++i)
{
const T val1 = T( t0_distr(t0_engine) );
const T val2 = T( t0_distr(t0_engine) );
mem[i] = std::complex<T>(val1, val2);
}
}
#else
{
arma_rng::randn< std::complex<T> >::fill_simple(mem, N);
}
#endif
}
inline
static
void
@@ -888,7 +775,7 @@ struct arma_rng::randg
inline
static
void
fill_simple(eT* mem, const uword N, const double a, const double b)
fill(eT* mem, const uword N, const double a, const double b)
{
#if defined(ARMA_USE_EXTERN_RNG)
{
@@ -909,64 +796,6 @@ struct arma_rng::randg
}
#endif
}
inline
static
void
fill(eT* mem, const uword N, const double a, const double b)
{
#if defined(ARMA_USE_OPENMP)
{
if((N < 512) || omp_in_parallel()) { arma_rng::randg<eT>::fill_simple(mem, N, a, b); return; }
typedef std::mt19937_64 motor_type;
typedef std::mt19937_64::result_type ovum_type;
typedef std::gamma_distribution<double> distr_type;
const uword n_threads = uword( mp_thread_limit::get() );
std::vector<motor_type> g_motor(n_threads);
std::vector<distr_type> g_distr(n_threads);
const distr_type g_distr_base(a,b);
for(uword t=0; t < n_threads; ++t)
{
motor_type& g_motor_t = g_motor[t];
distr_type& g_distr_t = g_distr[t];
g_motor_t.seed( ovum_type(t) + ovum_type(arma_rng::randi<ovum_type>()) );
g_distr_t.param( g_distr_base.param() );
}
const uword chunk_size = N / n_threads;
#pragma omp parallel for schedule(static) num_threads(int(n_threads))
for(uword t=0; t < n_threads; ++t)
{
const uword start = (t+0) * chunk_size;
const uword endp1 = (t+1) * chunk_size;
motor_type& g_motor_t = g_motor[t];
distr_type& g_distr_t = g_distr[t];
for(uword i=start; i < endp1; ++i) { mem[i] = eT( g_distr_t(g_motor_t)); }
}
motor_type& g_motor_0 = g_motor[0];
distr_type& g_distr_0 = g_distr[0];
for(uword i=(n_threads*chunk_size); i < N; ++i) { mem[i] = eT( g_distr_0(g_motor_0)); }
}
#else
{
arma_rng::randg<eT>::fill_simple(mem, N, a, b);
}
#endif
}
};