AMReX_Random.H

Go to the documentation of this file.

#ifndef AMREX_RAND_H
#define AMREX_RAND_H
#include <AMReX_Config.H>
 
#include <AMReX.H>
#include <AMReX_GpuQualifiers.H>
#include <AMReX_ParallelDescriptor.H>
#include <AMReX_RandomEngine.H>
#include <limits>
#include <cmath>
#include <cstdint>
 
namespace amrex
{
    Real Random ();
 
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    Real Random (RandomEngine const& random_engine)
    {
#if defined (__SYCL_DEVICE_ONLY__)
        mkl::rng::device::uniform<Real> distr;
        return mkl::rng::device::generate(distr, *random_engine.engine);
#else
#ifdef BL_USE_FLOAT
        AMREX_IF_ON_DEVICE((
                AMREX_HIP_OR_CUDA(
                        return 1.0f - hiprand_uniform(random_engine.rand_state); ,
                        return 1.0f - curand_uniform(random_engine.rand_state);
                )
        ))
#else
        AMREX_IF_ON_DEVICE((
                AMREX_HIP_OR_CUDA(
                        return 1.0 - hiprand_uniform_double(random_engine.rand_state); ,
                        return 1.0 - curand_uniform_double(random_engine.rand_state);
                )
        ))
#endif
        AMREX_IF_ON_HOST((
                amrex::ignore_unused(random_engine);
                return Random();
        ))
#endif
    }
 
    Real RandomNormal (Real mean, Real stddev);
 
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    Real RandomNormal (Real mean, Real stddev, RandomEngine const& random_engine)
    {
#if defined (__SYCL_DEVICE_ONLY__)
        mkl::rng::device::gaussian<Real> distr(mean, stddev);
        return mkl::rng::device::generate(distr, *random_engine.engine);
#else
#ifdef BL_USE_FLOAT
        AMREX_IF_ON_DEVICE((
                AMREX_HIP_OR_CUDA(
                        return stddev * hiprand_normal(random_engine.rand_state) + mean; ,
                        return stddev * curand_normal(random_engine.rand_state) + mean;
                )
        ))
#else
        AMREX_IF_ON_DEVICE((
                AMREX_HIP_OR_CUDA(
                        return stddev * hiprand_normal_double(random_engine.rand_state) + mean; ,
                        return stddev * curand_normal_double(random_engine.rand_state) + mean;
                )
        ))
#endif
        AMREX_IF_ON_HOST((
                amrex::ignore_unused(random_engine);
                return RandomNormal(mean, stddev);
        ))
#endif
    }
 
    unsigned int RandomPoisson (Real lambda);
 
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    unsigned int RandomPoisson (Real lambda, RandomEngine const& random_engine)
    {
#if defined (__SYCL_DEVICE_ONLY__)
        mkl::rng::device::poisson<unsigned int> distr(lambda);
        return mkl::rng::device::generate(distr, *random_engine.engine);
#else
        AMREX_IF_ON_DEVICE((
                AMREX_HIP_OR_CUDA(
                        return hiprand_poisson(random_engine.rand_state, lambda); ,
                        return curand_poisson(random_engine.rand_state, lambda);
                )
        ))
        AMREX_IF_ON_HOST((
                amrex::ignore_unused(random_engine);
                return RandomPoisson(lambda);
        ))
#endif
    }
 
    Real RandomGamma (Real alpha, Real beta);
 
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    Real RandomGamma (Real alpha, Real beta, RandomEngine const& random_engine)
    {
        AMREX_ASSERT(alpha > 0);
        AMREX_ASSERT(beta > 0);
 
        AMREX_IF_ON_DEVICE((
        if (alpha < 1)
        {
            Real u = amrex::Random(random_engine);
            return RandomGamma(1.0_rt + alpha, beta, random_engine) * std::pow(u, 1.0_rt / alpha);
        }
 
        {
            Real x, v, u;
            Real d = alpha - 1.0_rt / 3.0_rt;
            Real c = (1.0_rt / 3.0_rt) / std::sqrt(d);
 
            while (1) {
                do {
                    x = amrex::RandomNormal(0.0_rt, 1.0_rt, random_engine);
                    v = 1.0_rt + c * x;
                } while (v <= 0.0_rt);
 
                v = v * v * v;
                u = amrex::Random(random_engine);
 
                if (u < 1.0_rt - 0.0331_rt * x * x * x * x) {
                    break;
                }
 
                if (std::log(u) < 0.5_rt * x * x + d * (1.0_rt - v + std::log(v))) {
                    break;
                }
            }
            return beta * d * v;
        }
        ))
 
        AMREX_IF_ON_HOST((
                amrex::ignore_unused(random_engine);
                return RandomGamma(alpha, beta);
        ))
 
    }
 
    unsigned int Random_int (unsigned int n); // [0,n-1]
 
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    unsigned int Random_int (unsigned int n, RandomEngine const& random_engine)
    {
#if defined(__SYCL_DEVICE_ONLY__)
        mkl::rng::device::uniform<unsigned int> distr(0,n);
        return mkl::rng::device::generate(distr, *random_engine.engine);
#else
        AMREX_IF_ON_DEVICE((
                unsigned int rand;
                constexpr unsigned int RAND_M = 4294967295; // 2**32-1
                do {
                AMREX_HIP_OR_CUDA( rand = hiprand(random_engine.rand_state);,
                                rand =  curand(random_engine.rand_state) );
                } while (rand > (RAND_M - RAND_M % n));
                return rand % n;
        ))
        AMREX_IF_ON_HOST((
                amrex::ignore_unused(random_engine);
                return Random_int(n);
        ))
#endif
    }
 
    ULong Random_long (ULong n); // [0,n-1]
 
    void FillRandom (Real* p, Long N);
 
    void FillRandomNormal (Real* p, Long N, Real mean, Real stddev);
 
    namespace detail {
        inline ULong DefaultGpuSeed () {
            return ParallelDescriptor::MyProc()*1234567ULL + 12345ULL;
        }
    }
 
    void InitRandom (ULong cpu_seed, int nprocs=ParallelDescriptor::NProcs(),
                     ULong gpu_seed = detail::DefaultGpuSeed());
 
    void ResetRandomSeed (ULong cpu_seed, ULong gpu_seed = detail::DefaultGpuSeed());
 
    void SaveRandomState (std::ostream& os);
 
    void RestoreRandomState (std::istream& is, int nthreads_old, int nstep_old);
 
    void UniqueRandomSubset (Vector<int> &uSet, int setSize, int poolSize,
                             bool printSet = false);
 
    void DeallocateRandomSeedDevArray ();
}
 
#endif