| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 //
28 // G4UniformRandPool implementation
29 //
30 // Author: A.Dotti (SLAC)
31 // -------------------------------------------
32
33 #include "G4UniformRandPool.hh"
34
35 #include "G4AutoDelete.hh"
36 #include "G4Threading.hh"
37 #include "globals.hh"
38
39 #include <algorithm>
40 #include <climits>
41 #include <cstdlib>
42 #include <cstring>
43
44 // Not aligned memory
45 //
46 void create_pool(G4double*& buffer, G4int ps)
47
48 void destroy_pool(G4double*& buffer) { delete[
49
50 #if defined(WIN32) || defined(__MINGW32__)
51 // No bother with WIN
52 void create_pool_align(G4double*& buffer, G4in
53 void destroy_pool_align(G4double*& buffer) { d
54
55 #else
56
57 // Align memory pools
58 // Assumption is: static_assert(sizeof(G4doubl
59 //
60 void create_pool_align(G4double*& buffer, G4in
61 {
62 // POSIX standard way
63 G4int errcode = posix_memalign((void**) &buf
64 ps * sizeof(G
65 if(errcode != 0)
66 {
67 G4Exception("G4UniformRandPool::create_poo
68 FatalException, "Cannot alloca
69 return;
70 }
71 return;
72 }
73
74 void destroy_pool_align(G4double*& buffer) { f
75 #endif
76
77 G4UniformRandPool::G4UniformRandPool()
78 {
79 if(sizeof(G4double) * CHAR_BIT == 64)
80 {
81 create_pool_align(buffer, size);
82 }
83 else
84 {
85 create_pool(buffer, size);
86 }
87 Fill(size);
88 }
89
90 G4UniformRandPool::G4UniformRandPool(G4int siz
91 : size(siz)
92 {
93 if(sizeof(G4double) * CHAR_BIT == 64)
94 {
95 create_pool_align(buffer, size);
96 }
97 else
98 {
99 create_pool(buffer, size);
100 }
101 Fill(size);
102 }
103
104 G4UniformRandPool::~G4UniformRandPool()
105 {
106 if(sizeof(G4double) * CHAR_BIT == 64)
107 {
108 destroy_pool_align(buffer);
109 }
110 else
111 {
112 destroy_pool(buffer);
113 }
114 }
115
116 void G4UniformRandPool::Resize(/*PoolSize_t*/
117 {
118 if(newSize != size)
119 {
120 destroy_pool(buffer);
121 create_pool(buffer, newSize);
122 size = newSize;
123 currentIdx = 0;
124 }
125 currentIdx = 0;
126 }
127
128 void G4UniformRandPool::Fill(G4int howmany)
129 {
130 assert(howmany > 0 && howmany <= size);
131
132 // Fill buffer with random numbers
133 //
134 G4Random::getTheEngine()->flatArray(howmany,
135 currentIdx = 0;
136 }
137
138 void G4UniformRandPool::GetMany(G4double* rnds
139 {
140 assert(rnds != 0 && howmany > 0);
141
142 // if ( howmany <= 0 ) return;
143 // We generate at max "size" numbers at once
144 // We do not want to use recursive calls (ex
145 // We need to deal with the case howmany>si
146 // So:
147 // how many times I need to get "size" numbe
148
149 const G4int maxcycles = howmany / size;
150
151 // This is the rest
152 //
153 const G4int peel = howmany % size;
154 assert(peel < size);
155
156 // Ok from now on I will get random numbers
157 // Note that if howmany<size maxcycles == 0
158 //
159 G4int cycle = 0;
160
161 // Consider the case howmany>size, then maxc
162 // and we will request at least "size" rng,
163 // let's start with a fresh buffer of number
164 //
165 if(maxcycles > 0 && currentIdx > 0)
166 {
167 assert(currentIdx <= size);
168 Fill(currentIdx); //<size?currentIdx:size
169 }
170 for(; cycle < maxcycles; ++cycle)
171 {
172 // We can use memcpy of std::copy, it turn
173 // performance-wise equivalent (expected),
174 // little bit faster, I use that
175 //
176 memcpy(rnds + (cycle * size), buffer, size
177 // std::copy(buffer,buffer+size,rnds+(cycl
178
179 // Get a new set of numbers
180 //
181 Fill(size); // Now currentIdx is 0 again
182 }
183
184 // If maxcycles>0 last think we did was to c
185 // so currentIdx == 0
186 // and it is guaranteed that peel<size, we h
187 // but if maxcycles==0 currentIdx can be wha
188 // enough fresh numbers
189 //
190 if(currentIdx + peel >= size)
191 {
192 Fill(currentIdx < size ? currentIdx : size
193 }
194 memcpy(rnds + (cycle * size), buffer + curre
195 // std::copy(buffer+currentIdx,buffer+(curre
196
197 // Advance index, we are done
198 //
199 currentIdx += peel;
200 assert(currentIdx <= size);
201 }
202
203 // Static interfaces implementing CLHEP method
204
205 namespace
206 {
207 G4ThreadLocal G4UniformRandPool* rndpool = n
208 }
209
210 G4double G4UniformRandPool::flat()
211 {
212 if(rndpool == nullptr)
213 {
214 rndpool = new G4UniformRandPool;
215 G4AutoDelete::Register(rndpool);
216 }
217 return rndpool->GetOne();
218 }
219
220 void G4UniformRandPool::flatArray(G4int howman
221 {
222 if(rndpool == nullptr)
223 {
224 rndpool = new G4UniformRandPool;
225 G4AutoDelete::Register(rndpool);
226 }
227 rndpool->GetMany(rnds, (unsigned int) howman
228 }
229