Bridge++  Ver. 1.1.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
fopr_Clover_SF.cpp
Go to the documentation of this file.
1 
14 #include "fopr_Clover_SF.h"
15 
16 #ifdef USE_PARAMETERS_FACTORY
17 #include "parameters_factory.h"
18 #endif
19 
20 using std::valarray;
21 using std::string;
22 
23 //- parameter entries
24 namespace {
25  void append_entry(Parameters& param)
26  {
27  param.Register_string("gamma_matrix_type", "NULL");
28  param.Register_double("hopping_parameter", 0.0);
29  param.Register_double("clover_coefficient", 0.0);
30  param.Register_int_vector("boundary_condition", std::valarray<int>());
31  param.Register_double_vector("phi", std::valarray<double>());
32  param.Register_double_vector("phipr", std::valarray<double>());
33 
34  param.Register_string("verbose_level", "NULL");
35  }
36 
37 
38 #ifdef USE_PARAMETERS_FACTORY
39  bool init_param = ParametersFactory::Register("Fopr.Clover_SF", append_entry);
40 #endif
41 }
42 //- end
43 
44 //- parameters class
45 Parameters_Fopr_Clover_SF::Parameters_Fopr_Clover_SF() { append_entry(*this); }
46 //- end
47 
48 //====================================================================
49 void Fopr_Clover_SF::set_parameters(const Parameters& params)
50 {
51  const string str_vlevel = params.get_string("verbose_level");
52 
53  m_vl = vout.set_verbose_level(str_vlevel);
54 
55  //- fetch and check input parameters
56  string str_gmset_type;
57  double kappa, cSW;
58  valarray<int> bc;
59  valarray<double> phi, phipr;
60 
61  int err = 0;
62  int err_optional = 0;
63  err_optional += params.fetch_string("gamma_matrix_type", str_gmset_type);
64  err += params.fetch_double("hopping_parameter", kappa);
65  err += params.fetch_double("clover_coefficient", cSW);
66  err += params.fetch_int_vector("boundary_condition", bc);
67  err += params.fetch_double_vector("phi", phi);
68  err += params.fetch_double_vector("phipr", phipr);
69 
70  if (err) {
71  vout.crucial(m_vl, "Fopr_Clover_SF: fetch error, input parameter not found.\n");
72  abort();
73  }
74 
75 
76  m_repr = str_gmset_type;
77 
78  set_parameters(kappa, cSW, bc, &phi[0], &phipr[0]);
79 }
80 
81 
82 //====================================================================
83 void Fopr_Clover_SF::set_parameters(double kappa, double cSW, valarray<int> bc,
84  double *phi, double *phipr)
85 {
86  //- print input parameters
87  vout.general(m_vl, "Parameters of Fopr_Clover_SF:\n");
88  vout.general(m_vl, " kappa = %8.4f\n", kappa);
89  vout.general(m_vl, " cSW = %8.4f\n", cSW);
90  for (int mu = 0; mu < m_Ndim; ++mu) {
91  vout.general(m_vl, " boundary[%d] = %2d\n", mu, bc[mu]);
92  }
93  vout.general(m_vl, " phi1 = %12.6f\n", phi[0]);
94  vout.general(m_vl, " phi2 = %12.6f\n", phi[1]);
95  vout.general(m_vl, " phi3 = %12.6f\n", phi[2]);
96  vout.general(m_vl, " phipr1= %12.6f\n", phipr[0]);
97  vout.general(m_vl, " phipr2= %12.6f\n", phipr[1]);
98  vout.general(m_vl, " phipr3= %12.6f\n", phipr[2]);
99 
100  //- range check
101  // NB. kappa,cSW == 0 is allowed.
102  // NB. phi,phipr == 0 is allowed.
103  assert(bc.size() == m_Ndim);
104 
105  //- store values
106  m_kappa = kappa;
107  m_cSW = cSW;
108 
109  for (int i = 0; i < 3; ++i) {
110  m_phi[i] = phi[i];
111  m_phipr[i] = phipr[i];
112  }
113 
114  for (int mu = 0; mu < m_Ndim; ++mu) {
115  m_boundary[mu] = bc[mu];
116  }
117 
118  //- propagate parameters
119  m_fopr_w->set_parameters(m_kappa, m_boundary);
120 }
121 
122 
123 //====================================================================
124 namespace {
125  inline double mult_uv_r(double *u, double *v)
126  {
127  return u[0] * v[0] - u[1] * v[1]
128  + u[2] * v[2] - u[3] * v[3]
129  + u[4] * v[4] - u[5] * v[5];
130  }
131 
132 
133  inline double mult_uv_i(double *u, double *v)
134  {
135  return u[0] * v[1] + u[1] * v[0]
136  + u[2] * v[3] + u[3] * v[2]
137  + u[4] * v[5] + u[5] * v[4];
138  }
139 }
140 
141 //====================================================================
142 void Fopr_Clover_SF::init(string repr)
143 {
144  m_Nvol = CommonParameters::Nvol();
145  m_Ndim = CommonParameters::Ndim();
146  m_Nc = CommonParameters::Nc();
147  m_Nd = CommonParameters::Nd();
148  m_NinF = 2 * m_Nc * m_Nd;
149 
150  m_U = 0;
151 
152  m_repr = repr;
153 
154  m_boundary.resize(m_Ndim);
155  m_GM.resize(m_Ndim + 1);
156  m_SG.resize(m_Ndim * m_Ndim);
157 
158  GammaMatrixSet *gmset = GammaMatrixSet::New(m_repr);
159 
160  m_GM[0] = gmset->get_GM(gmset->GAMMA1);
161  m_GM[1] = gmset->get_GM(gmset->GAMMA2);
162  m_GM[2] = gmset->get_GM(gmset->GAMMA3);
163  m_GM[3] = gmset->get_GM(gmset->GAMMA4);
164  m_GM[4] = gmset->get_GM(gmset->GAMMA5);
165 
166  m_SG[sg_index(0, 1)] = gmset->get_GM(gmset->SIGMA12);
167  m_SG[sg_index(1, 2)] = gmset->get_GM(gmset->SIGMA23);
168  m_SG[sg_index(2, 0)] = gmset->get_GM(gmset->SIGMA31);
169  m_SG[sg_index(3, 0)] = gmset->get_GM(gmset->SIGMA41);
170  m_SG[sg_index(3, 1)] = gmset->get_GM(gmset->SIGMA42);
171  m_SG[sg_index(3, 2)] = gmset->get_GM(gmset->SIGMA43);
172 
173  m_SG[sg_index(1, 0)] = m_SG[sg_index(0, 1)].mult(-1);
174  m_SG[sg_index(2, 1)] = m_SG[sg_index(1, 2)].mult(-1);
175  m_SG[sg_index(0, 2)] = m_SG[sg_index(2, 0)].mult(-1);
176  m_SG[sg_index(0, 3)] = m_SG[sg_index(3, 0)].mult(-1);
177  m_SG[sg_index(1, 3)] = m_SG[sg_index(3, 1)].mult(-1);
178  m_SG[sg_index(2, 3)] = m_SG[sg_index(3, 2)].mult(-1);
179 
180  m_SG[sg_index(0, 0)] = gmset->get_GM(gmset->UNITY);
181  m_SG[sg_index(1, 1)] = gmset->get_GM(gmset->UNITY);
182  m_SG[sg_index(2, 2)] = gmset->get_GM(gmset->UNITY);
183  m_SG[sg_index(3, 3)] = gmset->get_GM(gmset->UNITY);
184  // these 4 gamma matrices are actually not used.
185 
186 
187  // m_fopr_w = new Fopr_Wilson_SF(repr);
188  // Dirac reps. only!
189  m_fopr_w = new Fopr_Wilson_SF();
190 
191  if (m_repr == "Dirac") {
192  m_csw = &Fopr_Clover_SF::mult_csw_dirac;
193  // }else if(m_repr == "Chiral"){
194  // m_csw = &Fopr_Clover_SF::mult_csw_chiral;
195  }
196 
197  delete gmset;
198 }
199 
200 
201 //====================================================================
202 void Fopr_Clover_SF::tidyup()
203 {
204  delete m_fopr_w;
205 }
206 
207 
208 //====================================================================
209 void Fopr_Clover_SF::DdagD(Field& w, const Field& f)
210 {
211  assert(f.nex() == 1);
212  Field w2(f.nin(), f.nvol(), 1);
213 
214  D(w2, f);
215  mult_gm5(w, w2);
216  D(w2, w);
217  mult_gm5(w, w2);
218 }
219 
220 
221 //====================================================================
222 void Fopr_Clover_SF::Ddag(Field& w, const Field& f)
223 {
224  assert(f.nex() == 1);
225  Field w2(f.nin(), f.nvol(), 1);
226 
227  mult_gm5(w, f);
228  D(w2, w);
229  mult_gm5(w, w2);
230 }
231 
232 
233 //====================================================================
234 void Fopr_Clover_SF::H(Field& w, const Field& f)
235 {
236  assert(f.nex() == 1);
237  Field w2(f.nin(), f.nvol(), 1);
238 
239  D(w2, f);
240  mult_gm5(w, w2);
241 }
242 
243 
244 //====================================================================
245 void Fopr_Clover_SF::D(Field& w, const Field& f)
246 {
247  assert(f.nex() == 1);
248  int Nvol = f.nvol();
249 
250  Field_F w2(Nvol, 1);
251 
252  // Field_F w3(f);
253  // setzero.set_boundary_zero(w3);
254  // m_fopr_w->D(w,w3);
255  // mult_csw(w2,w3);
256 
257  m_fopr_w->D(w, f);
258  mult_csw(w2, (Field_F)f);
259  w -= (Field)w2;
260  setzero.set_boundary_zero(w);
261 }
262 
263 
264 //====================================================================
265 const Field_F Fopr_Clover_SF::DdagD(const Field_F& f)
266 {
267  Field_F w(f.nvol(), f.nex());
268 
269  DdagD(w, f);
270 
271  return w;
272 }
273 
274 
275 //====================================================================
276 const Field_F Fopr_Clover_SF::Ddag(const Field_F& f)
277 {
278  Field_F w(f.nvol(), f.nex());
279 
280  Ddag(w, f);
281 
282  return w;
283 }
284 
285 
286 //====================================================================
287 const Field_F Fopr_Clover_SF::H(const Field_F& f)
288 {
289  Field_F w(f.nvol(), f.nex());
290 
291  H(w, f);
292 
293  return w;
294 }
295 
296 
297 //====================================================================
298 const Field_F Fopr_Clover_SF::D(const Field_F& f)
299 {
300  Field_F w(f.nvol(), f.nex());
301 
302  D(w, f);
303 
304  return w;
305 }
306 
307 
308 //====================================================================
309 void Fopr_Clover_SF::mult_isigma(Field_F& v, const Field_F& w,
310  const int mu, const int nu)
311 {
312  assert(mu != nu);
313  v.mult_iGM(m_SG[sg_index(mu, nu)], w);
314 }
315 
316 
317 //====================================================================
318 void Fopr_Clover_SF::mult_csw(Field_F& v, const Field_F& w)
319 {
320  (this->*m_csw)(v, w);
321 }
322 
323 
324 //====================================================================
325 void Fopr_Clover_SF::mult_csw_chiral(Field_F& v, const Field_F& w)
326 {
327  assert(w.nex() == 1);
328 
329  int Nc = CommonParameters::Nc();
330  int Nvc = 2 * Nc;
331  int Ndf = 2 * Nc * Nc;
332  int Nd = CommonParameters::Nd();
333  int Nvol = w.nvol();
334 
335  int id1 = 0;
336  int id2 = Nvc;
337  int id3 = Nvc * 2;
338  int id4 = Nvc * 3;
339 
340  double *w2 = const_cast<Field_F *>(&w)->ptr(0);
341  double *v2 = v.ptr(0);
342 
343  double *Bx = m_Bx.ptr(0);
344  double *By = m_By.ptr(0);
345  double *Bz = m_Bz.ptr(0);
346  double *Ex = m_Ex.ptr(0);
347  double *Ey = m_Ey.ptr(0);
348  double *Ez = m_Ez.ptr(0);
349 
350  v = 0.0;
351 
352  for (int site = 0; site < Nvol; ++site) {
353  int iv = Nvc * Nd * site;
354  int ig = Ndf * site;
355 
356  for (int ic = 0; ic < Nc; ++ic) {
357  int icr = 2 * ic;
358  int ici = icr + 1;
359  int icg = ic * Nvc + ig;
360 
361  // isigma_23 * Bx
362  v2[icr + id1 + iv] -= mult_uv_i(&Bx[icg], &w2[id2 + iv]);
363  v2[ici + id1 + iv] += mult_uv_r(&Bx[icg], &w2[id2 + iv]);
364  v2[icr + id2 + iv] -= mult_uv_i(&Bx[icg], &w2[id1 + iv]);
365  v2[ici + id2 + iv] += mult_uv_r(&Bx[icg], &w2[id1 + iv]);
366 
367  v2[icr + id3 + iv] -= mult_uv_i(&Bx[icg], &w2[id4 + iv]);
368  v2[ici + id3 + iv] += mult_uv_r(&Bx[icg], &w2[id4 + iv]);
369  v2[icr + id4 + iv] -= mult_uv_i(&Bx[icg], &w2[id3 + iv]);
370  v2[ici + id4 + iv] += mult_uv_r(&Bx[icg], &w2[id3 + iv]);
371 
372  // isigma_31 * By
373  v2[icr + id1 + iv] += mult_uv_r(&By[icg], &w2[id2 + iv]);
374  v2[ici + id1 + iv] += mult_uv_i(&By[icg], &w2[id2 + iv]);
375  v2[icr + id2 + iv] -= mult_uv_r(&By[icg], &w2[id1 + iv]);
376  v2[ici + id2 + iv] -= mult_uv_i(&By[icg], &w2[id1 + iv]);
377 
378  v2[icr + id3 + iv] += mult_uv_r(&By[icg], &w2[id4 + iv]);
379  v2[ici + id3 + iv] += mult_uv_i(&By[icg], &w2[id4 + iv]);
380  v2[icr + id4 + iv] -= mult_uv_r(&By[icg], &w2[id3 + iv]);
381  v2[ici + id4 + iv] -= mult_uv_i(&By[icg], &w2[id3 + iv]);
382 
383  // isigma_12 * Bz
384  v2[icr + id1 + iv] -= mult_uv_i(&Bz[icg], &w2[id1 + iv]);
385  v2[ici + id1 + iv] += mult_uv_r(&Bz[icg], &w2[id1 + iv]);
386  v2[icr + id2 + iv] += mult_uv_i(&Bz[icg], &w2[id2 + iv]);
387  v2[ici + id2 + iv] -= mult_uv_r(&Bz[icg], &w2[id2 + iv]);
388 
389  v2[icr + id3 + iv] -= mult_uv_i(&Bz[icg], &w2[id3 + iv]);
390  v2[ici + id3 + iv] += mult_uv_r(&Bz[icg], &w2[id3 + iv]);
391  v2[icr + id4 + iv] += mult_uv_i(&Bz[icg], &w2[id4 + iv]);
392  v2[ici + id4 + iv] -= mult_uv_r(&Bz[icg], &w2[id4 + iv]);
393 
394  // isigma_41 * Ex
395  v2[icr + id1 + iv] += mult_uv_i(&Ex[icg], &w2[id2 + iv]);
396  v2[ici + id1 + iv] -= mult_uv_r(&Ex[icg], &w2[id2 + iv]);
397  v2[icr + id2 + iv] += mult_uv_i(&Ex[icg], &w2[id1 + iv]);
398  v2[ici + id2 + iv] -= mult_uv_r(&Ex[icg], &w2[id1 + iv]);
399 
400  v2[icr + id3 + iv] -= mult_uv_i(&Ex[icg], &w2[id4 + iv]);
401  v2[ici + id3 + iv] += mult_uv_r(&Ex[icg], &w2[id4 + iv]);
402  v2[icr + id4 + iv] -= mult_uv_i(&Ex[icg], &w2[id3 + iv]);
403  v2[ici + id4 + iv] += mult_uv_r(&Ex[icg], &w2[id3 + iv]);
404 
405  // isigma_42 * Ey
406  v2[icr + id1 + iv] -= mult_uv_r(&Ey[icg], &w2[id2 + iv]);
407  v2[ici + id1 + iv] -= mult_uv_i(&Ey[icg], &w2[id2 + iv]);
408  v2[icr + id2 + iv] += mult_uv_r(&Ey[icg], &w2[id1 + iv]);
409  v2[ici + id2 + iv] += mult_uv_i(&Ey[icg], &w2[id1 + iv]);
410 
411  v2[icr + id3 + iv] += mult_uv_r(&Ey[icg], &w2[id4 + iv]);
412  v2[ici + id3 + iv] += mult_uv_i(&Ey[icg], &w2[id4 + iv]);
413  v2[icr + id4 + iv] -= mult_uv_r(&Ey[icg], &w2[id3 + iv]);
414  v2[ici + id4 + iv] -= mult_uv_i(&Ey[icg], &w2[id3 + iv]);
415 
416  // isigma_43 * Ez
417  v2[icr + id1 + iv] += mult_uv_i(&Ez[icg], &w2[id1 + iv]);
418  v2[ici + id1 + iv] -= mult_uv_r(&Ez[icg], &w2[id1 + iv]);
419  v2[icr + id2 + iv] -= mult_uv_i(&Ez[icg], &w2[id2 + iv]);
420  v2[ici + id2 + iv] += mult_uv_r(&Ez[icg], &w2[id2 + iv]);
421 
422  v2[icr + id3 + iv] -= mult_uv_i(&Ez[icg], &w2[id3 + iv]);
423  v2[ici + id3 + iv] += mult_uv_r(&Ez[icg], &w2[id3 + iv]);
424  v2[icr + id4 + iv] += mult_uv_i(&Ez[icg], &w2[id4 + iv]);
425  v2[ici + id4 + iv] -= mult_uv_r(&Ez[icg], &w2[id4 + iv]);
426  }
427  }
428 
429  v *= m_kappa * m_cSW;
430 }
431 
432 
433 //====================================================================
434 void Fopr_Clover_SF::mult_csw_dirac(Field_F& v, const Field_F& w)
435 {
436  assert(w.nex() == 1);
437 
438  int Nc = CommonParameters::Nc();
439  int Nvc = 2 * Nc;
440  int Ndf = 2 * Nc * Nc;
441  int Nd = CommonParameters::Nd();
442  int Nvol = w.nvol();
443 
444  int id1 = 0;
445  int id2 = Nvc;
446  int id3 = Nvc * 2;
447  int id4 = Nvc * 3;
448 
449  double *w2 = const_cast<Field_F *>(&w)->ptr(0);
450  double *v2 = v.ptr(0);
451 
452  double *Bx = m_Bx.ptr(0);
453  double *By = m_By.ptr(0);
454  double *Bz = m_Bz.ptr(0);
455  double *Ex = m_Ex.ptr(0);
456  double *Ey = m_Ey.ptr(0);
457  double *Ez = m_Ez.ptr(0);
458 
459  v = 0.0;
460 
461  for (int site = 0; site < Nvol; ++site) {
462  int iv = Nvc * Nd * site;
463  int ig = Ndf * site;
464 
465  for (int ic = 0; ic < Nc; ++ic) {
466  int icr = 2 * ic;
467  int ici = icr + 1;
468  int icg = ic * Nvc + ig;
469 
470  // isigma_23 * Bx
471  v2[icr + id1 + iv] -= mult_uv_i(&Bx[icg], &w2[id2 + iv]);
472  v2[ici + id1 + iv] += mult_uv_r(&Bx[icg], &w2[id2 + iv]);
473  v2[icr + id2 + iv] -= mult_uv_i(&Bx[icg], &w2[id1 + iv]);
474  v2[ici + id2 + iv] += mult_uv_r(&Bx[icg], &w2[id1 + iv]);
475 
476  v2[icr + id3 + iv] -= mult_uv_i(&Bx[icg], &w2[id4 + iv]);
477  v2[ici + id3 + iv] += mult_uv_r(&Bx[icg], &w2[id4 + iv]);
478  v2[icr + id4 + iv] -= mult_uv_i(&Bx[icg], &w2[id3 + iv]);
479  v2[ici + id4 + iv] += mult_uv_r(&Bx[icg], &w2[id3 + iv]);
480 
481  // isigma_31 * By
482  v2[icr + id1 + iv] += mult_uv_r(&By[icg], &w2[id2 + iv]);
483  v2[ici + id1 + iv] += mult_uv_i(&By[icg], &w2[id2 + iv]);
484  v2[icr + id2 + iv] -= mult_uv_r(&By[icg], &w2[id1 + iv]);
485  v2[ici + id2 + iv] -= mult_uv_i(&By[icg], &w2[id1 + iv]);
486 
487  v2[icr + id3 + iv] += mult_uv_r(&By[icg], &w2[id4 + iv]);
488  v2[ici + id3 + iv] += mult_uv_i(&By[icg], &w2[id4 + iv]);
489  v2[icr + id4 + iv] -= mult_uv_r(&By[icg], &w2[id3 + iv]);
490  v2[ici + id4 + iv] -= mult_uv_i(&By[icg], &w2[id3 + iv]);
491 
492  // isigma_12 * Bz
493  v2[icr + id1 + iv] -= mult_uv_i(&Bz[icg], &w2[id1 + iv]);
494  v2[ici + id1 + iv] += mult_uv_r(&Bz[icg], &w2[id1 + iv]);
495  v2[icr + id2 + iv] += mult_uv_i(&Bz[icg], &w2[id2 + iv]);
496  v2[ici + id2 + iv] -= mult_uv_r(&Bz[icg], &w2[id2 + iv]);
497 
498  v2[icr + id3 + iv] -= mult_uv_i(&Bz[icg], &w2[id3 + iv]);
499  v2[ici + id3 + iv] += mult_uv_r(&Bz[icg], &w2[id3 + iv]);
500  v2[icr + id4 + iv] += mult_uv_i(&Bz[icg], &w2[id4 + iv]);
501  v2[ici + id4 + iv] -= mult_uv_r(&Bz[icg], &w2[id4 + iv]);
502 
503  // isigma_41 * Ex
504  v2[icr + id1 + iv] += mult_uv_i(&Ex[icg], &w2[id4 + iv]);
505  v2[ici + id1 + iv] -= mult_uv_r(&Ex[icg], &w2[id4 + iv]);
506  v2[icr + id2 + iv] += mult_uv_i(&Ex[icg], &w2[id3 + iv]);
507  v2[ici + id2 + iv] -= mult_uv_r(&Ex[icg], &w2[id3 + iv]);
508 
509  v2[icr + id3 + iv] += mult_uv_i(&Ex[icg], &w2[id2 + iv]);
510  v2[ici + id3 + iv] -= mult_uv_r(&Ex[icg], &w2[id2 + iv]);
511  v2[icr + id4 + iv] += mult_uv_i(&Ex[icg], &w2[id1 + iv]);
512  v2[ici + id4 + iv] -= mult_uv_r(&Ex[icg], &w2[id1 + iv]);
513 
514  // isigma_42 * Ey
515  v2[icr + id1 + iv] -= mult_uv_r(&Ey[icg], &w2[id4 + iv]);
516  v2[ici + id1 + iv] -= mult_uv_i(&Ey[icg], &w2[id4 + iv]);
517  v2[icr + id2 + iv] += mult_uv_r(&Ey[icg], &w2[id3 + iv]);
518  v2[ici + id2 + iv] += mult_uv_i(&Ey[icg], &w2[id3 + iv]);
519 
520  v2[icr + id3 + iv] -= mult_uv_r(&Ey[icg], &w2[id2 + iv]);
521  v2[ici + id3 + iv] -= mult_uv_i(&Ey[icg], &w2[id2 + iv]);
522  v2[icr + id4 + iv] += mult_uv_r(&Ey[icg], &w2[id1 + iv]);
523  v2[ici + id4 + iv] += mult_uv_i(&Ey[icg], &w2[id1 + iv]);
524 
525  // isigma_43 * Ez
526  v2[icr + id1 + iv] += mult_uv_i(&Ez[icg], &w2[id3 + iv]);
527  v2[ici + id1 + iv] -= mult_uv_r(&Ez[icg], &w2[id3 + iv]);
528  v2[icr + id2 + iv] -= mult_uv_i(&Ez[icg], &w2[id4 + iv]);
529  v2[ici + id2 + iv] += mult_uv_r(&Ez[icg], &w2[id4 + iv]);
530 
531  v2[icr + id3 + iv] += mult_uv_i(&Ez[icg], &w2[id1 + iv]);
532  v2[ici + id3 + iv] -= mult_uv_r(&Ez[icg], &w2[id1 + iv]);
533  v2[icr + id4 + iv] -= mult_uv_i(&Ez[icg], &w2[id2 + iv]);
534  v2[ici + id4 + iv] += mult_uv_r(&Ez[icg], &w2[id2 + iv]);
535  }
536  }
537 
538  v *= m_kappa * m_cSW;
539 }
540 
541 
542 //====================================================================
543 void Fopr_Clover_SF::set_csw()
544 {
545  set_fieldstrength(m_Bx, 1, 2);
546  set_fieldstrength(m_By, 2, 0);
547  set_fieldstrength(m_Bz, 0, 1);
548  set_fieldstrength(m_Ex, 3, 0);
549  set_fieldstrength(m_Ey, 3, 1);
550  set_fieldstrength(m_Ez, 3, 2);
551 }
552 
553 
561 //====================================================================
562 void Fopr_Clover_SF::set_fieldstrength(Field_G& Fst,
563  const int mu, const int nu)
564 {
565  int Nvol = CommonParameters::Nvol();
566 
567  Staples_SF staple;
568 
569  staple.set_parameters(m_phi, m_phipr);
570 
571  Field_G Cup(Nvol, 1), Cdn(Nvol, 1);
572  Field_G Umu(Nvol, 1);
573  Field_G v(Nvol, 1), v2(Nvol, 1);
574 
575  Cup = staple.upper(*m_U, mu, nu);
576  Cdn = staple.lower(*m_U, mu, nu);
577  Umu.setpart_ex(0, *m_U, mu);
578 
579  Fst.mult_Field_Gnd(0, Umu, 0, Cup, 0);
580  Fst.multadd_Field_Gnd(0, Umu, 0, Cdn, 0, -1.0);
581 
582  v.mult_Field_Gdn(0, Cup, 0, Umu, 0);
583  v.multadd_Field_Gdn(0, Cdn, 0, Umu, 0, -1.0);
584 
585  m_shift.forward(v2, v, mu);
586 
587  Fst += v2;
588 
589  Fst.ah_Field_G(0);
590  Fst *= 0.25;
591 }
592 
593 
594 //====================================================================
595 //============================================================END=====