Bridge++  Ver. 1.1.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
randomNumbers.cpp
Go to the documentation of this file.
1 
15 #include "randomNumbers.h"
16 
17 #include <iostream>
18 #include <fstream>
19 #include <cassert>
20 
21 #include "communicator.h"
22 #include "index_lex.h"
23 #include "index_eo.h"
24 
25 using Bridge::vout;
26 
27 namespace {
28  const double sq2r = 1.0 / sqrt(2.0);
29 // const double pi = 3.141592653589793;
30  const double pi = 4.0 * atan(1.0);
31  const double pi2 = 2.0 * pi;
32 }
33 
34 //====================================================================
35 void RandomNumbers::gauss(double& rn1, double& rn2)
36 {
37  // Two Gaussian random number with deviation 1/\sqrt(2).
38  double r1 = get();
39  double r2 = get();
40 
41  // double slg1 = sqrt( -log(r1)*2.0 )*sq2r;
42  double slg1 = sqrt(-log(1 - r1) * 2.0) / sqrt(2.0);
43  double ang1 = pi2 * r2;
44 
45  rn1 = slg1 * cos(ang1);
46  rn2 = slg1 * sin(ang1);
47 }
48 
49 
50 //====================================================================
51 void RandomNumbers::uniform(Field& f)
52 {
53  int Nin = f.nin();
54  int Nvol = f.nvol();
55  int Nex = f.nex();
56 
57  for (int ex = 0; ex < Nex; ++ex) {
58  for (int site = 0; site < Nvol; ++site) {
59  for (int in = 0; in < Nin; ++in) {
60  double rn = get();
61  f.set(in, site, ex, rn);
62  }
63  }
64  }
65 }
66 
67 
68 //====================================================================
69 void RandomNumbers::gauss(Field& f)
70 {
71  int Nin = f.nin();
72  int Nvol = f.nvol();
73  int Nex = f.nex();
74 
75  double rn1, rn2;
76 
77  if ((Nin % 2) == 0) { // Nin is even
78  int Nin2 = Nin / 2;
79  for (int ex = 0; ex < Nex; ++ex) {
80  for (int site = 0; site < Nvol; ++site) {
81  for (int in = 0; in < Nin2; ++in) {
82  //assert(2*in+1 < Nin);
83  gauss(rn1, rn2);
84  f.set(2 * in, site, ex, rn1);
85  f.set(2 * in + 1, site, ex, rn2);
86  }
87  }
88  }
89  } else {
90  for (int ex = 0; ex < Nex; ++ex) {
91  for (int site = 0; site < Nvol; ++site) {
92  for (int in = 0; in < Nin; ++in) {
93  gauss(rn1, rn2);
94  f.set(in, site, ex, rn1);
95  // rn2 is not used
96  }
97  }
98  }
99  }
100 }
101 
102 
103 //====================================================================
104 void RandomNumbers::gauss_lex_global(Field& f)
105 {
106  Index_lex idx;
107 
108  int Nin = f.nin();
109  int Nvol = f.nvol();
110  int Nex = f.nex();
111 
112  int Lvol = CommonParameters::Lvol();
113 
114  Field f_tmp(Nin, Lvol, Nex);
115 
116  gauss(f_tmp);
117 
118  int Lx = CommonParameters::Lx();
119  int Ly = CommonParameters::Ly();
120  int Lz = CommonParameters::Lz();
121  int Lt = CommonParameters::Lt();
122  int Nx = CommonParameters::Nx();
123  int Ny = CommonParameters::Ny();
124  int Nz = CommonParameters::Nz();
125  int Nt = CommonParameters::Nt();
126  int Ipx = Communicator::ipe(0);
127  int Ipy = Communicator::ipe(1);
128  int Ipz = Communicator::ipe(2);
129  int Ipt = Communicator::ipe(3);
130 
131  for (int ex = 0; ex < Nex; ++ex) {
132  for (int t = 0; t < Nt; ++t) {
133  int t2 = t + Ipt * Nt;
134  for (int z = 0; z < Nz; ++z) {
135  int z2 = z + Ipz * Nz;
136  for (int y = 0; y < Ny; ++y) {
137  int y2 = y + Ipy * Ny;
138  for (int x = 0; x < Nx; ++x) {
139  int x2 = x + Ipx * Nx;
140  int site = idx.site(x, y, z, t);
141  int site2 = x2 + Lx * (y2 + Ly * (z2 + Lz * t2));
142  assert(site < Nvol);
143  assert(site2 < Lvol);
144  for (int in = 0; in < Nin; ++in) {
145  f.set(in, site, ex, f_tmp.cmp(in, site2, ex));
146  }
147  }
148  }
149  }
150  }
151  }
152 }
153 
154 
155 //====================================================================
156 void RandomNumbers::uniform_lex_global(Field& f)
157 {
158  Index_lex idx;
159 
160  int Nin = f.nin();
161  int Nvol = f.nvol();
162  int Nex = f.nex();
163 
164  int Lvol = CommonParameters::Lvol();
165 
166  Field f_tmp(Nin, Lvol, Nex);
167 
168  uniform(f_tmp);
169 
170  int Lx = CommonParameters::Lx();
171  int Ly = CommonParameters::Ly();
172  int Lz = CommonParameters::Lz();
173  int Lt = CommonParameters::Lt();
174  int Nx = CommonParameters::Nx();
175  int Ny = CommonParameters::Ny();
176  int Nz = CommonParameters::Nz();
177  int Nt = CommonParameters::Nt();
178  int Ipx = Communicator::ipe(0);
179  int Ipy = Communicator::ipe(1);
180  int Ipz = Communicator::ipe(2);
181  int Ipt = Communicator::ipe(3);
182 
183  for (int ex = 0; ex < Nex; ++ex) {
184  for (int t = 0; t < Nt; ++t) {
185  int t2 = t + Ipt * Nt;
186  for (int z = 0; z < Nz; ++z) {
187  int z2 = z + Ipz * Nz;
188  for (int y = 0; y < Ny; ++y) {
189  int y2 = y + Ipy * Ny;
190  for (int x = 0; x < Nx; ++x) {
191  int x2 = x + Ipx * Nx;
192  int site = idx.site(x, y, z, t);
193  int site2 = x2 + Lx * (y2 + Ly * (z2 + Lz * t2));
194  assert(site < Nvol);
195  assert(site2 < Lvol);
196  for (int in = 0; in < Nin; ++in) {
197  f.set(in, site, ex, f_tmp.cmp(in, site2, ex));
198  }
199  }
200  }
201  }
202  }
203  }
204 }
205 
206 
207 //====================================================================
208 void RandomNumbers::gauss_eo_global(Field& f)
209 {
210  Index_eo idx;
211 
212  int Nin = f.nin();
213  int Nvol = f.nvol();
214  int Nex = f.nex();
215 
216  int Lvol2 = CommonParameters::Lvol() / 2;
217 
218  Field f_tmp(Nin, Lvol2, Nex);
219 
220  gauss(f_tmp);
221 
222  int Lx2 = CommonParameters::Lx() / 2;
223  int Ly = CommonParameters::Ly();
224  int Lz = CommonParameters::Lz();
225  int Lt = CommonParameters::Lt();
226  int Nx2 = CommonParameters::Nx() / 2;
227  int Ny = CommonParameters::Ny();
228  int Nz = CommonParameters::Nz();
229  int Nt = CommonParameters::Nt();
230  int Ipx = Communicator::ipe(0);
231  int Ipy = Communicator::ipe(1);
232  int Ipz = Communicator::ipe(2);
233  int Ipt = Communicator::ipe(3);
234 
235  for (int ex = 0; ex < Nex; ++ex) {
236  for (int t = 0; t < Nt; ++t) {
237  int t2 = t + Ipt * Nt;
238  for (int z = 0; z < Nz; ++z) {
239  int z2 = z + Ipz * Nz;
240  for (int y = 0; y < Ny; ++y) {
241  int y2 = y + Ipy * Ny;
242  for (int x = 0; x < Nx2; ++x) {
243  int x2 = x + Ipx * Nx2;
244  int site = idx.site(x, y, z, t, 0);
245  int site2 = x2 + Lx2 * (y2 + Ly * (z2 + Lz * t2));
246  assert(site < Nvol);
247  assert(site2 < Lvol2);
248  for (int in = 0; in < Nin; ++in) {
249  f.set(in, site, ex, f_tmp.cmp(in, site2, ex));
250  }
251  }
252  }
253  }
254  }
255  }
256 }
257 
258 
259 //====================================================================
260 //============================================================END=====