alpaka
Abstraction Library for Parallel Kernel Acceleration
RandDefault.hpp
Go to the documentation of this file.
1 /* Copyright 2022 Jeffrey Kelling, Jan Stephan, Bernhard Manfred Gruber
2  * SPDX-License-Identifier: MPL-2.0
3  */
4 
5 #pragma once
6 
7 #include "alpaka/core/Common.hpp"
8 #include "alpaka/math/Traits.hpp"
10 #include "alpaka/rand/Traits.hpp"
11 
12 #include <algorithm>
13 #include <limits>
14 #include <type_traits>
15 
16 namespace alpaka::rand
17 {
18  class RandDefault : public concepts::Implements<ConceptRand, RandDefault>
19  {
20  };
21 
22  namespace distribution::gpu
23  {
24  namespace detail
25  {
26  template<typename TFloat>
27  struct BitsType;
28 
29  template<>
30  struct BitsType<float>
31  {
32  using type = std::uint32_t;
33  };
34 
35  template<>
36  struct BitsType<double>
37  {
38  using type = std::uint64_t;
39  };
40  } // namespace detail
41 
42  //! The GPU random number normal distribution.
43  template<typename T>
45  {
46  static_assert(std::is_integral_v<T>, "Return type of UniformUint must be integral.");
47 
48  public:
49  UniformUint() = default;
50 
51  template<typename TEngine>
52  ALPAKA_FN_HOST_ACC auto operator()(TEngine& engine) -> T
53  {
54  using BitsT = typename TEngine::result_type;
55  T ret = 0;
56  constexpr auto N = sizeof(T) / sizeof(BitsT);
57  for(unsigned int a = 0; a < N; ++a)
58  {
59  ret
60  ^= (static_cast<T>(engine())
61  << (sizeof(BitsT) * std::numeric_limits<unsigned char>::digits * a));
62  }
63  return ret;
64  }
65  };
66 
67  //! The GPU random number uniform distribution.
68  template<typename T>
70  {
71  static_assert(std::is_floating_point_v<T>, "Return type of UniformReal must be floating point.");
72 
73  using BitsT = typename detail::BitsType<T>::type;
74 
75  public:
76  UniformReal() = default;
77 
78  template<typename TEngine>
79  ALPAKA_FN_HOST_ACC auto operator()(TEngine& engine) -> T
80  {
81  constexpr BitsT limit = static_cast<BitsT>(1) << std::numeric_limits<T>::digits;
82  BitsT const b = UniformUint<BitsT>()(engine);
83  auto const ret = static_cast<T>(b & (limit - 1)) / limit;
84  return ret;
85  }
86  };
87 
88  /*! The GPU random number normal distribution.
89  *
90  * \note
91  * This type contains state and is not thread-safe: To be used
92  * per thread, not shared.
93  *
94  * \note When reproducibility is a concern, each instance of
95  * this class should be used with only on random engine
96  * instance, or two consecutive number should be generated with
97  * each engine used. This is due to the implicit caching of one
98  * Gaussian random number.
99  */
100  template<typename Acc, typename T>
102  {
103  static_assert(std::is_floating_point_v<T>, "Return type of NormalReal must be floating point.");
104 
105  Acc const* m_acc;
106  T m_cache = std::numeric_limits<T>::quiet_NaN();
107 
108  public:
109  /*! \warning Retains a reference to \p acc, thus must not outlive it.
110  */
111  ALPAKA_FN_HOST_ACC constexpr NormalReal(Acc const& acc) : m_acc(&acc)
112  {
113  }
114 
115  // All copy operations (and thus also move since we don't declare those and they fall back to copy) do NOT
116  // copy m_cache. This way we can ensure that the following holds:
117  // NormalReal<Acc> a(acc), b(acc);
118  // Engine<Acc> e(acc);
119  // assert(a(e) != b(e)); // because of two engine invocations
120  // b = a;
121  // assert(a(e) != b(e)); // because of two engine invocations
122 
123  ALPAKA_FN_HOST_ACC constexpr NormalReal(NormalReal const& other) : m_acc(other.m_acc)
124  {
125  }
126 
127  ALPAKA_FN_HOST_ACC constexpr auto operator=(NormalReal const& other) -> NormalReal&
128  {
129  m_acc = other.m_acc;
130  return *this;
131  }
132 
133  template<typename TEngine>
134  ALPAKA_FN_HOST_ACC auto operator()(TEngine& engine) -> T
135  {
136  constexpr auto sigma = T{1};
137  constexpr auto mu = T{0};
138  if(math::isnan(*m_acc, m_cache))
139  {
140  UniformReal<T> uni;
141 
142  T u1, u2;
143  do
144  {
145  u1 = uni(engine);
146  u2 = uni(engine);
147  } while(u1 <= std::numeric_limits<T>::epsilon());
148 
149  // compute z0 and z1
150  T const mag = sigma * math::sqrt(*m_acc, static_cast<T>(-2.) * math::log(*m_acc, u1));
151  constexpr T twoPi = static_cast<T>(2. * math::constants::pi);
152  // getting two normal number out of this, store one for later
153  m_cache = mag * static_cast<T>(math::cos(*m_acc, twoPi * u2)) + mu;
154 
155  return mag * static_cast<T>(math::sin(*m_acc, twoPi * u2)) + mu;
156  }
157 
158  T const ret = m_cache;
159  m_cache = std::numeric_limits<T>::quiet_NaN();
160  return ret;
161  }
162  };
163  } // namespace distribution::gpu
164 
165  namespace distribution::trait
166  {
167  //! The GPU device random number float normal distribution get trait specialization.
168  template<typename T>
169  struct CreateNormalReal<RandDefault, T, std::enable_if_t<std::is_floating_point_v<T>>>
170  {
171  template<typename TAcc>
173  {
174  return {acc};
175  }
176  };
177 
178  //! The GPU device random number float uniform distribution get trait specialization.
179  template<typename T>
180  struct CreateUniformReal<RandDefault, T, std::enable_if_t<std::is_floating_point_v<T>>>
181  {
183  {
184  return {};
185  }
186  };
187 
188  //! The GPU device random number integer uniform distribution get trait specialization.
189  template<typename T>
190  struct CreateUniformUint<RandDefault, T, std::enable_if_t<std::is_integral_v<T>>>
191  {
193  {
194  return {};
195  }
196  };
197  } // namespace distribution::trait
198 
199  namespace engine::trait
200  {
201  //! The GPU device random number default generator get trait specialization.
202  template<>
204  {
205  template<typename TAcc>
207  TAcc const& /* acc */,
208  std::uint32_t const& seed,
209  std::uint32_t const& subsequence,
210  std::uint32_t const& offset) -> Philox4x32x10
211  {
212  return {seed, subsequence, offset};
213  }
214  };
215  } // namespace engine::trait
216 } // namespace alpaka::rand
ALPAKA_FN_HOST_ACC auto operator()(TEngine &engine) -> T
constexpr ALPAKA_FN_HOST_ACC NormalReal(NormalReal const &other)
constexpr ALPAKA_FN_HOST_ACC auto operator=(NormalReal const &other) -> NormalReal &
constexpr ALPAKA_FN_HOST_ACC NormalReal(Acc const &acc)
The GPU random number uniform distribution.
Definition: RandDefault.hpp:70
ALPAKA_FN_HOST_ACC auto operator()(TEngine &engine) -> T
Definition: RandDefault.hpp:79
The GPU random number normal distribution.
Definition: RandDefault.hpp:45
ALPAKA_FN_HOST_ACC auto operator()(TEngine &engine) -> T
Definition: RandDefault.hpp:52
#define ALPAKA_FN_HOST_ACC
Definition: Common.hpp:39
constexpr double pi
Definition: Traits.hpp:61
ALPAKA_NO_HOST_ACC_WARNING ALPAKA_FN_HOST_ACC auto sqrt(T const &sqrt_ctx, TArg const &arg)
Computes the square root of arg.
Definition: Traits.hpp:1441
ALPAKA_NO_HOST_ACC_WARNING ALPAKA_FN_HOST_ACC auto cos(T const &cos_ctx, TArg const &arg)
Computes the cosine (measured in radians).
Definition: Traits.hpp:1060
ALPAKA_NO_HOST_ACC_WARNING ALPAKA_FN_HOST_ACC auto isnan(T const &ctx, TArg const &arg)
Checks if given value is NaN.
Definition: Traits.hpp:1192
ALPAKA_NO_HOST_ACC_WARNING ALPAKA_FN_HOST_ACC auto sin(T const &sin_ctx, TArg const &arg)
Computes the sine (measured in radians).
Definition: Traits.hpp:1393
ALPAKA_NO_HOST_ACC_WARNING ALPAKA_FN_HOST_ACC auto log(T const &log_ctx, TArg const &arg)
Computes the the natural (base e) logarithm of arg.
Definition: Traits.hpp:1210
constexpr auto offset
Definition: Extent.hpp:34
typename trait::AccType< T >::type Acc
The accelerator type trait alias template to remove the ::type.
Definition: Traits.hpp:58
Tag used in class inheritance hierarchies that describes that a specific concept (TConcept) is implem...
Definition: Concepts.hpp:15
static ALPAKA_FN_HOST_ACC auto createNormalReal(TAcc const &acc) -> gpu::NormalReal< TAcc, T >
The random number float normal distribution get trait.
Definition: Traits.hpp:27
static ALPAKA_FN_HOST_ACC auto createUniformReal(RandDefault const &) -> gpu::UniformReal< T >
The random number float uniform distribution get trait.
Definition: Traits.hpp:31
static ALPAKA_FN_HOST_ACC auto createUniformUint(RandDefault const &) -> gpu::UniformUint< T >
The random number integer uniform distribution get trait.
Definition: Traits.hpp:35
static ALPAKA_FN_HOST_ACC auto createDefault(TAcc const &, std::uint32_t const &seed, std::uint32_t const &subsequence, std::uint32_t const &offset) -> Philox4x32x10
The random number default generator engine get trait.
Definition: Traits.hpp:82