Bridge++  Ver. 2.0.2
afopr_Domainwall-tmpl.h
Go to the documentation of this file.
1 
10 #include "lib/Fopr/afopr_Domainwall.h"
11 
12 #include <stdio.h>
13 #include <stdlib.h>
14 #include <assert.h>
15 using namespace std;
16 
17 #include "lib/ResourceManager/threadManager.h"
18 #include "lib/Parameters/commonParameters.h"
19 #include "lib/Communicator/communicator.h"
20 
21 
22 template<typename AFIELD>
23 const std::string AFopr_Domainwall<AFIELD>::class_name
24  = "AFopr_Domainwall";
25 //====================================================================
26 template<typename AFIELD>
27 void AFopr_Domainwall<AFIELD>::init(const Parameters& params)
28 {
29  ThreadManager::assert_single_thread(class_name);
30 
31  m_vl = CommonParameters::Vlevel();
32 
33  vout.general(m_vl, "%s: construction\n", class_name.c_str());
34  vout.increase_indent();
35 
36  int Nc = CommonParameters::Nc();
37  int Nd = CommonParameters::Nd();
38  m_NinF = 2 * Nc * Nd;
39 
40  m_Nvol = CommonParameters::Nvol();
41  m_Ndim = CommonParameters::Ndim();
42 
43  // setup verbose level
44  string vlevel = params.get_string("verbose_level");
45  m_vl = vout.set_verbose_level(vlevel);
46 
47  int err = 0;
48 
49  // setup kernel operator
50  err += params.fetch_string("kernel_type", m_kernel_type);
51  if (err > 0) {
52  vout.crucial(m_vl, "%s: Error: kernel_type is not specified.\n",
53  class_name.c_str());
54  exit(EXIT_FAILURE);
55  }
56 
57  Parameters params_kernel = params;
58  double M0;
59  err += params.fetch_double("domain_wall_height", M0);
60  if (err > 0) {
61  vout.crucial(m_vl, "Error at %s: domain_wall_height is not specified.\n",
62  class_name.c_str());
63  exit(EXIT_FAILURE);
64  }
65  m_M0 = real_t(M0);
66 
67  double kappa = 1.0 / (8.0 - 2.0 * M0);
68  params_kernel.set_double("hopping_parameter", kappa);
69 
70  // Factory is assumed to work
71  m_foprw = AFopr<AFIELD>::New(m_kernel_type, params_kernel);
72  m_kernel_created = true;
73 
74  m_foprw->set_mode("D");
75 
76  m_Ns = 0;
77 
78  set_parameters(params);
79 
80  m_w4.reset(m_NinF, m_Nvol, 1);
81  m_v4.reset(m_NinF, m_Nvol, 1);
82  m_t4.reset(m_NinF, m_Nvol, 1);
83  m_y4.reset(m_NinF, m_Nvol, 1);
84 
85  if (needs_convert()) {
86  m_w4lex.reset(m_NinF, m_Nvol, 1);
87  m_v4lex.reset(m_NinF, m_Nvol, 1);
88  }
89 
90  vout.decrease_indent();
91  vout.general(m_vl, "%s: construction finished.\n",
92  class_name.c_str());
93 }
94 
95 
96 //====================================================================
97 template<typename AFIELD>
98 void AFopr_Domainwall<AFIELD>::init(AFOPR *fopr,
99  const Parameters& params)
100 {
101  ThreadManager::assert_single_thread(class_name);
102 
103  m_vl = CommonParameters::Vlevel();
104 
105  vout.general(m_vl, "%s: construction\n", class_name.c_str());
106  vout.increase_indent();
107 
108  int Nc = CommonParameters::Nc();
109  int Nd = CommonParameters::Nd();
110  m_NinF = 2 * Nc * Nd;
111 
112  m_Nvol = CommonParameters::Nvol();
113  m_Ndim = CommonParameters::Ndim();
114 
115  // setup verbose level
116  string vlevel = params.get_string("verbose_level");
117  m_vl = vout.set_verbose_level(vlevel);
118 
119  vout.general(m_vl, "%s: Initialization start\n", class_name.c_str());
120  int err = 0;
121 
122  m_foprw = fopr;
123  m_kernel_type = "unknown";
124  m_kernel_created = false;
125 
126  m_foprw->set_mode("D");
127 
128  m_Ns = 0;
129 
130  set_parameters(params);
131 
132  m_w4.reset(m_NinF, m_Nvol, 1);
133  m_v4.reset(m_NinF, m_Nvol, 1);
134  m_t4.reset(m_NinF, m_Nvol, 1);
135  m_y4.reset(m_NinF, m_Nvol, 1);
136 
137  if (needs_convert()) {
138  m_w4lex.reset(m_NinF, m_Nvol, 1);
139  m_v4lex.reset(m_NinF, m_Nvol, 1);
140  }
141 
142  vout.decrease_indent();
143  vout.general(m_vl, "%s: construction finished.\n",
144  class_name.c_str());
145 }
146 
147 
148 //====================================================================
149 template<typename AFIELD>
150 void AFopr_Domainwall<AFIELD>::init(AFOPR *fopr)
151 {
152  ThreadManager::assert_single_thread(class_name);
153 
154  m_vl = CommonParameters::Vlevel();
155 
156  vout.general(m_vl, "%s: construction\n", class_name.c_str());
157  vout.increase_indent();
158 
159  int Nc = CommonParameters::Nc();
160  int Nd = CommonParameters::Nd();
161  m_NinF = 2 * Nc * Nd;
162 
163  m_Nvol = CommonParameters::Nvol();
164  m_Ndim = CommonParameters::Ndim();
165 
166  vout.general(m_vl, "%s: Initialization start\n", class_name.c_str());
167  int err = 0;
168 
169  m_foprw = fopr;
170  m_kernel_type = "unknown";
171  m_kernel_created = false;
172 
173  m_foprw->set_mode("D");
174 
175  m_Ns = 0;
176 
177  m_w4.reset(m_NinF, m_Nvol, 1);
178  m_v4.reset(m_NinF, m_Nvol, 1);
179  m_t4.reset(m_NinF, m_Nvol, 1);
180  m_y4.reset(m_NinF, m_Nvol, 1);
181 
182  if (needs_convert()) {
183  m_w4lex.reset(m_NinF, m_Nvol, 1);
184  m_v4lex.reset(m_NinF, m_Nvol, 1);
185  }
186 
187  vout.decrease_indent();
188  vout.general(m_vl, "%s: construction finished.\n",
189  class_name.c_str());
190 }
191 
192 
193 //====================================================================
194 template<typename AFIELD>
195 void AFopr_Domainwall<AFIELD>::tidyup()
196 {
197  if (m_kernel_created == true) delete m_foprw;
198 }
199 
200 
201 //====================================================================
202 template<typename AFIELD>
203 void AFopr_Domainwall<AFIELD>::set_parameters(const Parameters& params)
204 {
205  std::string vlevel;
206  if (!params.fetch_string("verbose_level", vlevel)) {
207  m_vl = vout.set_verbose_level(vlevel);
208  }
209 
210  //- fetch and check input parameters
211  string gmset_type;
212  double mq, M0;
213  int Ns;
214  std::vector<int> bc;
215  double b, c;
216 
217  int err_optional = 0;
218  err_optional += params.fetch_string("gamma_matrix_type", gmset_type);
219 
220  int err = 0;
221  err += params.fetch_double("quark_mass", mq);
222  err += params.fetch_double("domain_wall_height", M0);
223  err += params.fetch_int("extent_of_5th_dimension", Ns);
224  err += params.fetch_int_vector("boundary_condition", bc);
225 
226  if (err) {
227  vout.crucial(m_vl, "Error at %s: input parameter not found.\n",
228  class_name.c_str());
229  exit(EXIT_FAILURE);
230  }
231 
232  int err2 = 0;
233  err2 += params.fetch_double("coefficient_b", b);
234  err2 += params.fetch_double("coefficient_c", c);
235 
236  if (err2) {
237  vout.general(m_vl, " coefficients b, c are not provided:"
238  " set to Shamir's form.\n");
239  b = 1.0;
240  c = 0.0;
241  }
242 
243  set_parameters(real_t(mq), real_t(M0), Ns, bc,
244  real_t(b), real_t(c));
245 
246  if (real_t(M0) != m_M0) set_kernel_parameters(params);
247 }
248 
249 
250 //====================================================================
251 template<typename AFIELD>
252 void AFopr_Domainwall<AFIELD>::get_parameters(Parameters& params) const
253 {
254  params.set_string("kernel_type", m_kernel_type);
255  params.set_double("quark_mass", double(m_mq));
256  params.set_double("domain_wall_height", double(m_M0));
257  params.set_int("extent_of_5th_dimension", m_Ns);
258  params.set_int_vector("boundary_condition", m_boundary);
259  params.set_double("coefficient_b", double(m_b[0]));
260  params.set_double("coefficient_c", double(m_c[0]));
261  params.set_string("gamma_matrix_type", m_repr);
262 
263  params.set_string("verbose_level", vout.get_verbose_level(m_vl));
264 }
265 
266 
267 //====================================================================
268 template<typename AFIELD>
269 void AFopr_Domainwall<AFIELD>::set_parameters(
270  const real_t mq,
271  const real_t M0,
272  const int Ns,
273  const std::vector<int> bc,
274  const real_t b,
275  const real_t c)
276 {
277  int ith = ThreadManager::get_thread_id();
278 
279  if (ith == 0) {
280  m_M0 = M0;
281  m_mq = mq;
282  m_Ns = Ns;
283 
284  assert(bc.size() == m_Ndim);
285  if (m_boundary.size() != m_Ndim) m_boundary.resize(m_Ndim);
286  for (int mu = 0; mu < m_Ndim; ++mu) {
287  m_boundary[mu] = bc[mu];
288  }
289 
290  if (m_b.size() != m_Ns) {
291  m_b.resize(m_Ns);
292  m_c.resize(m_Ns);
293  }
294  for (int is = 0; is < m_Ns; ++is) {
295  m_b[is] = real_t(b);
296  m_c[is] = real_t(c);
297  }
298  }
299 
300  vout.general(m_vl, "%s: input parameters\n", class_name.c_str());
301  vout.general(m_vl, " mq = %8.4f\n", m_mq);
302  vout.general(m_vl, " M0 = %8.4f\n", m_M0);
303  vout.general(m_vl, " Ns = %4d\n", m_Ns);
304  for (int mu = 0; mu < m_Ndim; ++mu) {
305  vout.general(m_vl, " boundary[%d] = %2d\n", mu, m_boundary[mu]);
306  }
307  vout.general(m_vl, " coefficients:\n");
308  for (int is = 0; is < m_Ns; ++is) {
309  vout.general(m_vl, " b[%2d] = %16.10f c[%2d] = %16.10f\n",
310  is, m_b[is], is, m_c[is]);
311  }
312 
313  set_precond_parameters();
314 
315  // working 5d vectors.
316  if (m_w1.nex() != Ns) {
317  m_w1.reset(m_NinF, m_Nvol, m_Ns);
318  m_v1.reset(m_NinF, m_Nvol, m_Ns);
319  m_v2.reset(m_NinF, m_Nvol, m_Ns);
320  }
321 }
322 
323 
324 //====================================================================
325 template<typename AFIELD>
326 void AFopr_Domainwall<AFIELD>::set_parameters(
327  const real_t mq,
328  const real_t M0,
329  const int Ns,
330  const std::vector<int> bc,
331  const std::vector<real_t> vec_b,
332  const std::vector<real_t> vec_c)
333 {
334  int ith = ThreadManager::get_thread_id();
335 
336  if (ith == 0) {
337  m_M0 = M0;
338  m_mq = mq;
339  m_Ns = Ns;
340 
341  assert(bc.size() == m_Ndim);
342  if (m_boundary.size() != m_Ndim) m_boundary.resize(m_Ndim);
343  for (int mu = 0; mu < m_Ndim; ++mu) {
344  m_boundary[mu] = bc[mu];
345  }
346 
347  if (m_b.size() != m_Ns) {
348  m_b.resize(m_Ns);
349  m_c.resize(m_Ns);
350  }
351  for (int is = 0; is < m_Ns; ++is) {
352  m_b[is] = real_t(vec_b[is]);
353  m_c[is] = real_t(vec_c[is]);
354  }
355  }
356 
357 
358  vout.general(m_vl, "%s: parameters\n", class_name.c_str());
359  vout.general(m_vl, " mq = %8.4f\n", m_mq);
360  vout.general(m_vl, " M0 = %8.4f\n", m_M0);
361  vout.general(m_vl, " Ns = %4d\n", m_Ns);
362  for (int mu = 0; mu < m_Ndim; ++mu) {
363  vout.general(m_vl, " boundary[%d] = %2d\n", mu, m_boundary[mu]);
364  }
365  vout.general(m_vl, " coefficients:\n");
366  for (int is = 0; is < m_Ns; ++is) {
367  vout.general(m_vl, " b[%2d] = %16.10f c[%2d] = %16.10f\n",
368  is, m_b[is], is, m_c[is]);
369  }
370 
371  set_precond_parameters();
372 
373  // working 5d vectors.
374  if (m_w1.nex() != Ns) {
375  m_w1.reset(m_NinF, m_Nvol, m_Ns);
376  m_v1.reset(m_NinF, m_Nvol, m_Ns);
377  m_v2.reset(m_NinF, m_Nvol, m_Ns);
378  }
379 
380 }
381 
382 
383 //====================================================================
384 template<typename AFIELD>
385 void AFopr_Domainwall<AFIELD>::set_coefficients(
386  const std::vector<real_t> vec_b,
387  const std::vector<real_t> vec_c)
388 {
389  if ((vec_b.size() != m_Ns) || (vec_c.size() != m_Ns)) {
390  vout.crucial(m_vl, "%s: size of coefficient vectors incorrect.\n",
391  class_name.c_str());
392  }
393 
394  vout.general(m_vl, "%s: coefficient vectors are set:\n",
395  class_name.c_str());
396 
397  for (int is = 0; is < m_Ns; ++is) {
398  m_b[is] = real_t(vec_b[is]);
399  m_c[is] = real_t(vec_c[is]);
400  vout.general(m_vl, "b[%2d] = %16.10f c[%2d] = %16.10f\n",
401  is, m_b[is], is, m_c[is]);
402  }
403 
404  set_precond_parameters();
405 }
406 
407 
408 //====================================================================
409 template<typename AFIELD>
410 void AFopr_Domainwall<AFIELD>::set_kernel_parameters(
411  const Parameters& params)
412 {
413  Parameters params_kernel = params;
414 
415  double M0;
416  params.fetch_double("domain_wall_height", M0);
417 
418  double kappa = 1.0 / (8.0 - 2.0 * M0);
419  params_kernel.set_double("hopping_parameter", kappa);
420 
421  m_foprw->set_parameters(params_kernel);
422 }
423 
424 
425 //====================================================================
426 template<typename AFIELD>
427 void AFopr_Domainwall<AFIELD>::set_precond_parameters()
428 {
429  if (m_dp.size() != m_Ns) {
430  m_dp.resize(m_Ns);
431  m_dm.resize(m_Ns);
432  m_e.resize(m_Ns - 1);
433  m_f.resize(m_Ns - 1);
434  }
435 
436  for (int is = 0; is < m_Ns; ++is) {
437  m_dp[is] = 1.0 + m_b[is] * (4.0 - m_M0);
438  m_dm[is] = 1.0 - m_c[is] * (4.0 - m_M0);
439  }
440 
441  m_e[0] = m_mq * m_dm[m_Ns - 1] / m_dp[0];
442  m_f[0] = m_mq * m_dm[0];
443  for (int is = 1; is < m_Ns - 1; ++is) {
444  m_e[is] = m_e[is - 1] * m_dm[is - 1] / m_dp[is];
445  m_f[is] = m_f[is - 1] * m_dm[is] / m_dp[is - 1];
446  }
447 
448  m_g = m_e[m_Ns - 2] * m_dm[m_Ns - 2];
449 }
450 
451 
452 //====================================================================
453 template<typename AFIELD>
454 void AFopr_Domainwall<AFIELD>::convert(AFIELD& v, const Field& w)
455 {
456  if (!needs_convert()) {
457  vout.crucial(m_vl, "%s: convert is not necessary.\n",
458  class_name.c_str());
459  exit(EXIT_FAILURE);
460  }
461 
462 #pragma omp barrier
463 
464  int Nex = w.nex();
465  for (int ex = 0; ex < Nex; ++ex) {
466  copy(m_w4lex, 0, w, ex);
467  m_foprw->convert(m_v4lex, m_w4lex);
468  copy(v, ex, m_v4lex, 0);
469  }
470 
471 #pragma omp barrier
472 }
473 
474 
475 //====================================================================
476 template<typename AFIELD>
477 void AFopr_Domainwall<AFIELD>::reverse(Field& v, const AFIELD& w)
478 {
479  if (!needs_convert()) {
480  vout.crucial(m_vl, "%s: convert is not necessary.\n",
481  class_name.c_str());
482  exit(EXIT_FAILURE);
483  }
484 
485 #pragma omp barrier
486 
487  int Nex = w.nex();
488  for (int ex = 0; ex < Nex; ++ex) {
489  copy(m_v4lex, 0, w, ex);
490  m_foprw->reverse(m_w4lex, m_v4lex);
491  copy(v, ex, m_w4lex, 0);
492  }
493 
494 #pragma omp barrier
495 }
496 
497 
498 //====================================================================
499 template<typename AFIELD>
500 void AFopr_Domainwall<AFIELD>::set_mode(std::string mode)
501 {
502 #pragma omp barrier
503 
504  int ith = ThreadManager::get_thread_id();
505  if (ith == 0) m_mode = mode;
506 
507 #pragma omp barrier
508 }
509 
510 
511 //====================================================================
512 template<typename AFIELD>
513 void AFopr_Domainwall<AFIELD>::mult(AFIELD& v, const AFIELD& w)
514 {
515  if (m_mode == "D") {
516  D(v, w);
517  } else if (m_mode == "Ddag") {
518  Ddag(v, w);
519  } else if (m_mode == "DdagD") {
520  DdagD(v, w);
521  } else if (m_mode == "DDdag") {
522  DDdag(v, w);
523  } else if (m_mode == "H") {
524  H(v, w);
525  } else if (m_mode == "Hdag") {
526  Hdag(v, w);
527  } else if (m_mode == "D_prec") {
528  D_prec(v, w);
529  } else if (m_mode == "Ddag_prec") {
530  Ddag_prec(v, w);
531  } else if (m_mode == "DdagD_prec") {
532  DdagD_prec(v, w);
533  } else if (m_mode == "Prec") {
534  Prec(v, w);
535  } else {
536  vout.crucial(m_vl, "mode undeifined in %s.\n", class_name.c_str());
537  exit(EXIT_FAILURE);
538  }
539 }
540 
541 
542 //====================================================================
543 template<typename AFIELD>
544 void AFopr_Domainwall<AFIELD>::mult_dag(AFIELD& v, const AFIELD& w)
545 {
546  if (m_mode == "D") {
547  Ddag(v, w);
548  } else if (m_mode == "Ddag") {
549  D(v, w);
550  } else if (m_mode == "DdagD") {
551  DdagD(v, w);
552  } else if (m_mode == "DDdag") {
553  DDdag(v, w);
554  } else if (m_mode == "H") {
555  Hdag(v, w);
556  } else if (m_mode == "Hdag") {
557  H(v, w);
558  } else if (m_mode == "D_prec") {
559  Ddag_prec(v, w);
560  } else if (m_mode == "Ddag_prec") {
561  D_prec(v, w);
562  } else if (m_mode == "DdagD_prec") {
563  DdagD_prec(v, w);
564  } else if (m_mode == "Prec") {
565  Precdag(v, w);
566  } else {
567  vout.crucial(m_vl, "mode undeifined in %s.\n", class_name.c_str());
568  exit(EXIT_FAILURE);
569  }
570 }
571 
572 
573 //====================================================================
574 template<typename AFIELD>
575 void AFopr_Domainwall<AFIELD>::mult(AFIELD& v, const AFIELD& w,
576  std::string mode)
577 {
578  assert(w.check_size(m_NinF, m_Nvol, m_Ns));
579  assert(v.check_size(m_NinF, m_Nvol, m_Ns));
580 
581  if (mode == "Prec") {
582  Prec(v, w);
583  } else if (mode == "Precdag") {
584  Precdag(v, w);
585  } else if (mode == "D") {
586  D(v, w);
587  } else if (mode == "Ddag") {
588  Ddag(v, w);
589  } else if (mode == "DdagD") {
590  DdagD(v, w);
591  } else if (mode == "DDdag") {
592  DDdag(v, w);
593  } else if (mode == "D_prec") {
594  D_prec(v, w);
595  } else if (mode == "Ddag_prec") {
596  Ddag_prec(v, w);
597  } else if (mode == "DdagD_prec") {
598  DdagD_prec(v, w);
599  } else {
600  vout.crucial(m_vl, "%s: undefined mode = %s\n",
601  class_name.c_str(), mode.c_str());
602  exit(EXIT_FAILURE);
603  }
604 }
605 
606 
607 //====================================================================
608 template<typename AFIELD>
609 void AFopr_Domainwall<AFIELD>::mult_dag(AFIELD& v, const AFIELD& w,
610  std::string mode)
611 {
612  assert(w.check_size(m_NinF, m_Nvol, m_Ns));
613  assert(v.check_size(m_NinF, m_Nvol, m_Ns));
614 
615  if (mode == "Prec") {
616  Precdag(v, w);
617  } else if (mode == "Precdag") {
618  Prec(v, w);
619  } else if (mode == "D") {
620  Ddag(v, w);
621  } else if (mode == "Ddag") {
622  D(v, w);
623  } else if (mode == "DdagD") {
624  DdagD(v, w);
625  } else if (mode == "DDdag") {
626  DDdag(v, w);
627  } else if (mode == "D_prec") {
628  Ddag_prec(v, w);
629  } else if (mode == "Ddag_prec") {
630  D_prec(v, w);
631  } else if (mode == "DdagD_prec") {
632  DdagD_prec(v, w);
633  } else {
634  std::cout << "mode undeifined in AFopr_Domainwall.\n";
635  abort();
636  }
637 }
638 
639 
640 //====================================================================
641 template<typename AFIELD>
642 void AFopr_Domainwall<AFIELD>::mult_gm5(AFIELD& v, const AFIELD& w)
643 {
644  int Nex = w.nex();
645  assert(w.check_size(m_NinF, m_Nvol, Nex));
646  assert(v.check_size(m_NinF, m_Nvol, Nex));
647 
648  // omp barrier at the beginning and end of m_foprw->mult_gm5 is
649  // assumed.
650  if (Nex == 1) {
651  m_foprw->mult_gm5(v, w);
652  } else {
653 #pragma omp barrier
654  for (int ex = 0; ex < Nex; ++ex) {
655  copy(m_w4, 0, w, ex);
656  m_foprw->mult_gm5(m_v4, m_w4);
657  copy(v, ex, m_v4, 0);
658 #pragma omp barrier
659  }
660  }
661 }
662 
663 
664 //====================================================================
665 template<typename AFIELD>
666 void AFopr_Domainwall<AFIELD>::mult_chproj_4d(AFIELD& v,
667  const AFIELD& w,
668  const int ipm)
669 {
670 #pragma omp barrier
671 
672  copy(v, w);
673 
674  m_foprw->mult_gm5(m_w4, w);
675 
676  axpy(v, real_t(ipm), m_w4);
677  scal(v, real_t(0.5));
678 
679 #pragma omp barrier
680 }
681 
682 
683 //====================================================================
684 template<typename AFIELD>
685 void AFopr_Domainwall<AFIELD>::mult_gm5_4d(AFIELD& v, const AFIELD& w)
686 {
687  m_foprw->mult_gm5(v, w);
688 }
689 
690 
691 //====================================================================
692 template<typename AFIELD>
693 void AFopr_Domainwall<AFIELD>::DdagD(AFIELD& v, const AFIELD& w)
694 {
695  assert(w.check_size(m_NinF, m_Nvol, m_Ns));
696  assert(v.check_size(m_NinF, m_Nvol, m_Ns));
697 
698  D(m_v1, w);
699  Ddag(v, m_v1);
700 }
701 
702 
703 //====================================================================
704 template<typename AFIELD>
705 void AFopr_Domainwall<AFIELD>::DDdag(AFIELD& v, const AFIELD& w)
706 {
707  assert(w.check_size(m_NinF, m_Nvol, m_Ns));
708  assert(v.check_size(m_NinF, m_Nvol, m_Ns));
709 
710  Ddag(m_v1, w);
711  D(v, m_v1);
712 }
713 
714 
715 //====================================================================
716 template<typename AFIELD>
717 void AFopr_Domainwall<AFIELD>::DdagD_prec(AFIELD& v, const AFIELD& w)
718 {
719  assert(w.check_size(m_NinF, m_Nvol, m_Ns));
720  assert(v.check_size(m_NinF, m_Nvol, m_Ns));
721 
722  D_prec(m_v1, w);
723  Ddag_prec(v, m_v1);
724 }
725 
726 
727 //====================================================================
728 template<typename AFIELD>
729 void AFopr_Domainwall<AFIELD>::D_prec(AFIELD& v, const AFIELD& w)
730 {
731  L_inv(v, w);
732  U_inv(m_v2, v);
733  D(v, m_v2);
734 }
735 
736 
737 //====================================================================
738 template<typename AFIELD>
739 void AFopr_Domainwall<AFIELD>::Ddag_prec(AFIELD& v, const AFIELD& w)
740 {
741  Ddag(v, w);
742  Udag_inv(m_v2, v);
743  Ldag_inv(v, m_v2);
744 }
745 
746 
747 //====================================================================
748 template<typename AFIELD>
749 void AFopr_Domainwall<AFIELD>::Prec(AFIELD& v, const AFIELD& w)
750 {
751  L_inv(m_v2, w);
752  U_inv(v, m_v2);
753 }
754 
755 
756 //====================================================================
757 template<typename AFIELD>
758 void AFopr_Domainwall<AFIELD>::Precdag(AFIELD& v, const AFIELD& w)
759 {
760  Udag_inv(m_v2, w);
761  Ldag_inv(v, m_v2);
762 }
763 
764 
765 //====================================================================
766 template<typename AFIELD>
767 void AFopr_Domainwall<AFIELD>::H(AFIELD& v, const AFIELD& w)
768 {
769  D(m_v2, w);
770  mult_gm5R(v, m_v2);
771 }
772 
773 
774 //====================================================================
775 template<typename AFIELD>
776 void AFopr_Domainwall<AFIELD>::Hdag(AFIELD& v, const AFIELD& w)
777 {
778  mult_gm5R(m_v2, w);
779  Ddag(v, m_v2);
780 }
781 
782 
783 //====================================================================
784 template<typename AFIELD>
785 void AFopr_Domainwall<AFIELD>::mult_gm5R(AFIELD& v, const AFIELD& w)
786 {
787 #pragma omp barrier
788 
789  for (int is = 0; is < m_Ns; ++is) {
790  copy(m_w4, 0, w, is);
791  mult_gm5_4d(m_v4, m_w4);
792  copy(v, m_Ns - 1 - is, m_v4, 0);
793  }
794 
795 #pragma omp barrier
796 }
797 
798 
799 //====================================================================
800 template<typename AFIELD>
801 void AFopr_Domainwall<AFIELD>::mult_R(AFIELD& v, const AFIELD& w)
802 {
803  assert(w.check_size(m_NinF, m_Nvol, m_Ns));
804  assert(v.check_size(m_NinF, m_Nvol, m_Ns));
805 
806 #pragma omp barrier
807 
808  for (int is = 0; is < m_Ns; ++is) {
809  copy(v, m_Ns - 1 - is, w, is);
810  }
811 
812 #pragma omp barrier
813 }
814 
815 
816 //====================================================================
817 template<typename AFIELD>
818 void AFopr_Domainwall<AFIELD>::D(AFIELD& v, const AFIELD& w)
819 {
820 #pragma omp barrier
821 
822  for (int is = 0; is < m_Ns; ++is) {
823  m_y4.set(0.0);
824  int is_up = (is + 1) % m_Ns;
825  real_t Fup = -0.5;
826  if (is == m_Ns - 1) Fup = 0.5 * m_mq;
827  copy(m_v4, 0, w, is_up);
828  mult_gm5_4d(m_t4, m_v4);
829  axpy(m_v4, real_t(-1.0), m_t4);
830  axpy(m_y4, 0, Fup, m_v4, 0);
831 
832  int is_dn = (is - 1 + m_Ns) % m_Ns;
833  real_t Fdn = -0.5;
834  if (is == 0) Fdn = 0.5 * m_mq;
835  copy(m_v4, 0, w, is_dn);
836  mult_gm5_4d(m_t4, m_v4);
837  axpy(m_v4, real_t(1.0), m_t4);
838  axpy(m_y4, 0, Fdn, m_v4, 0);
839 
840  copy(m_w4, 0, w, is);
841 
842  copy(v, is, m_w4, 0);
843  axpy(v, is, real_t(1.0), m_y4, 0);
844 
845  scal(m_w4, m_b[is]);
846  axpy(m_w4, -m_c[is], m_y4);
847 
848  m_foprw->mult(m_v4, m_w4);
849 
850  axpy(v, is, real_t(4.0) - m_M0, m_v4, 0);
851  }
852 
853 #pragma omp barrier
854 }
855 
856 
857 //====================================================================
858 template<typename AFIELD>
859 void AFopr_Domainwall<AFIELD>::Ddag(AFIELD& v, const AFIELD& w)
860 {
861 #pragma omp barrier
862 
863  v.set(0.0);
864 #pragma omp barrier
865 
866  for (int is = 0; is < m_Ns; ++is) {
867  copy(m_w4, 0, w, is);
868  copy(m_y4, m_w4);
869  m_foprw->mult_dag(m_v4, m_w4);
870 
871  axpy(m_w4, m_b[is] * (real_t(4.0) - m_M0), m_v4);
872  axpy(v, is, real_t(1.0), m_w4, 0);
873 
874  axpy(m_y4, -m_c[is] * (real_t(4.0) - m_M0), m_v4);
875 
876  int is_up = (is + 1) % m_Ns;
877  real_t Fup = -0.5;
878  if (is_up == 0) Fup = 0.5 * m_mq;
879  mult_gm5_4d(m_t4, m_y4);
880  axpy(m_t4, real_t(-1.0), m_y4); // m_t4 = - (1 - gm5) * m_y4
881  axpy(v, is_up, -Fup, m_t4, 0);
882 
883  int is_dn = (is - 1 + m_Ns) % m_Ns;
884  real_t Fdn = -0.5;
885  if (is_dn == m_Ns - 1) Fdn = 0.5 * m_mq;
886  mult_gm5_4d(m_t4, m_y4);
887  axpy(m_t4, real_t(1.0), m_y4); // m_t4 = (1 + gm5) * m_y4
888  axpy(v, is_dn, Fdn, m_t4, 0);
889  }
890 
891 #pragma omp barrier
892 }
893 
894 
895 //====================================================================
896 template<typename AFIELD>
897 void AFopr_Domainwall<AFIELD>::L_inv(AFIELD& v, const AFIELD& w)
898 {
899 #pragma omp barrier
900 
901  copy(v, 0, w, 0);
902  copy(m_y4, 0, w, 0);
903  scal(m_y4, m_e[0]);
904 
905  for (int is = 1; is < m_Ns - 1; ++is) {
906  copy(v, is, w, is);
907 
908  copy(m_v4, 0, v, is - 1);
909  mult_gm5_4d(m_t4, m_v4);
910  axpy(m_v4, real_t(1.0), m_t4);
911  scal(m_v4, real_t(0.5) * m_dm[is] / m_dp[is - 1]);
912 
913  axpy(v, is, real_t(1.0), m_v4, 0);
914  copy(m_t4, 0, v, is);
915  axpy(m_y4, m_e[is], m_t4);
916  }
917 
918  int is = m_Ns - 1;
919  copy(v, is, w, is);
920  copy(m_v4, 0, v, is - 1);
921  mult_gm5_4d(m_t4, m_v4);
922  axpy(m_v4, real_t(1.0), m_t4);
923  scal(m_v4, real_t(0.5) * m_dm[is] / m_dp[is - 1]);
924 
925  axpy(v, is, real_t(1.0), m_v4, 0);
926 
927  mult_gm5_4d(m_t4, m_y4);
928  axpy(m_y4, real_t(-1.0), m_t4);
929  scal(m_y4, real_t(-0.5));
930  axpy(v, is, real_t(1.0), m_y4, 0);
931 
932 #pragma omp barrier
933 }
934 
935 
936 //====================================================================
937 template<typename AFIELD>
938 void AFopr_Domainwall<AFIELD>::U_inv(AFIELD& v, const AFIELD& w)
939 {
940  v.set(0.0);
941 
942 #pragma omp barrier
943 
944  int is = m_Ns - 1;
945  copy(m_y4, 0, w, is);
946  scal(m_y4, real_t(1.0) / (m_dp[is] + m_g));
947  copy(v, is, m_y4, 0);
948 
949  mult_gm5_4d(m_w4, m_y4);
950  axpy(m_y4, real_t(1.0), m_w4);
951  scal(m_y4, real_t(0.5));
952 
953  for (int is = m_Ns - 2; is >= 0; --is) {
954  copy(m_v4, 0, w, is);
955 
956  copy(m_w4, 0, v, is + 1);
957  mult_gm5_4d(m_t4, m_w4);
958  axpy(m_w4, real_t(-1.0), m_t4);
959 
960  axpy(m_v4, real_t(0.5) * m_dm[is], m_w4);
961 
962  axpy(m_v4, -m_f[is], m_y4);
963 
964  scal(m_v4, real_t(1.0) / m_dp[is]);
965 
966  copy(v, is, m_v4, 0);
967  }
968 
969 #pragma omp barrier
970 }
971 
972 
973 //====================================================================
974 template<typename AFIELD>
975 void AFopr_Domainwall<AFIELD>::Udag_inv(AFIELD& v, const AFIELD& w)
976 {
977 #pragma omp barrier
978 
979  copy(m_v4, 0, w, 0);
980  scal(m_v4, real_t(1.0) / m_dp[0]);
981  copy(v, 0, m_v4, 0);
982 
983  copy(m_y4, m_v4);
984  scal(m_y4, m_f[0]);
985 
986  for (int is = 1; is < m_Ns - 1; ++is) {
987  copy(m_t4, 0, w, is);
988 
989  copy(m_v4, 0, v, is - 1);
990  mult_gm5_4d(m_w4, m_v4);
991  axpy(m_v4, real_t(-1.0), m_w4);
992  axpy(m_t4, real_t(0.5) * m_dm[is - 1], m_v4);
993 
994  scal(m_t4, real_t(1.0) / m_dp[is]);
995  copy(v, is, m_t4, 0);
996 
997  axpy(m_y4, m_f[is], m_t4);
998  }
999 
1000  int is = m_Ns - 1;
1001 
1002  copy(m_t4, 0, w, is);
1003 
1004  copy(m_v4, 0, v, is - 1);
1005  mult_gm5_4d(m_w4, m_v4);
1006  axpy(m_v4, real_t(-1.0), m_w4);
1007  axpy(m_t4, real_t(0.5) * m_dm[is - 1], m_v4);
1008 
1009  mult_gm5_4d(m_w4, m_y4);
1010  axpy(m_y4, real_t(1.0), m_w4);
1011  scal(m_y4, real_t(0.5));
1012 
1013  axpy(m_t4, real_t(-1.0), m_y4);
1014 
1015  scal(m_t4, real_t(1.0) / (m_dp[is] + m_g));
1016  copy(v, is, m_t4, 0);
1017 
1018 #pragma omp barrier
1019 }
1020 
1021 
1022 //====================================================================
1023 template<typename AFIELD>
1024 void AFopr_Domainwall<AFIELD>::Ldag_inv(AFIELD& v, const AFIELD& w)
1025 {
1026  int is = m_Ns - 1;
1027 
1028 #pragma omp barrier
1029 
1030  copy(v, is, w, is);
1031 
1032  copy(m_y4, 0, w, is);
1033  mult_gm5_4d(m_w4, m_y4);
1034  axpy(m_y4, real_t(-1.0), m_w4);
1035  scal(m_y4, real_t(0.5));
1036 
1037  for (int is = m_Ns - 2; is >= 0; --is) {
1038  copy(v, is, w, is);
1039 
1040  copy(m_v4, 0, v, is + 1);
1041  mult_gm5_4d(m_w4, m_v4);
1042  axpy(m_v4, real_t(1.0), m_w4);
1043  scal(m_v4, real_t(0.5) * m_dm[is + 1] / m_dp[is]);
1044 
1045  axpy(m_v4, -m_e[is], m_y4);
1046 
1047  axpy(v, is, real_t(1.0), m_v4, 0);
1048  }
1049 
1050 #pragma omp barrier
1051 }
1052 
1053 
1054 //====================================================================
1055 template<typename AFIELD>
1056 double AFopr_Domainwall<AFIELD>::flop_count(std::string mode)
1057 {
1058  int Lvol = CommonParameters::Lvol();
1059  double vsite = static_cast<double>(Lvol);
1060  double vNs = static_cast<double>(m_Ns);
1061 
1062  // double flop_Wilson = m_foprw->flop_count("D");
1063  double flop_Wilson = m_foprw->flop_count();
1064 
1065  double axpy1 = static_cast<double>(2 * m_NinF);
1066  double scal1 = static_cast<double>(1 * m_NinF);
1067 
1068  double flop_DW = vNs * (flop_Wilson + vsite * (6 * axpy1 + 2 * scal1));
1069  // In Ddag case, flop_Wilson + 7 axpy which equals flop_DW.
1070 
1071  double flop_LU_inv = 2.0 * vsite *
1072  ((3.0 * axpy1 + scal1) * (vNs - 1.0) + axpy1 + 2.0 * scal1);
1073 
1074  double flop = 0.0;
1075  if (mode == "Prec") {
1076  flop = flop_LU_inv;
1077  } else if ((mode == "D") || (mode == "Ddag")) {
1078  flop = flop_DW;
1079  } else if (mode == "DdagD") {
1080  flop = 2.0 * flop_DW;
1081  } else if ((mode == "D_prec") || (mode == "Ddag_prec")) {
1082  flop = flop_LU_inv + flop_DW;
1083  } else if (mode == "DdagD_prec") {
1084  flop = 2.0 * (flop_LU_inv + flop_DW);
1085  } else {
1086  vout.crucial(m_vl, "Error at %s: input repr is undefined.\n",
1087  class_name.c_str());
1088  exit(EXIT_FAILURE);
1089  }
1090 
1091  return flop;
1092 }
1093 
1094 
1095 //============================================================END=====
AFopr_Domainwall::Precdag
void Precdag(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall-tmpl.h:758
AFopr_Domainwall::mult_gm5
void mult_gm5(AFIELD &, const AFIELD &)
multiplies gamma_5 matrix.
Definition: afopr_Domainwall-tmpl.h:642
CommonParameters::Lvol
static long_t Lvol()
Definition: commonParameters.h:95
AFopr_Domainwall::DDdag
void DDdag(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall-tmpl.h:705
Parameters::set_string
void set_string(const string &key, const string &value)
Definition: parameters.cpp:39
AFopr
Definition: afopr.h:48
AFopr_Domainwall::U_inv
void U_inv(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall-tmpl.h:938
CommonParameters::Ndim
static int Ndim()
Definition: commonParameters.h:117
Field::set
void set(const int jin, const int site, const int jex, double v)
Definition: field.h:175
Parameters
Class for parameters.
Definition: parameters.h:46
Parameters::set_double
void set_double(const string &key, const double value)
Definition: parameters.cpp:33
AFopr_Domainwall::Udag_inv
void Udag_inv(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall-tmpl.h:975
Bridge::BridgeIO::decrease_indent
void decrease_indent()
Definition: bridgeIO.h:86
AFopr_Domainwall::DdagD
void DdagD(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall-tmpl.h:693
AFopr_Domainwall::reverse
void reverse(Field &, const AFIELD &)
reverse AField to Field.
Definition: afopr_Domainwall-tmpl.h:477
Bridge::BridgeIO::increase_indent
void increase_indent()
Definition: bridgeIO.h:85
Field::nex
int nex() const
Definition: field.h:128
Field::check_size
bool check_size(const int nin, const int nvol, const int nex) const
checking size parameters. [23 May 2016 H.Matsufuru]
Definition: field.h:135
CommonParameters::Nvol
static int Nvol()
Definition: commonParameters.h:109
AFopr_Domainwall
Domain-wall fermion operator.
Definition: afopr_Domainwall.h:38
axpy
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
Definition: field.cpp:380
real_t
double real_t
Definition: bridgeQXS_Clover_coarse_double.cpp:16
AFopr_Domainwall::tidyup
void tidyup()
final tidyup.
Definition: afopr_Domainwall-tmpl.h:195
AFopr::set_mode
virtual void set_mode(std::string mode)
setting the mode of multiplication if necessary. Default implementation here is just to avoid irrelev...
Definition: afopr.h:81
AFopr_Domainwall::mult
void mult(AFIELD &v, const AFIELD &w)
multiplies fermion operator to a given field.
Definition: afopr_Domainwall-tmpl.h:513
AFopr_Domainwall::Prec
void Prec(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall-tmpl.h:749
AFopr_Domainwall::DdagD_prec
void DdagD_prec(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall-tmpl.h:717
copy
void copy(Field &y, const Field &x)
copy(y, x): y = x
Definition: field.cpp:212
AFopr_Domainwall::Ddag_prec
void Ddag_prec(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall-tmpl.h:739
AFopr_Domainwall::flop_count
double flop_count()
this returns the number of floating point number operations.
Definition: afopr_Domainwall.h:155
AFopr_Domainwall::set_kernel_parameters
void set_kernel_parameters(const Parameters &params)
set parameters of kernel operaotr.
Definition: afopr_Domainwall-tmpl.h:410
AFopr_Domainwall::convert
void convert(AFIELD &, const Field &)
convert Field to AField for this class.
Definition: afopr_Domainwall-tmpl.h:454
CommonParameters::Nc
static int Nc()
Definition: commonParameters.h:115
AFopr_Domainwall::mult_gm5_4d
void mult_gm5_4d(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall-tmpl.h:685
AFopr_Domainwall::Ldag_inv
void Ldag_inv(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall-tmpl.h:1024
AFopr_Domainwall::get_parameters
void get_parameters(Parameters &params) const
gets parameters by a Parameter object: to be implemented in a subclass.
Definition: afopr_Domainwall-tmpl.h:252
AFopr_Domainwall::set_precond_parameters
void set_precond_parameters()
set parameters for preconditioning.
Definition: afopr_Domainwall-tmpl.h:427
AFopr_Domainwall::D_prec
void D_prec(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall-tmpl.h:729
Parameters::fetch_int_vector
int fetch_int_vector(const string &key, vector< int > &value) const
Definition: parameters.cpp:429
afopr_Domainwall.h
threadManager.h
Parameters::set_int_vector
void set_int_vector(const string &key, const vector< int > &value)
Definition: parameters.cpp:45
AFopr_Domainwall::set_coefficients
void set_coefficients(const std::vector< real_t > b, const std::vector< real_t > c)
set coefficients if they depend in s.
Definition: afopr_Domainwall-tmpl.h:385
AFopr_Domainwall::mult_dag
void mult_dag(AFIELD &v, const AFIELD &w)
hermitian conjugate of mult.
Definition: afopr_Domainwall-tmpl.h:544
CommonParameters::Nd
static int Nd()
Definition: commonParameters.h:116
CommonParameters::Vlevel
static Bridge::VerboseLevel Vlevel()
Definition: commonParameters.h:122
Bridge::BridgeIO::set_verbose_level
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:133
AFopr_Domainwall::set_parameters
void set_parameters(const Parameters &params)
sets parameters by a Parameter object: to be implemented in a subclass.
Definition: afopr_Domainwall-tmpl.h:203
Parameters::set_int
void set_int(const string &key, const int value)
Definition: parameters.cpp:36
scal
void scal(Field &x, const double a)
scal(x, a): x = a * x
Definition: field.cpp:261
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
commonParameters.h
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
AFopr_Domainwall::set_mode
void set_mode(std::string mode)
setting the mode of multiplication if necessary. Default implementation here is just to avoid irrelev...
Definition: afopr_Domainwall-tmpl.h:500
Field
Container of Field-type object.
Definition: field.h:46
AFopr_Domainwall::Hdag
void Hdag(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall-tmpl.h:776
communicator.h
AFopr_Domainwall::mult_R
void mult_R(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall-tmpl.h:801
ThreadManager::get_thread_id
static int get_thread_id()
returns thread id.
Definition: threadManager.cpp:253
AFopr_Domainwall::mult_gm5R
void mult_gm5R(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall-tmpl.h:785
Parameters::fetch_int
int fetch_int(const string &key, int &value) const
Definition: parameters.cpp:346
AFopr_Domainwall::D
void D(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall-tmpl.h:818
AFopr_Domainwall::Ddag
void Ddag(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall-tmpl.h:859
Bridge::BridgeIO::general
void general(const char *format,...)
Definition: bridgeIO.cpp:200
AFopr_Domainwall::H
void H(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall-tmpl.h:767
AFopr_Domainwall::L_inv
void L_inv(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall-tmpl.h:897
AFopr_Domainwall::real_t
AFIELD::real_t real_t
Definition: afopr_Domainwall.h:42
AFopr_Domainwall::mult_chproj_4d
void mult_chproj_4d(AFIELD &, const AFIELD &, const int ipm)
Definition: afopr_Domainwall-tmpl.h:666
ThreadManager::assert_single_thread
static void assert_single_thread(const std::string &class_name)
assert currently running on single thread.
Definition: threadManager.cpp:372
AFopr_Domainwall::init
void init(const Parameters &params)
initial setup.
Definition: afopr_Domainwall-tmpl.h:27
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