CASToR  2.0
Tomographic Reconstruction (PET/SPECT/CT)
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
sRandomNumberGenerator.cc
Go to the documentation of this file.
1 /*
2 This file is part of CASToR.
3 
4  CASToR is free software: you can redistribute it and/or modify it under the
5  terms of the GNU General Public License as published by the Free Software
6  Foundation, either version 3 of the License, or (at your option) any later
7  version.
8 
9  CASToR is distributed in the hope that it will be useful, but WITHOUT ANY
10  WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12  details.
13 
14  You should have received a copy of the GNU General Public License along with
15  CASToR (in file GNU_GPL.TXT). If not, see <http://www.gnu.org/licenses/>.
16 
17 Copyright 2017-2018 all CASToR contributors listed below:
18 
19  --> current contributors: Thibaut MERLIN, Simon STUTE, Didier BENOIT, Claude COMTAT, Marina FILIPOVIC, Mael MILLARDET
20  --> past contributors: Valentin VIELZEUF
21 
22 This is CASToR version 2.0.
23 */
24 
34 
36 
37 // =====================================================================
38 // ---------------------------------------------------------------------
39 // ---------------------------------------------------------------------
40 // =====================================================================
45 {
46  mp_Instance = NULL;
47  m_verbose = 5;
48  mp_Engines.clear();
49  mp_extraEngines.clear();
50  m_seed = -1;
51 }
52 
53 
54 
55 // =====================================================================
56 // ---------------------------------------------------------------------
57 // ---------------------------------------------------------------------
58 // =====================================================================
59 /*
60  \brief Destructor of sRandomNumberGenerator. Do nothing by default
61 */
63 {
64 }
65 
66 
67 
68 // =====================================================================
69 // ---------------------------------------------------------------------
70 // ---------------------------------------------------------------------
71 // =====================================================================
72 
73 int sRandomNumberGenerator::Initialize(int a_nbThreads, int a_nbExtra)
74 {
75  if(m_verbose >=3) Cout("sRandomNumberGenerator::Initialize ..."<< endl);
76 
77  int64_t seed = -1;
78 
79  // call a real random generator, only once because this operating system
80  // random source might be limited in size and time
81  std::random_device rd;
82 
83  // save the initial seed
84  m_seed = rd();
85 
86  #ifdef CASTOR_MPI
87  int mpi_rank_temp = 0;
88  int mpi_size_temp = 0;
89  MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank_temp);
90  MPI_Comm_size(MPI_COMM_WORLD, &mpi_size_temp);
91 
92  // helper variable
93  int64_t temp_seed = -1;
94 
95  // if more than one MPI instance, generate seeds for each instance
96  // and dispatch them, otherwise keep the initial seed
97  if (mpi_size_temp>1)
98  {
99  if (mpi_rank_temp==0)
100  {
101  Engine mpi_generator(m_seed);
102  for (int p=0; p<mpi_size_temp; p++)
103  {
104  m_seed = mpi_generator();
105 
106  Cout("sRandomNumberGenerator::Seed for rank " << p << " is " << m_seed << endl);
107 
108  if (p==0)
109  temp_seed = m_seed;
110  else
111  MPI_Send(&m_seed, 1, MPI_LONG, p, 0, MPI_COMM_WORLD);
112  }
113  m_seed = temp_seed;
114  }
115  else
116  {
117  MPI_Status* status = NULL;
118  MPI_Recv(&m_seed, 1, MPI_LONG, 0, 0, MPI_COMM_WORLD, status);
119  }
120  // wait for all the processes to have their seeds
121  MPI_Barrier(MPI_COMM_WORLD);
122  }
123  #endif
124 
125  // seed the initial temporary PRNG, which will generate seeds for actual PRNGs (threaded and non-threaded)
126  Engine initialGenerator(m_seed);
127 
128  for (int th=0; th<a_nbThreads; th++)
129  {
130  seed = initialGenerator();
131  if (m_verbose>=3) Cout("sRandomNumberGenerator::Initialize() -> Seed for thread " << th << " : " << seed << endl);
132  mp_Engines.push_back(Engine(seed));
133  }
134 
135  for(int ex=0; ex<a_nbExtra; ex++)
136  {
137  seed = initialGenerator();
138  if(m_verbose>=3) Cout("sRandomNumberGenerator::Initialize() - >Seed for additional nonthreaded generator "<<ex<<" : " << seed << endl);
139  mp_extraEngines.push_back(Engine(seed));
140  }
141 
142  return 0;
143 }
144 
145 
146 
147 // =====================================================================
148 // ---------------------------------------------------------------------
149 // ---------------------------------------------------------------------
150 // =====================================================================
151 
152 int sRandomNumberGenerator::Initialize(int64_t a_seed, int a_nbThreads, int a_nbExtra)
153 {
154  if(m_verbose >=3) Cout("sRandomNumberGenerator::Initialize with provided seed "<<a_seed<< endl);
155 
156  // temporary variable
157  int64_t seed = -1;
158 
159  if (a_seed<0)
160  {
161  Cout("***** sRandomNumberGenerator::Initialize()-> Error : seed for RNG should be >=0 !" << endl);
162  return 1;
163  }
164 
165  // save the initial seed
166  m_seed = a_seed;
167 
168  #ifdef CASTOR_MPI
169  int mpi_rank_temp = 0;
170  int mpi_size_temp = 0;
171  MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank_temp);
172  MPI_Comm_size(MPI_COMM_WORLD, &mpi_size_temp);
173 
174  // helper variable
175  int64_t temp_seed = -1;
176 
177  // if more than one MPI instance, generate seeds for each instance
178  // and dispatch them, otherwise keep the initial seed
179  if (mpi_size_temp>1)
180  {
181  if (mpi_rank_temp==0)
182  {
183  Engine mpi_generator(a_seed);
184  for (int p=0; p<mpi_size_temp; p++)
185  {
186  m_seed = mpi_generator();
187 
188  Cout("sRandomNumberGenerator::Seed for rank "<<p<<" is "<<m_seed<<""<< endl);
189 
190  if (p==0)
191  temp_seed = m_seed;
192  else
193  MPI_Send(&m_seed, 1, MPI_LONG, p, 0, MPI_COMM_WORLD);
194  }
195  m_seed = temp_seed;
196  }
197  else
198  {
199  MPI_Status* status = NULL;
200  MPI_Recv(&m_seed, 1, MPI_LONG, 0, 0, MPI_COMM_WORLD, status);
201  }
202  MPI_Barrier(MPI_COMM_WORLD);
203  }
204  #endif
205 
206  // seed the initial temporary PRNG, which will generate seeds for actual PRNGs (threaded and non-threaded)
207  Engine initialGenerator(m_seed);
208 
209  for(int th=0; th<a_nbThreads; th++)
210  {
211  seed = initialGenerator();
212  mp_Engines.push_back(Engine(seed));
213  if (m_verbose>=3) Cout("sRandomNumberGenerator::Initialize() -> Seed for thread " << th << " : " << seed << endl);
214  }
215 
216  for(int ex=0; ex<a_nbExtra; ex++)
217  {
218  seed = initialGenerator();
219  mp_extraEngines.push_back(Engine(seed));
220  if (m_verbose>=3) Cout("sRandomNumberGenerator::Initialize() -> Seed for additional nonthreaded random generator "<<ex<< " : " << seed << endl);
221  }
222 
223  return 0;
224 }
225 
226 
227 
228 // =====================================================================
229 // ---------------------------------------------------------------------
230 // ---------------------------------------------------------------------
231 // =====================================================================
232 
234 {
235  #ifdef CASTOR_VERBOSE
236  if (m_verbose >=4) Cout("sRandomNumberGenerator::GenerateRdmNber..." << endl);
237  #endif
238 
239  int id=0;
240  #ifdef CASTOR_OMP
241  id = omp_get_thread_num();
242  #endif
243  return mp_Distribution(mp_Engines[id]);
244 }
245 
246 
247 
248 // =====================================================================
249 // ---------------------------------------------------------------------
250 // ---------------------------------------------------------------------
251 // =====================================================================
252 
254 {
255  #ifdef CASTOR_VERBOSE
256  if(m_verbose >=4) Cout("sRandomNumberGenerator::GenerateExtraRdmNber..."<< endl);
257  #endif
258 
259  return mp_Distribution(mp_extraEngines[a_nb]);
260 }
261 
262 
263 // =====================================================================
264 // ---------------------------------------------------------------------
265 // ---------------------------------------------------------------------
266 // =====================================================================
267 
269 {
270  #ifdef CASTOR_VERBOSE
271  if(m_verbose >=4) Cout("sRandomNumberGenerator::GetExtraGenerator..."<< endl);
272  #endif
273 
274  return mp_extraEngines[a_nb];
275 }
#define HPFLTNB
Definition: gVariables.hh:83
Engine & GetExtraGenerator(int a_nb=0)
Get the not thread safe additional random generator, for use various sequential parts of an otherwise...
HPFLTNB GenerateExtraRdmNber(int a_nb=0)
Generate a random number using the specified additional not thread safe random generator, for use in sequential parts of an otherwise multithreaded code.
int Initialize(int a_nbThreads, int a_nbExtra)
Instantiate pseudo random number generators, one per thread by default, and additional extra ones if ...
sRandomNumberGenerator()
Constructor of sRandomNumberGenerator. Do nothing by default as it is a singleton clasee...
HPFLTNB GenerateRdmNber()
Generate a random number for the thread which index is recovered from the OpenMP function.
Singleton class that generate a thread-safe random generator number for openMP As singleton...
Declaration of class sRandomNumberGenerator.
#define Cout(MESSAGE)
static sRandomNumberGenerator * mp_Instance
~sRandomNumberGenerator()
Destructor of sRandomNumberGenerator. Do nothing by default.