su2hmc
Loading...
Searching...
No Matches
random.c
Go to the documentation of this file.
1
6#ifdef __USE_MKL__
7#include <mkl.h>
8//Bad practice? Yes but it is convenient
9#endif
10#include "random.h"
11#include <time.h>
12
13//Declaring external variables
14#ifdef __RANLUX__
16gsl_rng *ranlux_instd;
18unsigned long seed;
19#else
21long seed;
22#endif
23#ifndef M_PI
25#define M_PI acos(-1)
26#endif
27
28#ifdef __RANLUX__
29inline int ranset(unsigned long *seed)
30#else
31inline int ranset(long *seed)
32#endif
33{
34#ifdef __RANLUX__
35 ranlux_instd=gsl_rng_alloc(gsl_rng_ranlxd2);
36 gsl_rng_set(ranlux_instd,*seed);
37 return 0;
38#else
39 return 0;
40#endif
41}
42int Par_ranread(char *filename, double *ranval){
43 const char funcname[] = "Par_psread";
44 FILE *dest;
45 if(!rank){
46 if(!(dest = fopen(filename, "rb"))){
47 fprintf(stderr, "Error %i in %s: Failed to open %s.\nExiting...\n\n", OPENERROR, funcname, filename);
48#if(nproc>1)
49 MPI_Abort(comm,OPENERROR);
50#else
51 exit(OPENERROR);
52#endif
53
54 }
55 fread(&ranval, sizeof(ranval), 1, dest);
56 fclose(dest);
57 }
58#if(nproc>1)
59 Par_dcopy(ranval);
60#endif
61 return 0;
62}
63#ifdef __RANLUX__
64int Par_ranset(unsigned long *seed,int iread)
65#else
66int Par_ranset(long *seed,int iread)
67#endif
68{
69 const char funcname[] = "Par_ranset";
70 //If we're not using the master thread, we need to change the seed
71#ifdef _DEBUG
72 printf("Master seed: %lu\t",*seed);
73#endif
74 if(rank)
75 *seed *= 1.0f+8.0f*(float)rank/(float)(size-1);
76#ifdef _DEBUG
77 printf("Rank: %i\tSeed %lu\n",rank, *seed);
78#endif
79 //Next we set the seed using ranset
80 //This is one of the really weird FORTRAN 66-esque functions with ENTRY points, so good luck!
81#ifdef __RANLUX__
82 return ranset(seed);
83#else
84 return 0;
85#endif
86}
87double Par_granf(){
88 const char funcname[] = "Par_granf";
89 double ran_val=0;
90 if(!rank){
91#ifdef __RANLUX__
92 ran_val = gsl_rng_uniform(ranlux_instd);
93 #else
94 ran_val = ran2(&seed);
95#endif
96 }
97#if(nproc>1)
98 Par_dcopy(&ran_val);
99#endif
100 return ran_val;
101}
102int Gauss_z(Complex *ps, unsigned int n, const Complex mu, const double sigma){
103 const char funcname[] = "Gauss_z";
104 if(n<=0){
105 fprintf(stderr, "Error %i in %s: Array cannot have length %i.\nExiting...\n\n",
106 ARRAYLEN, funcname, n);
107#if(nproc>1)
108 MPI_Abort(comm,ARRAYLEN);
109#else
110 exit(ARRAYLEN);
111#endif
112 }
113#pragma unroll
114 for(int i=0;i<n;i++){
115 /* Marsaglia Method for fun
116 do{
117 u=sfmt_genrand_real1(sfmt);
118 v=sfmt_genrand_real1(sfmt);
119 r=u*u+v*v;
120 }while(0<r & r<1);
121 r=sqrt(r);
122 r=sqrt(-2.0*log(r)/r)*sigma;
123 ps[i] = mu+u*r + I*(mu+v*r);
124 */
125#ifdef __RANLUX__
126 double r =sigma*sqrt(-2*log(gsl_rng_uniform(ranlux_instd)));
127 double theta=2.0*M_PI*gsl_rng_uniform(ranlux_instd);
128#else
129 double r =sigma*sqrt(-2*log(ran2(&seed)));
130 double theta=2.0*M_PI*ran2(&seed);
131#endif
132 ps[i]=r*(cos(theta)+sin(theta)*I)+mu;
133 }
134 return 0;
135}
136int Gauss_c(Complex_f *ps, unsigned int n, const Complex_f mu, const float sigma){
137 const char funcname[] = "Gauss_z";
138 if(n<=0){
139 fprintf(stderr, "Error %i in %s: Array cannot have length %i.\nExiting...\n\n",
140 ARRAYLEN, funcname, n);
141#if(nproc>1)
142 MPI_Abort(comm,ARRAYLEN);
143#else
144 exit(ARRAYLEN);
145#endif
146 }
147#pragma unroll
148 for(int i=0;i<n;i++){
149 /* Marsaglia Method for fun
150 do{
151 u=sfmt_genrand_real1(sfmt);
152 v=sfmt_genrand_real1(sfmt);
153 r=u*u+v*v;
154 }while(0<r & r<1);
155 r=sqrt(r);
156 r=sqrt(-2.0*log(r)/r)*sigma;
157 ps[i] = mu+u*r + I*(mu+v*r);
158 */
159#ifdef __RANLUX__
160 float r =sigma*sqrt(-2*log(gsl_rng_uniform(ranlux_instd)));
161 float theta=2.0*M_PI*gsl_rng_uniform(ranlux_instd);
162#else
163 float r =sigma*sqrt(-2*log(ran2(&seed)));
164 float theta=2.0*M_PI*ran2(&seed);
165#endif
166 ps[i]=r*(cos(theta)+mu+sin(theta)*I)+mu;
167 }
168 return 0;
169}
170int Gauss_d(double *ps, unsigned int n, const double mu, const double sigma){
171 const char funcname[] = "Gauss_z";
172 //The FORTRAN Code had two different Gauss Routines. gaussp having unit
173 //mean and variance and gauss0 where the variance would appear to be 1/sqrt(2)
174 //(Since we multiply by sqrt(-ln(r)) instead of sqrt(-2ln(r)) )
175 if(n<=0){
176 fprintf(stderr, "Error %i in %s: Array cannot have length %i.\nExiting...\n\n",
177 ARRAYLEN, funcname, n);
178#if(nproc>1)
179 MPI_Abort(comm,ARRAYLEN);
180#else
181 exit(ARRAYLEN);
182#endif
183 }
184 int i;
185 double r, u, v;
186 //If n is odd we calculate the last index seperately and the rest in pairs
187 if(n%2==1){
188 n--;
189#ifdef __RANLUX__
190 r=2.0*M_PI*gsl_rng_uniform(ranlux_instd);
191 ps[n]=sqrt(-2*log(gsl_rng_uniform(ranlux_instd)))*cos(r);
192#else
193 r=2.0*M_PI*ran2(&seed);
194 ps[n]=sqrt(-2*log(ran2(&seed)))*cos(r);
195#endif
196 }
197 for(i=0;i<n;i+=2){
198 /* Marsaglia Method for fun
199 do{
200 u=sfmt_genrand_real1(sfmt);
201 v=sfmt_genrand_real1(sfmt);
202 r=u*u+v*v;
203 }while(0<r & r<1);
204 r=sqrt(r);
205 r=sqrt(-2.0*log(r)/r)*sigma;
206 ps[i] = mu+u*r;
207 ps[i+1]=mu+v*r;
208 */
209#ifdef __RANLUX__
210 u=sqrt(-2*log(gsl_rng_uniform(ranlux_instd)))*sigma;
211 r=2.0*M_PI*gsl_rng_uniform(ranlux_instd);
212#else
213 u=sqrt(-2*log(ran2(&seed)))*sigma;
214 r=2.0*M_PI*ran2(&seed);
215#endif
216 ps[i]=u*cos(r)+mu;
217 ps[i+1]=u*sin(r)+mu;
218 }
219 return 0;
220}
221int Gauss_f(float *ps, unsigned int n, const float mu, const float sigma){
222 const char funcname[] = "Gauss_z";
223 //The FORTRAN Code had two different Gauss Routines. gaussp having unit
224 //mean and variance and gauss0 where the variance would appear to be 1/sqrt(2)
225 //(Since we multiply by sqrt(-ln(r)) instead of sqrt(-2ln(r)) )
226 if(n<=0){
227 fprintf(stderr, "Error %i in %s: Array cannot have length %i.\nExiting...\n\n",
228 ARRAYLEN, funcname, n);
229#if(nproc>1)
230 MPI_Abort(comm,ARRAYLEN);
231#else
232 exit(ARRAYLEN);
233#endif
234 }
235 int i;
236 float r, u, v;
237 //If n is odd we calculate the last index seperately and the rest in pairs
238 if(n%2==1){
239 n--;
240#ifdef __RANLUX__
241 r=2.0*M_PI*gsl_rng_uniform(ranlux_instd);
242 ps[n]=sqrt(-2*log(gsl_rng_uniform(ranlux_instd)))*cos(r);
243#else
244 r=2.0*M_PI*ran2(&seed);
245 ps[n]=sqrt(-2*log(ran2(&seed)))*cos(r);
246#endif
247 }
248#ifdef __RANLUX__
249 r=2.0*M_PI*gsl_rng_uniform(ranlux_instd);
250 ps[n]=sqrt(-2*log(gsl_rng_uniform(ranlux_instd)))*cos(r);
251#else
252 r=2.0*M_PI*ran2(&seed);
253 ps[n]=sqrt(-2*log(ran2(&seed)))*cos(r);
254#endif
255 for(i=0;i<n;i+=2){
256 /* Marsaglia Method for fun
257 do{
258 u=sfmt_genrand_real1(sfmt);
259 v=sfmt_genrand_real1(sfmt);
260 r=u*u+v*v;
261 }while(0<r & r<1);
262 r=sqrt(r);
263 r=sqrt(-2.0*log(r)/r)*sigma;
264 ps[i] = mu+u*r;
265 ps[i+1]=mu+v*r;
266 */
267#ifdef __RANLUX__
268 u=sqrt(-2*log(gsl_rng_uniform(ranlux_instd)))*sigma;
269 r=2.0*M_PI*gsl_rng_uniform(ranlux_instd);
270#else
271 u=sqrt(-2*log(ran2(&seed)))*sigma;
272 r=2.0*M_PI*ran2(&seed);
273#endif
274 ps[i]=u*cos(r)+mu;
275 ps[i+1]=u*sin(r)+mu;
276 }
277 return 0;
278}
279#ifndef __RANLUX__
280double ran2(long *idum) {
281 long k;
282 int j;
283 static long idum2=123456789;
284 static long iy=0;
285 static long iv[NTAB];
287#pragma omp threadprivate(idum2, iy, iv)
288 //No worries
289 double temp;
290
291 if (*idum <= 0) {
292 if (-(*idum) < 1) *idum=1;
293 else *idum = -(*idum);
294 {
295 idum2=(*idum);
296
297 for(j=NTAB+7;j>=0;j--) {
298 k=(*idum)/IQ1;
299 *idum=IA1*(*idum-k*IQ1)-k*IR1;
300 if (*idum < 0) *idum += IM1;
301 if (j < NTAB)
302 iv[j] = *idum;
303 }
304 iy=iv[0];
305 }
306 }
307 k=(*idum)/IQ1;
308 *idum=IA1*(*idum-k*IQ1)-k*IR1;
309 if (*idum < 0) *idum += IM1;
310 k=idum2/IQ2;
311 idum2=IA2*(idum2-k*IQ2)-k*IR2;
312
313 if (idum2 < 0) idum2 += IM2; j=iy/NDIV;
314 iy=iv[j]-idum2;
315 iv[j] = *idum;
316 if (iy < 1) iy += IMM1;
317 if ((temp=AM*iy) > RNMX)
318 return RNMX;
319
320 else return temp;
321
322}
323#endif
324
ARRAYLEN
#define ARRAYLEN
Impossible value for array length.
Definition errorcodes.h:64
OPENERROR
#define OPENERROR
Error opening file.
Definition errorcodes.h:32
Par_dcopy
int Par_dcopy(double *dval)
Broadcasts a double to the other processes.
Gauss_f
int Gauss_f(float *ps, unsigned int n, const float mu, const float sigma)
Generates a vector of normally distributed random single precision numbers using the Box-Muller Metho...
Definition random.c:221
Par_ranread
int Par_ranread(char *filename, double *ranval)
Reads ps from a file Since this function is very similar to Par_sread, I'm not really going to commen...
Definition random.c:42
Gauss_z
int Gauss_z(Complex *ps, unsigned int n, const Complex mu, const double sigma)
Generates a vector of normally distributed random double precision complex numbers using the Box-Mull...
Definition random.c:102
Gauss_c
int Gauss_c(Complex_f *ps, unsigned int n, const Complex_f mu, const float sigma)
Generates a vector of normally distributed random single precision complex numbers using the Box-Mull...
Definition random.c:136
Gauss_d
int Gauss_d(double *ps, unsigned int n, const double mu, const double sigma)
Generates a vector of normally distributed random double precision numbers using the Box-Muller Metho...
Definition random.c:170
Par_ranset
int Par_ranset(unsigned long *seed, int iread)
Uses the rank to get a new seed. Copying from the FORTRAN description here c create new seeds in rang...
Definition random.c:64
Par_granf
double Par_granf()
Generates a random double which is then sent to the other ranks.
Definition random.c:87
ranset
int ranset(unsigned long *seed)
Seed the ranlux generator from GSL.
Definition random.c:29
size
int size
The number of MPI ranks in total.
Definition par_mpi.c:20
rank
int rank
The MPI rank.
Definition par_mpi.c:20
M_PI
#define M_PI
if not defined elsewhere
Definition random.c:25
random.h
Header for random number configuration.
seed
unsigned long seed
RANLUX seed.
Definition random.c:18
ranlux_instd
gsl_rng * ranlux_instd
RANLUX instance.
Definition random.c:16
Complex
#define Complex
Double precision complex number.
Definition sizes.h:64
Complex_f
#define Complex_f
Single precision complex number.
Definition sizes.h:62
I
#define I
Define I in double precision using C standard notation.
Definition thrust_complex.h:28