HILA
Loading...
Searching...
No Matches
random.cpp
1#include <cmath>
2#include "hila.h"
3#include "plumbing/random.h"
4
5/////////////////////////////////////////////////////////////////////////
6
7#include <random>
8
9// #ifndef OPENMP
10
11// static variable which holds the random state
12// Use 64-bit mersenne twister
13static std::mt19937_64 mersenne_twister_gen;
14
15// random numbers are in interval [0,1)
16static std::uniform_real_distribution<double> real_rnd_dist(0.0, 1.0);
17
18// #endif
19
20
21// In GPU code hila::random() defined in hila_gpu.cpp
22#if !defined(CUDA) && !defined(HIP)
23double hila::random() {
24 return real_rnd_dist(mersenne_twister_gen);
25}
26
27#endif
28
29// Generate random number in non-kernel (non-loop) code. Not meant to
30// be used in "user code"
31double hila::host_random() {
32 return real_rnd_dist(mersenne_twister_gen);
33}
34
35/////////////////////////////////////////////////////////////////////////
36
37static bool rng_is_initialized = false;
38
39namespace hila {
40
41/**
42 * @brief Random shuffling of rng seed for MPI nodes.
43 * @details Do it in a manner makes it difficult to give the same seed by mistake
44 * and also avoids giving the same seed for 2 nodes
45 * For single MPI node seed remains unchanged
46 */
47uint64_t shuffle_rng_seed(uint64_t seed) {
48
49 uint64_t n = hila::myrank();
50 if (hila::partitions.number() > 1)
51 n += hila::partitions.mylattice() * hila::number_of_nodes();
52
53 return (seed + n) ^ (n << 31);
54}
55
56
57/**
58 *@param seed unsigned int_64
59 *@note On MPI, this function shuffles different seed values for all MPI ranks.
60 */
61void initialize_host_rng(uint64_t seed) {
62
63 seed = hila::shuffle_rng_seed(seed);
64
65 // #if !defined(OPENMP)
66 mersenne_twister_gen.seed(seed);
67 // warm it up
68 for (int i = 0; i < 9000; i++)
69 mersenne_twister_gen();
70 // #endif
71}
72
73} // namespace hila
74
75
76/**
77 *@details The optional 2nd argument indicates whether to initialize the RNG on GPU device:
78 * `hila::device_rng_on` (default) or `hila::device_rng_off`. This argument does nothing if no GPU
79 * platform. If `hila::device_rng_off` is used, `onsites()` -loops cannot contain random number calls
80 * (Runtime error will be flagged and program exits).
81 *
82 * Seed is shuffled so that different nodes
83 * get different rng seeds. If `seed == 0`,
84 * seed is generated through using the `time()` -function.
85 */
86void hila::seed_random(uint64_t seed, bool device_init) {
87
88 rng_is_initialized = true;
89
90 if (!lattice.is_initialized()) {
91 hila::error("lattice.setup() must be called before hila::seed_random()");
92 }
93
94
95 uint64_t n = hila::myrank();
96 if (hila::partitions.number() > 1)
97 n += hila::partitions.mylattice() * hila::number_of_nodes();
98
99 if (seed == 0) {
100 // get seed from time
101 if (hila::myrank() == 0) {
102 struct timespec tp;
103
104 clock_gettime(CLOCK_MONOTONIC, &tp);
105 seed = tp.tv_sec;
106 seed = (seed << 30) ^ tp.tv_nsec;
107 hila::out0 << "Random seed from time: " << seed << '\n';
108 }
109 hila::broadcast(seed);
110 }
111
112 if (hila::partitions.number() > 1)
113 seed = seed ^ ((static_cast<uint64_t>(hila::partitions.mylattice())) << 28);
114
115#ifndef SITERAND
116
117 hila::out0 << "Using node random numbers, seed for node 0: " << seed << std::endl;
118
120
121#if defined(CUDA) || defined(HIP)
122
123 // we can use the same seed, the generator is different
124 if (device_init) {
126 } else {
127 hila::out0 << "Not initializing GPU random numbers\n";
128 }
129
130#endif
131
132 // taus_initialize();
133
134#else
135
136 // TODO: SITERAND is not yet implemented!
137 // This is usually used only for occasional benchmarking, where identical output
138 // independent of the node number is desired
139
140 hila::out0 << "*** SITERAND is in use!\n";
141
142 random_seed_arr = (unsigned short(*)[3])memalloc(3 * node.sites * sizeof(unsigned short));
143 forallsites(i) {
144 random_seed_arr[i][0] = (unsigned short)(seed + site[i].index);
145 random_seed_arr[i][1] = (unsigned short)(seed + 2 * site[i].index);
146 random_seed_arr[i][2] = (unsigned short)(seed + 3 * site[i].index);
147 }
148
149 random_seed_ptr = random_seed_arr[0];
150
151#endif
152}
153
154////////////////////////////////////////////////////////////////////
155// Def here gpu rng functions for non-gpu
156////////////////////////////////////////////////////////////////////
157
158#if !(defined(CUDA) || defined(HIP))
159
160/**
161 *@details `hila::random()` does not work inside `onsites()` after this,
162 *unless seeded again using `initialize_device_rng()`. Frees the memory RNG takes on the device.
163 */
164void hila::free_device_rng() {}
165
166/**
167 *@details Returns `true` on non-GPU archs.
168 */
169
170/**
171 *@details Returns `true` on non-GPU archs.
172 */
174 return true;
175}
176
177/**
178 *@details This function shuffles the seed for different MPI ranks on MPI. Called by `seed_random()` unless its 2nd
179 *argument is `hila::device_rng_off`. This can reinitialize device RNG free'd by `free_device_rng()`.
180 */
181void hila::initialize_device_rng(uint64_t seed) {}
182
183#endif
184
185
186#define VARIANCE 1.0
187
188double hila::gaussrand2(double &out2) {
189
190 double phi, urnd, r;
191 phi = 2.0 * M_PI * hila::random();
192
193 // this should not really trigger
194 do {
195 urnd = 1.0 - hila::random();
196 } while (urnd == 0.0);
197
198 r = sqrt(-::log(urnd) * (2.0 * VARIANCE));
199 out2 = r * cos(phi);
200 return r * sin(phi);
201}
202
203#if !defined(CUDA) && !defined(HIP)
204
205/**
206 * @details By default these gives random numbers with variance \f$1.0\f$ and expectation value \f$0.0\f$, i.e.
207 * \f[
208 * e^{-(\frac{x^{2}}{2})}
209 * \f]
210 * with variance
211 * \f[
212 * < x^{2}-0.0> = 1
213 * \f]
214 *
215 * If you want random numbers with variance \f$ \sigma^{2} \f$, multiply the
216 * result by \f$ \sqrt{\sigma^{2}} \f$ i.e.,
217 * \code {.cpp}
218 * sqrt(sigma * sigma) * gaussrand();
219 * \endcode
220 *
221 * @return double
222 */
223double hila::gaussrand() {
224 static double second;
225 static bool draw_new = true;
226 if (draw_new) {
227 draw_new = false;
228 return hila::gaussrand2(second);
229 }
230 draw_new = true;
231 return second;
232}
233
234#else
235
236// Cuda and other stuff which does not accept
237// static variables - just throw away another gaussian number.
238
239// #pragma hila loop function contains rng
241 double second;
242 return hila::gaussrand2(second);
243}
244
245#endif
246
247/**
248 *@returns bool
249 */
251 return rng_is_initialized;
252}
253
254
255///////////////////////////////////////////////////////////////
256// RNG initialization check - emitted on loops
257///////////////////////////////////////////////////////////////
258
259/**
260 *@details program quits with error message if RNG is not initialized.
261 *It also quit with error messages if the device RNG is not initialized.
262 */
264 bool isinit;
265
266 if (!rng_is_initialized) {
267 hila::out0 << "ERROR: trying to use random numbers without initializing the generator"
268 << std::endl;
270 }
271#if defined(CUDA) || defined(HIP)
272 if (!hila::is_device_rng_on()) {
273 hila::out0 << "ERROR: GPU random number generator is not initialized and onsites()-loop is "
274 "using random numbers"
275 << std::endl;
277 }
278#endif
279}
280
281
282
283
284
285
Implement hila::swap for gauge fields.
Definition array.h:982
void initialize_device_rng(uint64_t seed)
Initialize device random number generator on GPUs, if application run on GPU platform....
Definition hila_gpu.cpp:71
void free_device_rng()
Free GPU RNG state, does nothing on non-GPU archs.
Definition hila_gpu.cpp:107
double random()
Real valued uniform random number generator.
Definition hila_gpu.cpp:121
int myrank()
rank of this node
Definition com_mpi.cpp:237
int number_of_nodes()
how many nodes there are
Definition com_mpi.cpp:248
bool is_device_rng_on()
Check if the RNG on GPU is allocated and ready to use.
Definition hila_gpu.cpp:58
uint64_t shuffle_rng_seed(uint64_t seed)
Random shuffling of rng seed for MPI nodes.
Definition random.cpp:47
double gaussrand()
Gaussian random generation routine.
Definition random.cpp:240
std::ostream out0
This writes output only from main process (node 0)
void seed_random(uint64_t seed, bool device_rng=true)
Seed random generators with 64-bit unsigned value. On MPI shuffles the seed so that different MPI ran...
Definition random.cpp:86
T broadcast(T &var, int rank=0)
Broadcast the value of var to all MPI ranks from rank (default=0).
Definition com_mpi.h:170
bool is_rng_seeded()
Check if RNG is seeded already.
Definition random.cpp:250
void initialize_host_rng(uint64_t seed)
Initialize host (CPU) random number generator separately, done implicitly by seed_random()
Definition random.cpp:61
double gaussrand2(double &out2)
hila::gaussrand2 returns 2 Gaussian distributed random numbers with variance .
Definition random.cpp:188
void check_that_rng_is_initialized()
Check if RNG is initialized, do what the name says.
Definition random.cpp:263
void terminate(int status)