su2hmc
Loading...
Searching...
No Matches
PRNG

Functions

int ranset (unsigned long *seed)
 Seed the ranlux generator from GSL.
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 range seed to 9*seed c having a range of 0*seed gave an unfortunate pattern c in the underlying value of ds(1) (it was always 10 times bigger c on the last processor). This does not appear to happen with 9.
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-Muller Method.
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 Method.
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-Muller Method.
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 Method.
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 comment it check there if you are confused about things.
double Par_granf ()
 Generates a random double which is then sent to the other ranks.

Detailed Description

Function Documentation

◆ 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-Muller Method.

Parameters
[out]psThe output array
[in]nThe array length
[in]mumean
[in]sigmavariance
Returns
Zero on success integer error code otherwise

Definition at line 136 of file random.c.

136 {
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}
ARRAYLEN
#define ARRAYLEN
Impossible value for array length.
Definition errorcodes.h:64
M_PI
#define M_PI
if not defined elsewhere
Definition random.c:25
seed
unsigned long seed
RANLUX seed.
Definition random.c:18
ranlux_instd
gsl_rng * ranlux_instd
RANLUX instance.
Definition random.c:16
I
#define I
Define I in double precision using C standard notation.
Definition thrust_complex.h:28

References ARRAYLEN, Complex_f, I, M_PI, ranlux_instd, and seed.

Here is the caller graph for this function:

◆ 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 Method.

Parameters
[out]psThe output array
[in]nThe array length
[in]mumean
[in]sigmavariance
Returns
Zero on success integer error code otherwise

Definition at line 170 of file random.c.

170 {
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}

References ARRAYLEN, M_PI, ranlux_instd, and seed.

Here is the caller graph for this function:

◆ 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 Method.

Parameters
[out]psThe output array
[in]nThe array length
[in]mumean
[in]sigmavariance
Returns
Zero on success integer error code otherwise

Definition at line 221 of file random.c.

221 {
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}

References ARRAYLEN, M_PI, ranlux_instd, and seed.

◆ 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-Muller Method.

Parameters
[out]psThe output array
[in]nThe array length
[in]mumean
[in]sigmavariance
Returns
Zero on success integer error code otherwise

Definition at line 102 of file random.c.

102 {
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}

References ARRAYLEN, Complex, I, M_PI, ranlux_instd, and seed.

◆ Par_granf()

double Par_granf ( )

Generates a random double which is then sent to the other ranks.

Returns
the random number generated

Definition at line 87 of file random.c.

87 {
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}
Par_dcopy
int Par_dcopy(double *dval)
Broadcasts a double to the other processes.
rank
int rank
The MPI rank.
Definition par_mpi.c:20

References Par_dcopy(), rank, ranlux_instd, and seed.

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 comment it check there if you are confused about things.

Parameters
[in]filenameThe name of the file we're reading from
[in]ranvalThe destination for the file's contents
Returns
Zero on success, integer error code otherwise

Definition at line 42 of file random.c.

42 {
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}
OPENERROR
#define OPENERROR
Error opening file.
Definition errorcodes.h:32

References OPENERROR, Par_dcopy(), and rank.

Here is the call graph for this function:

◆ 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 range seed to 9*seed c having a range of 0*seed gave an unfortunate pattern c in the underlying value of ds(1) (it was always 10 times bigger c on the last processor). This does not appear to happen with 9.

Parameters
[in]seedThe seed from the rank in question.
[in]ireadDo we read from file or not. Don't remember why it's here as it's not used
Returns
Zero on success, integer error code otherwise

Definition at line 64 of file random.c.

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}
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

References rank, ranset(), seed, and size.

Here is the call graph for this function:
Here is the caller graph for this function:

◆ ranset()

int ranset ( unsigned long * seed)
inline

Seed the ranlux generator from GSL.

Parameters
[in]seedpointer to seed
Returns
0

Definition at line 29 of file random.c.

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}

References ranlux_instd, and seed.

Here is the caller graph for this function: