Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4UniformRandPool.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id:$
28 
29 // ------------------------------------------------------------
30 
31 #include "G4UniformRandPool.hh"
32 #include "globals.hh"
33 
34 #include <climits>
35 #include <stdlib.h>
36 #include <algorithm>
37 #include <cstring>
38 
39 // Not aligned memory
40 //
42 {
43  buffer = new G4double[ps];
44 }
45 
47 {
48  delete[] buffer;
49 }
50 
51 
52 #if defined(WIN32)
53 // No bother with WIN
55 {
56  create_pool(buffer,ps);
57 }
59 {
60  destroy_pool(buffer);
61 }
62 
63 #else
64 
65 // Align memory pools
66 // Assumption is: static_assert(sizeof(G4double)*CHAR_BIT==64)
67 //
69 {
70 //#if defined(__ICC) || (__INTEL_COMPILER)
71 // //INTEL optimized way
72 // buffer = (G4doulbe*)mm_allign(ps*sizeof(G4double),sizeof(G4double)*CHAR_BIT);
73 //#else
74 
75  // POSIX standard way
76  G4int errcode = posix_memalign( (void**) &buffer ,
77  sizeof(G4double)*CHAR_BIT,
78  ps*sizeof(G4double));
79  if ( errcode != 0 )
80  {
81  G4Exception("G4UniformRandPool::create_pool_align()",
82  "InvalidCondition", FatalException,
83  "Cannot allocate aligned buffer");
84  return;
85  }
86 //#endif
87 
88  return;
89 }
90 
92 {
93  //#if defined(__ICC) || (__INTEL_COMPILER)
94  //mm_free(buffer);
95  //#else
96 
97  free(buffer);
98 
99  //#endif
100 }
101 #endif
102 
104  : size(G4UNIFORMRANDPOOL_DEFAULT_POOLSIZE), buffer(0), currentIdx(0)
105 {
106  if ( sizeof(G4double)*CHAR_BIT==64 )
107  {
108  create_pool_align(buffer,size);
109  }
110  else
111  {
112  create_pool(buffer,size);
113  }
114  Fill(size);
115 }
116 
118  : size(siz), buffer(0), currentIdx(0)
119 {
120  if ( sizeof(G4double)*CHAR_BIT==64 )
121  {
122  create_pool_align(buffer,size);
123  }
124  else
125  {
126  create_pool(buffer,size);
127  }
128  Fill(size);
129 
130 }
131 
133 {
134  if ( sizeof(G4double)*CHAR_BIT==64 )
135  {
136  destroy_pool_align(buffer);
137  }
138  else
139  {
140  destroy_pool(buffer);
141  }
142 }
143 
144 void G4UniformRandPool::Resize(/*PoolSize_t*/ G4int newSize )
145 {
146  if ( newSize != size )
147  {
148  destroy_pool(buffer);
149  create_pool(buffer,newSize);
150  size=newSize;
151  currentIdx = 0;
152  }
153  currentIdx = 0;
154 }
155 
157 {
158  assert(rnds!=0 && howmany>0);
159 
160  // if ( howmany <= 0 ) return;
161  // We generate at max "size" numbers at once, and
162  // We do not want to use recursive calls (expensive).
163  // We need to deal with the case howmany>size
164  // So:
165  // how many times I need to get "size" numbers?
166 
167  const G4int maxcycles = howmany/size;
168 
169  // This is the rest
170  //
171  const G4int peel = howmany%size;
172  assert(peel<size);
173 
174  // Ok from now I will get random numbers in group of "size"
175  // Note that if howmany<size maxcycles == 0
176  //
177  G4int cycle = 0;
178 
179  // Consider the case howmany>size, then maxcycles>=1
180  // and we will request at least "size" rng, so
181  // let's start with a fresh buffer of numbers if needed
182  //
183  if ( maxcycles>0 && currentIdx>0 )
184  {
185  assert(currentIdx<=size);
186  Fill(currentIdx);//<size?currentIdx:size);
187  }
188  for ( ; cycle < maxcycles ; ++cycle )
189  {
190  // We can use memcpy of std::copy, it turns out that the two are basically
191  // performance-wise equivalent (expected), since in my tests memcpy is a
192  // little bit faster, I use that
193  // std::copy(buffer,buffer+size,rnds+(cycle*size));
194  //
195  memcpy(rnds+(cycle*size),buffer,sizeof(G4double)*size );
196 
197  // Get a new set of numbers
198  //
199  Fill(size); // Now currentIdx is 0 again
200  }
201 
202  // If maxcycles>0 last think we did was to call Fill(size)
203  // so currentIdx == 0
204  // and it is guranteed that peel<size, we have enough fresh random numbers
205  // but if maxcycles==0 currentIdx can be whatever, let's make sure we have
206  // enough fresh numbers
207  //
208  if (currentIdx + peel >= size)
209  {
210  Fill(currentIdx<size?currentIdx:size);
211  }
212  memcpy(rnds+(cycle*size) , buffer+currentIdx , sizeof(G4double)*peel );
213  // std::copy(buffer+currentIdx,buffer+(currentIdx+peel), rnds+(cycle*size));
214 
215  // Advance index, we are done
216  //
217  currentIdx+=peel;
218  assert(currentIdx<=size);
219 }
220 
221 // Static interfaces implementing CLHEP methods
222 
223 #include "G4Threading.hh"
224 #include "G4AutoDelete.hh"
225 
226 namespace
227 {
228  G4ThreadLocal G4UniformRandPool* rndpool = 0;
229 }
230 
232 {
233  if ( rndpool == 0 )
234  {
235  rndpool = new G4UniformRandPool;
236  G4AutoDelete::Register(rndpool);
237  }
238  return rndpool->GetOne();
239 }
240 
242 {
243  if ( rndpool == 0 )
244  {
245  rndpool = new G4UniformRandPool;
246  G4AutoDelete::Register(rndpool);
247  }
248  rndpool->GetMany(rnds,(unsigned int)howmany);
249 }
void destroy_pool(G4double *&buffer)
#define G4UNIFORMRANDPOOL_DEFAULT_POOLSIZE
void create_pool_align(G4double *&buffer, G4int ps)
void destroy_pool_align(G4double *&buffer)
#define buffer
Definition: xmlparse.cc:628
static G4double flat()
#define G4ThreadLocal
Definition: tls.hh:89
int G4int
Definition: G4Types.hh:78
void Register(T *inst)
Definition: G4AutoDelete.hh:65
static void flatArray(G4int howmany, G4double *rnds)
void GetMany(G4double *rnds, G4int howMany)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
void create_pool(G4double *&buffer, G4int ps)
void Resize(G4int newSize)
double G4double
Definition: G4Types.hh:76
static constexpr double ps
Definition: G4SIunits.hh:172