Bridge++  Ver. 1.1.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
field_F_imp.cpp
Go to the documentation of this file.
1 
14 #include "field_F.h"
15 
16 #include "bridgeIO.h"
17 using Bridge::vout;
18 
19 using std::valarray;
20 
21 // This implementation only applies to SU(3) group and Nd=4 case.
22 #define NC 3
23 #define NC2 6
24 #define NDF 18
25 #define ND 4
26 #define NCD 12
27 #define NCD2 24
28 
29 //====================================================================
30 void Field_F::check()
31 {
32  assert(NC == CommonParameters::Nc());
33  assert(ND == CommonParameters::Nd());
34 }
35 
36 
37 //====================================================================
38 void Field_F::mult_GM(const GammaMatrix& gm, const Field_F& x)
39 {
40  assert(x.nex() == m_Nex);
41  assert(x.nvol() == m_Nvol);
42 
43  int id[ND];
44  int idc_r[ND];
45  int idc_i[ND];
46  double gv_r[ND];
47  double gv_i[ND];
48 
49  for (int s = 0; s < ND; ++s) {
50  id[s] = gm.index(s);
51  gv_r[s] = gm.value_r(s);
52  gv_i[s] = gm.value_i(s);
53  idc_r[s] = gm.index_c(s);
54  idc_i[s] = 1 - idc_r[s];
55  }
56 
57  double *w;
58  w = const_cast<Field_F *>(&x)->ptr(0);
59 
60  for (int ex = 0; ex < m_Nex; ++ex) {
61  for (int site = 0; site < m_Nvol; ++site) {
62  int iv = NCD2 * (site + m_Nvol * ex);
63  for (int s = 0; s < ND; ++s) {
64  int iv2 = s * NC2 + iv;
65  int iw2 = id[s] * NC2 + iv;
66  for (int ic = 0; ic < NC; ++ic) {
67  field[2 * ic + iv2] = gv_r[s] * w[2 * ic + idc_r[s] + iw2];
68  field[2 * ic + 1 + iv2] = gv_i[s] * w[2 * ic + idc_i[s] + iw2];
69  }
70  }
71  }
72  }
73 }
74 
75 
76 //====================================================================
77 void Field_F::mult_iGM(const GammaMatrix& gm, const Field_F& x)
78 {
79  assert(x.nex() == m_Nex);
80  assert(x.nvol() == m_Nvol);
81 
82  int id[ND];
83  int idc_r[ND];
84  int idc_i[ND];
85  double gv_r[ND];
86  double gv_i[ND];
87 
88  for (int s = 0; s < ND; ++s) {
89  id[s] = gm.index(s);
90  gv_r[s] = -gm.value_i(s);
91  gv_i[s] = gm.value_r(s);
92  idc_r[s] = 1 - gm.index_c(s);
93  idc_i[s] = 1 - idc_r[s];
94  }
95 
96  double *w;
97  w = const_cast<Field_F *>(&x)->ptr(0);
98 
99  for (int ex = 0; ex < m_Nex; ++ex) {
100  for (int site = 0; site < m_Nvol; ++site) {
101  int iv = NCD2 * (site + m_Nvol * ex);
102  for (int s = 0; s < ND; ++s) {
103  int iv2 = s * NC2 + iv;
104  int iw2 = id[s] * NC2 + iv;
105  for (int ic = 0; ic < NC; ++ic) {
106  field[2 * ic + iv2] = gv_r[s] * w[2 * ic + idc_r[s] + iw2];
107  field[2 * ic + 1 + iv2] = gv_i[s] * w[2 * ic + idc_i[s] + iw2];
108  }
109  }
110  }
111  }
112 }
113 
114 
115 //====================================================================
116 void Field_F::mult_GMproj(const int pm, const GammaMatrix& gm,
117  const Field_F& w)
118 {
119  assert(w.nvol() == m_Nvol);
120  assert(w.nex() == m_Nex);
121 
122  Field_F v(w.nvol(), w.nex());
123  v.mult_GM(gm, w);
124 
125  int size = m_Nin * m_Nvol * m_Nex;
126  if (pm == 1) {
127  for (int i = 0; i < size; ++i) {
128  field[i] = 0.5 * (w.field[i] + v.field[i]);
129  }
130  } else if (pm == -1) {
131  for (int i = 0; i < size; ++i) {
132  field[i] = 0.5 * (w.field[i] - v.field[i]);
133  }
134  } else {
135  vout.general(m_vl, "Field_F: wrong pm = %d\n", pm);
136  abort();
137  }
138 }
139 
140 
141 //====================================================================
142 void Field_F::mult_GMproj2(const int pm, const GammaMatrix& gm,
143  const Field_F& w)
144 {
145  assert(w.nvol() == m_Nvol);
146  assert(w.nex() == m_Nex);
147 
148  Field_F v(w.nvol(), w.nex());
149  v.mult_GM(gm, w);
150 
151  int size = m_Nin * m_Nvol * m_Nex;
152  if (pm == 1) {
153  for (int i = 0; i < size; ++i) {
154  field[i] = w.field[i] + v.field[i];
155  }
156  } else if (pm == -1) {
157  for (int i = 0; i < size; ++i) {
158  field[i] = w.field[i] - v.field[i];
159  }
160  } else {
161  vout.general(m_vl, "Field_F: wrong pm = %d\n", pm);
162  abort();
163  }
164 }
165 
166 
167 //====================================================================
168 void Field_F::mult_Field_Gn(int ex, const Field_G& U, int ex1,
169  const Field_F& x, int ex2)
170 {
171  assert(ex <= m_Nex);
172  assert(ex1 <= U.nex());
173  assert(ex2 <= x.nex());
174  assert(U.nvol() == m_Nvol);
175  assert(x.nvol() == m_Nvol);
176 
177  double *w;
178  double *g;
179  w = const_cast<Field_F *>(&x)->ptr(0);
180  g = const_cast<Field_G *>(&U)->ptr(0);
181 
182  for (int site = 0; site < m_Nvol; ++site) {
183  int ig = NDF * (site + m_Nvol * ex1);
184  int iw = NCD2 * (site + m_Nvol * ex2);
185  int iv = NCD2 * (site + m_Nvol * ex);
186  for (int s = 0; s < ND; ++s) {
187  for (int ic = 0; ic < NC; ++ic) {
188  int ig2 = ic * NC2 + ig;
189  int iw2 = s * NC2 + iw;
190  int iv2 = s * NC2 + iv;
191  field[2 * ic + iv2] =
192  g[0 + ig2] * w[0 + iw2] - g[1 + ig2] * w[1 + iw2]
193  + g[2 + ig2] * w[2 + iw2] - g[3 + ig2] * w[3 + iw2]
194  + g[4 + ig2] * w[4 + iw2] - g[5 + ig2] * w[5 + iw2];
195  field[2 * ic + 1 + iv2] =
196  g[0 + ig2] * w[1 + iw2] + g[1 + ig2] * w[0 + iw2]
197  + g[2 + ig2] * w[3 + iw2] + g[3 + ig2] * w[2 + iw2]
198  + g[4 + ig2] * w[5 + iw2] + g[5 + ig2] * w[4 + iw2];
199  }
200  }
201  }
202 }
203 
204 
205 //====================================================================
206 void Field_F::mult_Field_Gd(int ex, const Field_G& U, int ex1,
207  const Field_F& x, int ex2)
208 {
209  assert(ex <= m_Nex);
210  assert(ex1 <= U.nex());
211  assert(ex2 <= x.nex());
212  assert(U.nvol() == m_Nvol);
213  assert(x.nvol() == m_Nvol);
214 
215  double *w;
216  double *g;
217  w = const_cast<Field_F *>(&x)->ptr(0);
218  g = const_cast<Field_G *>(&U)->ptr(0);
219 
220  for (int site = 0; site < m_Nvol; ++site) {
221  int ig = NDF * (site + m_Nvol * ex1);
222  int iw = NCD2 * (site + m_Nvol * ex2);
223  int iv = NCD2 * (site + m_Nvol * ex);
224  for (int s = 0; s < ND; ++s) {
225  for (int ic = 0; ic < NC; ++ic) {
226  int ig2 = ic * 2 + ig;
227  int iw2 = s * NC2 + iw;
228  int iv2 = s * NC2 + iv;
229  field[2 * ic + iv2] =
230  g[0 + ig2] * w[0 + iw2] + g[1 + ig2] * w[1 + iw2]
231  + g[6 + ig2] * w[2 + iw2] + g[7 + ig2] * w[3 + iw2]
232  + g[12 + ig2] * w[4 + iw2] + g[13 + ig2] * w[5 + iw2];
233  field[2 * ic + 1 + iv2] =
234  g[0 + ig2] * w[1 + iw2] - g[1 + ig2] * w[0 + iw2]
235  + g[6 + ig2] * w[3 + iw2] - g[7 + ig2] * w[2 + iw2]
236  + g[12 + ig2] * w[5 + iw2] - g[13 + ig2] * w[4 + iw2];
237  }
238  }
239  }
240 }
241 
242 
243 //====================================================================
244 void Field_F::multadd_Field_Gn(int ex, const Field_G& U, int ex1,
245  const Field_F& x, int ex2, double a)
246 {
247  assert(ex <= m_Nex);
248  assert(ex1 <= U.nex());
249  assert(ex2 <= x.nex());
250  assert(U.nvol() == m_Nvol);
251  assert(x.nvol() == m_Nvol);
252 
253  double *w;
254  double *g;
255  w = const_cast<Field_F *>(&x)->ptr(0);
256  g = const_cast<Field_G *>(&U)->ptr(0);
257 
258  for (int site = 0; site < m_Nvol; ++site) {
259  int ig = NDF * (site + m_Nvol * ex1);
260  int iw = NCD2 * (site + m_Nvol * ex2);
261  int iv = NCD2 * (site + m_Nvol * ex);
262  for (int s = 0; s < ND; ++s) {
263  for (int ic = 0; ic < NC; ++ic) {
264  int ig2 = ic * NC2 + ig;
265  int iw2 = s * NC2 + iw;
266  int iv2 = s * NC2 + iv;
267  field[2 * ic + iv2] += a * (
268  g[0 + ig2] * w[0 + iw2] - g[1 + ig2] * w[1 + iw2]
269  + g[2 + ig2] * w[2 + iw2] - g[3 + ig2] * w[3 + iw2]
270  + g[4 + ig2] * w[4 + iw2] - g[5 + ig2] * w[5 + iw2]);
271  field[2 * ic + 1 + iv2] += a * (
272  g[0 + ig2] * w[1 + iw2] + g[1 + ig2] * w[0 + iw2]
273  + g[2 + ig2] * w[3 + iw2] + g[3 + ig2] * w[2 + iw2]
274  + g[4 + ig2] * w[5 + iw2] + g[5 + ig2] * w[4 + iw2]);
275  }
276  }
277  }
278 }
279 
280 
281 //====================================================================
282 void Field_F::multadd_Field_Gd(int ex, const Field_G& U, int ex1,
283  const Field_F& x, int ex2, double a)
284 {
285  assert(ex <= m_Nex);
286  assert(ex1 <= U.nex());
287  assert(ex2 <= x.nex());
288  assert(U.nvol() == m_Nvol);
289  assert(x.nvol() == m_Nvol);
290 
291  double *w;
292  double *g;
293  w = const_cast<Field_F *>(&x)->ptr(0);
294  g = const_cast<Field_G *>(&U)->ptr(0);
295 
296  for (int site = 0; site < m_Nvol; ++site) {
297  int ig = NDF * (site + m_Nvol * ex1);
298  int iw = NCD2 * (site + m_Nvol * ex2);
299  int iv = NCD2 * (site + m_Nvol * ex);
300  for (int s = 0; s < ND; ++s) {
301  for (int ic = 0; ic < NC; ++ic) {
302  int ig2 = ic * 2 + ig;
303  int iw2 = s * NC2 + iw;
304  int iv2 = s * NC2 + iv;
305  field[2 * ic + iv2] += a * (
306  g[0 + ig2] * w[0 + iw2] + g[1 + ig2] * w[1 + iw2]
307  + g[6 + ig2] * w[2 + iw2] + g[7 + ig2] * w[3 + iw2]
308  + g[12 + ig2] * w[4 + iw2] + g[13 + ig2] * w[5 + iw2]);
309  field[2 * ic + 1 + iv2] += a * (
310  g[0 + ig2] * w[1 + iw2] - g[1 + ig2] * w[0 + iw2]
311  + g[6 + ig2] * w[3 + iw2] - g[7 + ig2] * w[2 + iw2]
312  + g[12 + ig2] * w[5 + iw2] - g[13 + ig2] * w[4 + iw2]);
313  }
314  }
315  }
316 }
317 
318 
319 //====================================================================
320 //============================================================END=====