|
| 1 | +// |
| 2 | +// AMReX Interface to Mersenne Twistor |
| 3 | +// |
| 4 | + |
| 5 | +/* A C-program for MT19937: Real number version (1999/10/28) */ |
| 6 | +/* genrand() generates one pseudorandom real number (double) */ |
| 7 | +/* which is uniformly distributed on [0,1]-interval, for each */ |
| 8 | +/* call. sgenrand(seed) sets initial values to the working area */ |
| 9 | +/* of 624 words. Before genrand(), sgenrand(seed) must be */ |
| 10 | +/* called once. (seed is any 32-bit integer.) */ |
| 11 | +/* Integer generator is obtained by modifying two lines. */ |
| 12 | +/* Coded by Takuji Nishimura, considering the suggestions by */ |
| 13 | +/* Topher Cooper and Marc Rieffel in July-Aug. 1997. */ |
| 14 | + |
| 15 | +/* This library is free software under the Artistic license: */ |
| 16 | +/* see the file COPYING distributed together with this code. */ |
| 17 | +/* For the verification of the code, its output sequence file */ |
| 18 | +/* mt19937-1.out is attached (2001/4/2) */ |
| 19 | + |
| 20 | +/* Copyright (C) 1997, 1999 Makoto Matsumoto and Takuji Nishimura. */ |
| 21 | +/* Any feedback is very welcome. For any question, comments, */ |
| 22 | +/* see http://www.math.keio.ac.jp/matumoto/emt.html or email */ |
| 23 | +/* matumoto@math.keio.ac.jp */ |
| 24 | + |
| 25 | +/* REFERENCE */ |
| 26 | +/* M. Matsumoto and T. Nishimura, */ |
| 27 | +/* "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform */ |
| 28 | +/* Pseudo-Random Number Generator", */ |
| 29 | +/* ACM Transactions on Modeling and Computer Simulation, */ |
| 30 | +/* Vol. 8, No. 1, January 1998, pp 3--30. */ |
| 31 | + |
| 32 | +#ifndef MERSENNETWISTER_H_ |
| 33 | +#define MERSENNETWISTER_H_ |
| 34 | + |
| 35 | +#include <AMReX_Vector.H> |
| 36 | + |
| 37 | +namespace pele::physics::turbforcing::mersenne_twister { |
| 38 | +class mt19937 |
| 39 | +{ |
| 40 | +public: |
| 41 | + using seed_type = unsigned long; |
| 42 | + explicit mt19937(seed_type seed = 4357UL); |
| 43 | + mt19937(seed_type array[], int array_len); |
| 44 | + void rewind(); |
| 45 | + void reset(unsigned long seed); |
| 46 | + |
| 47 | + double d_value(); // [0,1] random numbers |
| 48 | + double d1_value(); // [0,1) random numbers |
| 49 | + double d2_value(); // (0,1) random numbers |
| 50 | + |
| 51 | + long l_value(); // [0,2^31-1] random numbers |
| 52 | + unsigned long u_value(); // [0,2^32-1] random numbers |
| 53 | + |
| 54 | + void save(amrex::Vector<unsigned long>& state) const; |
| 55 | + int RNGstatesize() const; |
| 56 | + void restore(const amrex::Vector<unsigned long>& state); |
| 57 | + |
| 58 | +private: |
| 59 | + enum { N = 624 }; |
| 60 | + static unsigned long init_seed; |
| 61 | + static unsigned long mt[N]; // the array for the state vector |
| 62 | + static int mti; // mti==N+1 means mt[N] is not initialized |
| 63 | +#ifdef _OPENMP |
| 64 | +#pragma omp threadprivate(init_seed, mt, mti) |
| 65 | +#endif |
| 66 | + void sgenrand(unsigned long seed); |
| 67 | + unsigned long igenrand(); |
| 68 | + void reload(); |
| 69 | +}; |
| 70 | + |
| 71 | +/* |
| 72 | + Generates one pseudorandom real number (double) which is |
| 73 | + uniformly distributed on [0,1]-interval for each call. |
| 74 | +
|
| 75 | + Accepts any 32-bit integer as a seed -- uses 4357 as the default. |
| 76 | +
|
| 77 | + Has a period of 2**19937. |
| 78 | +
|
| 79 | + Mersenne Twister Home Page: http://www.math.keio.ac.jp/matumoto/emt.html |
| 80 | +
|
| 81 | + There is also an entry point for Fortran callable as: |
| 82 | +
|
| 83 | + REAL_T rn |
| 84 | + call blutilrand(rn) |
| 85 | +
|
| 86 | + Internally, blutilrand() calls a static Mersenne Twister object (the |
| 87 | + same one used by AMReX::Random()) to get a value in [0,1] and then |
| 88 | + sets "rn" to that value. |
| 89 | +*/ |
| 90 | +double Random(); // [0,1] |
| 91 | +double Random1(); // [0,1) |
| 92 | +double Random2(); // (0,1) |
| 93 | +unsigned long Random_int(unsigned long n); // [0,n-1], where n<=2^32-1 |
| 94 | +/* Set the seed of the random number generator. |
| 95 | +
|
| 96 | + There is also an entry point for Fortran callable as: |
| 97 | +
|
| 98 | + INTEGER seed |
| 99 | + call blutilinitrand(seed) |
| 100 | +
|
| 101 | + or |
| 102 | +
|
| 103 | + INTEGER seed |
| 104 | + call blinitrand(seed) |
| 105 | +*/ |
| 106 | +void InitRandom(unsigned long seed); |
| 107 | + |
| 108 | +void ResetRandomSeed(unsigned long seed); |
| 109 | +// |
| 110 | +// Save and restore random state. |
| 111 | +// |
| 112 | +// state.size() == 626 on return from Save & on entry to Restore. |
| 113 | +// |
| 114 | +void SaveRandomState(amrex::Vector<unsigned long>& state); |
| 115 | + |
| 116 | +int sizeofRandomState(); |
| 117 | + |
| 118 | +void RestoreRandomState(const amrex::Vector<unsigned long>& state); |
| 119 | +// |
| 120 | +// Create a unique subset of random numbers from a pool |
| 121 | +// of integers in the range [0, poolSize - 1] |
| 122 | +// the set will be in the order they are found |
| 123 | +// setSize must be <= poolSize |
| 124 | +// uSet will be resized to setSize |
| 125 | +// if you want all processors to have the same set, |
| 126 | +// call this on one processor and broadcast the array |
| 127 | +// |
| 128 | +void UniqueRandomSubset( |
| 129 | + amrex::Vector<int>& uSet, int setSize, int poolSize, bool printSet = false); |
| 130 | +} // namespace pele::physics::turbforcing::mersenne_twister |
| 131 | + |
| 132 | +#endif |
0 commit comments