14 template<
typename AFIELD>
16 =
"AFopr_Domainwall_5din_eo";
22 #define START_TIMER(var_timer) var_timer->start()
23 #define STOP_TIMER(var_timer) var_timer->stop()
25 #define START_TIMER(var_timer)
26 #define STOP_TIMER(var_timer)
30 template<
typename AFIELD>
42 vout.
general(m_vl,
"%s: construction\n", class_name.c_str());
49 if (repr !=
"Dirac") {
50 vout.
crucial(
"Error at %s: unsupported gamma-matrix type: %s\n",
51 class_name.c_str(), repr.c_str());
65 timer_mult_Deo.reset(
new Timer(
"afopr_Domainwall_5din_eo: Deo & dagger "));
66 timer_mult_Dee_inv.reset(
new Timer(
"afopr_Domainwall_5din_eo: Dee_inv & dag "));
67 timer_pack.reset(
new Timer(
"afopr_Domainwall_5din_eo: pack "));
68 timer_bulk.reset(
new Timer(
"afopr_Domainwall_5din_eo: bulk "));
69 timer_boundary.reset(
new Timer(
"afopr_Domainwall_5din_eo: boundary "));
70 timer_comm.reset(
new Timer(
"afopr_Domainwall_5din_eo: comm "));
71 timer_comm_recv_wait.reset(
new Timer(
"afopr_Domainwall_5din_eo: comm_recv_wait "));
72 timer_comm_send_wait.reset(
new Timer(
"afopr_Domainwall_5din_eo: comm_send_wait "));
73 timer_comm_recv_start.reset(
new Timer(
"afopr_Domainwall_5din_eo: comm_recv_start"));
74 timer_comm_send_start.reset(
new Timer(
"afopr_Domainwall_5din_eo: comm_send_start"));
76 vout.
detailed(m_vl,
"%s: detailed timer was initialized\n", class_name.c_str());
81 m_Nvcd = 2 * Nc * m_Nd;
97 m_Ieo_origin = (ipe1 * m_Ny + ipe2 * m_Nz + ipe3 * m_Nt) % 2;
114 if (!params.
fetch_int(
"require_communication", req_comm)) {
115 vout.
general(m_vl,
"req_comm = %d (input)\n", req_comm);
117 vout.
general(m_vl,
"req_comm = %d (default)\n", req_comm);
121 for (
int mu = 0; mu < m_Ndim; ++mu) {
124 do_comm_any += do_comm[mu];
125 vout.
general(
"do_comm[%d] = %d\n", mu, do_comm[mu]);
129 set_parameters(params);
131 m_Nbdsize.resize(m_Ndim);
132 int Nbdin = (m_Nvcd / 2) * m_Ns;
133 m_Nbdsize[0] = Nbdin * ceil_nwp((m_Ny * m_Nz * m_Nt + 1)/2);
134 m_Nbdsize[1] = Nbdin * ceil_nwp(m_Nx2 * m_Nz * m_Nt);
135 m_Nbdsize[2] = Nbdin * ceil_nwp(m_Nx2 * m_Ny * m_Nt);
136 m_Nbdsize[3] = Nbdin * ceil_nwp(m_Nx2 * m_Ny * m_Nz);
141 int Nst_pad2 = 2 * ceil_nwp(m_Nst2);
142 m_Ueo.reset(m_Ndf, Nst_pad2, m_Ndim);
150 template<
typename AFIELD>
156 for(
int mu = 0; mu < m_Ndim; ++mu){
177 timer_mult_Deo->report();
178 timer_mult_Dee_inv->report();
179 timer_pack->report();
180 timer_bulk->report();
181 timer_boundary->report();
182 timer_comm->report();
183 timer_comm_recv_wait->report();
184 timer_comm_send_wait->report();
185 timer_comm_recv_start->report();
186 timer_comm_send_start->report();
193 template<
typename AFIELD>
198 chsend_up.resize(m_Ndim);
199 chrecv_up.resize(m_Ndim);
200 chsend_dn.resize(m_Ndim);
201 chrecv_dn.resize(m_Ndim);
203 for (
int mu = 0; mu < m_Ndim; ++mu) {
204 size_t Nvsize = m_Nbdsize[mu] *
sizeof(
real_t);
206 chsend_dn[mu].send_init(Nvsize, mu, -1);
207 chsend_up[mu].send_init(Nvsize, mu, 1);
209 chrecv_up[mu].recv_init(Nvsize, mu, 1);
210 chrecv_dn[mu].recv_init(Nvsize, mu, -1);
212 void *buf_up = (
void *)chsend_dn[mu].ptr();
213 chrecv_up[mu].recv_init(Nvsize, mu, 1, buf_up);
214 void *buf_dn = (
void *)chsend_up[mu].ptr();
215 chrecv_dn[mu].recv_init(Nvsize, mu, -1, buf_dn);
218 if (do_comm[mu] == 1) {
219 chset_send.append(chsend_up[mu]);
220 chset_send.append(chsend_dn[mu]);
221 chset_recv.append(chrecv_up[mu]);
222 chset_recv.append(chrecv_dn[mu]);
245 template<
typename AFIELD>
260 int err_optional = 0;
261 err_optional += params.
fetch_string(
"gamma_matrix_type", m_repr);
263 vout.
crucial(m_vl,
" gamma_matrix_type is not given\n");
267 err_optional += params.
fetch_string(
"code_implementation", m_impl);
269 vout.
crucial(m_vl,
" code_implementation is not given\n");
276 err += params.
fetch_int(
"extent_of_5th_dimension", Ns);
280 vout.
crucial(m_vl,
"Error at %s: input parameter not found.\n",
290 vout.
general(m_vl,
" coefficients b, c are not provided:"
291 " set to Shamir's form.\n");
299 vout.
general(m_vl,
" parameter alpha is not provided: set to 1.0.\n");
309 template<
typename AFIELD>
313 params.
set_string(
"kernel_type", m_kernel_type);
314 params.
set_string(
"gamma_matrix_type", m_repr);
315 params.
set_string(
"code_implementation", m_impl);
316 params.
set_double(
"quark_mass",
double(m_mq));
317 params.
set_double(
"domain_wall_height",
double(m_M0));
318 params.
set_int(
"extent_of_5th_dimension", m_Ns);
320 params.
set_double(
"coefficient_b",
double(m_b[0]));
321 params.
set_double(
"coefficient_c",
double(m_c[0]));
322 params.
set_double(
"parameter_alpha",
double(m_alpha));
323 params.
set_string(
"gamma_matrix_type", m_repr);
330 template<
typename AFIELD>
335 const std::vector<int> bc,
348 m_NinF = m_Nvcd * m_Ns;
351 assert(bc.size() == m_Ndim);
352 if (m_boundary.size() != m_Ndim) m_boundary.resize(m_Ndim);
354 for (
int mu = 0; mu < m_Ndim; ++mu) {
355 m_boundary[mu] = bc[mu];
362 m_bc2[mu] = m_boundary[mu];
366 if (m_b.size() != m_Ns) {
370 for (
int is = 0; is < m_Ns; ++is) {
378 vout.
general(m_vl,
"%s: input parameters\n", class_name.c_str());
379 vout.
general(m_vl,
" gamma matrix repr.: %s\n", m_repr.c_str());
380 vout.
general(m_vl,
" code implementation: %s\n", m_impl.c_str());
384 for (
int mu = 0; mu < m_Ndim; ++mu) {
385 vout.
general(m_vl,
" boundary[%d] = %2d\n", mu, m_boundary[mu]);
388 for (
int is = 0; is < m_Ns; ++is) {
389 vout.
general(m_vl,
" b[%2d] = %16.10f c[%2d] = %16.10f\n",
390 is, m_b[is], is, m_c[is]);
395 if (m_w1.nex() != Ns) {
396 m_w1.reset(m_NinF, m_Nst2, 1);
397 m_v1.reset(m_NinF, m_Nst2, 1);
398 m_v2.reset(m_NinF, m_Nst2, 1);
401 set_precond_parameters();
408 template<
typename AFIELD>
410 const std::vector<real_t> vec_b,
411 const std::vector<real_t> vec_c)
415 if ((vec_b.size() != m_Ns) || (vec_c.size() != m_Ns)) {
416 vout.
crucial(m_vl,
"%s: size of coefficient vectors incorrect.\n",
420 vout.
general(m_vl,
"%s: coefficient vectors are set:\n",
425 for (
int is = 0; is < m_Ns; ++is) {
428 vout.
general(m_vl,
"b[%2d] = %16.10f c[%2d] = %16.10f\n",
429 is, m_b[is], is, m_c[is]);
433 set_precond_parameters();
439 template<
typename AFIELD>
445 if (m_dp.size() != m_Ns) {
448 m_dpinv.resize(m_Ns);
449 m_e.resize(m_Ns - 1);
450 m_f.resize(m_Ns - 1);
453 for (
int is = 0; is < m_Ns; ++is) {
454 m_dp[is] = m_alpha * (1.0 + m_b[is] * (4.0 - m_M0));
455 m_dm[is] = m_alpha * (1.0 - m_c[is] * (4.0 - m_M0));
458 m_e[0] = m_mq * m_dm[m_Ns - 1] / m_dp[0];
459 m_f[0] = m_mq * m_dm[0]/m_alpha;
460 for (
int is = 1; is < m_Ns - 1; ++is) {
461 m_e[is] = m_e[is - 1] * m_dm[is - 1] / m_dp[is];
462 m_f[is] = m_f[is - 1] * m_dm[is] / m_dp[is - 1];
465 m_g = m_e[m_Ns - 2] * m_dm[m_Ns - 2];
467 for (
int is = 0; is < m_Ns - 1; ++is) {
468 m_dpinv[is] = 1.0 / m_dp[is];
470 m_dpinv[m_Ns - 1] = 1.0 / (m_dp[m_Ns - 1] + m_g);
475 set_matrix5d_inverse();
480 template<
typename AFIELD>
485 int mat_size = m_Nd2 * m_Ns;
486 m_mat_inv.resize(mat_size * mat_size);
488 for(
int isb = 0; isb < m_Ns; ++isb){
489 for(
int idb = 0; idb < m_Nd2; ++idb){
491 int inb = 2 * Nc * (idb * 2 + m_Nd * isb);
493 int idxb = m_index_eo.idxh(inb, m_NinF, site, 0);
496 D_ee_inv(m_v2, m_w1, 0);
498 for(
int isa = 0; isa < m_Ns; ++isa){
499 for(
int ida = 0; ida < m_Nd2; ++ida){
500 int ina = 2 * Nc * (ida * 2 + m_Nd * isa);
501 int idxa = m_index_eo.idxh(ina, m_NinF, site, 0);
502 m_mat_inv[mat_index(ida, isa, idb, isb)] = m_v2.cmp(idxa);
513 template<
typename AFIELD>
518 vout.
detailed(m_vl,
"%s: set_config is called: num_threads = %d\n",
519 class_name.c_str(), nth);
527 vout.
detailed(m_vl,
"%s: set_config finished\n", class_name.c_str());
531 template<
typename AFIELD>
543 template<
typename AFIELD>
554 template<
typename AFIELD>
561 vout.
crucial(m_vl,
"%s: convert: field size irrelevant\n");
565 int ith, nth, isite, nsite;
566 set_threadtask_afopr(ith, nth, isite, nsite, Nst);
570 for (
int site = isite; site < nsite; ++site) {
571 for (
int is = 0; is < m_Ns; ++is) {
572 for (
int ivcd = 0; ivcd <
NVCD; ++ivcd) {
573 int in_alt = ivcd +
NVCD * is;
575 v.set_host(index.idx(in_alt, m_NinF, site, 0), vt);
586 template<
typename AFIELD>
593 vout.
crucial(m_vl,
"%s: convert: field size irrelevant\n");
597 int ith, nth, isite, nsite;
598 set_threadtask_afopr(ith, nth, isite, nsite, Nst);
604 for (
int site = isite; site < nsite; ++site) {
605 for (
int is = 0; is < m_Ns; ++is) {
606 for (
int ivcd = 0; ivcd <
NVCD; ++ivcd) {
607 int in_alt = ivcd +
NVCD * is;
608 double vt = double(w.cmp_host(index.idx(in_alt, m_NinF, site, 0)));
609 v.
set(ivcd, site, is, vt);
619 template<
typename AFIELD>
625 if (ith == 0) m_mode = mode;
633 template<
typename AFIELD>
639 }
else if (m_mode ==
"Ddag") {
642 }
else if (m_mode ==
"DdagD") {
650 }
else if (m_mode ==
"Deo") {
652 }
else if (m_mode ==
"Doe") {
654 }
else if (m_mode ==
"Dee") {
656 }
else if (m_mode ==
"Doo") {
658 }
else if (m_mode ==
"Dee_inv") {
661 }
else if (m_mode ==
"Doo_inv") {
665 vout.
crucial(m_vl,
"mode undeifined in %s.\n", class_name.c_str());
666 vout.
crucial(m_vl,
"in mult, mode=%s.\n", m_mode.c_str());
673 template<
typename AFIELD>
679 }
else if (m_mode ==
"Ddag") {
682 }
else if (m_mode ==
"DdagD") {
687 }
else if (m_mode ==
"Deo") {
689 }
else if (m_mode ==
"Doe") {
691 }
else if (m_mode ==
"Dee") {
693 }
else if (m_mode ==
"Doo") {
695 }
else if (m_mode ==
"Dee_inv") {
696 Ddag_ee_inv(v, w, 0);
698 }
else if (m_mode ==
"Doo_inv") {
699 Ddag_ee_inv(v, w, 1);
702 vout.
crucial(m_vl,
"mode undeifined in %s.\n", class_name.c_str());
703 vout.
crucial(m_vl,
"in mult_dag, mode=%s.\n", m_mode.c_str());
710 template<
typename AFIELD>
719 }
else if (mode ==
"Doe") {
721 }
else if (mode ==
"Dee") {
723 }
else if (mode ==
"Doo") {
725 }
else if (mode ==
"Dee_inv") {
728 }
else if (mode ==
"Doo_inv") {
732 vout.
crucial(m_vl,
"mode undeifined in %s.\n", class_name.c_str());
733 vout.
crucial(m_vl,
"in mult, mode=%s.\n", m_mode.c_str());
740 template<
typename AFIELD>
750 }
else if (mode ==
"Doe") {
752 }
else if (mode ==
"Dee") {
754 }
else if (mode ==
"Doo") {
756 }
else if (mode ==
"Dee_inv") {
757 Ddag_ee_inv(v, w, 0);
759 }
else if (mode ==
"Doo_inv") {
760 Ddag_ee_inv(v, w, 1);
763 vout.
crucial(m_vl,
"mode undeifined in %s.\n", class_name.c_str());
764 vout.
crucial(m_vl,
"in mult, mode=%s.\n", m_mode.c_str());
770 template<
typename AFIELD>
792 template<
typename AFIELD>
805 template<
typename AFIELD>
814 Ddag_eo(m_v2, m_v1, 1);
815 LUdag_inv(m_v1, m_v2);
816 Ddag_eo(m_v2, m_v1, 0);
825 template<
typename AFIELD>
829 D_ee_inv_alt(m_v2, m_v1, 1);
831 D_ee_inv_alt(m_v2, m_v1, 0);
836 Ddag_ee_inv_alt(v, m_v1, 0);
838 Ddag_ee_inv_alt(v, m_v2, 1);
846 template<
typename AFIELD>
850 D_ee_inv_alt(m_v2, m_v1, 1);
852 D_ee_inv_alt(m_v2, m_v1, 0);
860 template<
typename AFIELD>
868 Ddag_ee_inv_alt(m_v1, w, 0);
869 Ddag_eo(m_v2, m_v1, 1);
870 Ddag_ee_inv_alt(m_v1, m_v2, 1);
871 Ddag_eo(m_v2, m_v1, 0);
880 template<
typename AFIELD>
893 vp, wp, m_Ns, m_Nsize);
899 template<
typename AFIELD>
921 template<
typename AFIELD>
934 vp, wp, m_Ns, m_Nsize);
941 template<
typename AFIELD>
953 if(nth > 1) ith_kernel = 1;
958 real_t *up = m_Ueo.ptr(0);
959 int jeo = (m_Ieo_origin + ieo) % 2;
961 if(do_comm_any > 0 && ith == 0){
967 if (ith == ith_kernel){
970 yp, wp, m_mq, m_M0, m_Ns,
971 &m_b[0], &m_c[0], m_alpha, m_Nsize);
973 if (do_comm_any > 0) {
985 buf1_xp, buf1_xm, buf1_yp, buf1_ym,
986 buf1_zp, buf1_zm, buf1_tp, buf1_tm,
988 m_Ns, m_bc, m_Nsize, do_comm, ieo, jeo, 0);
995 if(do_comm_any > 0 && ith == 0){
1002 if (ith == ith_kernel){
1006 vp, up, yp, m_Ns, m_bc2,
1007 m_Nsize, do_comm, ieo, jeo, 0);
1009 #ifdef USE_DOMAINWALL_5DIN_4D_KERNEL
1011 vp, up, yp, m_Ns, m_bc2,
1012 m_Nsize, do_comm, ieo, jeo, 0);
1014 vout.
crucial(m_vl,
"%s: 4D kernel not compiled\n",
1015 class_name.c_str());
1022 if(do_comm_any > 0 && ith == 0){
1030 if(do_comm_any > 0 && ith == 0){
1034 if(do_comm_any > 0 && ith == ith_kernel){
1048 buf2_xp, buf2_xm, buf2_yp, buf2_ym,
1049 buf2_zp, buf2_zm, buf2_tp, buf2_tm,
1051 m_Nsize, do_comm, ieo, jeo);
1055 if(do_comm_any > 0 && ith == 0){
1068 template<
typename AFIELD>
1081 if(nth > 1) ith_kernel = 1;
1086 real_t *up = m_Ueo.ptr(0);
1087 int jeo = (m_Ieo_origin + ieo) % 2;
1089 if(do_comm_any > 0 && ith == 0){
1095 if (ith == ith_kernel){
1097 if (do_comm_any > 0) {
1110 buf1_xp, buf1_xm, buf1_yp, buf1_ym,
1111 buf1_zp, buf1_zm, buf1_tp, buf1_tm,
1113 m_Ns, m_bc, m_Nsize, do_comm, ieo, jeo, 1);
1119 if(do_comm_any > 0 && ith == 0){
1126 if(ith == ith_kernel){
1130 yp, up, wp, m_Ns, m_bc2,
1131 m_Nsize, do_comm, ieo, jeo, 1);
1133 #ifdef USE_DOMAINWALL_5DIN_4D_KERNEL
1135 yp, up, wp, m_Ns, m_bc2,
1136 m_Nsize, do_comm, ieo, jeo, 1);
1138 vout.
crucial(m_vl,
"%s: 4D kernel not compiled\n",
1139 class_name.c_str());
1146 if(do_comm_any > 0 && ith == 0){
1154 if(do_comm_any > 0 && ith == 0){
1158 if (ith == ith_kernel) {
1159 if (do_comm_any > 0) {
1173 buf2_xp, buf2_xm, buf2_yp, buf2_ym,
1174 buf2_zp, buf2_zm, buf2_tp, buf2_tm,
1176 m_Nsize, do_comm, ieo, jeo);
1181 vp, yp, m_mq, m_M0, m_Ns,
1182 &m_b[0], &m_c[0], m_alpha, m_Nsize);
1185 if(do_comm_any > 0 && ith == 0){
1198 template<
typename AFIELD>
1210 vp, wp, m_mq, m_M0, m_Ns,
1211 &m_b[0], &m_c[0], m_alpha, m_Nsize);
1220 template<
typename AFIELD>
1232 vp, wp, m_mq, m_M0, m_Ns,
1233 &m_b[0], &m_c[0], m_alpha, m_Nsize);
1242 template<
typename AFIELD>
1253 #ifdef USE_DOMAINWALL_5DIN_EE_MATINV_KERNEL
1260 vp, wp, 1, m_Ns, &m_mat_inv[0], m_Nsize);
1263 vp, wp, 1, m_Ns, &m_mat_inv[0], m_Nsize);
1267 vout.
crucial(m_vl,
"%s: Dee inverse matrix mult not compiled\n",
1268 class_name.c_str());
1279 template<
typename AFIELD>
1290 #ifdef USE_DOMAINWALL_5DIN_EE_MATINV_KERNEL
1297 vp, wp, -1, m_Ns, &m_mat_inv[0], m_Nsize);
1300 vp, wp, -1, m_Ns, &m_mat_inv[0], m_Nsize);
1304 vout.
crucial(m_vl,
"%s: Dee inverse matrix mult not compiled\n",
1305 class_name.c_str());
1316 template<
typename AFIELD>
1329 template<
typename AFIELD>
1342 template<
typename AFIELD>
1356 vp, wp, m_Ns, m_Nsize,
1357 &m_e[0], &m_f[0], &m_dpinv[0], &m_dm[0], m_alpha);
1368 template<
typename AFIELD>
1381 vp, wp, m_Ns, m_Nsize,
1382 &m_e[0], &m_f[0], &m_dpinv[0], &m_dm[0], m_alpha);
1393 template<
typename AFIELD>
1399 double vsite =
static_cast<double>(Lvol);
1400 double vNs =
static_cast<double>(m_Ns);
1402 double flop_Wilson_hop;
1407 if (m_repr ==
"Dirac") {
1408 flop_Wilson_hop =
static_cast<double>(
1409 vNs * Nc * Nd * (6 * (4 * Nc + 2) + 2 * (4 * Nc + 1))) * vsite;
1410 flop_pre =
static_cast<double>(vNs * Nc * Nd * 14) * vsite;
1420 flop_LU_inv =
static_cast<double>(Nc * Nd * (22 * vNs -16)) * vsite;
1421 }
else if (m_repr ==
"Chiral") {
1422 flop_Wilson_hop =
static_cast<double>(
1423 vNs * Nc * Nd * (8 * (4 * Nc + 2))) * vsite;
1424 flop_pre =
static_cast<double>(vNs * Nc * Nd * 6) * vsite;
1431 flop_LU_inv =
static_cast<double>(Nc * Nd * (10 * vNs -8)) * vsite;
1433 double flop_axpy =
static_cast<double>(vNs * 2 * Nc * Nd) * vsite;
1439 double flop_Deo = flop_Wilson_hop + flop_pre;
1446 if ((mode ==
"D") || (mode ==
"Ddag")) {
1447 flop = flop_Deo + flop_LU_inv + 0.5 * flop_axpy;
1448 }
else if (mode ==
"DdagD") {
1449 flop = 2 * flop_Deo + 2 * flop_LU_inv + flop_axpy;
1450 }
else if ((mode ==
"Dee_inv") || (mode ==
"Doo_inv")) {
1451 flop = 0.5 * flop_LU_inv;
1452 }
else if ((mode ==
"Dee") || (mode ==
"Doo")) {
1453 flop = 0.5 * flop_pre;
1454 }
else if ((mode ==
"Doe") || (mode ==
"Doe")) {
1455 flop = 0.5 * flop_Deo;
1457 vout.
crucial(m_vl,
"Error at %s: input mode is undefined: %s.\n",
1458 class_name.c_str(), mode.c_str());