Bridge++  Ver. 2.0.2
fopr_NonRelativistic.cpp
Go to the documentation of this file.
1 
10 #include "lib/Fopr/fopr_NonRelativistic.h"
11 
12 #include "lib/ResourceManager/threadManager.h"
13 #include "lib/Measurements/Gauge/fieldStrength.h"
14 
15 #ifdef USE_FACTORY_AUTOREGISTER
16 namespace {
17  bool init = Fopr_NonRelativistic::register_factory();
18 }
19 #endif
20 
21 const std::string Fopr_NonRelativistic::class_name
22  = "Fopr_NonRelativistic";
23 
24 //====================================================================
25 namespace { // inline functions
26  inline void set_threadtask(int& ith, int& nth, int& is, int& ns,
27  const int size)
28  {
29  nth = ThreadManager_OpenMP::get_num_threads();
30  ith = ThreadManager_OpenMP::get_thread_id();
31  is = size * ith / nth;
32  ns = size * (ith + 1) / nth;
33  }
34 
35 
36  double mult_Gn_r(const double *g, const double *w, int Nc)
37  {
38  return g[0] * w[0] - g[1] * w[1]
39  + g[2] * w[2] - g[3] * w[3]
40  + g[4] * w[4] - g[5] * w[5];
41  }
42 
43 
44  double mult_Gn_i(const double *g, const double *w, int Nc)
45  {
46  return g[0] * w[1] + g[1] * w[0]
47  + g[2] * w[3] + g[3] * w[2]
48  + g[4] * w[5] + g[5] * w[4];
49  }
50 
51 
52  double mult_Gd_r(const double *g, const double *w, int Nc)
53  {
54  return g[0] * w[0] + g[1] * w[1]
55  + g[6] * w[2] + g[7] * w[3]
56  + g[12] * w[4] + g[13] * w[5];
57  }
58 
59 
60  double mult_Gd_i(const double *g, const double *w, int Nc)
61  {
62  return g[0] * w[1] - g[1] * w[0]
63  + g[6] * w[3] - g[7] * w[2]
64  + g[12] * w[5] - g[13] * w[4];
65  }
66 } // end of nameless namespace
67 
68 //====================================================================
69 void Fopr_NonRelativistic::init(const Parameters& params)
70 {
71  ThreadManager::assert_single_thread(class_name);
72 
73  m_vl = CommonParameters::Vlevel();
74  vout.general(m_vl, "%s: construction\n", class_name.c_str());
75  vout.increase_indent();
76 
77  m_Nc = CommonParameters::Nc();
78  m_Nvc = 2 * m_Nc;
79  m_Nd = CommonParameters::Nd();
80  m_Nd2 = m_Nd / 2;
81  m_Nx = CommonParameters::Nx();
82  m_Ny = CommonParameters::Ny();
83  m_Nz = CommonParameters::Nz();
84  m_Nt = CommonParameters::Nt();
85  m_Ltime = CommonParameters::Lt();
86  m_Nvol = CommonParameters::Nvol();
87  m_Ndim = CommonParameters::Ndim();
88  m_Nspc = m_Nx * m_Ny * m_Nz;
89 
90  m_boundary.resize(m_Ndim);
91 
92  m_U.reset(m_Nvol, m_Ndim);
93 
94  m_repr = "Dirac"; // only Dirac representation is available.
95 
96  // number of correction terms in this implementation
97  m_num_correct = 6;
98 
99  set_parameters(params);
100 
101  m_Fstr.resize(6);
102  for (int j = 0; j < 6; ++j) {
103  m_Fstr[j].reset(m_Nvol, 1);
104  }
105 
106  int NinF2 = 2 * m_Nc * m_Nd2;
107 
108  m_Xt.resize(m_Nt);
109  for (int it = 0; it < m_Nt; ++it) {
110  m_Xt[it].reset(NinF2, m_Nspc, 1);
111  m_Xt[it].set(0.0);
112  }
113 
114  m_Zt1.reset(NinF2, m_Nspc, 1);
115  m_Zt2.reset(NinF2, m_Nspc, 1);
116  m_Yt1.reset(NinF2, m_Nspc, 1);
117  m_Yt2.reset(NinF2, m_Nspc, 1);
118  m_Yt3.reset(NinF2, m_Nspc, 1);
119  m_Wt1.reset(NinF2, m_Nspc, 1);
120  m_Wt2.reset(NinF2, m_Nspc, 1);
121  m_Wt3.reset(NinF2, m_Nspc, 1);
122 
123  int Nbuf_x = NinF2 * m_Nspc / m_Nx;
124  buf1_x.resize(Nbuf_x);
125  buf2_x.resize(Nbuf_x);
126 
127  int Nbuf_y = NinF2 * m_Nspc / m_Ny;
128  buf1_y.resize(Nbuf_y);
129  buf2_y.resize(Nbuf_y);
130 
131  int Nbuf_z = NinF2 * m_Nspc / m_Nz;
132  buf1_z.resize(Nbuf_z);
133  buf2_z.resize(Nbuf_z);
134 
135  vout.decrease_indent();
136  vout.general(m_vl, "%s: construction finished.\n",
137  class_name.c_str());
138 }
139 
140 
141 //====================================================================
142 void Fopr_NonRelativistic::tidyup()
143 {
144  // do nothing.
145 }
146 
147 
148 //====================================================================
149 void Fopr_NonRelativistic::set_parameters(const Parameters& params)
150 {
151  const std::string str_vlevel = params.get_string("verbose_level");
152 
153  vout.general(m_vl, "%s: parameter setup:\n", class_name.c_str());
154 
155  m_vl = vout.set_verbose_level(str_vlevel);
156 
157  //- fetch and check input parameters
158  double MQ;
159  int nstab;
160  double u0;
161  std::vector<int> bc;
162  std::vector<double> coeff;
163  std::string evolution_type;
164  std::string correction_terms;
165 
166  int err = 0;
167  err += params.fetch_double("quark_mass", MQ);
168  err += params.fetch_int("stabilization_parameter", nstab);
169  err += params.fetch_double("mean_field_value", u0);
170  err += params.fetch_int_vector("boundary_condition", bc);
171  err += params.fetch_double_vector("coefficients", coeff);
172  err += params.fetch_string("evolution_type", evolution_type);
173  err += params.fetch_string("correction_terms", correction_terms);
174 
175  if (err) {
176  vout.crucial(m_vl, "Error at %s: input parameter not found.\n",
177  class_name.c_str());
178  exit(EXIT_FAILURE);
179  }
180 
181  set_parameters(MQ, nstab, u0, coeff, bc,
182  evolution_type, correction_terms);
183 
184  vout.general(m_vl, "%s: parameter setup finished.\n",
185  class_name.c_str());
186 }
187 
188 
189 //====================================================================
190 void Fopr_NonRelativistic::set_parameters(
191  const double MQ,
192  const int nstab,
193  const double u0,
194  const std::vector<double> coeff,
195  const std::vector<int> bc,
196  const std::string evolution_type,
197  const std::string correction_terms)
198 {
199  assert(bc.size() == m_Ndim);
200 
201  m_MQ = MQ;
202  m_u0 = u0;
203  m_nstab = nstab;
204  m_coeff = coeff;
205  m_boundary = bc;
206  m_evolution_type = evolution_type;
207  m_correction_terms = correction_terms;
208 
209  if (m_num_correct != m_coeff.size()) {
210  vout.crucial(m_vl, "Error at %s: inconsistent number of coefficients.\n",
211  class_name.c_str());
212  exit(EXIT_FAILURE);
213  }
214 
215  vout.general(m_vl, "%s: parameters:\n", class_name.c_str());
216  vout.general(m_vl, " quark mass = %12.8f\n", m_MQ);
217  vout.general(m_vl, " stabilization parameter = %2d\n", m_nstab);
218  vout.general(m_vl, " mean-field value = %12.8f\n", m_u0);
219  vout.general(m_vl, " evolution_type = %s\n",
220  m_evolution_type.c_str());
221  vout.general(m_vl, " correction_terms = %s\n",
222  m_correction_terms.c_str());
223  for (int i = 0; i < m_num_correct; ++i) {
224  vout.general(m_vl, " coeff[%d] = %12.8f\n", i, m_coeff[i]);
225  }
226  for (int mu = 0; mu < m_Ndim; ++mu) {
227  vout.general(m_vl, " boundary[%d] = %2d\n", mu, m_boundary[mu]);
228  }
229 }
230 
231 
232 //====================================================================
233 void Fopr_NonRelativistic::get_parameters(Parameters& params) const
234 {
235  params.set_double("quark_mass", m_MQ);
236  params.set_int("stabilization_parameter", m_nstab);
237  params.set_double("mean_field_value", m_u0);
238  params.set_int_vector("boundary_condition", m_boundary);
239  params.set_double_vector("coefficients", m_coeff);
240  params.set_string("evolution_type", m_evolution_type);
241  params.set_string("correction_terms", m_correction_terms);
242  params.set_string("verbose_level", vout.get_verbose_level(m_vl));
243 }
244 
245 
246 //====================================================================
247 void Fopr_NonRelativistic::set_config(Field *U)
248 {
249  copy(m_U, *U);
250  scal(m_U, 1.0 / m_u0); // applying mean-field improvement.
251 
252  // Note: the following definition of field strength is not traceless.
253  FieldStrength fstr;
254  fstr.construct_Fmunu_1x1(m_Fstr[0], 2, 1, m_U);
255  fstr.construct_Fmunu_1x1(m_Fstr[1], 0, 2, m_U);
256  fstr.construct_Fmunu_1x1(m_Fstr[2], 1, 0, m_U);
257  fstr.construct_Fmunu_1x1(m_Fstr[3], 0, 3, m_U);
258  fstr.construct_Fmunu_1x1(m_Fstr[4], 1, 3, m_U);
259  fstr.construct_Fmunu_1x1(m_Fstr[5], 2, 3, m_U);
260 }
261 
262 
263 //====================================================================
264 void Fopr_NonRelativistic::set_mode(std::string mode)
265 {
266 #pragma omp barrier
267 
268  int ith = ThreadManager::get_thread_id();
269  if (ith == 0) m_mode = mode;
270 
271 #pragma omp barrier
272 }
273 
274 
275 //====================================================================
276 void Fopr_NonRelativistic::mult(Field& v, const Field& w)
277 {
278  if (m_mode == "Evolve") {
279  evolve(v, w);
280  } else if (m_mode == "Rotation") {
281  rotation(v, w, 1);
282  } else {
283  vout.crucial(m_vl, "Error at %s: undefined mode = %s.\n",
284  class_name.c_str(), m_mode.c_str());
285  exit(EXIT_FAILURE);
286  }
287 }
288 
289 
290 //====================================================================
291 void Fopr_NonRelativistic::mult_dag(Field& v, const Field& w)
292 {
293  if (m_mode == "Rotation") {
294  rotation(v, w, -1);
295  } else {
296  vout.crucial(m_vl, "Error at %s: undefined mode = %s.\n",
297  class_name.c_str(), m_mode.c_str());
298  exit(EXIT_FAILURE);
299  }
300 }
301 
302 
303 //====================================================================
304 void Fopr_NonRelativistic::mult(Field& v, const Field& w,
305  const std::string mode)
306 {
307  if (mode == "Evolve") {
308  evolve(v, w);
309  } else if (mode == "Rotation") {
310  rotation(v, w, 1);
311  } else {
312  vout.crucial(m_vl, "Error at %s: undefined mode = %s.\n",
313  class_name.c_str(), m_mode.c_str());
314  exit(EXIT_FAILURE);
315  }
316 }
317 
318 
319 //====================================================================
320 void Fopr_NonRelativistic::mult_dag(Field& v, const Field& w,
321  const std::string mode)
322 {
323  if (mode == "Rotation") {
324  rotation(v, w, -1);
325  } else {
326  vout.crucial(m_vl, "Error at %s: undefined mode = %s.\n",
327  class_name.c_str(), m_mode.c_str());
328  exit(EXIT_FAILURE);
329  }
330 }
331 
332 
333 //====================================================================
334 void Fopr_NonRelativistic::mult_gm5(Field& v, const Field& w)
335 {
336  int ith, nth, is, ns;
337  set_threadtask(ith, nth, is, ns, m_Nvol);
338 
339 #pragma omp barrier
340 
341  for (int site = is; site < ns; ++site) {
342  for (int id = 0; id < m_Nd2; ++id) {
343  for (int ivc = 0; ivc < m_Nvc; ++ivc) {
344  int ivcd = ivc + m_Nvc * id;
345  int ivcd2 = ivc + m_Nvc * (2 + id);
346  double w1 = w.cmp(ivcd, site, 0);
347  double w2 = w.cmp(ivcd2, site, 0);
348  v.set(ivcd2, site, 0, w1);
349  v.set(ivcd, site, 0, w2);
350  }
351  }
352  }
353 #pragma omp barrier
354 }
355 
356 
357 //====================================================================
358 void Fopr_NonRelativistic::mult_up(const int mu, Field& v, const Field& w)
359 {
360  vout.crucial(m_vl, "Error at %s: mult_up undefined.\n",
361  class_name.c_str());
362  exit(EXIT_FAILURE);
363 }
364 
365 
366 //====================================================================
367 void Fopr_NonRelativistic::mult_dn(const int mu, Field& v, const Field& w)
368 {
369  vout.crucial(m_vl, "Error at %s: mult_dn undefined.\n",
370  class_name.c_str());
371  exit(EXIT_FAILURE);
372 }
373 
374 
375 //====================================================================
376 void Fopr_NonRelativistic::evolve(Field& v, const Field& w)
377 {
378  vout.detailed(m_vl, "%s: evolution of heavy quark field start.\n",
379  class_name.c_str());
380 
381  int IPEt = Communicator::ipe(3);
382 
383  // setting source time slice and source vector
384  set_source(m_Yt1, w);
385 
386  int it_src = m_itime_src % m_Nt;
387  int ipet_src = m_itime_src / m_Nt;
388  if (ipet_src == IPEt) copy(m_Xt[it_src], m_Yt1);
389 #pragma omp barrier
390 
391  // time evolution
392  for (int itime = 1; itime < m_Ltime - 1; ++itime) {
393  int itime1 = (itime - 1 + m_itime_src) % m_Ltime;
394  int it1 = itime1 % m_Nt;
395  int itime2 = (itime + m_itime_src) % m_Ltime;
396  int it2 = itime2 % m_Nt;
397 
398  evolve_impl(m_Xt[it2], m_Xt[it1], itime2);
399 
400  if (itime2 / m_Nt == Communicator::ipe(3)) {
401  copy(m_Zt1, m_Xt[it2]);
402  } else {
403  m_Zt1.set(0.0);
404  }
405 #pragma omp barrier
406 
407  double xnorm2 = m_Zt1.norm2();
408  vout.detailed(m_vl, " itime = %d norm2 = %18.14e\n",
409  itime2, xnorm2);
410  }
411 
412  // set 4D spinor field
413 
414  int ith, nth, is, ns;
415  set_threadtask(ith, nth, is, ns, m_Nspc);
416 
417  for (int it = 0; it < m_Nt; ++it) {
418  for (int ispc = is; ispc < ns; ++ispc) {
419  for (int id = 0; id < m_Nd2; ++id) {
420  for (int ivc = 0; ivc < m_Nvc; ++ivc) {
421  int ivcd = ivc + m_Nvc * id;
422  int ivcd2 = ivc + m_Nvc * (2 + id);
423  double Xt1 = m_Xt[it].cmp(ivcd, ispc, 0);
424  v.set(ivcd, ispc + m_Nspc * it, 0, Xt1);
425  v.set(ivcd2, ispc + m_Nspc * it, 0, 0.0);
426  }
427  }
428  }
429  }
430 #pragma omp barrier
431 }
432 
433 
434 //====================================================================
435 void Fopr_NonRelativistic::rotation(Field& v, const Field& w,
436  const int jd)
437 {
438 #pragma omp barrier
439 
440  int ith, nth, is, ns;
441  set_threadtask(ith, nth, is, ns, m_Nspc);
442 
443  for (int it = 0; it < m_Nt; ++it) {
444  for (int ispc = is; ispc < ns; ++ispc) {
445  for (int id = 0; id < m_Nd2; ++id) {
446  for (int ivc = 0; ivc < m_Nvc; ++ivc) {
447  int ivcd = ivc + m_Nvc * id;
448  int ivcd2 = ivc + m_Nvc * (2 + id);
449  double w1 = w.cmp(ivcd, ispc + m_Nspc * it, 0);
450  double w2 = w.cmp(ivcd2, ispc + m_Nspc * it, 0);
451  m_Zt1.set(ivcd, ispc, 0, w1);
452  m_Zt2.set(ivcd, ispc, 0, w2);
453  }
454  }
455  }
456 #pragma omp barrier
457 
458  int IPEt = Communicator::ipe(3);
459  int itime = it + m_Nt * IPEt;
460 
461  // upper components of v:
462  copy(m_Xt[0], m_Zt1);
463 #pragma omp barrier
464 
465  add_R2(m_Xt[0], m_Zt2, -jd, itime);
466 
467  if (m_correction_terms == "next-leading") {
468  add_R3(m_Xt[0], m_Zt1, itime);
469  add_R4(m_Xt[0], m_Zt1, itime);
470  add_R5(m_Xt[0], m_Zt2, itime);
471  }
472 
473  for (int ispc = is; ispc < ns; ++ispc) {
474  for (int id = 0; id < m_Nd2; ++id) {
475  for (int ivc = 0; ivc < m_Nvc; ++ivc) {
476  int ivcd = ivc + m_Nvc * id;
477  double v1 = m_Xt[0].cmp(ivcd, ispc, 0);
478  v.set(ivcd, ispc + m_Nspc * it, 0, v1);
479  }
480  }
481  }
482 #pragma omp barrier
483 
484  // lower components of v:
485  copy(m_Xt[0], m_Zt2);
486 #pragma omp barrier
487 
488  add_R2(m_Xt[0], m_Zt1, jd, itime);
489 
490  if (m_correction_terms == "next-leading") {
491  add_R3(m_Xt[0], m_Zt2, itime);
492  add_R4(m_Xt[0], m_Zt2, itime);
493  add_R5(m_Xt[0], m_Zt1, itime);
494  }
495 
496  for (int ispc = is; ispc < ns; ++ispc) {
497  for (int id = 0; id < m_Nd2; ++id) {
498  for (int ivc = 0; ivc < m_Nvc; ++ivc) {
499  int ivcd = ivc + m_Nvc * id;
500  int ivcd2 = ivc + m_Nvc * (2 + id);
501  double v2 = m_Xt[0].cmp(ivcd, ispc, 0);
502  v.set(ivcd2, ispc + m_Nspc * it, 0, v2);
503  }
504  }
505  }
506 #pragma omp barrier
507  }
508 }
509 
510 
511 //====================================================================
512 void Fopr_NonRelativistic::add_R2(Field& Xt, const Field& Yt,
513  const int jd, const int itime)
514 {
515  double fac = -1.0 / (2.0 * m_MQ);
516 
517  m_Yt3.set(0.0);
518 
519  calc_Delta1(m_Yt1, Yt, 0, itime);
520  mult_sigma1(m_Yt2, m_Yt1);
521  axpy(m_Yt3, fac, m_Yt2);
522 #pragma omp barrier
523 
524  calc_Delta1(m_Yt1, Yt, 1, itime);
525  mult_sigma2(m_Yt2, m_Yt1);
526  axpy(m_Yt3, fac, m_Yt2);
527 #pragma omp barrier
528 
529  calc_Delta1(m_Yt1, Yt, 2, itime);
530  mult_sigma3(m_Yt2, m_Yt1);
531  axpy(m_Yt3, fac, m_Yt2);
532 #pragma omp barrier
533 
534  m_Yt3.xI();
535  axpy(Xt, double(jd), m_Yt3);
536 }
537 
538 
539 //====================================================================
540 void Fopr_NonRelativistic::add_R3(Field& Xt, const Field& Yt,
541  const int itime)
542 {
543  double fac = 1.0 / (8.0 * m_MQ * m_MQ);
544 
545  calc_Delta2(m_Yt1, Yt, itime);
546  axpy(Xt, fac, m_Yt1);
547 #pragma omp barrier
548 }
549 
550 
551 //====================================================================
552 void Fopr_NonRelativistic::add_R4(Field& Xt, const Field& Yt,
553  const int itime)
554 {
555  double fac = 1.0 / (8.0 * m_MQ * m_MQ);
556 
557  m_Yt3.set(0.0);
558 #pragma omp barrier
559 
560  mult_F(m_Yt1, Yt, 0, itime);
561  mult_sigma1(m_Yt2, m_Yt1);
562  axpy(m_Yt3, fac, m_Yt2);
563 #pragma omp barrier
564 
565  mult_F(m_Yt1, Yt, 1, itime);
566  mult_sigma2(m_Yt2, m_Yt1);
567  axpy(m_Yt3, fac, m_Yt2);
568 #pragma omp barrier
569 
570  mult_F(m_Yt1, Yt, 2, itime);
571  mult_sigma3(m_Yt2, m_Yt1);
572  axpy(m_Yt3, fac, m_Yt2);
573 #pragma omp barrier
574 
575  axpy(Xt, 1.0, m_Yt3);
576 #pragma omp barrier
577 }
578 
579 
580 //====================================================================
581 void Fopr_NonRelativistic::add_R5(Field& Xt, const Field& Yt,
582  const int itime)
583 {
584  double fac = 1.0 / (4.0 * m_MQ * m_MQ);
585 
586  m_Yt3.set(0.0);
587 #pragma omp barrier
588 
589  mult_F(m_Yt1, Yt, 3, itime);
590  mult_sigma1(m_Yt2, m_Yt1);
591  axpy(m_Yt3, fac, m_Yt2);
592 #pragma omp barrier
593 
594  mult_F(m_Yt1, Yt, 4, itime);
595  mult_sigma2(m_Yt2, m_Yt1);
596  axpy(m_Yt3, fac, m_Yt2);
597 #pragma omp barrier
598 
599  mult_F(m_Yt1, Yt, 5, itime);
600  mult_sigma3(m_Yt2, m_Yt1);
601  axpy(m_Yt3, fac, m_Yt2);
602 #pragma omp barrier
603 
604  axpy(Xt, -1.0, m_Yt3);
605 #pragma omp barrier
606 }
607 
608 
609 //====================================================================
610 void Fopr_NonRelativistic::evolve_impl(Field& Xt, const Field& Yt,
611  const int itime)
612 {
613  if (m_evolution_type == "evolve-A") {
614  evolve_typeA(Xt, Yt, itime);
615  } else if (m_evolution_type == "evolve-B") {
616  evolve_typeB(Xt, Yt, itime);
617  } else {
618  vout.crucial(m_vl, "Error in %s: unsupported evolution_type: %s\n",
619  class_name.c_str(), m_evolution_type.c_str());
620  exit(EXIT_FAILURE);
621  }
622 
623  ThreadManager_OpenMP::sync_barrier_all();
624 }
625 
626 
627 //====================================================================
628 void Fopr_NonRelativistic::evolve_typeA(Field& Xt, const Field& Yt,
629  const int itime)
630 {
631 #pragma omp barrier
632 
633  vout.paranoiac(m_vl, " evolve-A: itime = %d\n", itime);
634 
635  int itime1 = (itime - 1 + m_Ltime) % m_Ltime;
636 
637  copy(m_Zt1, Yt);
638 #pragma omp barrier
639 
640  evolve_H0(m_Zt2, m_Zt1, itime1);
641 
642  evolve_deltaH(m_Zt1, m_Zt2, itime1);
643 
644  evolve_U4(m_Zt2, m_Zt1, itime1);
645 
646  evolve_deltaH(m_Zt1, m_Zt2, itime);
647 
648  evolve_H0(m_Zt2, m_Zt1, itime);
649 
650  int ipet = itime / m_Nt;
651  int IPEt = Communicator::ipe(3);
652 
653  if (ipet == IPEt) copy(Xt, m_Zt2);
654 #pragma omp barrier
655 }
656 
657 
658 //====================================================================
659 void Fopr_NonRelativistic::evolve_typeB(Field& Xt, const Field& Yt,
660  const int itime)
661 {
662 #pragma omp barrier
663 
664  vout.paranoiac(m_vl, " evolve-B: itime = %d\n", itime);
665 
666  int itime1 = (itime - 1 + m_Ltime) % m_Ltime;
667 
668  copy(m_Zt1, Yt);
669 #pragma omp barrier
670 
671  if (itime1 != m_itime_src) {
672  evolve_deltaH(m_Zt2, m_Zt1, itime1);
673  m_Zt1.set(0.0);
674  copy(m_Zt1, m_Zt2);
675  }
676 #pragma omp barrier
677 
678  evolve_H0(m_Zt2, m_Zt1, itime1);
679 
680  evolve_U4(m_Zt1, m_Zt2, itime1);
681 
682  evolve_H0(m_Zt2, m_Zt1, itime);
683 
684  int ipet = itime / m_Nt;
685  int IPEt = Communicator::ipe(3);
686 
687  if (ipet == IPEt) copy(Xt, m_Zt2);
688 #pragma omp barrier
689 }
690 
691 
692 //====================================================================
693 void Fopr_NonRelativistic::evolve_H0(Field& Xt, const Field& Yt,
694  const int itime)
695 {
696 #pragma omp barrier
697 
698  double fac = 1.0 / (4.0 * m_MQ * double(m_nstab));
699 
700  copy(m_Yt1, Yt);
701 #pragma omp barrier
702 
703  for (int i = 0; i < m_nstab; ++i) {
704  calc_Delta2(Xt, m_Yt1, itime);
705  aypx(fac, Xt, m_Yt1);
706 #pragma omp barrier
707 
708  if (i < m_nstab - 1) {
709  copy(m_Yt1, Xt);
710 #pragma omp barrier
711  }
712  }
713 }
714 
715 
716 //====================================================================
717 void Fopr_NonRelativistic::evolve_U4(Field& Xt, const Field& Yt,
718  const int itime)
719 {
720  int NPEt = Communicator::npe(3);
721  int IPEt = Communicator::ipe(3);
722 
723  int it = itime % m_Nt;
724  int ipet = itime / m_Nt;
725 
726  ThreadManager_OpenMP::sync_barrier_all();
727 
728  if (ipet == IPEt) mult_Gd(m_Yt1, Yt, 3, itime);
729 
730  ThreadManager_OpenMP::sync_barrier_all();
731 
732  if (it < m_Nt - 1) {
733  if (ipet == IPEt) {
734  copy(Xt, m_Yt1);
735 #pragma omp barrier
736  }
737  } else {
738 #pragma omp master
739  {
740  int size = m_Nvc * m_Nd2 * m_Nspc;
741  double *yp1 = m_Yt1.ptr(0);
742  double *yp2 = m_Yt2.ptr(0);
743  Communicator::exchange(size, yp2, yp1, 3, -1, 7);
744  }
745 #pragma omp barrier
746 
747  int ipet_up = (ipet + 1) % NPEt;
748  if (ipet_up == IPEt) copy(Xt, m_Yt2);
749  }
750 
751  ThreadManager_OpenMP::sync_barrier_all();
752 }
753 
754 
755 //====================================================================
756 void Fopr_NonRelativistic::evolve_deltaH(Field& Xt, const Field& Yt,
757  const int itime)
758 {
759  // Be sure that the number of correction terms implemented here
760  // coincides with the variable m_num_correct (currently 6).
761 
762 #pragma omp barrier
763 
764  Xt.set(0.0);
765  add_deltaH1(Xt, Yt, itime);
766 
767  if (m_correction_terms == "next-leading") {
768  add_deltaH2(Xt, Yt, itime);
769  add_deltaH3(Xt, Yt, itime);
770  add_deltaH5(Xt, Yt, itime);
771  add_deltaH46(Xt, Yt, itime);
772  }
773 
774  double fac;
775  if (m_evolution_type == "evolve-A") {
776  fac = -0.5;
777  } else if (m_evolution_type == "evolve-B") {
778  fac = -1.0;
779  } else {
780  vout.crucial(m_vl, "Error in %s: unsupported evolution_type: %s\n",
781  class_name.c_str(), m_evolution_type.c_str());
782  exit(EXIT_FAILURE);
783  }
784 
785  aypx(fac, Xt, Yt);
786 #pragma omp barrier
787 }
788 
789 
790 //====================================================================
791 void Fopr_NonRelativistic::add_deltaH1(Field& Xt, const Field& Yt,
792  const int itime)
793 {
794  double fac = -m_coeff[0] / (2.0 * m_MQ);
795 
796  mult_F(m_Yt1, Yt, 0, itime);
797  mult_sigma1(m_Yt2, m_Yt1);
798  axpy(Xt, fac, m_Yt2);
799 #pragma omp barrier
800 
801  mult_F(m_Yt1, Yt, 1, itime);
802  mult_sigma2(m_Yt2, m_Yt1);
803  axpy(Xt, fac, m_Yt2);
804 #pragma omp barrier
805 
806  mult_F(m_Yt1, Yt, 2, itime);
807  mult_sigma3(m_Yt2, m_Yt1);
808  axpy(Xt, fac, m_Yt2);
809 #pragma omp barrier
810 }
811 
812 
813 //====================================================================
814 void Fopr_NonRelativistic::add_deltaH2(Field& Xt, const Field& Yt,
815  const int itime)
816 {
817  double fac = m_coeff[1] / (8.0 * m_MQ * m_MQ);
818 
819  m_Yt3.set(0.0);
820 
821  mult_F(m_Yt1, Yt, 3, itime);
822  calc_Delta1(m_Yt2, m_Yt1, 0, itime);
823  axpy(m_Yt3, 1.0, m_Yt2);
824 #pragma omp barrier
825 
826  mult_F(m_Yt1, Yt, 4, itime);
827  calc_Delta1(m_Yt2, m_Yt1, 1, itime);
828  axpy(m_Yt3, 1.0, m_Yt2);
829 #pragma omp barrier
830 
831  mult_F(m_Yt1, Yt, 5, itime);
832  calc_Delta1(m_Yt2, m_Yt1, 2, itime);
833  axpy(m_Yt3, 1.0, m_Yt2);
834 #pragma omp barrier
835 
836  calc_Delta1(m_Yt1, Yt, 0, itime);
837  mult_F(m_Yt2, m_Yt1, 3, itime);
838  axpy(m_Yt3, -1.0, m_Yt2);
839 #pragma omp barrier
840 
841  calc_Delta1(m_Yt1, Yt, 1, itime);
842  mult_F(m_Yt2, m_Yt1, 4, itime);
843  axpy(m_Yt3, -1.0, m_Yt2);
844 #pragma omp barrier
845 
846  calc_Delta1(m_Yt1, Yt, 2, itime);
847  mult_F(m_Yt2, m_Yt1, 5, itime);
848  axpy(m_Yt3, -1.0, m_Yt2);
849 #pragma omp barrier
850 
851  m_Yt3.xI();
852  axpy(Xt, fac, m_Yt3);
853 #pragma omp barrier
854 }
855 
856 
857 //====================================================================
858 void Fopr_NonRelativistic::add_deltaH3(Field& Xt, const Field& Yt,
859  const int itime)
860 {
861  double fac = -m_coeff[2] / (8.0 * m_MQ * m_MQ);
862 
863  // x-component
864  m_Yt3.set(0.0);
865 #pragma omp barrier
866 
867  mult_F(m_Yt1, Yt, 5, itime); // E_z
868  calc_Delta1(m_Yt2, m_Yt1, 1, itime); // D_y
869  axpy(m_Yt3, 1.0, m_Yt2); // Yt3 += D_y * E_z
870 #pragma omp barrier
871 
872  mult_F(m_Yt1, Yt, 4, itime); // E_y
873  calc_Delta1(m_Yt2, m_Yt1, 2, itime); // D_z
874  axpy(m_Yt3, -1.0, m_Yt2); // Yt3 -= D_z * E_y
875 #pragma omp barrier
876 
877  calc_Delta1(m_Yt1, Yt, 2, itime); // D_z
878  mult_F(m_Yt2, m_Yt1, 4, itime); // E_y
879  axpy(m_Yt3, -1.0, m_Yt2); // Yt3 += E_y * D_z
880 #pragma omp barrier
881 
882  calc_Delta1(m_Yt1, Yt, 1, itime); // D_y
883  mult_F(m_Yt2, m_Yt1, 5, itime); // E_z
884  axpy(m_Yt3, 1.0, m_Yt2); // Yt3 -= E_z * D_y
885 #pragma omp barrier
886 
887  mult_sigma1(m_Yt1, m_Yt3);
888  axpy(Xt, fac, m_Yt1);
889 #pragma omp barrier
890 
891  // y-component
892  m_Yt3.set(0.0);
893 #pragma omp barrier
894 
895  mult_F(m_Yt1, Yt, 3, itime); // E_x
896  calc_Delta1(m_Yt2, m_Yt1, 2, itime); // D_z
897  axpy(m_Yt3, 1.0, m_Yt2); // Yt3 += D_z * E_x
898 #pragma omp barrier
899 
900  mult_F(m_Yt1, Yt, 5, itime); // E_z
901  calc_Delta1(m_Yt2, m_Yt1, 0, itime); // D_x
902  axpy(m_Yt3, -1.0, m_Yt2); // Yt3 -= D_x * E_z
903 #pragma omp barrier
904 
905  calc_Delta1(m_Yt1, Yt, 0, itime); // D_x
906  mult_F(m_Yt2, m_Yt1, 5, itime); // E_z
907  axpy(m_Yt3, -1.0, m_Yt2); // Yt3 += E_z * D_x
908 #pragma omp barrier
909 
910  calc_Delta1(m_Yt1, Yt, 2, itime); // D_z
911  mult_F(m_Yt2, m_Yt1, 3, itime); // E_x
912  axpy(m_Yt3, 1.0, m_Yt2); // Yt3 -= E_x * D_z
913 #pragma omp barrier
914 
915  mult_sigma2(m_Yt1, m_Yt3);
916  axpy(Xt, fac, m_Yt1);
917 #pragma omp barrier
918 
919  // z-component
920  m_Yt3.set(0.0);
921 #pragma omp barrier
922 
923  mult_F(m_Yt1, Yt, 4, itime); // E_y
924  calc_Delta1(m_Yt2, m_Yt1, 0, itime); // D_x
925  axpy(m_Yt3, 1.0, m_Yt2); // Yt3 += D_x * E_y
926 #pragma omp barrier
927 
928  mult_F(m_Yt1, Yt, 3, itime); // E_x
929  calc_Delta1(m_Yt2, m_Yt1, 1, itime); // D_y
930  axpy(m_Yt3, -1.0, m_Yt2); // Yt3 -= D_y * E_x
931 #pragma omp barrier
932 
933  calc_Delta1(m_Yt1, Yt, 1, itime); // D_y
934  mult_F(m_Yt2, m_Yt1, 3, itime); // E_x
935  axpy(m_Yt3, -1.0, m_Yt2); // Yt3 += E_x * D_y
936 #pragma omp barrier
937 
938  calc_Delta1(m_Yt1, Yt, 0, itime); // D_x
939  mult_F(m_Yt2, m_Yt1, 4, itime); // E_y
940  axpy(m_Yt3, 1.0, m_Yt2); // Yt3 -= E_y * D_x
941 #pragma omp barrier
942 
943  mult_sigma3(m_Yt1, m_Yt3);
944  axpy(Xt, fac, m_Yt1);
945 #pragma omp barrier
946 }
947 
948 
949 //====================================================================
950 void Fopr_NonRelativistic::add_deltaH5(Field& Xt, const Field& Yt,
951  const int itime)
952 {
953  double fac = m_coeff[4] / (24.0 * m_MQ);
954 
955  calc_Delta4(m_Yt1, Yt, itime);
956  axpy(Xt, fac, m_Yt1);
957 #pragma omp barrier
958 }
959 
960 
961 //====================================================================
962 void Fopr_NonRelativistic::add_deltaH46(Field& Xt, const Field& Yt,
963  const int itime)
964 {
965  double fac4 = -m_coeff[3] / (8.0 * m_MQ * m_MQ * m_MQ);
966  double fac6 = -m_coeff[5] / (16.0 * double(m_nstab) * m_MQ * m_MQ);
967  double fac = fac4 + fac6;
968 
969  calc_Delta2(m_Yt1, Yt, itime);
970  calc_Delta2(m_Yt2, m_Yt1, itime);
971  axpy(Xt, fac, m_Yt2);
972 #pragma omp barrier
973 }
974 
975 
976 //====================================================================
977 void Fopr_NonRelativistic::calc_Delta2(Field& Xt, const Field& Yt,
978  const int itime)
979 {
980  int ipet = itime / m_Nt;
981  int IPEt = Communicator::ipe(3);
982  if (ipet != IPEt) return;
983 
984 #pragma omp barrier
985 
986  Xt.set(0.0);
987  axpy(Xt, -6.0, Yt);
988 #pragma omp barrier
989 
990  for (int idir = 0; idir < m_Ndim - 1; ++idir) {
991  shift_backward(m_Wt1, Yt, idir, itime);
992  mult_Gn(m_Wt2, m_Wt1, idir, itime);
993  axpy(Xt, 1.0, m_Wt2);
994 #pragma omp barrier
995 
996  mult_Gd(m_Wt1, Yt, idir, itime);
997  shift_forward(m_Wt2, m_Wt1, idir, itime);
998  axpy(Xt, 1.0, m_Wt2);
999 #pragma omp barrier
1000  }
1001 }
1002 
1003 
1004 //====================================================================
1005 void Fopr_NonRelativistic::calc_Delta4(Field& Xt, const Field& Yt,
1006  const int itime)
1007 {
1008  int ipet = itime / m_Nt;
1009  int IPEt = Communicator::ipe(3);
1010  if (ipet != IPEt) return;
1011 
1012 #pragma omp barrier
1013 
1014  Xt.set(0.0);
1015 #pragma omp barrier
1016 
1017  for (int idir = 0; idir < m_Ndim - 1; ++idir) {
1018  m_Wt3.set(0.0);
1019 
1020  shift_backward(m_Wt1, Yt, idir, itime);
1021  mult_Gn(m_Wt2, m_Wt1, idir, itime);
1022  axpy(m_Wt3, 1.0, m_Wt2);
1023 #pragma omp barrier
1024 
1025  mult_Gd(m_Wt1, Yt, idir, itime);
1026  shift_forward(m_Wt2, m_Wt1, idir, itime);
1027  axpy(m_Wt3, 1.0, m_Wt2);
1028 
1029  axpy(m_Wt3, -2.0, Yt);
1030 #pragma omp barrier
1031 
1032  shift_backward(m_Wt1, m_Wt3, idir, itime);
1033  mult_Gn(m_Wt2, m_Wt1, idir, itime);
1034  axpy(Xt, 1.0, m_Wt2);
1035 #pragma omp barrier
1036 
1037  mult_Gd(m_Wt1, m_Wt3, idir, itime);
1038  shift_forward(m_Wt2, m_Wt1, idir, itime);
1039  axpy(Xt, 1.0, m_Wt2);
1040 
1041  axpy(Xt, -2.0, m_Wt3);
1042 #pragma omp barrier
1043  }
1044 }
1045 
1046 
1047 //====================================================================
1048 void Fopr_NonRelativistic::calc_Delta1(Field& Xt, const Field& Yt,
1049  const int idir, const int itime)
1050 {
1051  int ipet = itime / m_Nt;
1052  int IPEt = Communicator::ipe(3);
1053  if (ipet != IPEt) return;
1054 
1055 #pragma omp barrier
1056 
1057  Xt.set(0.0);
1058 #pragma omp barrier
1059 
1060  shift_backward(m_Wt1, Yt, idir, itime);
1061  mult_Gn(m_Wt2, m_Wt1, idir, itime);
1062  axpy(Xt, 0.5, m_Wt2);
1063 #pragma omp barrier
1064 
1065  mult_Gd(m_Wt1, Yt, idir, itime);
1066  shift_forward(m_Wt2, m_Wt1, idir, itime);
1067  axpy(Xt, -0.5, m_Wt2);
1068 #pragma omp barrier
1069 }
1070 
1071 
1072 //====================================================================
1073 void Fopr_NonRelativistic::set_source(Field& Xt, const Field& w)
1074 {
1075  // here check of Field size is necessary.
1076 
1077  vout.detailed(m_vl, " source norm2 = %e.\n", w.norm2());
1078  // vout.paranoiac(m_vl, " source norm2 = %e.\n", w.norm2());
1079 
1080 #pragma omp master
1081  {
1082  int IPEt = Communicator::ipe(m_Ndim - 1);
1083  int Nvcd2 = m_Nvc * m_Nd2; // only upper spinor components matter
1084  double epsilon = 1.0e-10;
1085 
1086  // extracting source time slice
1087  int ndef_src = 0;
1088  for (int itime = 0; itime < m_Ltime; ++itime) {
1089  int it = itime % m_Nt;
1090  int ipet = itime / m_Nt;
1091 
1092  double norm2t = 0.0;
1093  if (ipet == IPEt) {
1094  for (int ispc = 0; ispc < m_Nspc; ++ispc) {
1095  for (int ivcd = 0; ivcd < Nvcd2; ++ivcd) {
1096  double wt = w.cmp(ivcd, ispc + m_Nspc * it, 0);
1097  norm2t += wt * wt;
1098  }
1099  }
1100  }
1101  norm2t = Communicator::reduce_sum(norm2t);
1102 
1103  if (norm2t > epsilon) {
1104  m_itime_src = itime;
1105  ndef_src += 1;
1106  }
1107  }
1108 
1109  if (ndef_src == 0) {
1110  vout.crucial(m_vl, "Error in %s: source is not set.\n",
1111  class_name.c_str());
1112  exit(EXIT_FAILURE);
1113  }
1114 
1115  if (ndef_src > 1) {
1116  vout.crucial(m_vl, "Error in %s: source is multiplly set.\n",
1117  class_name.c_str());
1118  exit(EXIT_FAILURE);
1119  }
1120 
1121  // setting source vector on single time slice
1122  int it = m_itime_src % m_Nt;
1123  int ipet = m_itime_src / m_Nt;
1124  if (ipet == IPEt) {
1125  for (int ispc = 0; ispc < m_Nspc; ++ispc) {
1126  for (int id = 0; id < m_Nd2; ++id) {
1127  for (int ivc = 0; ivc < m_Nvc; ++ivc) {
1128  int ivcd = ivc + m_Nvc * id;
1129  double wt = w.cmp(ivcd, ispc + m_Nspc * it, 0);
1130  Xt.set(ivcd, ispc, 0, wt);
1131  }
1132  }
1133  }
1134  } else {
1135  Xt.set(0.0);
1136  }
1137  } // #pragma omp master
1138 
1139  vout.detailed(m_vl, " source is set on itime = %d.\n", m_itime_src);
1140  vout.detailed(m_vl, " source norm2 = %e.\n", Xt.norm2());
1141 
1142 #pragma omp barrier
1143 }
1144 
1145 
1146 //====================================================================
1147 void Fopr_NonRelativistic::shift_forward(Field& Xt, const Field& Yt,
1148  const int idir, const int itime)
1149 { // Xt(x) = Yt(x-\hat{idir})
1150  if (idir == 0) {
1151  shift_xdn(Xt, Yt, itime);
1152  } else if (idir == 1) {
1153  shift_ydn(Xt, Yt, itime);
1154  } else if (idir == 2) {
1155  shift_zdn(Xt, Yt, itime);
1156  } else {
1157  vout.crucial(m_vl, "%s: irrelevant direction = %d\n",
1158  class_name.c_str(), idir);
1159  exit(EXIT_FAILURE);
1160  }
1161 }
1162 
1163 
1164 //====================================================================
1165 void Fopr_NonRelativistic::shift_backward(Field& Xt, const Field& Yt,
1166  const int idir, const int itime)
1167 { // Xt(x) = Yt(x+\hat{idir})
1168  if (idir == 0) {
1169  shift_xup(Xt, Yt, itime);
1170  } else if (idir == 1) {
1171  shift_yup(Xt, Yt, itime);
1172  } else if (idir == 2) {
1173  shift_zup(Xt, Yt, itime);
1174  } else {
1175  vout.crucial(m_vl, "%s: irrelevant direction = %d\n",
1176  class_name.c_str(), idir);
1177  exit(EXIT_FAILURE);
1178  }
1179 }
1180 
1181 
1182 //====================================================================
1183 void Fopr_NonRelativistic::shift_xup(Field& Xt, const Field& Yt,
1184  const int itime)
1185 {
1186  int idir = 0;
1187 
1188  int ipet = itime / m_Nt;
1189  int IPEt = Communicator::ipe(3);
1190  if (ipet != IPEt) return;
1191 
1192  double *xp = Xt.ptr(0);
1193  const double *yp = Yt.ptr(0);
1194 
1195  int NinF2 = 2 * m_Nc * m_Nd2;
1196  int Nyz = m_Ny * m_Nz;
1197 
1198  double bc2 = 1.0;
1199  if (Communicator::ipe(idir) == 0) bc2 = m_boundary[idir];
1200 
1201 #pragma omp barrier
1202 
1203  int ith, nth, is, ns;
1204  set_threadtask(ith, nth, is, ns, Nyz);
1205 
1206  // boundary
1207  int ix = 0;
1208  for (int iyz = is; iyz < ns; ++iyz) {
1209  int nei = ix + m_Nx * iyz;
1210  for (int in = 0; in < NinF2; ++in) {
1211  buf1_x[in + NinF2 * iyz] = bc2 * yp[in + NinF2 * nei];
1212  }
1213  }
1214 #pragma omp barrier
1215 
1216 #pragma omp master
1217  {
1218  const int size = NinF2 * Nyz;
1219  Communicator::exchange(size, &buf2_x[0], &buf1_x[0], 0, 1, 0);
1220  }
1221 #pragma omp barrier
1222 
1223  set_threadtask(ith, nth, is, ns, m_Nspc);
1224 
1225  for (int site = is; site < ns; ++site) {
1226  int ix = site % m_Nx;
1227  int iyz = site / m_Nx;
1228  if (ix < m_Nx - 1) { // bulk
1229  int nei = ix + 1 + m_Nx * iyz;
1230  for (int in = 0; in < NinF2; ++in) {
1231  xp[in + NinF2 * site] = yp[in + NinF2 * nei];
1232  }
1233  } else { // boundary
1234  for (int in = 0; in < NinF2; ++in) {
1235  xp[in + NinF2 * site] = buf2_x[in + NinF2 * iyz];
1236  }
1237  }
1238  }
1239 
1240 #pragma omp barrier
1241 }
1242 
1243 
1244 //====================================================================
1245 void Fopr_NonRelativistic::shift_xdn(Field& Xt, const Field& Yt,
1246  const int itime)
1247 {
1248  int idir = 0;
1249 
1250  int ipet = itime / m_Nt;
1251  int IPEt = Communicator::ipe(3);
1252  if (ipet != IPEt) return;
1253 
1254  double *xp = Xt.ptr(0);
1255  const double *yp = Yt.ptr(0);
1256 
1257  int NinF2 = 2 * m_Nc * m_Nd2;
1258  int Nyz = m_Ny * m_Nz;
1259 
1260  double bc2 = 1.0;
1261  if (Communicator::ipe(idir) == 0) bc2 = m_boundary[idir];
1262 
1263 #pragma omp barrier
1264 
1265  int ith, nth, is, ns;
1266  set_threadtask(ith, nth, is, ns, Nyz);
1267 
1268  // boundary
1269  int ix = m_Nx - 1;
1270  for (int iyz = is; iyz < ns; ++iyz) {
1271  int nei = ix + m_Nx * iyz;
1272  for (int in = 0; in < NinF2; ++in) {
1273  buf1_x[in + NinF2 * iyz] = yp[in + NinF2 * nei];
1274  }
1275  }
1276 #pragma omp barrier
1277 
1278 #pragma omp master
1279  {
1280  const int size = NinF2 * Nyz;
1281  Communicator::exchange(size, &buf2_x[0], &buf1_x[0], 0, -1, 4);
1282  }
1283 #pragma omp barrier
1284 
1285  set_threadtask(ith, nth, is, ns, m_Nspc);
1286 
1287  for (int site = is; site < ns; ++site) {
1288  int ix = site % m_Nx;
1289  int iyz = site / m_Nx;
1290  if (ix > 0) { // bulk
1291  int nei = ix - 1 + m_Nx * iyz;
1292  for (int in = 0; in < NinF2; ++in) {
1293  xp[in + NinF2 * site] = yp[in + NinF2 * nei];
1294  }
1295  } else { // boundary
1296  for (int in = 0; in < NinF2; ++in) {
1297  xp[in + NinF2 * site] = bc2 * buf2_x[in + NinF2 * iyz];
1298  }
1299  }
1300  }
1301 #pragma omp barrier
1302 }
1303 
1304 
1305 //====================================================================
1306 void Fopr_NonRelativistic::shift_yup(Field& Xt, const Field& Yt,
1307  const int itime)
1308 {
1309  int idir = 1;
1310 
1311  int ipet = itime / m_Nt;
1312  int IPEt = Communicator::ipe(3);
1313  if (ipet != IPEt) return;
1314 
1315  double *xp = Xt.ptr(0);
1316  const double *yp = Yt.ptr(0);
1317 
1318  int NinF2 = 2 * m_Nc * m_Nd2;
1319  int Nxz = m_Nx * m_Nz;
1320 
1321  double bc2 = 1.0;
1322  if (Communicator::ipe(idir) == 0) bc2 = m_boundary[idir];
1323 
1324 #pragma omp barrier
1325 
1326  int ith, nth, is, ns;
1327  set_threadtask(ith, nth, is, ns, Nxz);
1328 
1329  // boundary
1330  int iy = 0;
1331  for (int ixz = is; ixz < ns; ++ixz) {
1332  int ix = ixz % m_Nx;
1333  int iz = ixz / m_Nx;
1334  int nei = ix + m_Nx * (iy + m_Ny * iz);
1335  for (int in = 0; in < NinF2; ++in) {
1336  buf1_y[in + NinF2 * ixz] = bc2 * yp[in + NinF2 * nei];
1337  }
1338  }
1339 #pragma omp barrier
1340 
1341 #pragma omp master
1342  {
1343  const int size = NinF2 * Nxz;
1344  Communicator::exchange(size, &buf2_y[0], &buf1_y[0], 1, 1, 1);
1345  }
1346 #pragma omp barrier
1347 
1348  set_threadtask(ith, nth, is, ns, m_Nspc);
1349 
1350  for (int site = is; site < ns; ++site) {
1351  int ix = site % m_Nx;
1352  int iy = (site / m_Nx) % m_Ny;
1353  int iz = site / (m_Nx * m_Ny);
1354  if (iy < m_Ny - 1) { // bulk
1355  int nei = ix + m_Nx * (iy + 1 + m_Ny * iz);
1356  for (int in = 0; in < NinF2; ++in) {
1357  xp[in + NinF2 * site] = yp[in + NinF2 * nei];
1358  }
1359  } else { // boundary
1360  int ixz = ix + m_Nx * iz;
1361  for (int in = 0; in < NinF2; ++in) {
1362  xp[in + NinF2 * site] = buf2_y[in + NinF2 * ixz];
1363  }
1364  }
1365  }
1366 #pragma omp barrier
1367 }
1368 
1369 
1370 //====================================================================
1371 void Fopr_NonRelativistic::shift_ydn(Field& Xt, const Field& Yt,
1372  const int itime)
1373 {
1374  int idir = 1;
1375 
1376  int ipet = itime / m_Nt;
1377  int IPEt = Communicator::ipe(3);
1378  if (ipet != IPEt) return;
1379 
1380  double *xp = Xt.ptr(0);
1381  const double *yp = Yt.ptr(0);
1382 
1383  int NinF2 = 2 * m_Nc * m_Nd2;
1384  int Nxz = m_Nx * m_Nz;
1385 
1386  double bc2 = 1.0;
1387  if (Communicator::ipe(idir) == 0) bc2 = m_boundary[idir];
1388 
1389 #pragma omp barrier
1390 
1391  int ith, nth, is, ns;
1392  set_threadtask(ith, nth, is, ns, Nxz);
1393 
1394  // boundary
1395  int iy = m_Ny - 1;
1396  for (int ixz = is; ixz < ns; ++ixz) {
1397  int ix = ixz % m_Nx;
1398  int iz = ixz / m_Nx;
1399  int nei = ix + m_Nx * (iy + m_Ny * iz);
1400  for (int in = 0; in < NinF2; ++in) {
1401  buf1_y[in + NinF2 * ixz] = yp[in + NinF2 * nei];
1402  }
1403  }
1404 #pragma omp barrier
1405 
1406 #pragma omp master
1407  {
1408  const int size = NinF2 * Nxz;
1409  Communicator::exchange(size, &buf2_y[0], &buf1_y[0], 1, -1, 5);
1410  }
1411 #pragma omp barrier
1412 
1413  set_threadtask(ith, nth, is, ns, m_Nspc);
1414 
1415  for (int site = is; site < ns; ++site) {
1416  int ix = site % m_Nx;
1417  int iy = (site / m_Nx) % m_Ny;
1418  int iz = site / (m_Nx * m_Ny);
1419  if (iy > 0) { // bulk
1420  int nei = ix + m_Nx * (iy - 1 + m_Ny * iz);
1421  for (int in = 0; in < NinF2; ++in) {
1422  xp[in + NinF2 * site] = yp[in + NinF2 * nei];
1423  }
1424  } else { // boundary
1425  int ixz = ix + m_Nx * iz;
1426  for (int in = 0; in < NinF2; ++in) {
1427  xp[in + NinF2 * site] = bc2 * buf2_y[in + NinF2 * ixz];
1428  }
1429  }
1430  }
1431 
1432 #pragma omp barrier
1433 }
1434 
1435 
1436 //====================================================================
1437 void Fopr_NonRelativistic::shift_zup(Field& Xt, const Field& Yt,
1438  const int itime)
1439 {
1440  int idir = 2;
1441 
1442  int ipet = itime / m_Nt;
1443  int IPEt = Communicator::ipe(3);
1444  if (ipet != IPEt) return;
1445 
1446  double *xp = Xt.ptr(0);
1447  const double *yp = Yt.ptr(0);
1448 
1449  int NinF2 = 2 * m_Nc * m_Nd2;
1450  int Nxy = m_Nx * m_Ny;
1451  double bc2 = 1.0;
1452  if (Communicator::ipe(idir) == 0) bc2 = m_boundary[idir];
1453 
1454 #pragma omp barrier
1455 
1456  int ith, nth, is, ns;
1457  set_threadtask(ith, nth, is, ns, Nxy);
1458 
1459  // boundary
1460  int iz = 0;
1461  for (int ixy = is; ixy < ns; ++ixy) {
1462  int nei = ixy + Nxy * iz;
1463  for (int in = 0; in < NinF2; ++in) {
1464  buf1_z[in + NinF2 * ixy] = bc2 * yp[in + NinF2 * nei];
1465  }
1466  }
1467 #pragma omp barrier
1468 
1469 #pragma omp master
1470  {
1471  const int size = NinF2 * Nxy;
1472  Communicator::exchange(size, &buf2_z[0], &buf1_z[0], 2, 1, 2);
1473  }
1474 #pragma omp barrier
1475 
1476  set_threadtask(ith, nth, is, ns, m_Nspc);
1477 
1478  for (int site = is; site < ns; ++site) {
1479  int ixy = site % Nxy;
1480  int iz = site / Nxy;
1481  if (iz < m_Nz - 1) { // bulk
1482  int nei = ixy + Nxy * (iz + 1);
1483  for (int in = 0; in < NinF2; ++in) {
1484  xp[in + NinF2 * site] = yp[in + NinF2 * nei];
1485  }
1486  } else { // boundary
1487  for (int in = 0; in < NinF2; ++in) {
1488  xp[in + NinF2 * site] = buf2_z[in + NinF2 * ixy];
1489  }
1490  }
1491  }
1492 #pragma omp barrier
1493 }
1494 
1495 
1496 //====================================================================
1497 void Fopr_NonRelativistic::shift_zdn(Field& Xt, const Field& Yt,
1498  const int itime)
1499 {
1500  int idir = 2;
1501 
1502  int ipet = itime / m_Nt;
1503  int IPEt = Communicator::ipe(3);
1504  if (ipet != IPEt) return;
1505 
1506  double *xp = Xt.ptr(0);
1507  const double *yp = Yt.ptr(0);
1508 
1509  int NinF2 = 2 * m_Nc * m_Nd2;
1510  int Nxy = m_Nx * m_Ny;
1511  double bc2 = 1.0;
1512  if (Communicator::ipe(idir) == 0) bc2 = m_boundary[idir];
1513 
1514 #pragma omp barrier
1515 
1516  int ith, nth, is, ns;
1517  set_threadtask(ith, nth, is, ns, Nxy);
1518 
1519  // boundary
1520  int iz = m_Nz - 1;
1521  for (int ixy = is; ixy < ns; ++ixy) {
1522  int nei = ixy + Nxy * iz;
1523  for (int in = 0; in < NinF2; ++in) {
1524  buf1_z[in + NinF2 * ixy] = yp[in + NinF2 * nei];
1525  }
1526  }
1527 #pragma omp barrier
1528 
1529 #pragma omp master
1530  {
1531  const int size = NinF2 * Nxy;
1532  Communicator::exchange(size, &buf2_z[0], &buf1_z[0], 2, -1, 6);
1533  }
1534 #pragma omp barrier
1535 
1536  set_threadtask(ith, nth, is, ns, m_Nspc);
1537 
1538  // bulk
1539  for (int site = is; site < ns; ++site) {
1540  int ixy = site % Nxy;
1541  int iz = site / Nxy;
1542  if (iz > 0) {
1543  int nei = ixy + Nxy * (iz - 1);
1544  for (int in = 0; in < NinF2; ++in) {
1545  xp[in + NinF2 * site] = yp[in + NinF2 * nei];
1546  }
1547  } else { // boundary
1548  for (int in = 0; in < NinF2; ++in) {
1549  xp[in + NinF2 * site] = bc2 * buf2_z[in + NinF2 * ixy];
1550  }
1551  }
1552  }
1553 #pragma omp barrier
1554 }
1555 
1556 
1557 //====================================================================
1558 void Fopr_NonRelativistic::mult_sigma1(Field& Xt, const Field& Yt)
1559 {
1560  int Nc2 = 2 * m_Nc;
1561  int Ncd = Nc2 * m_Nd2;
1562 
1563  double *xp = Xt.ptr(0);
1564  const double *yp = Yt.ptr(0);
1565 
1566  int ith, nth, is, ns;
1567  set_threadtask(ith, nth, is, ns, m_Nspc);
1568 
1569 #pragma omp barrier
1570 
1571  for (int site = is; site < ns; ++site) {
1572  int iv = Ncd * site;
1573  for (int ic = 0; ic < m_Nc; ++ic) {
1574  xp[2 * ic + Nc2 * 0 + iv] = yp[2 * ic + Nc2 * 1 + iv];
1575  xp[2 * ic + 1 + Nc2 * 0 + iv] = yp[2 * ic + 1 + Nc2 * 1 + iv];
1576  xp[2 * ic + Nc2 * 1 + iv] = yp[2 * ic + Nc2 * 0 + iv];
1577  xp[2 * ic + 1 + Nc2 * 1 + iv] = yp[2 * ic + 1 + Nc2 * 0 + iv];
1578  }
1579  }
1580 #pragma omp barrier
1581 }
1582 
1583 
1584 //====================================================================
1585 void Fopr_NonRelativistic::mult_sigma2(Field& Xt, const Field& Yt)
1586 {
1587  int Nc2 = 2 * m_Nc;
1588  int Ncd = Nc2 * m_Nd2;
1589 
1590  double *xp = Xt.ptr(0);
1591  const double *yp = Yt.ptr(0);
1592 
1593  int ith, nth, is, ns;
1594  set_threadtask(ith, nth, is, ns, m_Nspc);
1595 
1596 #pragma omp barrier
1597 
1598  for (int site = is; site < ns; ++site) {
1599  int iv = Ncd * site;
1600  for (int ic = 0; ic < m_Nc; ++ic) {
1601  xp[2 * ic + Nc2 * 0 + iv] = yp[2 * ic + 1 + Nc2 * 1 + iv];
1602  xp[2 * ic + 1 + Nc2 * 0 + iv] = -yp[2 * ic + Nc2 * 1 + iv];
1603  xp[2 * ic + Nc2 * 1 + iv] = -yp[2 * ic + 1 + Nc2 * 0 + iv];
1604  xp[2 * ic + 1 + Nc2 * 1 + iv] = yp[2 * ic + Nc2 * 0 + iv];
1605  }
1606  }
1607 #pragma omp barrier
1608 }
1609 
1610 
1611 //====================================================================
1612 void Fopr_NonRelativistic::mult_sigma3(Field& Xt, const Field& Yt)
1613 {
1614  int Nc2 = 2 * m_Nc;
1615  int Ncd = Nc2 * m_Nd2;
1616 
1617  double *xp = Xt.ptr(0);
1618  const double *yp = Yt.ptr(0);
1619 
1620  int ith, nth, is, ns;
1621  set_threadtask(ith, nth, is, ns, m_Nspc);
1622 
1623 #pragma omp barrier
1624 
1625  for (int site = is; site < ns; ++site) {
1626  int iv = Ncd * site;
1627  for (int ic = 0; ic < m_Nc; ++ic) {
1628  xp[2 * ic + Nc2 * 0 + iv] = yp[2 * ic + Nc2 * 0 + iv];
1629  xp[2 * ic + 1 + Nc2 * 0 + iv] = yp[2 * ic + 1 + Nc2 * 0 + iv];
1630  xp[2 * ic + Nc2 * 1 + iv] = -yp[2 * ic + Nc2 * 1 + iv];
1631  xp[2 * ic + 1 + Nc2 * 1 + iv] = -yp[2 * ic + 1 + Nc2 * 1 + iv];
1632  }
1633  }
1634 #pragma omp barrier
1635 }
1636 
1637 
1638 //====================================================================
1639 void Fopr_NonRelativistic::mult_Gn(Field& Xt, const Field& Yt,
1640  const int idir, const int itime)
1641 {
1642  int ipet = itime / m_Nt;
1643  int IPEt = Communicator::ipe(3);
1644  if (ipet != IPEt) return;
1645 
1646  int it = itime % m_Nt;
1647  const double *up = m_U.ptr(0, 0, idir);
1648 
1649  double *xp = Xt.ptr(0);
1650  const double *yp = Yt.ptr(0);
1651 
1652  int Nc2 = 2 * m_Nc;
1653  int Ndf = 2 * m_Nc * m_Nc;
1654  int Ncd = 2 * m_Nc * m_Nd2;
1655 
1656  int ith, nth, is, ns;
1657  set_threadtask(ith, nth, is, ns, m_Nspc);
1658 
1659 #pragma omp barrier
1660 
1661  for (int site = is; site < ns; ++site) {
1662  int ig = Ndf * (site + m_Nspc * it);
1663  int iv = Ncd * site;
1664  for (int id = 0; id < m_Nd2; ++id) {
1665  for (int ic = 0; ic < m_Nc; ++ic) {
1666  int ig2 = Nc2 * ic + ig;
1667  int iv2 = Nc2 * id + iv;
1668  xp[2 * ic + iv2] = mult_Gn_r(&up[ig2], &yp[iv2], m_Nc);
1669  xp[2 * ic + 1 + iv2] = mult_Gn_i(&up[ig2], &yp[iv2], m_Nc);
1670  }
1671  }
1672  }
1673 #pragma omp barrier
1674 }
1675 
1676 
1677 //====================================================================
1678 void Fopr_NonRelativistic::mult_F(Field& Xt, const Field& Yt,
1679  const int icomp, const int itime)
1680 {
1681  int ipet = itime / m_Nt;
1682  int IPEt = Communicator::ipe(3);
1683  if (ipet != IPEt) return;
1684 
1685  int it = itime % m_Nt;
1686  const double *up = m_Fstr[icomp].ptr(0);
1687 
1688  double *xp = Xt.ptr(0);
1689  const double *yp = Yt.ptr(0);
1690 
1691  int Nc2 = 2 * m_Nc;
1692  int Ndf = 2 * m_Nc * m_Nc;
1693  int Ncd = 2 * m_Nc * m_Nd2;
1694 
1695  int ith, nth, is, ns;
1696  set_threadtask(ith, nth, is, ns, m_Nspc);
1697 
1698 #pragma omp barrier
1699 
1700  for (int site = is; site < ns; ++site) {
1701  int ig = Ndf * (site + m_Nspc * it);
1702  int iv = Ncd * site;
1703  for (int id = 0; id < m_Nd2; ++id) {
1704  for (int ic = 0; ic < m_Nc; ++ic) {
1705  int ig2 = Nc2 * ic + ig;
1706  int iv2 = Nc2 * id + iv;
1707  xp[2 * ic + iv2] = mult_Gn_r(&up[ig2], &yp[iv2], m_Nc);
1708  xp[2 * ic + 1 + iv2] = mult_Gn_i(&up[ig2], &yp[iv2], m_Nc);
1709  }
1710  }
1711  }
1712 #pragma omp barrier
1713 }
1714 
1715 
1716 //====================================================================
1717 void Fopr_NonRelativistic::mult_Gd(Field& Xt, const Field& Yt,
1718  const int idir, const int itime)
1719 {
1720  int ipet = itime / m_Nt;
1721  int IPEt = Communicator::ipe(3);
1722  if (ipet != IPEt) return;
1723 
1724  int it = itime % m_Nt;
1725  const double *up = m_U.ptr(0, 0, idir);
1726 
1727  double *xp = Xt.ptr(0);
1728  const double *yp = Yt.ptr(0);
1729 
1730  int Nc2 = 2 * m_Nc;
1731  int Ndf = 2 * m_Nc * m_Nc;
1732  int Ncd = 2 * m_Nc * m_Nd2;
1733 
1734  int ith, nth, is, ns;
1735  set_threadtask(ith, nth, is, ns, m_Nspc);
1736 
1737 #pragma omp barrier
1738 
1739  for (int site = is; site < ns; ++site) {
1740  int ig = Ndf * (site + m_Nspc * it);
1741  int iv = Ncd * site;
1742  for (int id = 0; id < m_Nd2; ++id) {
1743  for (int ic = 0; ic < m_Nc; ++ic) {
1744  int ig2 = 2 * ic + ig;
1745  int iv2 = Nc2 * id + iv;
1746  xp[2 * ic + iv2] = mult_Gd_r(&up[ig2], &yp[iv2], m_Nc);
1747  xp[2 * ic + 1 + iv2] = mult_Gd_i(&up[ig2], &yp[iv2], m_Nc);
1748  }
1749  }
1750  }
1751 
1752 #pragma omp barrier
1753 }
1754 
1755 
1756 //====================================================================
1757 double Fopr_NonRelativistic::flop_count()
1758 {
1759  return flop_count("Dirac");
1760 }
1761 
1762 
1763 //====================================================================
1764 double Fopr_NonRelativistic::flop_count(const std::string repr)
1765 {
1766  // not implemented yet
1767 
1768  const double gflop = 0.0;
1769 
1770  return gflop;
1771 }
1772 
1773 
1774 //============================================================END=====
Fopr_NonRelativistic::mult
void mult(Field &v, const Field &w)
mult with m_mode
Definition: fopr_NonRelativistic.cpp:276
CommonParameters::Ny
static int Ny()
Definition: commonParameters.h:106
CommonParameters::Nz
static int Nz()
Definition: commonParameters.h:107
Fopr_NonRelativistic::m_Zt1
Field m_Zt1
Definition: fopr_NonRelativistic.h:72
Parameters::set_string
void set_string(const string &key, const string &value)
Definition: parameters.cpp:39
Fopr_NonRelativistic::m_Ny
int m_Ny
Definition: fopr_NonRelativistic.h:54
Fopr_NonRelativistic::init
void init(const Parameters &params)
initial setup.
Definition: fopr_NonRelativistic.cpp:69
Fopr_NonRelativistic::buf1_z
std::vector< double > buf1_z
Definition: fopr_NonRelativistic.h:67
Fopr_NonRelativistic::evolve_deltaH
void evolve_deltaH(Field &, const Field &, const int itime)
evolution with correction terms
Definition: fopr_NonRelativistic.cpp:756
ThreadManager::get_num_threads
static int get_num_threads()
returns available number of threads.
Definition: threadManager.cpp:246
CommonParameters::Ndim
static int Ndim()
Definition: commonParameters.h:117
fopr_NonRelativistic.h
Fopr_NonRelativistic::add_deltaH1
void add_deltaH1(Field &, const Field &, const int itime)
correction term deltaH(1): spin-magnetic interaction
Definition: fopr_NonRelativistic.cpp:791
Field::set
void set(const int jin, const int site, const int jex, double v)
Definition: field.h:175
Fopr_NonRelativistic::mult_dag
void mult_dag(Field &v, const Field &w)
not available
Definition: fopr_NonRelativistic.cpp:291
Parameters
Class for parameters.
Definition: parameters.h:46
Fopr_NonRelativistic::m_nstab
int m_nstab
stabilization parameter
Definition: fopr_NonRelativistic.h:38
Fopr_NonRelativistic::m_repr
std::string m_repr
gamma matrix repr. (set to Dirac)
Definition: fopr_NonRelativistic.h:49
Fopr_NonRelativistic::shift_zdn
void shift_zdn(Field &Xt, const Field &Yt, const int itime)
Definition: fopr_NonRelativistic.cpp:1497
Parameters::set_double
void set_double(const string &key, const double value)
Definition: parameters.cpp:33
Bridge::BridgeIO::decrease_indent
void decrease_indent()
Definition: bridgeIO.h:86
Fopr_NonRelativistic::m_Nd
int m_Nd
Definition: fopr_NonRelativistic.h:53
Fopr_NonRelativistic::m_Fstr
std::vector< Field_G > m_Fstr
Field strenth (0-2: B, 3-5: E)
Definition: fopr_NonRelativistic.h:63
Bridge::BridgeIO::increase_indent
void increase_indent()
Definition: bridgeIO.h:85
Fopr_NonRelativistic::m_u0
double m_u0
mean-field value
Definition: fopr_NonRelativistic.h:40
Bridge::BridgeIO::detailed
void detailed(const char *format,...)
Definition: bridgeIO.cpp:219
Fopr_NonRelativistic::m_boundary
std::vector< int > m_boundary
boundary conditions
Definition: fopr_NonRelativistic.h:45
Fopr_NonRelativistic::add_R2
void add_R2(Field &v, const Field &w, const int jd, const int itime)
Definition: fopr_NonRelativistic.cpp:512
fieldStrength.h
Fopr_NonRelativistic::mult_sigma2
void mult_sigma2(Field &Xt, const Field &Yt)
multiplying Pauli matrix $\sigma_2$
Definition: fopr_NonRelativistic.cpp:1585
Fopr_NonRelativistic::add_R5
void add_R5(Field &v, const Field &w, const int itime)
Definition: fopr_NonRelativistic.cpp:581
Fopr_NonRelativistic::add_R3
void add_R3(Field &v, const Field &w, const int itime)
Definition: fopr_NonRelativistic.cpp:540
CommonParameters::Nvol
static int Nvol()
Definition: commonParameters.h:109
Fopr_NonRelativistic::m_Yt2
Field m_Yt2
Definition: fopr_NonRelativistic.h:71
aypx
void aypx(const double a, Field &y, const Field &x)
aypx(y, a, x): y := a * y + x
Definition: field.cpp:509
Fopr_NonRelativistic::buf2_x
std::vector< double > buf2_x
communication buffer in x-dir.
Definition: fopr_NonRelativistic.h:65
Fopr_NonRelativistic::m_Ndim
int m_Ndim
Definition: fopr_NonRelativistic.h:55
Fopr_NonRelativistic::m_Nc
int m_Nc
Definition: fopr_NonRelativistic.h:53
Field::xI
void xI()
Definition: field.cpp:138
Fopr_NonRelativistic::shift_ydn
void shift_ydn(Field &Xt, const Field &Yt, const int itime)
Definition: fopr_NonRelativistic.cpp:1371
Fopr_NonRelativistic::mult_dn
void mult_dn(const int mu, Field &v, const Field &w)
transporter in lower direction : not available
Definition: fopr_NonRelativistic.cpp:367
axpy
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
Definition: field.cpp:380
Fopr_NonRelativistic::mult_Gn
void mult_Gn(Field &Xt, const Field &Yt, const int idir, const int itime)
nultiplication of gauge field
Definition: fopr_NonRelativistic.cpp:1639
Fopr_NonRelativistic::m_coeff
std::vector< double > m_coeff
coefficients of correction terms
Definition: fopr_NonRelativistic.h:44
Parameters::set_double_vector
void set_double_vector(const string &key, const vector< double > &value)
Definition: parameters.cpp:42
Fopr_NonRelativistic::buf1_y
std::vector< double > buf1_y
Definition: fopr_NonRelativistic.h:66
Fopr_NonRelativistic::m_Wt1
Field m_Wt1
Definition: fopr_NonRelativistic.h:70
Fopr_NonRelativistic::buf2_y
std::vector< double > buf2_y
communication buffer in y-dir.
Definition: fopr_NonRelativistic.h:66
copy
void copy(Field &y, const Field &x)
copy(y, x): y = x
Definition: field.cpp:212
Fopr_NonRelativistic::calc_Delta1
void calc_Delta1(Field &, const Field &, const int idir, const int itime)
first order covariant derivative
Definition: fopr_NonRelativistic.cpp:1048
Fopr_NonRelativistic::add_deltaH3
void add_deltaH3(Field &, const Field &, const int itime)
correction term deltaH(3)
Definition: fopr_NonRelativistic.cpp:858
Bridge::BridgeIO::paranoiac
void paranoiac(const char *format,...)
Definition: bridgeIO.cpp:238
Fopr_NonRelativistic::m_Nspc
int m_Nspc
spatial volume (m_Nx * m_Ny * m_Nz)
Definition: fopr_NonRelativistic.h:56
Fopr_NonRelativistic::get_parameters
void get_parameters(Parameters &params) const
gets parameters by a Parameter object: to be implemented in a subclass.
Definition: fopr_NonRelativistic.cpp:233
Fopr_NonRelativistic::flop_count
double flop_count()
number of floating point operations: not implemented
Definition: fopr_NonRelativistic.cpp:1757
FieldStrength
field strength construction.
Definition: fieldStrength.h:42
Fopr_NonRelativistic::set_source
void set_source(Field &Xt, const Field &b)
extract source time slice and 3D source vector from 4D source.
Definition: fopr_NonRelativistic.cpp:1073
Fopr_NonRelativistic::shift_forward
void shift_forward(Field &Xt, const Field &Yt, const int idir, const int itime)
shifting 3D field forward in idir-direction
Definition: fopr_NonRelativistic.cpp:1147
Fopr_NonRelativistic::m_Zt2
Field m_Zt2
working field on each time slice.
Definition: fopr_NonRelativistic.h:72
Field::norm2
double norm2() const
Definition: field.cpp:113
CommonParameters::Nx
static int Nx()
Definition: commonParameters.h:105
Fopr_NonRelativistic::m_vl
Bridge::VerboseLevel m_vl
verbose level
Definition: fopr_NonRelativistic.h:42
Fopr_NonRelativistic::class_name
static const std::string class_name
Definition: fopr_NonRelativistic.h:34
Fopr_NonRelativistic::mult_gm5
void mult_gm5(Field &v, const Field &w)
multiply $\gamma_5$ : not available
Definition: fopr_NonRelativistic.cpp:334
CommonParameters::Lt
static int Lt()
Definition: commonParameters.h:94
Communicator::reduce_sum
static int reduce_sum(int count, dcomplex *recv_buf, dcomplex *send_buf, int pattern=0)
make a global sum of an array of dcomplex over the communicator. pattern specifies the dimensions to ...
Definition: communicator.cpp:263
Fopr_NonRelativistic::m_evolution_type
std::string m_evolution_type
type of evolution equation
Definition: fopr_NonRelativistic.h:46
CommonParameters::Nc
static int Nc()
Definition: commonParameters.h:115
Fopr_NonRelativistic::m_correction_terms
std::string m_correction_terms
order of correction terms
Definition: fopr_NonRelativistic.h:47
Fopr_NonRelativistic::m_Yt1
Field m_Yt1
Definition: fopr_NonRelativistic.h:71
Fopr_NonRelativistic::calc_Delta2
void calc_Delta2(Field &, const Field &, const int itime)
second order covariant derivative
Definition: fopr_NonRelativistic.cpp:977
Fopr_NonRelativistic::add_R4
void add_R4(Field &v, const Field &w, const int itime)
Definition: fopr_NonRelativistic.cpp:552
Fopr_NonRelativistic::add_deltaH5
void add_deltaH5(Field &, const Field &, const int itime)
correction term deltaH(5)
Definition: fopr_NonRelativistic.cpp:950
Fopr_NonRelativistic::m_Nz
int m_Nz
Definition: fopr_NonRelativistic.h:54
CommonParameters::Nt
static int Nt()
Definition: commonParameters.h:108
Fopr_NonRelativistic::rotation
void rotation(Field &v, const Field &w, const int jd)
field rotation (jd = 1: normal, -1: conjugate)
Definition: fopr_NonRelativistic.cpp:435
Parameters::fetch_int_vector
int fetch_int_vector(const string &key, vector< int > &value) const
Definition: parameters.cpp:429
Fopr_NonRelativistic::shift_xdn
void shift_xdn(Field &Xt, const Field &Yt, const int itime)
Definition: fopr_NonRelativistic.cpp:1245
Fopr_NonRelativistic::evolve_H0
void evolve_H0(Field &, const Field &, const int itime)
evolution with kinetic term
Definition: fopr_NonRelativistic.cpp:693
Fopr_NonRelativistic::m_Yt3
Field m_Yt3
working field on each time slice.
Definition: fopr_NonRelativistic.h:71
ThreadManager::sync_barrier_all
static void sync_barrier_all()
barrier among all the threads and nodes.
Definition: threadManager.cpp:276
Communicator::npe
static int npe(const int dir)
logical grid extent
Definition: communicator.cpp:112
Fopr_NonRelativistic::evolve_typeA
void evolve_typeA(Field &, const Field &, const int itime)
evolution equation according to [1]
Definition: fopr_NonRelativistic.cpp:628
Fopr_NonRelativistic::evolve_typeB
void evolve_typeB(Field &, const Field &, const int itime)
evolution equation according to [2]
Definition: fopr_NonRelativistic.cpp:659
threadManager.h
Fopr_NonRelativistic::set_config
void set_config(Field *U)
setting gauge configuration and field strength
Definition: fopr_NonRelativistic.cpp:247
Field_G::reset
void reset(const int Nvol, const int Nex)
Definition: field_G.h:79
Fopr_NonRelativistic::m_Nx
int m_Nx
Definition: fopr_NonRelativistic.h:54
Fopr_NonRelativistic::m_itime_src
int m_itime_src
source time slice
Definition: fopr_NonRelativistic.h:59
Fopr_NonRelativistic::m_MQ
double m_MQ
heavy quark mass
Definition: fopr_NonRelativistic.h:37
Fopr_NonRelativistic::calc_Delta4
void calc_Delta4(Field &, const Field &, const int itime)
fourth order covariant derivative
Definition: fopr_NonRelativistic.cpp:1005
Field::reset
void reset(const int Nin, const int Nvol, const int Nex, const element_type cmpl=Element_type::COMPLEX)
Definition: field.h:95
Parameters::set_int_vector
void set_int_vector(const string &key, const vector< int > &value)
Definition: parameters.cpp:45
Fopr_NonRelativistic::m_Xt
std::vector< Field > m_Xt
heavy quark field.
Definition: fopr_NonRelativistic.h:69
Fopr_NonRelativistic::m_num_correct
int m_num_correct
number of correction terms
Definition: fopr_NonRelativistic.h:39
Field::cmp
double cmp(const int jin, const int site, const int jex) const
Definition: field.h:143
Field::ptr
const double * ptr(const int jin, const int site, const int jex) const
Definition: field.h:153
CommonParameters::Nd
static int Nd()
Definition: commonParameters.h:116
Fopr_NonRelativistic::m_U
Field_G m_U
gauge configuration
Definition: fopr_NonRelativistic.h:61
CommonParameters::Vlevel
static Bridge::VerboseLevel Vlevel()
Definition: commonParameters.h:122
Fopr_NonRelativistic::m_Wt3
Field m_Wt3
working field on each time slice.
Definition: fopr_NonRelativistic.h:70
Bridge::BridgeIO::set_verbose_level
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:133
Fopr_NonRelativistic::shift_yup
void shift_yup(Field &Xt, const Field &Yt, const int itime)
Definition: fopr_NonRelativistic.cpp:1306
Fopr_NonRelativistic::m_Ltime
int m_Ltime
global temporal extent
Definition: fopr_NonRelativistic.h:57
Fopr_NonRelativistic::set_parameters
void set_parameters(const Parameters &params)
seting parameters with a Parameter object
Definition: fopr_NonRelativistic.cpp:149
Parameters::set_int
void set_int(const string &key, const int value)
Definition: parameters.cpp:36
Fopr_NonRelativistic::buf1_x
std::vector< double > buf1_x
Definition: fopr_NonRelativistic.h:65
scal
void scal(Field &x, const double a)
scal(x, a): x = a * x
Definition: field.cpp:261
Communicator::ipe
static int ipe(const int dir)
logical coordinate of current proc.
Definition: communicator.cpp:105
Fopr_NonRelativistic::mult_Gd
void mult_Gd(Field &Xt, const Field &Yt, const int idir, const int itime)
nultiplication of gauge field (hermitian conjugate)
Definition: fopr_NonRelativistic.cpp:1717
Parameters::fetch_string
int fetch_string(const string &key, string &value) const
Definition: parameters.cpp:378
Parameters::fetch_double
int fetch_double(const string &key, double &value) const
Definition: parameters.cpp:327
Fopr_NonRelativistic::add_deltaH2
void add_deltaH2(Field &, const Field &, const int itime)
correction term deltaH(2)
Definition: fopr_NonRelativistic.cpp:814
Fopr_NonRelativistic::mult_sigma3
void mult_sigma3(Field &Xt, const Field &Yt)
multiplying Pauli matrix $\sigma_3$
Definition: fopr_NonRelativistic.cpp:1612
FieldStrength::construct_Fmunu_1x1
void construct_Fmunu_1x1(Field_G &Fmunu, const int mu, const int nu, const Field_G &U)
Constructs the anti-Hermitian field strength with four 1x1 plquette clover leaves.
Definition: fieldStrength.cpp:19
Fopr_NonRelativistic::mult_up
void mult_up(const int mu, Field &v, const Field &w)
transporter in upper direction : not available
Definition: fopr_NonRelativistic.cpp:358
Fopr_NonRelativistic::m_Nvol
int m_Nvol
Definition: fopr_NonRelativistic.h:55
Parameters::get_string
string get_string(const string &key) const
Definition: parameters.cpp:221
Bridge::BridgeIO::crucial
void crucial(const char *format,...)
Definition: bridgeIO.cpp:180
Fopr_NonRelativistic::m_Nvc
int m_Nvc
Definition: fopr_NonRelativistic.h:53
Fopr_NonRelativistic::m_mode
std::string m_mode
mode of multiplication
Definition: fopr_NonRelativistic.h:51
Fopr_NonRelativistic::set_mode
void set_mode(const std::string mode)
setting mult mode: 'Evolve' and 'Rotation'
Definition: fopr_NonRelativistic.cpp:264
Field
Container of Field-type object.
Definition: field.h:46
Communicator::exchange
static int exchange(int count, dcomplex *recv_buf, dcomplex *send_buf, int idir, int ipm, int tag)
receive array of dcomplex from upstream specified by idir and ipm, and send array to downstream.
Definition: communicator.cpp:207
Fopr_NonRelativistic::m_Wt2
Field m_Wt2
Definition: fopr_NonRelativistic.h:70
Fopr_NonRelativistic::evolve_U4
void evolve_U4(Field &, const Field &, const int itime)
evolution with one time slice ahead
Definition: fopr_NonRelativistic.cpp:717
ThreadManager::get_thread_id
static int get_thread_id()
returns thread id.
Definition: threadManager.cpp:253
Fopr_NonRelativistic::shift_backward
void shift_backward(Field &Xt, const Field &Yt, const int idir, const int itime)
shifting 3D field backward in idir-direction
Definition: fopr_NonRelativistic.cpp:1165
Parameters::fetch_double_vector
int fetch_double_vector(const string &key, vector< double > &value) const
Definition: parameters.cpp:410
Fopr_NonRelativistic::buf2_z
std::vector< double > buf2_z
communication buffer in z-dir.
Definition: fopr_NonRelativistic.h:67
Fopr_NonRelativistic::m_Nd2
int m_Nd2
Definition: fopr_NonRelativistic.h:53
Fopr_NonRelativistic::shift_xup
void shift_xup(Field &Xt, const Field &Yt, const int itime)
implemetation of shifting field
Definition: fopr_NonRelativistic.cpp:1183
Fopr_NonRelativistic::shift_zup
void shift_zup(Field &Xt, const Field &Yt, const int itime)
Definition: fopr_NonRelativistic.cpp:1437
Parameters::fetch_int
int fetch_int(const string &key, int &value) const
Definition: parameters.cpp:346
Bridge::BridgeIO::general
void general(const char *format,...)
Definition: bridgeIO.cpp:200
Fopr_NonRelativistic::evolve_impl
void evolve_impl(Field &, const Field &, const int itime)
called from evolve() and switchs evlution equation
Definition: fopr_NonRelativistic.cpp:610
Fopr_NonRelativistic::evolve
void evolve(Field &, const Field &)
evolution equation (facade).
Definition: fopr_NonRelativistic.cpp:376
Fopr_NonRelativistic::mult_sigma1
void mult_sigma1(Field &Xt, const Field &Yt)
multiplying Pauli matrix $\sigma_1$
Definition: fopr_NonRelativistic.cpp:1558
Fopr_NonRelativistic::m_Nt
int m_Nt
Definition: fopr_NonRelativistic.h:54
ThreadManager::assert_single_thread
static void assert_single_thread(const std::string &class_name)
assert currently running on single thread.
Definition: threadManager.cpp:372
Fopr_NonRelativistic::mult_F
void mult_F(Field &Xt, const Field &Yt, const int icomp, const int itime)
nultiplication of field strength
Definition: fopr_NonRelativistic.cpp:1678
Fopr_NonRelativistic::add_deltaH46
void add_deltaH46(Field &, const Field &, const int itime)
correction terms deltaH(4) + deltaH(6)
Definition: fopr_NonRelativistic.cpp:962
Bridge::vout
BridgeIO vout
Definition: bridgeIO.cpp:512
Bridge::BridgeIO::get_verbose_level
static std::string get_verbose_level(const VerboseLevel vl)
Definition: bridgeIO.cpp:154
Fopr_NonRelativistic::tidyup
void tidyup()
final clean up.
Definition: fopr_NonRelativistic.cpp:142