Randn.h

Go to the documentation of this file.

// Struct Randn
//   This struct can be used for sampling normal distribution function
//   on unbounded domain.
//
#ifndef IPPL_RANDN_H
#define IPPL_RANDN_H
 
#include "Random/Distribution.h"
#include "Random/Utility.h"
 
namespace ippl {
    namespace random {

        template <typename T, unsigned Dim>
        struct randn {
            using view_type =
                typename ippl::detail::ViewType<ippl::Vector<double, Dim>, 1>::view_type;
            using GeneratorPool = typename Kokkos::Random_XorShift64_Pool<>;
 
            // Output View for the random numbers
            view_type v;
 
            // The GeneratorPool
            GeneratorPool rand_pool;
 
            T mu[Dim];
            T sd[Dim];
            KOKKOS_INLINE_FUNCTION randn(view_type v_, GeneratorPool rand_pool_, T* mu_p, T* sd_p)
                : v(v_)
                , rand_pool(rand_pool_) {
                for (unsigned int i = 0; i < Dim; i++) {
                    mu[i] = mu_p[i];
                    sd[i] = sd_p[i];
                }
            }
 
            KOKKOS_INLINE_FUNCTION randn(view_type v_, GeneratorPool rand_pool_)
                : v(v_)
                , rand_pool(rand_pool_) {
                for (unsigned int i = 0; i < Dim; i++) {
                    mu[i] = 0.0;
                    sd[i] = 1.0;
                }
            }

            KOKKOS_INLINE_FUNCTION const T& getMu(unsigned int idx) const { return mu[idx]; }

            KOKKOS_INLINE_FUNCTION const T& getSd(unsigned int idx) const { return sd[idx]; }

            KOKKOS_INLINE_FUNCTION void operator()(const size_t i) const {
                // Get a random number state from the pool for the active thread
                typename GeneratorPool::generator_type rand_gen = rand_pool.get_state();
 
                for (unsigned d = 0; d < Dim; ++d) {
                    v(i)[d] = mu[d] + sd[d] * rand_gen.normal(0.0, 1.0);
                }
 
                // Give the state back, which will allow another thread to acquire it
                rand_pool.free_state(rand_gen);
            }
        };
    }  // namespace random
}  // namespace ippl
 
#endif