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