From de03e259fde547b82da18b8ec3ea1df597a6f85a Mon Sep 17 00:00:00 2001 From: conrad Date: Wed, 7 Dec 2022 13:45:31 +1000 Subject: [PATCH] ensure single-threaded and multi-threaded RNG results are the same --- include/armadillo_bits/arma_rng.hpp | 177 +--------------------------- 1 file changed, 3 insertions(+), 174 deletions(-) diff --git a/include/armadillo_bits/arma_rng.hpp b/include/armadillo_bits/arma_rng.hpp index 08d77d26..d5097db4 100644 --- a/include/armadillo_bits/arma_rng.hpp +++ b/include/armadillo_bits/arma_rng.hpp @@ -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::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 > 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()) ); - } - - 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& 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& 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::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 > inline static void - fill_simple(std::complex* mem, const uword N) + fill(std::complex* mem, const uword N) { #if defined(ARMA_RNG_ALT) { @@ -796,68 +745,6 @@ struct arma_rng::randn< std::complex > } - inline - static - void - fill(std::complex* mem, const uword N) - { - #if defined(ARMA_USE_OPENMP) - { - if((N < 512) || omp_in_parallel()) { arma_rng::randn< std::complex >::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 > 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()) ); - } - - 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& 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(val1, val2); - } - } - - std::mt19937_64& t0_engine = engine[0]; - std::normal_distribution& 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(val1, val2); - } - } - #else - { - arma_rng::randn< std::complex >::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::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 distr_type; - - const uword n_threads = uword( mp_thread_limit::get() ); - - std::vector g_motor(n_threads); - std::vector 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()) ); - - 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::fill_simple(mem, N, a, b); - } - #endif - } - };