/[escript]/trunk/esysUtils/src/EsysRandom.cpp
ViewVC logotype

Contents of /trunk/esysUtils/src/EsysRandom.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4696 - (show annotations)
Wed Feb 19 07:29:50 2014 UTC (5 years, 2 months ago) by jfenwick
File size: 6316 byte(s)
Correcting python libname for sage.
Fixing bug where the number of values in the shape was not considered for buffer and message size.

1 /*****************************************************************************
2 *
3 * Copyright (c) 2013-2014 by University of Queensland
4 * http://www.uq.edu.au
5 *
6 * Primary Business: Queensland, Australia
7 * Licensed under the Open Software License version 3.0
8 * http://www.opensource.org/licenses/osl-3.0.php
9 *
10 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 * Development 2012-2013 by School of Earth Sciences
12 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 *
14 *****************************************************************************/
15
16 #include <vector>
17 #include <boost/random/mersenne_twister.hpp>
18 #include <cstring>
19 #include "Esys_MPI.h"
20 #ifdef _OPENMP
21 #include <omp.h>
22 #endif
23 using namespace std;
24
25 namespace {
26
27 boost::mt19937 base; // used to seed all the other generators
28 vector<boost::mt19937*> gens;
29 vector<uint32_t> seeds;
30
31 void seedGens(long seed)
32 {
33 #ifdef _OPENMP
34 int numthreads=omp_get_max_threads();
35 #else
36 int numthreads=1;
37 #endif
38 if (gens.size()==0) // we haven't instantiated the generators yet
39 {
40 gens.resize(numthreads);
41 seeds.resize(numthreads);
42 }
43 if (seed!=0)
44 {
45 int i;
46 base.seed((uint32_t)seed); // without this cast, icc gets confused
47 for (int i=0;i<numthreads;++i)
48 {
49 uint32_t b=base();
50 seeds[i]=b; // initialise each generator with successive random values
51 }
52 #pragma omp parallel for private(i)
53 for (i=0;i<numthreads;++i)
54 {
55 gens[i]=new boost::mt19937(seeds[i]);
56 }
57 }
58 }
59
60
61 }
62
63 namespace esysUtils
64 {
65
66 // Put n random values into array
67 // Idea here is to create an array of seeds by feeding the original seed into the random generator
68 // The code at the beginning of the function to compute the seed if one is given is
69 // just supposed to introduce some variety (and ensure that multiple ranks don't get the same seed).
70 // I make no claim about how well these initial seeds are distributed
71 // uses openmp
72 // don't forget to call CHECK_FOR_EX_WRITE if using this on Data
73 void randomFillArray(long seed, double* array, size_t n)
74 {
75 static unsigned prevseed=0; // So if we create a bunch of objects we don't get the same start seed
76 if (seed==0) // for each one
77 {
78 if (prevseed==0)
79 {
80 time_t s=time(0);
81 seed=s;
82 }
83 else
84 {
85 seed=prevseed+419; // these numbers are arbitrary
86 if (seed>3040101) // I want to avoid overflow on 32bit systems
87 {
88 seed=((int)(seed)%0xABCD)+1;
89 }
90 }
91 }
92 // now we need to consider MPI since we don't want each rank to start with the same seed. Rank in COMM_WORLD will do
93 #ifdef ESYS_MPI
94 Esys_MPI_rank rank;
95 int mperr=MPI_Comm_rank(MPI_COMM_WORLD, &rank);
96 if (mperr!=MPI_SUCCESS) {
97 rank=0;
98 }
99 seed+=rank*5;
100 #endif
101 prevseed=seed;
102
103 boost::mt19937::result_type RMAX=base.max();
104 seedGens(seed);
105 long i;
106
107 #pragma omp parallel private(i)
108 {
109 int tnum=0;
110 #ifdef _OPENMP
111 tnum=omp_get_thread_num();
112 #endif
113 boost::mt19937& generator=*(gens[tnum]);
114
115 #pragma omp for schedule(static)
116 for (i=0;i<n;++i)
117 {
118 #ifdef _WIN32
119 array[i]=((double)generator())/RMAX;
120 #else
121 array[i]=((double)generator())/RMAX;
122 #endif
123 }
124 }
125 }
126
127 // see patternFillArray for details on parameters
128 void patternFillArray2D(size_t x, size_t y, double* array, size_t spacing, size_t basex, size_t basey, size_t numpoints)
129 {
130 memset(array, 0, x*y*sizeof(double)*numpoints);
131 size_t xoff=basex%spacing;
132 size_t yoff=basey%spacing;
133 for (int r=0;r<y;++r)
134 {
135 size_t step=((r+yoff)%spacing)?spacing:1;
136 for (int c=(spacing-xoff)%spacing;c<x;c+=step)
137 {
138 for (int p=0;p<numpoints;++p)
139 {
140 array[(c+r*x)*numpoints+p]=1+p;
141 }
142 }
143 }
144 }
145
146
147 // fill the array (which we assume is 3D with x by y by z points in it) with a pattern.
148 // The base? params give the coordinates (in # of elements) of the origin of _this_ rank
149 // used to ensure patterns are generated consistantly across multiple ranks
150 // This is only for internal debug so the patterns (or this function) may disappear
151 // without notice
152 void patternFillArray(int pattern, size_t x, size_t y, size_t z, double* array, size_t spacing, size_t basex, size_t basey, size_t basez, size_t numpoints)
153 {
154 if (pattern==0) // a cross pattern in the z=0 plane, repeated for each z layer
155 {
156 memset(array, 0, x*y*sizeof(double));
157 size_t xoff=basex%spacing;
158 size_t yoff=basey%spacing;
159 for (int r=0;r<y;++r)
160 {
161 size_t step=((r+yoff)%spacing)?spacing:1;
162 for (int c=(spacing-xoff)%spacing;c<x;c+=step)
163 {
164 for (int p=0;p<numpoints;++p)
165 {
166 array[(c+r*x)*numpoints+p]=p+1;
167 }
168 }
169 }
170 for (int l=1;l<z;++l)
171 {
172 memcpy(array+(x*y*l*numpoints), array, x*y*sizeof(double)*numpoints);
173 }
174 }
175 else // pattern 1. A grid in all 3 dimensions
176 {
177 if (z<2)
178 {
179 patternFillArray(0, x, y, z, array, spacing, basex, basey, basez, numpoints);
180 return; // this pattern needs a minimum of 2 layers
181 }
182 size_t xoff=basex%spacing;
183 size_t yoff=basey%spacing;
184 size_t zoff=basez%spacing;
185
186 double* buff1=new double[x*y*numpoints]; // stores the main cross pattern
187 double* buff2=new double[x*y*numpoints]; // stores the "verticals"
188 memset(buff1, 0, x*y*sizeof(double));
189 memset(buff2, 0, x*y*sizeof(double));
190 for (size_t r=0;r<y;++r)
191 {
192 size_t step=((r+yoff)%spacing)?spacing:1;
193 if (step==1)
194 {
195 for (size_t c=0;c<x;++c)
196 {
197 for (int p=0;p<numpoints;++p)
198 {
199 buff1[(c+r*x)*numpoints+p]=p+1;
200 }
201 }
202 }
203 else
204 {
205 for (size_t c=(spacing-xoff)%spacing;c<x;c+=step)
206 {
207 for (int p=0;p<numpoints;++p)
208 {
209 buff1[(c+r*x)*numpoints+p]=p+1;
210 }
211 }
212 }
213 }
214 for (size_t r=(spacing-yoff)%spacing;r<y;r+=spacing)
215 {
216 for (size_t c=(spacing-xoff)%spacing;c<x;c+=spacing)
217 {
218 for (int p=0;p<numpoints;++p)
219 {
220 buff2[(c+r*x)*numpoints+p]=p+1;
221 }
222 }
223 }
224 for (size_t l=0;l<z;++l)
225 {
226 if ((l+zoff)%spacing)
227 {
228 memcpy(array+(x*y*l*numpoints), buff2, x*y*sizeof(double)*numpoints);
229 }
230 else
231 {
232 memcpy(array+(x*y*l*numpoints), buff1, x*y*sizeof(double)*numpoints);
233 }
234 }
235 delete[] buff1;
236 delete[] buff2;
237 }
238
239
240
241 }
242
243 } // end namespace

  ViewVC Help
Powered by ViewVC 1.1.26