From 808f1e0e8c199204c7369fe2a033bc6041cbaa91 Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Tue, 15 Dec 2020 16:33:29 +0000 Subject: [PATCH 01/23] merge develop --- tests/solver/Test_zMADWF_prec.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/solver/Test_zMADWF_prec.cc b/tests/solver/Test_zMADWF_prec.cc index d1168764..f18e1d86 100644 --- a/tests/solver/Test_zMADWF_prec.cc +++ b/tests/solver/Test_zMADWF_prec.cc @@ -52,7 +52,7 @@ struct TestParams{ bool zmobius_inner; double lambda_max; //upper bound of H_T eigenvalue range required to generate zMobius approximation - TestParams(): load_config(true), config_file("ckpoint_lat.1000"), mass(0.01), + TestParams(): load_config(false), config_file("ckpoint_lat.1000"), mass(0.01), Ls_outer(24), b_plus_c_outer(2.0), resid_outer(1e-8), Ls_inner(12), b_plus_c_inner(1.0), resid_inner(1e-8), zmobius_inner(true), lambda_max(1.42), outer_precon("Standard"), inner_precon("Standard") {} @@ -246,7 +246,7 @@ void run(const TestParams ¶ms){ typename RunParamsInner::SchurSolverType SchurSolver_inner(CG_inner); ZeroGuesser Guess; - MADWF > madwf(D_outer, D_inner, PV_outer, SchurSolver_inner, Guess, params.resid_outer, 100, &update); + MADWF > madwf(D_outer, D_inner, PV_outer, SchurSolver_inner, Guess, params.resid_outer, 10000, &update); LatticeFermionD result_MADWF(FGrid_outer); result_MADWF = Zero(); From f36d6f3923b7632e169b6740c94cc39ecc6bc8a9 Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Thu, 17 Dec 2020 17:04:08 +0000 Subject: [PATCH 02/23] compiles on GPU. 3pt still wrong!!!! --- Grid/qcd/utils/A2Autils.h | 2 +- Grid/qcd/utils/BaryonUtils.h | 991 +++++++++++++++++++---------------- 2 files changed, 544 insertions(+), 449 deletions(-) diff --git a/Grid/qcd/utils/A2Autils.h b/Grid/qcd/utils/A2Autils.h index b63d8571..497927dd 100644 --- a/Grid/qcd/utils/A2Autils.h +++ b/Grid/qcd/utils/A2Autils.h @@ -1047,7 +1047,7 @@ A2Autils::ContractWWVV(std::vector &WWVV, { GridBase *grid = vs[0].Grid(); - int nd = grid->_ndimension; + //int nd = grid->_ndimension; int Nsimd = grid->Nsimd(); int N_t = WW_sd.dimensions()[0]; int N_s = WW_sd.dimensions()[1]; diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index 15516b56..25c71e3a 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -44,24 +44,24 @@ public: typedef typename ComplexField::vector_object vobj; typedef Lattice> SpinMatrixField; - typedef typename SpinMatrixField::vector_object sobj; + //typedef typename SpinMatrixField::vector_object sobj; - static const int epsilon[6][3] ; - static const Real epsilon_sgn[6]; + //static const int epsilon[6][3] ; + //static const Real epsilon_sgn[6]; private: - template + template accelerator_inline static void BaryonSite(const mobj &D1, - const mobj &D2, - const mobj &D3, - const Gamma GammaA_left, - const Gamma GammaB_left, - const Gamma GammaA_right, - const Gamma GammaB_right, - const int parity, - const bool * wick_contractions, - robj &result); - template + const mobj &D2, + const mobj &D3, + const Gamma GammaA_left, + const Gamma GammaB_left, + const Gamma GammaA_right, + const Gamma GammaB_right, + const int parity, + const bool * wick_contractions, + robj &result); + template accelerator_inline static void BaryonSiteMatrix(const mobj &D1, const mobj &D2, const mobj &D3, @@ -76,15 +76,15 @@ public: std::string qf, bool* wick_contractions); static void ContractBaryons(const PropagatorField &q1_left, - const PropagatorField &q2_left, - const PropagatorField &q3_left, - const Gamma GammaA_left, - const Gamma GammaB_left, - const Gamma GammaA_right, - const Gamma GammaB_right, - const bool* wick_contractions, - const int parity, - ComplexField &baryon_corr); + const PropagatorField &q2_left, + const PropagatorField &q3_left, + const Gamma GammaA_left, + const Gamma GammaB_left, + const Gamma GammaA_right, + const Gamma GammaB_right, + const bool* wick_contractions, + const int parity, + ComplexField &baryon_corr); static void ContractBaryonsMatrix(const PropagatorField &q1_left, const PropagatorField &q2_left, const PropagatorField &q3_left, @@ -96,16 +96,16 @@ public: SpinMatrixField &baryon_corr); template static void ContractBaryonsSliced(const mobj &D1, - const mobj &D2, - const mobj &D3, - const Gamma GammaA_left, - const Gamma GammaB_left, - const Gamma GammaA_right, - const Gamma GammaB_right, - const bool* wick_contractions, - const int parity, - const int nt, - robj &result); + const mobj &D2, + const mobj &D3, + const Gamma GammaA_left, + const Gamma GammaB_left, + const Gamma GammaA_right, + const Gamma GammaB_right, + const bool* wick_contractions, + const int parity, + const int nt, + robj &result); template static void ContractBaryonsSlicedMatrix(const mobj &D1, const mobj &D2, @@ -118,11 +118,11 @@ public: const int nt, robj &result); private: - template + template accelerator_inline static void BaryonGamma3ptGroup1Site( const mobj &Dq1_ti, const mobj2 &Dq2_spec, - const mobj2 &Dq3_spec, + // const mobj2 &Dq3_spec, const mobj &Dq4_tf, const Gamma GammaJ, const Gamma GammaBi, @@ -130,11 +130,11 @@ public: int wick_contraction, robj &result); - template + template accelerator_inline static void BaryonGamma3ptGroup2Site( const mobj2 &Dq1_spec, const mobj &Dq2_ti, - const mobj2 &Dq3_spec, + //const mobj2 &Dq3_spec, const mobj &Dq4_tf, const Gamma GammaJ, const Gamma GammaBi, @@ -142,10 +142,10 @@ public: int wick_contraction, robj &result); - template + template accelerator_inline static void BaryonGamma3ptGroup3Site( const mobj2 &Dq1_spec, - const mobj2 &Dq2_spec, + //const mobj2 &Dq2_spec, const mobj &Dq3_ti, const mobj &Dq4_tf, const Gamma GammaJ, @@ -167,86 +167,78 @@ public: const Gamma GammaBf, SpinMatrixField &stn_corr); private: - template + template accelerator_inline static void SigmaToNucleonQ1EyeSite(const mobj &Dq_loop, - const mobj2 &Du_spec, - const mobj &Dd_tf, - const mobj &Ds_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, - robj &result); - template + const mobj2 &Du_spec, + const mobj &Dd_tf, + const mobj &Ds_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + robj &result); + template accelerator_inline static void SigmaToNucleonQ1NonEyeSite(const mobj &Du_ti, - const mobj &Du_tf, - const mobj2 &Du_spec, - const mobj &Dd_tf, - const mobj &Ds_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, - robj &result); + const mobj &Du_tf, + const mobj2 &Du_spec, + const mobj &Dd_tf, + const mobj &Ds_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + robj &result); - template + template accelerator_inline static void SigmaToNucleonQ2EyeSite(const mobj &Dq_loop, - const mobj2 &Du_spec, - const mobj &Dd_tf, - const mobj &Ds_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, - robj &result); - template + const mobj2 &Du_spec, + const mobj &Dd_tf, + const mobj &Ds_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + robj &result); + template accelerator_inline static void SigmaToNucleonQ2NonEyeSite(const mobj &Du_ti, - const mobj &Du_tf, - const mobj2 &Du_spec, - const mobj &Dd_tf, - const mobj &Ds_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, - robj &result); + const mobj &Du_tf, + const mobj2 &Du_spec, + const mobj &Dd_tf, + const mobj &Ds_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + robj &result); public: template static void SigmaToNucleonEye(const PropagatorField &qq_loop, - const mobj &Du_spec, - const PropagatorField &qd_tf, - const PropagatorField &qs_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, - const std::string op, - SpinMatrixField &stn_corr); + const mobj &Du_spec, + const PropagatorField &qd_tf, + const PropagatorField &qs_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + const std::string op, + SpinMatrixField &stn_corr); template static void SigmaToNucleonNonEye(const PropagatorField &qq_ti, - const PropagatorField &qq_tf, - const mobj &Du_spec, - const PropagatorField &qd_tf, - const PropagatorField &qs_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, - const std::string op, - SpinMatrixField &stn_corr); + const PropagatorField &qq_tf, + const mobj &Du_spec, + const PropagatorField &qd_tf, + const PropagatorField &qs_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + const std::string op, + SpinMatrixField &stn_corr); }; - +/* template const int BaryonUtils::epsilon[6][3] = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}}; -/*template -const Complex BaryonUtils::epsilon_sgn[6] = {Complex(1), - Complex(1), - Complex(1), - Complex(-1), - Complex(-1), - Complex(-1)}; -*/ template const Real BaryonUtils::epsilon_sgn[6] = {1.,1.,1.,-1.,-1.,-1.}; - +*/ //This is the old version template -template +template accelerator_inline void BaryonUtils::BaryonSite(const mobj &D1, const mobj &D2, const mobj &D3, @@ -274,16 +266,20 @@ void BaryonUtils::BaryonSite(const mobj &D1, auto GBf_D3 = GammaB_f * D3; auto GAf_D3 = GammaA_f * D3; - for (int ie_f=0; ie_f < 6 ; ie_f++){ - int a_f = epsilon[ie_f][0]; //a - int b_f = epsilon[ie_f][1]; //b - int c_f = epsilon[ie_f][2]; //c - for (int ie_i=0; ie_i < 6 ; ie_i++){ - int a_i = epsilon[ie_i][0]; //a' - int b_i = epsilon[ie_i][1]; //b' - int c_i = epsilon[ie_i][2]; //c' + Real ee; - Real ee = epsilon_sgn[ie_f] * epsilon_sgn[ie_i]; + for (int ie_f=0; ie_f < 6 ; ie_f++){ + int a_f = (ie_f < 3 ? ie_f : (6-ie_f)%3 ); //epsilon[ie_n][0]; //a + int b_f = (ie_f < 3 ? (ie_f+1)%3 : (8-ie_f)%3 ); //epsilon[ie_n][1]; //b + int c_f = (ie_f < 3 ? (ie_f+2)%3 : (7-ie_f)%3 ); //epsilon[ie_n][2]; //c + int eSgn_f = (ie_f < 3 ? 1 : -1); + for (int ie_i=0; ie_i < 6 ; ie_i++){ + int a_i = (ie_i < 3 ? ie_i : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a' + int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b' + int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c' + int eSgn_i = (ie_i < 3 ? 1 : -1); + + ee = Real(eSgn_f * eSgn_i); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s]; //This is the \delta_{456}^{123} part if (wick_contraction[0]){ for (int rho=0; rho::BaryonSite(const mobj &D1, //New version without parity projection or trace template -template +template accelerator_inline void BaryonUtils::BaryonSiteMatrix(const mobj &D1, const mobj &D2, const mobj &D3, @@ -384,16 +380,21 @@ void BaryonUtils::BaryonSiteMatrix(const mobj &D1, auto GBf_D3 = GammaB_f * D3; auto GAf_D3 = GammaA_f * D3; - for (int ie_f=0; ie_f < 6 ; ie_f++){ - int a_f = epsilon[ie_f][0]; //a - int b_f = epsilon[ie_f][1]; //b - int c_f = epsilon[ie_f][2]; //c - for (int ie_i=0; ie_i < 6 ; ie_i++){ - int a_i = epsilon[ie_i][0]; //a' - int b_i = epsilon[ie_i][1]; //b' - int c_i = epsilon[ie_i][2]; //c' - Real ee = epsilon_sgn[ie_f] * epsilon_sgn[ie_i]; + Real ee; + + for (int ie_f=0; ie_f < 6 ; ie_f++){ + int a_f = (ie_f < 3 ? ie_f : (6-ie_f)%3 ); //epsilon[ie_n][0]; //a + int b_f = (ie_f < 3 ? (ie_f+1)%3 : (8-ie_f)%3 ); //epsilon[ie_n][1]; //b + int c_f = (ie_f < 3 ? (ie_f+2)%3 : (7-ie_f)%3 ); //epsilon[ie_n][2]; //c + int eSgn_f = (ie_f < 3 ? 1 : -1); + for (int ie_i=0; ie_i < 6 ; ie_i++){ + int a_i = (ie_i < 3 ? ie_i : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a' + int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b' + int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c' + int eSgn_i = (ie_i < 3 ? 1 : -1); + + ee = Real(eSgn_f * eSgn_i); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s]; //This is the \delta_{456}^{123} part if (wick_contraction[0]){ for (int rho_i=0; rho_i::WickContractions(std::string qi, std::string qf, bool* * Wick_Contractions function above */ template void BaryonUtils::ContractBaryons(const PropagatorField &q1_left, - const PropagatorField &q2_left, - const PropagatorField &q3_left, - const Gamma GammaA_left, - const Gamma GammaB_left, - const Gamma GammaA_right, - const Gamma GammaB_right, - const bool* wick_contractions, - const int parity, - ComplexField &baryon_corr) + const PropagatorField &q2_left, + const PropagatorField &q3_left, + const Gamma GammaA_left, + const Gamma GammaB_left, + const Gamma GammaA_right, + const Gamma GammaB_right, + const bool* wick_contractions, + const int parity, + ComplexField &baryon_corr) { assert(Ns==4 && "Baryon code only implemented for N_spin = 4"); @@ -519,10 +520,10 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_left, GridBase *grid = q1_left.Grid(); - autoView(vbaryon_corr, baryon_corr,CpuWrite); - autoView( v1 , q1_left, CpuRead); - autoView( v2 , q2_left, CpuRead); - autoView( v3 , q3_left, CpuRead); + autoView(vbaryon_corr , baryon_corr , AcceleratorWrite); + autoView( v1 , q1_left , AcceleratorRead); + autoView( v2 , q2_left , AcceleratorRead); + autoView( v3 , q3_left , AcceleratorRead); Real bytes =0.; bytes += grid->oSites() * (432.*sizeof(vComplex) + 126.*sizeof(int) + 36.*sizeof(Real)); @@ -538,12 +539,13 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_left, t =-usecond(); accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - auto D1 = v1[ss]; - auto D2 = v2[ss]; - auto D3 = v3[ss]; - vobj result=Zero(); + auto D1 = v1(ss); + auto D2 = v2(ss); + auto D3 = v3(ss); + typedef decltype(coalescedRead(vbaryon_corr[0])) cVec; + cVec result=Zero(); BaryonSite(D1,D2,D3,GammaA_left,GammaB_left,GammaA_right,GammaB_right,parity,wick_contractions,result); - vbaryon_corr[ss] = result; + coalescedWrite(vbaryon_corr[ss],result); } );//end loop over lattice sites t += usecond(); @@ -567,38 +569,21 @@ void BaryonUtils::ContractBaryonsMatrix(const PropagatorField &q1_left, assert(Nc==3 && "Baryon code only implemented for N_colour = 3"); GridBase *grid = q1_left.Grid(); - - autoView(vbaryon_corr, baryon_corr,CpuWrite); - autoView( v1 , q1_left, CpuRead); - autoView( v2 , q2_left, CpuRead); - autoView( v3 , q3_left, CpuRead); - // Real bytes =0.; - // bytes += grid->oSites() * (432.*sizeof(vComplex) + 126.*sizeof(int) + 36.*sizeof(Real)); - // for (int ie=0; ie < 6 ; ie++){ - // if(ie==0 or ie==3){ - // bytes += grid->oSites() * (4.*sizeof(int) + 4752.*sizeof(vComplex)) * wick_contractions[ie]; - // } - // else{ - // bytes += grid->oSites() * (64.*sizeof(int) + 5184.*sizeof(vComplex)) * wick_contractions[ie]; - // } - // } - // Real t=0.; - // t =-usecond(); + autoView(vbaryon_corr , baryon_corr , AcceleratorWrite); + autoView( v1 , q1_left , AcceleratorRead); + autoView( v2 , q2_left , AcceleratorRead); + autoView( v3 , q3_left , AcceleratorRead); accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - auto D1 = v1[ss]; - auto D2 = v2[ss]; - auto D3 = v3[ss]; - sobj result=Zero(); + auto D1 = v1(ss); + auto D2 = v2(ss); + auto D3 = v3(ss); + typedef decltype(coalescedRead(vbaryon_corr[0])) spinor; + spinor result=Zero(); BaryonSiteMatrix(D1,D2,D3,GammaA_left,GammaB_left,GammaA_right,GammaB_right,wick_contractions,result); - vbaryon_corr[ss] = result; + coalescedWrite(vbaryon_corr[ss],result); } );//end loop over lattice sites - - // t += usecond(); - - // std::cout << GridLogDebug << std::setw(10) << bytes/t*1.0e6/1024/1024/1024 << " GB/s " << std::endl; - } /* The array wick_contractions must be of length 6. The order * @@ -609,16 +594,16 @@ void BaryonUtils::ContractBaryonsMatrix(const PropagatorField &q1_left, template template void BaryonUtils::ContractBaryonsSliced(const mobj &D1, - const mobj &D2, - const mobj &D3, - const Gamma GammaA_left, - const Gamma GammaB_left, - const Gamma GammaA_right, - const Gamma GammaB_right, - const bool* wick_contractions, - const int parity, - const int nt, - robj &result) + const mobj &D2, + const mobj &D3, + const Gamma GammaA_left, + const Gamma GammaB_left, + const Gamma GammaA_right, + const Gamma GammaB_right, + const bool* wick_contractions, + const int parity, + const int nt, + robj &result) { assert(Ns==4 && "Baryon code only implemented for N_spin = 4"); @@ -664,11 +649,11 @@ void BaryonUtils::ContractBaryonsSlicedMatrix(const mobj &D1, * Dq3_spec is a quark line from t_i to t_f * Dq4_tf is a quark line from t_f to t_J */ template -template +template accelerator_inline void BaryonUtils::BaryonGamma3ptGroup1Site( const mobj &Dq1_ti, const mobj2 &Dq2_spec, - const mobj2 &Dq3_spec, + // const mobj2 &Dq3_spec, const mobj &Dq4_tf, const Gamma GammaJ, const Gamma GammaBi, @@ -678,41 +663,47 @@ void BaryonUtils::BaryonGamma3ptGroup1Site( { Gamma g5(Gamma::Algebra::Gamma5); - auto adjD4_g_D1 = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq1_ti; +// auto adjD4_g_D1 = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq1_ti; + auto adjD4 = g5 * adj(Dq4_tf) * g5 ; + auto adjD4_g_D1 = adjD4 * GammaJ * Dq1_ti; auto Gf_adjD4_g_D1 = GammaBf * adjD4_g_D1; auto D2_Gi = Dq2_spec * GammaBi; auto Gf_D2_Gi = GammaBf * D2_Gi; - auto Gf_D3 = GammaBf * Dq3_spec; - int a_f, b_f, c_f; - int a_i, b_i, c_i; +// auto Gf_D3 = GammaBf * Dq3_spec; // including a second mobj2 parameter leads to compilation error + auto Gf_D3 = GammaBf * Dq2_spec; //WRONG!!!!! - Real ee; - for (int ie_f=0; ie_f < 6 ; ie_f++){ - a_f = epsilon[ie_f][0]; //a - b_f = epsilon[ie_f][1]; //b - c_f = epsilon[ie_f][2]; //c - for (int ie_i=0; ie_i < 6 ; ie_i++){ - a_i = epsilon[ie_i][0]; //a' - b_i = epsilon[ie_i][1]; //b' - c_i = epsilon[ie_i][2]; //c' + Real ee; - ee = epsilon_sgn[ie_f] * epsilon_sgn[ie_i]; + for (int ie_f=0; ie_f < 6 ; ie_f++){ + int a_f = (ie_f < 3 ? ie_f : (6-ie_f)%3 ); //epsilon[ie_n][0]; //a + int b_f = (ie_f < 3 ? (ie_f+1)%3 : (8-ie_f)%3 ); //epsilon[ie_n][1]; //b + int c_f = (ie_f < 3 ? (ie_f+2)%3 : (7-ie_f)%3 ); //epsilon[ie_n][2]; //c + int eSgn_f = (ie_f < 3 ? 1 : -1); + for (int ie_i=0; ie_i < 6 ; ie_i++){ + int a_i = (ie_i < 3 ? ie_i : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a' + int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b' + int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c' + int eSgn_i = (ie_i < 3 ? 1 : -1); + + ee = Real(eSgn_f * eSgn_i); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s]; for (int alpha_f=0; alpha_f::BaryonGamma3ptGroup1Site( * Dq3_spec is a quark line from t_i to t_f * Dq4_tf is a quark line from t_f to t_J */ template -template +template accelerator_inline void BaryonUtils::BaryonGamma3ptGroup2Site( const mobj2 &Dq1_spec, const mobj &Dq2_ti, - const mobj2 &Dq3_spec, + // const mobj2 &Dq3_spec, const mobj &Dq4_tf, const Gamma GammaJ, const Gamma GammaBi, @@ -773,37 +764,40 @@ void BaryonUtils::BaryonGamma3ptGroup2Site( auto adjD4_g_D2_Gi = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq2_ti * GammaBi; auto Gf_adjD4_g_D2_Gi = GammaBf * adjD4_g_D2_Gi; auto Gf_D1 = GammaBf * Dq1_spec; - auto Gf_D3 = GammaBf * Dq3_spec; + //auto Gf_D3 = GammaBf * Dq3_spec; + auto Gf_D3 = GammaBf * Dq1_spec; // WRONG!!!!! - int a_f, b_f, c_f; - int a_i, b_i, c_i; - Real ee; + Real ee; - for (int ie_f=0; ie_f < 6 ; ie_f++){ - a_f = epsilon[ie_f][0]; //a - b_f = epsilon[ie_f][1]; //b - c_f = epsilon[ie_f][2]; //c - for (int ie_i=0; ie_i < 6 ; ie_i++){ - a_i = epsilon[ie_i][0]; //a' - b_i = epsilon[ie_i][1]; //b' - c_i = epsilon[ie_i][2]; //c' + for (int ie_f=0; ie_f < 6 ; ie_f++){ + int a_f = (ie_f < 3 ? ie_f : (6-ie_f)%3 ); //epsilon[ie_n][0]; //a + int b_f = (ie_f < 3 ? (ie_f+1)%3 : (8-ie_f)%3 ); //epsilon[ie_n][1]; //b + int c_f = (ie_f < 3 ? (ie_f+2)%3 : (7-ie_f)%3 ); //epsilon[ie_n][2]; //c + int eSgn_f = (ie_f < 3 ? 1 : -1); + for (int ie_i=0; ie_i < 6 ; ie_i++){ + int a_i = (ie_i < 3 ? ie_i : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a' + int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b' + int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c' + int eSgn_i = (ie_i < 3 ? 1 : -1); - ee = epsilon_sgn[ie_f] * epsilon_sgn[ie_i]; + ee = Real(eSgn_f * eSgn_i); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s]; for (int alpha_f=0; alpha_f::BaryonGamma3ptGroup2Site( * Dq3_ti is a quark line from t_i to t_J * Dq4_tf is a quark line from t_f to t_J */ template -template +template accelerator_inline void BaryonUtils::BaryonGamma3ptGroup3Site( const mobj2 &Dq1_spec, - const mobj2 &Dq2_spec, + // const mobj2 &Dq2_spec, const mobj &Dq3_ti, const mobj &Dq4_tf, const Gamma GammaJ, @@ -863,24 +857,25 @@ void BaryonUtils::BaryonGamma3ptGroup3Site( auto adjD4_g_D3 = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq3_ti; auto Gf_adjD4_g_D3 = GammaBf * adjD4_g_D3; auto Gf_D1 = GammaBf * Dq1_spec; - auto D2_Gi = Dq2_spec * GammaBi; + //auto D2_Gi = Dq2_spec * GammaBi; + auto D2_Gi = Dq1_spec * GammaBi; //WRONG!!!!!!!!!!!!!!!!! auto Gf_D2_Gi = GammaBf * D2_Gi; - int a_f, b_f, c_f; - int a_i, b_i, c_i; - Real ee; + Real ee; - for (int ie_f=0; ie_f < 6 ; ie_f++){ - a_f = epsilon[ie_f][0]; //a - b_f = epsilon[ie_f][1]; //b - c_f = epsilon[ie_f][2]; //c - for (int ie_i=0; ie_i < 6 ; ie_i++){ - a_i = epsilon[ie_i][0]; //a' - b_i = epsilon[ie_i][1]; //b' - c_i = epsilon[ie_i][2]; //c' + for (int ie_f=0; ie_f < 6 ; ie_f++){ + int a_f = (ie_f < 3 ? ie_f : (6-ie_f)%3 ); //epsilon[ie_n][0]; //a + int b_f = (ie_f < 3 ? (ie_f+1)%3 : (8-ie_f)%3 ); //epsilon[ie_n][1]; //b + int c_f = (ie_f < 3 ? (ie_f+2)%3 : (7-ie_f)%3 ); //epsilon[ie_n][2]; //c + int eSgn_f = (ie_f < 3 ? 1 : -1); + for (int ie_i=0; ie_i < 6 ; ie_i++){ + int a_i = (ie_i < 3 ? ie_i : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a' + int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b' + int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c' + int eSgn_i = (ie_i < 3 ? 1 : -1); - ee = epsilon_sgn[ie_f] * epsilon_sgn[ie_i]; + ee = Real(eSgn_f * eSgn_i); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s]; for (int alpha_f=0; alpha_f::BaryonGamma3pt( const Gamma GammaBf, SpinMatrixField &stn_corr) { + assert(Ns==4 && "Baryon code only implemented for N_spin = 4"); + assert(Nc==3 && "Baryon code only implemented for N_colour = 3"); + GridBase *grid = q_tf.Grid(); - autoView( vcorr, stn_corr, CpuWrite); - autoView( vq_ti , q_ti, CpuRead); - autoView( vq_tf , q_tf, CpuRead); + // autoView( vcorr, stn_corr, CpuWrite); + // autoView( vq_ti , q_ti, CpuRead); + // autoView( vq_tf , q_tf, CpuRead); + + // if (group == 1) { + // accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + // auto Dq_ti = vq_ti[ss]; + // auto Dq_tf = vq_tf[ss]; + // sobj result=Zero(); + // BaryonGamma3ptGroup1Site(Dq_ti,Dq_spec1,Dq_spec2,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); + // vcorr[ss] += result; + // });//end loop over lattice sites + // } else if (group == 2) { + // accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + // auto Dq_ti = vq_ti[ss]; + // auto Dq_tf = vq_tf[ss]; + // sobj result=Zero(); + // BaryonGamma3ptGroup2Site(Dq_spec1,Dq_ti,Dq_spec2,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); + // vcorr[ss] += result; + // });//end loop over lattice sites + // } else if (group == 3) { + // accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + // auto Dq_ti = vq_ti[ss]; + // auto Dq_tf = vq_tf[ss]; + // sobj result=Zero(); + // BaryonGamma3ptGroup3Site(Dq_spec1,Dq_spec2,Dq_ti,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); + + // vcorr[ss] += result; + // });//end loop over lattice sites + // } + + autoView( vcorr , stn_corr , AcceleratorWrite); + autoView( vq_ti , q_ti , AcceleratorRead); + autoView( vq_tf , q_tf , AcceleratorRead); if (group == 1) { accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - auto Dq_ti = vq_ti[ss]; - auto Dq_tf = vq_tf[ss]; - sobj result=Zero(); - BaryonGamma3ptGroup1Site(Dq_ti,Dq_spec1,Dq_spec2,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); - vcorr[ss] += result; + auto Dq_ti = vq_ti(ss); + auto Dq_tf = vq_tf(ss); + //sobj result=Zero(); + typedef decltype(coalescedRead(vcorr[0])) spinor; + spinor result=Zero(); + //BaryonGamma3ptGroup1Site(Dq_ti,Dq_spec1,Dq_spec2,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); + BaryonGamma3ptGroup1Site(Dq_ti,Dq_spec1,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); //WRONG + // vcorr[ss] += result; + coalescedWrite(vcorr[ss],result); });//end loop over lattice sites + } else if (group == 2) { accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - auto Dq_ti = vq_ti[ss]; - auto Dq_tf = vq_tf[ss]; - sobj result=Zero(); - BaryonGamma3ptGroup2Site(Dq_spec1,Dq_ti,Dq_spec2,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); - vcorr[ss] += result; + auto Dq_ti = vq_ti(ss); + auto Dq_tf = vq_tf(ss); + //sobj result=Zero(); + typedef decltype(coalescedRead(vcorr[0])) spinor; + spinor result=Zero(); + // BaryonGamma3ptGroup2Site(Dq_spec1,Dq_ti,Dq_spec2,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); + BaryonGamma3ptGroup2Site(Dq_spec1,Dq_ti,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); //WRONG + // vcorr[ss] += result; + coalescedWrite(vcorr[ss],result); });//end loop over lattice sites } else if (group == 3) { accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - auto Dq_ti = vq_ti[ss]; - auto Dq_tf = vq_tf[ss]; - sobj result=Zero(); - BaryonGamma3ptGroup3Site(Dq_spec1,Dq_spec2,Dq_ti,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); - - vcorr[ss] += result; + auto Dq_ti = vq_ti(ss); + auto Dq_tf = vq_tf(ss); + //sobj result=Zero(); + typedef decltype(coalescedRead(vcorr[0])) spinor; + spinor result=Zero(); + //BaryonGamma3ptGroup3Site(Dq_spec1,Dq_spec2,Dq_ti,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); + BaryonGamma3ptGroup3Site(Dq_spec1,Dq_ti,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); //WRONG + // vcorr[ss] += result; + coalescedWrite(vcorr[ss],result); });//end loop over lattice sites } + } /*********************************************************************** * End of BaryonGamma3pt-function code. * - * * + * * * The following code is for Sigma -> N rare hypeon decays * **********************************************************************/ @@ -997,49 +1039,60 @@ void BaryonUtils::BaryonGamma3pt( * Dd_tf is a quark line from t_f to t_H * Ds_ti is a quark line from t_i to t_H */ template -template +template accelerator_inline void BaryonUtils::SigmaToNucleonQ1EyeSite(const mobj &Dq_loop, - const mobj2 &Du_spec, - const mobj &Dd_tf, - const mobj &Ds_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, - robj &result) + const mobj2 &Du_spec, + const mobj &Dd_tf, + const mobj &Ds_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + robj &result) { Gamma g5(Gamma::Algebra::Gamma5); - auto DuG = Du_spec * GammaB_nucl; - // Gamma^B * Ds * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5) - auto GDsGDd = GammaB_sigma * Ds_ti * Gamma_H * g5 * adj(Dd_tf) * g5; - // Dq_loop * \gamma_\mu^L - auto DqG = Dq_loop * Gamma_H; + //auto Gn_adjDd_GH_Ds = GammaB_nucl * g5 * adj(Dd_tf) * g5 * Gamma_H * Ds_ti; + auto adjDd_GH_Ds = g5 * adj(Dd_tf) * g5 * Gamma_H * Ds_ti; + auto Gn_adjDd_GH_Ds = GammaB_nucl * adjDd_GH_Ds; + auto Du_Gs = Du_spec * GammaB_sigma; + auto Dq_GH = Dq_loop * Gamma_H; + auto Tr_Dq_GH = trace(Dq_GH)()()(); + + Real ee; for (int ie_n=0; ie_n < 6 ; ie_n++){ - int a_n = epsilon[ie_n][0]; //a - int b_n = epsilon[ie_n][1]; //b - int c_n = epsilon[ie_n][2]; //c - for (int ie_s=0; ie_s < 6 ; ie_s++){ - int a_s = epsilon[ie_s][0]; //a' - int b_s = epsilon[ie_s][1]; //b' - int c_s = epsilon[ie_s][2]; //c' - for (int alpha_s=0; alpha_s::SigmaToNucleonQ1EyeSite(const mobj &Dq_loop, * Dd_tf is a quark line from t_f to t_H * Ds_ti is a quark line from t_i to t_H */ template -template +template accelerator_inline void BaryonUtils::SigmaToNucleonQ1NonEyeSite(const mobj &Du_ti, - const mobj &Du_tf, - const mobj2 &Du_spec, - const mobj &Dd_tf, - const mobj &Ds_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, - robj &result) + const mobj &Du_tf, + const mobj2 &Du_spec, + const mobj &Dd_tf, + const mobj &Ds_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + robj &result) { Gamma g5(Gamma::Algebra::Gamma5); - auto DuG = Du_spec * GammaB_nucl; - auto adjDu = g5 * adj(Du_tf) * g5; - auto adjDuG = adjDu * GammaB_nucl; - // Gamma^B * Ds * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5) - auto GDsGDd = GammaB_sigma * Ds_ti * Gamma_H * g5 * adj(Dd_tf) * g5; - // Dq_loop * \gamma_\mu^L - auto DuGH = Du_ti * Gamma_H; + auto Du_Gs = Du_spec * GammaB_sigma; + //auto Gn_adjDd_GH_Ds = GammaB_nucl * g5 * adj(Dd_tf) * g5 * Gamma_H * Ds_ti; + auto adjDd_GH_Ds = g5 * adj(Dd_tf) * g5 * Gamma_H * Ds_ti; + auto Gn_adjDd_GH_Ds = GammaB_nucl * adjDd_GH_Ds; + auto adjDu_GH_Du = g5 * adj(Du_tf) * g5 * Gamma_H * Du_ti; + auto adjDu_GH_Du_Gs = adjDu_GH_Du * GammaB_sigma; + + Real ee; for (int ie_n=0; ie_n < 6 ; ie_n++){ - int a_n = epsilon[ie_n][0]; //a - int b_n = epsilon[ie_n][1]; //b - int c_n = epsilon[ie_n][2]; //c - for (int ie_s=0; ie_s < 6 ; ie_s++){ - int a_s = epsilon[ie_s][0]; //a' - int b_s = epsilon[ie_s][1]; //b' - int c_s = epsilon[ie_s][2]; //c' - for (int alpha_s=0; alpha_s::SigmaToNucleonQ1NonEyeSite(const mobj &Du_ti, * Dd_tf is a quark line from t_f to t_H * Ds_ti is a quark line from t_i to t_H */ template -template +template accelerator_inline void BaryonUtils::SigmaToNucleonQ2EyeSite(const mobj &Dq_loop, - const mobj2 &Du_spec, - const mobj &Dd_tf, - const mobj &Ds_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, - robj &result) + const mobj2 &Du_spec, + const mobj &Dd_tf, + const mobj &Ds_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + robj &result) { Gamma g5(Gamma::Algebra::Gamma5); - auto DuG = Du_spec * GammaB_nucl; - // Gamma^B * Ds * \gamma_\mu^L - auto GDsG = GammaB_sigma * Ds_ti * Gamma_H; - // Dq_loop * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5) - auto DqGDd = Dq_loop * Gamma_H * g5 * adj(Dd_tf) * g5; + //auto Gn_adjDd_GH_Duloop_GH_Ds = GammaB_nucl * g5 * adj(Dd_tf) * g5 * Gamma_H * Dq_loop * Gamma_H * Ds_ti; + auto adjDd_GH_Duloop_GH_Ds = g5 * adj(Dd_tf) * g5 * Gamma_H * Dq_loop * Gamma_H * Ds_ti; + auto Gn_adjDd_GH_Duloop_GH_Ds = GammaB_nucl * adjDd_GH_Duloop_GH_Ds; + auto Du_Gs = Du_spec * GammaB_sigma; + + Real ee; for (int ie_n=0; ie_n < 6 ; ie_n++){ - int a_n = epsilon[ie_n][0]; //a - int b_n = epsilon[ie_n][1]; //b - int c_n = epsilon[ie_n][2]; //c - for (int ie_s=0; ie_s < 6 ; ie_s++){ - int a_s = epsilon[ie_s][0]; //a' - int b_s = epsilon[ie_s][1]; //b' - int c_s = epsilon[ie_s][2]; //c' - for (int alpha_s=0; alpha_s::SigmaToNucleonQ2EyeSite(const mobj &Dq_loop, * Dd_tf is a quark line from t_f to t_H * Ds_ti is a quark line from t_i to t_H */ template -template +template accelerator_inline void BaryonUtils::SigmaToNucleonQ2NonEyeSite(const mobj &Du_ti, - const mobj &Du_tf, - const mobj2 &Du_spec, - const mobj &Dd_tf, - const mobj &Ds_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, - robj &result) + const mobj &Du_tf, + const mobj2 &Du_spec, + const mobj &Dd_tf, + const mobj &Ds_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + robj &result) { Gamma g5(Gamma::Algebra::Gamma5); - auto DuG = Du_spec * GammaB_nucl; - auto adjDu = g5 * adj(Du_tf) * g5; - auto adjDuG = adjDu * GammaB_nucl; - // Gamma^B * Ds * \gamma_\mu^L - auto GDsG = GammaB_sigma * Ds_ti * Gamma_H; - // Du * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5) - auto DuGDd = Du_ti * Gamma_H * g5 * adj(Dd_tf) * g5; + auto Du_Gs = Du_spec * GammaB_sigma; + auto adjDu_GH_Ds = g5 * adj(Du_tf) * g5 * Gamma_H * Ds_ti; + //auto Gn_adjDd_GH_Du = GammaB_nucl * g5 * adj(Dd_tf) * g5 * Gamma_H * Du_ti; + auto adjDd_GH_Du = g5 * adj(Dd_tf) * g5 * Gamma_H * Du_ti; + auto Gn_adjDd_GH_Du = GammaB_nucl * adjDd_GH_Du; // for some reason I needed to split this into two lines to avoid the compilation error 'error: identifier "Grid::Gamma::mul" is undefined in device code' + + auto Gn_adjDd_GH_Du_Gs = Gn_adjDd_GH_Du * GammaB_sigma; + + Real ee; for (int ie_n=0; ie_n < 6 ; ie_n++){ - int a_n = epsilon[ie_n][0]; //a - int b_n = epsilon[ie_n][1]; //b - int c_n = epsilon[ie_n][2]; //c - for (int ie_s=0; ie_s < 6 ; ie_s++){ - int a_s = epsilon[ie_s][0]; //a' - int b_s = epsilon[ie_s][1]; //b' - int c_s = epsilon[ie_s][2]; //c' - for (int alpha_s=0; alpha_s template void BaryonUtils::SigmaToNucleonEye(const PropagatorField &qq_loop, - const mobj &Du_spec, - const PropagatorField &qd_tf, - const PropagatorField &qs_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, - const std::string op, - SpinMatrixField &stn_corr) + const mobj &Du_spec, + const PropagatorField &qd_tf, + const PropagatorField &qs_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + const std::string op, + SpinMatrixField &stn_corr) { assert(Ns==4 && "Baryon code only implemented for N_spin = 4"); @@ -1229,39 +1316,43 @@ void BaryonUtils::SigmaToNucleonEye(const PropagatorField &qq_loop, GridBase *grid = qs_ti.Grid(); - autoView( vcorr, stn_corr, CpuWrite); - autoView( vq_loop , qq_loop, CpuRead); - autoView( vd_tf , qd_tf, CpuRead); - autoView( vs_ti , qs_ti, CpuRead); + autoView( vcorr , stn_corr , AcceleratorWrite); + autoView( vq_loop , qq_loop , AcceleratorRead); + autoView( vd_tf , qd_tf , AcceleratorRead); + autoView( vs_ti , qs_ti , AcceleratorRead); + + bool doQ1 = (op == "Q1"); + bool doQ2 = (op == "Q2"); accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - auto Dq_loop = vq_loop[ss]; - auto Dd_tf = vd_tf[ss]; - auto Ds_ti = vs_ti[ss]; - sobj result=Zero(); - if(op == "Q1"){ + auto Dq_loop = vq_loop(ss); + auto Dd_tf = vd_tf(ss); + auto Ds_ti = vs_ti(ss); + typedef decltype(coalescedRead(vcorr[0])) spinor; + spinor result=Zero(); + if(doQ1){ SigmaToNucleonQ1EyeSite(Dq_loop,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); - } else if(op == "Q2"){ + } else if(doQ2){ SigmaToNucleonQ2EyeSite(Dq_loop,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); } else { assert(0 && "Weak Operator not correctly specified"); } - vcorr[ss] = result; - } );//end loop over lattice sites + coalescedWrite(vcorr[ss],result); + });//end loop over lattice sites } template template void BaryonUtils::SigmaToNucleonNonEye(const PropagatorField &qq_ti, - const PropagatorField &qq_tf, - const mobj &Du_spec, - const PropagatorField &qd_tf, - const PropagatorField &qs_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, - const std::string op, - SpinMatrixField &stn_corr) + const PropagatorField &qq_tf, + const mobj &Du_spec, + const PropagatorField &qd_tf, + const PropagatorField &qs_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + const std::string op, + SpinMatrixField &stn_corr) { assert(Ns==4 && "Baryon code only implemented for N_spin = 4"); @@ -1269,27 +1360,31 @@ void BaryonUtils::SigmaToNucleonNonEye(const PropagatorField &qq_ti, GridBase *grid = qs_ti.Grid(); - autoView( vcorr , stn_corr, CpuWrite); - autoView( vq_ti , qq_ti, CpuRead); - autoView( vq_tf , qq_tf, CpuRead); - autoView( vd_tf , qd_tf, CpuRead); - autoView( vs_ti , qs_ti, CpuRead); - // accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - thread_for(ss,grid->oSites(),{ - auto Dq_ti = vq_ti[ss]; - auto Dq_tf = vq_tf[ss]; - auto Dd_tf = vd_tf[ss]; - auto Ds_ti = vs_ti[ss]; - sobj result=Zero(); - if(op == "Q1"){ + autoView( vcorr , stn_corr , AcceleratorWrite ); + autoView( vq_ti , qq_ti , AcceleratorRead ); + autoView( vq_tf , qq_tf , AcceleratorRead ); + autoView( vd_tf , qd_tf , AcceleratorRead ); + autoView( vs_ti , qs_ti , AcceleratorRead ); + + bool doQ1 = (op == "Q1"); + bool doQ2 = (op == "Q2"); + + accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + auto Dq_ti = vq_ti(ss); + auto Dq_tf = vq_tf(ss); + auto Dd_tf = vd_tf(ss); + auto Ds_ti = vs_ti(ss); + typedef decltype(coalescedRead(vcorr[0])) spinor; + spinor result=Zero(); + if(doQ1){ SigmaToNucleonQ1NonEyeSite(Dq_ti,Dq_tf,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); - } else if(op == "Q2"){ + } else if(doQ2){ SigmaToNucleonQ2NonEyeSite(Dq_ti,Dq_tf,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); } else { assert(0 && "Weak Operator not correctly specified"); } - vcorr[ss] = result; - } );//end loop over lattice sites + coalescedWrite(vcorr[ss],result); + });//end loop over lattice sites } NAMESPACE_END(Grid); From e759367d42c14c2de4eac0349b9800b780b63988 Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Fri, 8 Jan 2021 18:04:50 +0000 Subject: [PATCH 03/23] tested and working --- Grid/qcd/utils/BaryonUtils.h | 906 +++++++++++++++++------------------ 1 file changed, 427 insertions(+), 479 deletions(-) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index 25c71e3a..35358d05 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -40,14 +40,7 @@ public: typedef typename FImpl::FermionField FermionField; typedef typename FImpl::PropagatorField PropagatorField; - typedef typename FImpl::SitePropagator pobj; - typedef typename ComplexField::vector_object vobj; - typedef Lattice> SpinMatrixField; - //typedef typename SpinMatrixField::vector_object sobj; - - //static const int epsilon[6][3] ; - //static const Real epsilon_sgn[6]; private: template accelerator_inline @@ -122,7 +115,7 @@ public: static void BaryonGamma3ptGroup1Site( const mobj &Dq1_ti, const mobj2 &Dq2_spec, - // const mobj2 &Dq3_spec, + const mobj2 &Dq3_spec, const mobj &Dq4_tf, const Gamma GammaJ, const Gamma GammaBi, @@ -134,7 +127,7 @@ public: static void BaryonGamma3ptGroup2Site( const mobj2 &Dq1_spec, const mobj &Dq2_ti, - //const mobj2 &Dq3_spec, + const mobj2 &Dq3_spec, const mobj &Dq4_tf, const Gamma GammaJ, const Gamma GammaBi, @@ -145,7 +138,7 @@ public: template accelerator_inline static void BaryonGamma3ptGroup3Site( const mobj2 &Dq1_spec, - //const mobj2 &Dq2_spec, + const mobj2 &Dq2_spec, const mobj &Dq3_ti, const mobj &Dq4_tf, const Gamma GammaJ, @@ -230,13 +223,7 @@ public: const std::string op, SpinMatrixField &stn_corr); }; -/* -template -const int BaryonUtils::epsilon[6][3] = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}}; -template -const Real BaryonUtils::epsilon_sgn[6] = {1.,1.,1.,-1.,-1.,-1.}; -*/ -//This is the old version +//This computes a baryon contraction on a lattice site, including the spin-trace of the correlation matrix template template accelerator_inline void BaryonUtils::BaryonSite(const mobj &D1, @@ -251,20 +238,20 @@ void BaryonUtils::BaryonSite(const mobj &D1, robj &result) { - Gamma g4(Gamma::Algebra::GammaT); //needed for parity P_\pm = 0.5*(1 \pm \gamma_4) + Gamma g4(Gamma::Algebra::GammaT); //needed for parity P_\pm = 0.5*(1 \pm \gamma_4) - auto D1_GAi = D1 * GammaA_i; - auto D1_GAi_g4 = D1_GAi * g4; - auto D1_GAi_P = 0.5*(D1_GAi + (Real)parity * D1_GAi_g4); - auto GAf_D1_GAi_P = GammaA_f * D1_GAi_P; - auto GBf_D1_GAi_P = GammaB_f * D1_GAi_P; + auto D1_GAi = D1 * GammaA_i; + auto D1_GAi_g4 = D1_GAi * g4; + auto D1_GAi_P = 0.5*(D1_GAi + (Real)parity * D1_GAi_g4); + auto GAf_D1_GAi_P = GammaA_f * D1_GAi_P; + auto GBf_D1_GAi_P = GammaB_f * D1_GAi_P; - auto D2_GBi = D2 * GammaB_i; - auto GBf_D2_GBi = GammaB_f * D2_GBi; - auto GAf_D2_GBi = GammaA_f * D2_GBi; + auto D2_GBi = D2 * GammaB_i; + auto GBf_D2_GBi = GammaB_f * D2_GBi; + auto GAf_D2_GBi = GammaA_f * D2_GBi; - auto GBf_D3 = GammaB_f * D3; - auto GAf_D3 = GammaA_f * D3; + auto GBf_D3 = GammaB_f * D3; + auto GAf_D3 = GammaA_f * D3; Real ee; @@ -273,86 +260,87 @@ void BaryonUtils::BaryonSite(const mobj &D1, int b_f = (ie_f < 3 ? (ie_f+1)%3 : (8-ie_f)%3 ); //epsilon[ie_n][1]; //b int c_f = (ie_f < 3 ? (ie_f+2)%3 : (7-ie_f)%3 ); //epsilon[ie_n][2]; //c int eSgn_f = (ie_f < 3 ? 1 : -1); - for (int ie_i=0; ie_i < 6 ; ie_i++){ - int a_i = (ie_i < 3 ? ie_i : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a' - int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b' - int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c' - int eSgn_i = (ie_i < 3 ? 1 : -1); + for (int ie_i=0; ie_i < 6 ; ie_i++){ + int a_i = (ie_i < 3 ? ie_i : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a' + int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b' + int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c' + int eSgn_i = (ie_i < 3 ? 1 : -1); - ee = Real(eSgn_f * eSgn_i); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s]; - //This is the \delta_{456}^{123} part - if (wick_contraction[0]){ - for (int rho=0; rho::BaryonSiteMatrix(const mobj &D1, robj &result) { - auto D1_GAi = D1 * GammaA_i; - auto GAf_D1_GAi = GammaA_f * D1_GAi; - auto GBf_D1_GAi = GammaB_f * D1_GAi; + auto D1_GAi = D1 * GammaA_i; + auto GAf_D1_GAi = GammaA_f * D1_GAi; + auto GBf_D1_GAi = GammaB_f * D1_GAi; - auto D2_GBi = D2 * GammaB_i; - auto GBf_D2_GBi = GammaB_f * D2_GBi; - auto GAf_D2_GBi = GammaA_f * D2_GBi; - - auto GBf_D3 = GammaB_f * D3; - auto GAf_D3 = GammaA_f * D3; + auto D2_GBi = D2 * GammaB_i; + auto GBf_D2_GBi = GammaB_f * D2_GBi; + auto GAf_D2_GBi = GammaA_f * D2_GBi; + auto GBf_D3 = GammaB_f * D3; + auto GAf_D3 = GammaA_f * D3; Real ee; @@ -388,96 +375,101 @@ void BaryonUtils::BaryonSiteMatrix(const mobj &D1, int b_f = (ie_f < 3 ? (ie_f+1)%3 : (8-ie_f)%3 ); //epsilon[ie_n][1]; //b int c_f = (ie_f < 3 ? (ie_f+2)%3 : (7-ie_f)%3 ); //epsilon[ie_n][2]; //c int eSgn_f = (ie_f < 3 ? 1 : -1); - for (int ie_i=0; ie_i < 6 ; ie_i++){ - int a_i = (ie_i < 3 ? ie_i : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a' - int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b' - int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c' - int eSgn_i = (ie_i < 3 ? 1 : -1); + for (int ie_i=0; ie_i < 6 ; ie_i++){ + int a_i = (ie_i < 3 ? ie_i : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a' + int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b' + int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c' + int eSgn_i = (ie_i < 3 ? 1 : -1); - ee = Real(eSgn_f * eSgn_i); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s]; - //This is the \delta_{456}^{123} part - if (wick_contraction[0]){ - for (int rho_i=0; rho_i::ContractBaryons(const PropagatorField &q1_left, for (int ie=0; ie < 6 ; ie++){ if(ie==0 or ie==3){ bytes += grid->oSites() * (4.*sizeof(int) + 4752.*sizeof(vComplex)) * wick_contractions[ie]; - } - else{ + } else{ bytes += grid->oSites() * (64.*sizeof(int) + 5184.*sizeof(vComplex)) * wick_contractions[ie]; } } @@ -653,7 +644,7 @@ template accelerator_inline void BaryonUtils::BaryonGamma3ptGroup1Site( const mobj &Dq1_ti, const mobj2 &Dq2_spec, - // const mobj2 &Dq3_spec, + const mobj2 &Dq3_spec, const mobj &Dq4_tf, const Gamma GammaJ, const Gamma GammaBi, @@ -661,18 +652,14 @@ void BaryonUtils::BaryonGamma3ptGroup1Site( int wick_contraction, robj &result) { - Gamma g5(Gamma::Algebra::Gamma5); - -// auto adjD4_g_D1 = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq1_ti; - auto adjD4 = g5 * adj(Dq4_tf) * g5 ; - auto adjD4_g_D1 = adjD4 * GammaJ * Dq1_ti; - auto Gf_adjD4_g_D1 = GammaBf * adjD4_g_D1; - auto D2_Gi = Dq2_spec * GammaBi; - auto Gf_D2_Gi = GammaBf * D2_Gi; - -// auto Gf_D3 = GammaBf * Dq3_spec; // including a second mobj2 parameter leads to compilation error - auto Gf_D3 = GammaBf * Dq2_spec; //WRONG!!!!! + Gamma g5(Gamma::Algebra::Gamma5); + auto adjD4 = g5 * adj(Dq4_tf) * g5 ; + auto adjD4_g_D1 = adjD4 * GammaJ * Dq1_ti; + auto Gf_adjD4_g_D1 = GammaBf * adjD4_g_D1; + auto D2_Gi = Dq2_spec * GammaBi; + auto Gf_D2_Gi = GammaBf * D2_Gi; + auto Gf_D3 = GammaBf * Dq3_spec; Real ee; @@ -681,65 +668,65 @@ void BaryonUtils::BaryonGamma3ptGroup1Site( int b_f = (ie_f < 3 ? (ie_f+1)%3 : (8-ie_f)%3 ); //epsilon[ie_n][1]; //b int c_f = (ie_f < 3 ? (ie_f+2)%3 : (7-ie_f)%3 ); //epsilon[ie_n][2]; //c int eSgn_f = (ie_f < 3 ? 1 : -1); - for (int ie_i=0; ie_i < 6 ; ie_i++){ - int a_i = (ie_i < 3 ? ie_i : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a' - int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b' - int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c' - int eSgn_i = (ie_i < 3 ? 1 : -1); + for (int ie_i=0; ie_i < 6 ; ie_i++){ + int a_i = (ie_i < 3 ? ie_i : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a' + int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b' + int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c' + int eSgn_i = (ie_i < 3 ? 1 : -1); - ee = Real(eSgn_f * eSgn_i); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s]; + ee = Real(eSgn_f * eSgn_i); - for (int alpha_f=0; alpha_f accelerator_inline void BaryonUtils::BaryonGamma3ptGroup2Site( const mobj2 &Dq1_spec, const mobj &Dq2_ti, - // const mobj2 &Dq3_spec, + const mobj2 &Dq3_spec, const mobj &Dq4_tf, const Gamma GammaJ, const Gamma GammaBi, @@ -759,14 +746,12 @@ void BaryonUtils::BaryonGamma3ptGroup2Site( int wick_contraction, robj &result) { - Gamma g5(Gamma::Algebra::Gamma5); - - auto adjD4_g_D2_Gi = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq2_ti * GammaBi; - auto Gf_adjD4_g_D2_Gi = GammaBf * adjD4_g_D2_Gi; - auto Gf_D1 = GammaBf * Dq1_spec; - //auto Gf_D3 = GammaBf * Dq3_spec; - auto Gf_D3 = GammaBf * Dq1_spec; // WRONG!!!!! + Gamma g5(Gamma::Algebra::Gamma5); + auto adjD4_g_D2_Gi = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq2_ti * GammaBi; + auto Gf_adjD4_g_D2_Gi = GammaBf * adjD4_g_D2_Gi; + auto Gf_D1 = GammaBf * Dq1_spec; + auto Gf_D3 = GammaBf * Dq3_spec; Real ee; @@ -775,64 +760,64 @@ void BaryonUtils::BaryonGamma3ptGroup2Site( int b_f = (ie_f < 3 ? (ie_f+1)%3 : (8-ie_f)%3 ); //epsilon[ie_n][1]; //b int c_f = (ie_f < 3 ? (ie_f+2)%3 : (7-ie_f)%3 ); //epsilon[ie_n][2]; //c int eSgn_f = (ie_f < 3 ? 1 : -1); - for (int ie_i=0; ie_i < 6 ; ie_i++){ - int a_i = (ie_i < 3 ? ie_i : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a' - int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b' - int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c' - int eSgn_i = (ie_i < 3 ? 1 : -1); + for (int ie_i=0; ie_i < 6 ; ie_i++){ + int a_i = (ie_i < 3 ? ie_i : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a' + int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b' + int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c' + int eSgn_i = (ie_i < 3 ? 1 : -1); - ee = Real(eSgn_f * eSgn_i); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s]; + ee = Real(eSgn_f * eSgn_i); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s]; - for (int alpha_f=0; alpha_f template accelerator_inline void BaryonUtils::BaryonGamma3ptGroup3Site( const mobj2 &Dq1_spec, - // const mobj2 &Dq2_spec, + const mobj2 &Dq2_spec, const mobj &Dq3_ti, const mobj &Dq4_tf, const Gamma GammaJ, @@ -852,15 +837,13 @@ void BaryonUtils::BaryonGamma3ptGroup3Site( int wick_contraction, robj &result) { - Gamma g5(Gamma::Algebra::Gamma5); - - auto adjD4_g_D3 = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq3_ti; - auto Gf_adjD4_g_D3 = GammaBf * adjD4_g_D3; - auto Gf_D1 = GammaBf * Dq1_spec; - //auto D2_Gi = Dq2_spec * GammaBi; - auto D2_Gi = Dq1_spec * GammaBi; //WRONG!!!!!!!!!!!!!!!!! - auto Gf_D2_Gi = GammaBf * D2_Gi; + Gamma g5(Gamma::Algebra::Gamma5); + auto adjD4_g_D3 = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq3_ti; + auto Gf_adjD4_g_D3 = GammaBf * adjD4_g_D3; + auto Gf_D1 = GammaBf * Dq1_spec; + auto D2_Gi = Dq2_spec * GammaBi; + auto Gf_D2_Gi = GammaBf * D2_Gi; Real ee; @@ -869,62 +852,64 @@ void BaryonUtils::BaryonGamma3ptGroup3Site( int b_f = (ie_f < 3 ? (ie_f+1)%3 : (8-ie_f)%3 ); //epsilon[ie_n][1]; //b int c_f = (ie_f < 3 ? (ie_f+2)%3 : (7-ie_f)%3 ); //epsilon[ie_n][2]; //c int eSgn_f = (ie_f < 3 ? 1 : -1); - for (int ie_i=0; ie_i < 6 ; ie_i++){ - int a_i = (ie_i < 3 ? ie_i : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a' - int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b' - int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c' - int eSgn_i = (ie_i < 3 ? 1 : -1); + for (int ie_i=0; ie_i < 6 ; ie_i++){ + int a_i = (ie_i < 3 ? ie_i : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a' + int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b' + int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c' + int eSgn_i = (ie_i < 3 ? 1 : -1); - ee = Real(eSgn_f * eSgn_i); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s]; + ee = Real(eSgn_f * eSgn_i); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s]; - for (int alpha_f=0; alpha_f::BaryonGamma3pt( assert(Ns==4 && "Baryon code only implemented for N_spin = 4"); assert(Nc==3 && "Baryon code only implemented for N_colour = 3"); - GridBase *grid = q_tf.Grid(); + GridBase *grid = q_tf.Grid(); - // autoView( vcorr, stn_corr, CpuWrite); - // autoView( vq_ti , q_ti, CpuRead); - // autoView( vq_tf , q_tf, CpuRead); + autoView( vcorr , stn_corr , AcceleratorWrite); + autoView( vq_ti , q_ti , AcceleratorRead); + autoView( vq_tf , q_tf , AcceleratorRead); - // if (group == 1) { - // accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - // auto Dq_ti = vq_ti[ss]; - // auto Dq_tf = vq_tf[ss]; - // sobj result=Zero(); - // BaryonGamma3ptGroup1Site(Dq_ti,Dq_spec1,Dq_spec2,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); - // vcorr[ss] += result; - // });//end loop over lattice sites - // } else if (group == 2) { - // accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - // auto Dq_ti = vq_ti[ss]; - // auto Dq_tf = vq_tf[ss]; - // sobj result=Zero(); - // BaryonGamma3ptGroup2Site(Dq_spec1,Dq_ti,Dq_spec2,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); - // vcorr[ss] += result; - // });//end loop over lattice sites - // } else if (group == 3) { - // accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - // auto Dq_ti = vq_ti[ss]; - // auto Dq_tf = vq_tf[ss]; - // sobj result=Zero(); - // BaryonGamma3ptGroup3Site(Dq_spec1,Dq_spec2,Dq_ti,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); + Vector my_Dq_spec{Dq_spec1,Dq_spec2}; + mobj * Dq_spec_p = &my_Dq_spec[0]; - // vcorr[ss] += result; - // });//end loop over lattice sites - // } + if (group == 1) { + accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + auto Dq_ti = vq_ti(ss); + auto Dq_tf = vq_tf(ss); + typedef decltype(coalescedRead(vcorr[0])) spinor; + spinor result=Zero(); + BaryonGamma3ptGroup1Site(Dq_ti,Dq_spec_p[0],Dq_spec_p[1],Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); + coalescedWrite(vcorr[ss],result); + });//end loop over lattice sites - autoView( vcorr , stn_corr , AcceleratorWrite); - autoView( vq_ti , q_ti , AcceleratorRead); - autoView( vq_tf , q_tf , AcceleratorRead); - - if (group == 1) { - accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - auto Dq_ti = vq_ti(ss); - auto Dq_tf = vq_tf(ss); - //sobj result=Zero(); - typedef decltype(coalescedRead(vcorr[0])) spinor; - spinor result=Zero(); - //BaryonGamma3ptGroup1Site(Dq_ti,Dq_spec1,Dq_spec2,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); - BaryonGamma3ptGroup1Site(Dq_ti,Dq_spec1,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); //WRONG - // vcorr[ss] += result; - coalescedWrite(vcorr[ss],result); - });//end loop over lattice sites - - } else if (group == 2) { - accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - auto Dq_ti = vq_ti(ss); - auto Dq_tf = vq_tf(ss); - //sobj result=Zero(); - typedef decltype(coalescedRead(vcorr[0])) spinor; - spinor result=Zero(); - // BaryonGamma3ptGroup2Site(Dq_spec1,Dq_ti,Dq_spec2,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); - BaryonGamma3ptGroup2Site(Dq_spec1,Dq_ti,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); //WRONG - // vcorr[ss] += result; - coalescedWrite(vcorr[ss],result); - });//end loop over lattice sites - } else if (group == 3) { - accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - auto Dq_ti = vq_ti(ss); - auto Dq_tf = vq_tf(ss); - //sobj result=Zero(); - typedef decltype(coalescedRead(vcorr[0])) spinor; - spinor result=Zero(); - //BaryonGamma3ptGroup3Site(Dq_spec1,Dq_spec2,Dq_ti,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); - BaryonGamma3ptGroup3Site(Dq_spec1,Dq_ti,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); //WRONG - // vcorr[ss] += result; - coalescedWrite(vcorr[ss],result); - });//end loop over lattice sites - } + } else if (group == 2) { + accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + auto Dq_ti = vq_ti(ss); + auto Dq_tf = vq_tf(ss); + typedef decltype(coalescedRead(vcorr[0])) spinor; + spinor result=Zero(); + BaryonGamma3ptGroup2Site(Dq_spec_p[0],Dq_ti,Dq_spec_p[1],Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); + coalescedWrite(vcorr[ss],result); + });//end loop over lattice sites + } else if (group == 3) { + accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + auto Dq_ti = vq_ti(ss); + auto Dq_tf = vq_tf(ss); + typedef decltype(coalescedRead(vcorr[0])) spinor; + spinor result=Zero(); + BaryonGamma3ptGroup3Site(Dq_spec_p[0],Dq_spec_p[1],Dq_ti,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); + coalescedWrite(vcorr[ss],result); + });//end loop over lattice sites + } } @@ -1052,7 +1000,6 @@ void BaryonUtils::SigmaToNucleonQ1EyeSite(const mobj &Dq_loop, Gamma g5(Gamma::Algebra::Gamma5); - //auto Gn_adjDd_GH_Ds = GammaB_nucl * g5 * adj(Dd_tf) * g5 * Gamma_H * Ds_ti; auto adjDd_GH_Ds = g5 * adj(Dd_tf) * g5 * Gamma_H * Ds_ti; auto Gn_adjDd_GH_Ds = GammaB_nucl * adjDd_GH_Ds; auto Du_Gs = Du_spec * GammaB_sigma; @@ -1066,33 +1013,33 @@ void BaryonUtils::SigmaToNucleonQ1EyeSite(const mobj &Dq_loop, int b_n = (ie_n < 3 ? (ie_n+1)%3 : (8-ie_n)%3 ); //epsilon[ie_n][1]; //b int c_n = (ie_n < 3 ? (ie_n+2)%3 : (7-ie_n)%3 ); //epsilon[ie_n][2]; //c int eSgn_n = (ie_n < 3 ? 1 : -1); - for (int ie_s=0; ie_s < 6 ; ie_s++){ - int a_s = (ie_s < 3 ? ie_s : (6-ie_s)%3 ); //epsilon[ie_s][0]; //a' - int b_s = (ie_s < 3 ? (ie_s+1)%3 : (8-ie_s)%3 ); //epsilon[ie_s][1]; //b' - int c_s = (ie_s < 3 ? (ie_s+2)%3 : (7-ie_s)%3 ); //epsilon[ie_s][2]; //c' - int eSgn_s = (ie_s < 3 ? 1 : -1); + for (int ie_s=0; ie_s < 6 ; ie_s++){ + int a_s = (ie_s < 3 ? ie_s : (6-ie_s)%3 ); //epsilon[ie_s][0]; //a' + int b_s = (ie_s < 3 ? (ie_s+1)%3 : (8-ie_s)%3 ); //epsilon[ie_s][1]; //b' + int c_s = (ie_s < 3 ? (ie_s+2)%3 : (7-ie_s)%3 ); //epsilon[ie_s][2]; //c' + int eSgn_s = (ie_s < 3 ? 1 : -1); - ee = Real(eSgn_n * eSgn_s); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s]; + ee = Real(eSgn_n * eSgn_s); + for (int alpha_n=0; alpha_n::SigmaToNucleonQ1NonEyeSite(const mobj &Du_ti, Gamma g5(Gamma::Algebra::Gamma5); auto Du_Gs = Du_spec * GammaB_sigma; - //auto Gn_adjDd_GH_Ds = GammaB_nucl * g5 * adj(Dd_tf) * g5 * Gamma_H * Ds_ti; auto adjDd_GH_Ds = g5 * adj(Dd_tf) * g5 * Gamma_H * Ds_ti; auto Gn_adjDd_GH_Ds = GammaB_nucl * adjDd_GH_Ds; auto adjDu_GH_Du = g5 * adj(Du_tf) * g5 * Gamma_H * Du_ti; @@ -1129,40 +1075,41 @@ void BaryonUtils::SigmaToNucleonQ1NonEyeSite(const mobj &Du_ti, int b_n = (ie_n < 3 ? (ie_n+1)%3 : (8-ie_n)%3 ); //epsilon[ie_n][1]; //b int c_n = (ie_n < 3 ? (ie_n+2)%3 : (7-ie_n)%3 ); //epsilon[ie_n][2]; //c int eSgn_n = (ie_n < 3 ? 1 : -1); - for (int ie_s=0; ie_s < 6 ; ie_s++){ - int a_s = (ie_s < 3 ? ie_s : (6-ie_s)%3 ); //epsilon[ie_s][0]; //a' - int b_s = (ie_s < 3 ? (ie_s+1)%3 : (8-ie_s)%3 ); //epsilon[ie_s][1]; //b' - int c_s = (ie_s < 3 ? (ie_s+2)%3 : (7-ie_s)%3 ); //epsilon[ie_s][2]; //c' - int eSgn_s = (ie_s < 3 ? 1 : -1); + for (int ie_s=0; ie_s < 6 ; ie_s++){ + int a_s = (ie_s < 3 ? ie_s : (6-ie_s)%3 ); //epsilon[ie_s][0]; //a' + int b_s = (ie_s < 3 ? (ie_s+1)%3 : (8-ie_s)%3 ); //epsilon[ie_s][1]; //b' + int c_s = (ie_s < 3 ? (ie_s+2)%3 : (7-ie_s)%3 ); //epsilon[ie_s][2]; //c' + int eSgn_s = (ie_s < 3 ? 1 : -1); - ee = Real(eSgn_n * eSgn_s); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s]; + ee = Real(eSgn_n * eSgn_s); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s]; - for (int alpha_n=0; alpha_n::SigmaToNucleonQ2EyeSite(const mobj &Dq_loop, Gamma g5(Gamma::Algebra::Gamma5); - //auto Gn_adjDd_GH_Duloop_GH_Ds = GammaB_nucl * g5 * adj(Dd_tf) * g5 * Gamma_H * Dq_loop * Gamma_H * Ds_ti; auto adjDd_GH_Duloop_GH_Ds = g5 * adj(Dd_tf) * g5 * Gamma_H * Dq_loop * Gamma_H * Ds_ti; auto Gn_adjDd_GH_Duloop_GH_Ds = GammaB_nucl * adjDd_GH_Duloop_GH_Ds; auto Du_Gs = Du_spec * GammaB_sigma; @@ -1196,32 +1142,33 @@ void BaryonUtils::SigmaToNucleonQ2EyeSite(const mobj &Dq_loop, int b_n = (ie_n < 3 ? (ie_n+1)%3 : (8-ie_n)%3 ); //epsilon[ie_n][1]; //b int c_n = (ie_n < 3 ? (ie_n+2)%3 : (7-ie_n)%3 ); //epsilon[ie_n][2]; //c int eSgn_n = (ie_n < 3 ? 1 : -1); - for (int ie_s=0; ie_s < 6 ; ie_s++){ - int a_s = (ie_s < 3 ? ie_s : (6-ie_s)%3 ); //epsilon[ie_s][0]; //a' - int b_s = (ie_s < 3 ? (ie_s+1)%3 : (8-ie_s)%3 ); //epsilon[ie_s][1]; //b' - int c_s = (ie_s < 3 ? (ie_s+2)%3 : (7-ie_s)%3 ); //epsilon[ie_s][2]; //c' - int eSgn_s = (ie_s < 3 ? 1 : -1); + for (int ie_s=0; ie_s < 6 ; ie_s++){ + int a_s = (ie_s < 3 ? ie_s : (6-ie_s)%3 ); //epsilon[ie_s][0]; //a' + int b_s = (ie_s < 3 ? (ie_s+1)%3 : (8-ie_s)%3 ); //epsilon[ie_s][1]; //b' + int c_s = (ie_s < 3 ? (ie_s+2)%3 : (7-ie_s)%3 ); //epsilon[ie_s][2]; //c' + int eSgn_s = (ie_s < 3 ? 1 : -1); - ee = Real(eSgn_n * eSgn_s); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s]; + ee = Real(eSgn_n * eSgn_s); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s]; - for (int alpha_n=0; alpha_n::SigmaToNucleonQ2NonEyeSite(const mobj &Du_ti, auto Du_Gs = Du_spec * GammaB_sigma; auto adjDu_GH_Ds = g5 * adj(Du_tf) * g5 * Gamma_H * Ds_ti; - //auto Gn_adjDd_GH_Du = GammaB_nucl * g5 * adj(Dd_tf) * g5 * Gamma_H * Du_ti; auto adjDd_GH_Du = g5 * adj(Dd_tf) * g5 * Gamma_H * Du_ti; auto Gn_adjDd_GH_Du = GammaB_nucl * adjDd_GH_Du; // for some reason I needed to split this into two lines to avoid the compilation error 'error: identifier "Grid::Gamma::mul" is undefined in device code' @@ -1259,43 +1205,45 @@ void BaryonUtils::SigmaToNucleonQ2NonEyeSite(const mobj &Du_ti, int b_n = (ie_n < 3 ? (ie_n+1)%3 : (8-ie_n)%3 ); //epsilon[ie_n][1]; //b int c_n = (ie_n < 3 ? (ie_n+2)%3 : (7-ie_n)%3 ); //epsilon[ie_n][2]; //c int eSgn_n = (ie_n < 3 ? 1 : -1); - for (int ie_s=0; ie_s < 6 ; ie_s++){ - int a_s = (ie_s < 3 ? ie_s : (6-ie_s)%3 ); //epsilon[ie_s][0]; //a' - int b_s = (ie_s < 3 ? (ie_s+1)%3 : (8-ie_s)%3 ); //epsilon[ie_s][1]; //b' - int c_s = (ie_s < 3 ? (ie_s+2)%3 : (7-ie_s)%3 ); //epsilon[ie_s][2]; //c' - int eSgn_s = (ie_s < 3 ? 1 : -1); + for (int ie_s=0; ie_s < 6 ; ie_s++){ + int a_s = (ie_s < 3 ? ie_s : (6-ie_s)%3 ); //epsilon[ie_s][0]; //a' + int b_s = (ie_s < 3 ? (ie_s+1)%3 : (8-ie_s)%3 ); //epsilon[ie_s][1]; //b' + int c_s = (ie_s < 3 ? (ie_s+2)%3 : (7-ie_s)%3 ); //epsilon[ie_s][2]; //c' + int eSgn_s = (ie_s < 3 ? 1 : -1); - ee = Real(eSgn_n * eSgn_s); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s]; + ee = Real(eSgn_n * eSgn_s); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s]; - for (int alpha_n=0; alpha_n From 74de2d9742de7a85850cf40cbbe0215edc3405da Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Fri, 8 Jan 2021 18:28:36 +0000 Subject: [PATCH 04/23] whitespace changes --- Grid/qcd/utils/BaryonUtils.h | 156 +++++++++++++++++------------------ 1 file changed, 78 insertions(+), 78 deletions(-) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index 35358d05..80a80a76 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -53,7 +53,7 @@ public: const Gamma GammaB_right, const int parity, const bool * wick_contractions, - robj &result); + robj &result); template accelerator_inline static void BaryonSiteMatrix(const mobj &D1, const mobj &D2, @@ -63,7 +63,7 @@ public: const Gamma GammaA_right, const Gamma GammaB_right, const bool * wick_contractions, - robj &result); + robj &result); public: static void WickContractions(std::string qi, std::string qf, @@ -117,9 +117,9 @@ public: const mobj2 &Dq2_spec, const mobj2 &Dq3_spec, const mobj &Dq4_tf, - const Gamma GammaJ, - const Gamma GammaBi, - const Gamma GammaBf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, int wick_contraction, robj &result); @@ -129,9 +129,9 @@ public: const mobj &Dq2_ti, const mobj2 &Dq3_spec, const mobj &Dq4_tf, - const Gamma GammaJ, - const Gamma GammaBi, - const Gamma GammaBf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, int wick_contraction, robj &result); @@ -141,9 +141,9 @@ public: const mobj2 &Dq2_spec, const mobj &Dq3_ti, const mobj &Dq4_tf, - const Gamma GammaJ, - const Gamma GammaBi, - const Gamma GammaBf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, int wick_contraction, robj &result); public: @@ -155,9 +155,9 @@ public: const PropagatorField &q_tf, int group, int wick_contraction, - const Gamma GammaJ, - const Gamma GammaBi, - const Gamma GammaBf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, SpinMatrixField &stn_corr); private: template accelerator_inline @@ -165,9 +165,9 @@ public: const mobj2 &Du_spec, const mobj &Dd_tf, const mobj &Ds_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, robj &result); template accelerator_inline static void SigmaToNucleonQ1NonEyeSite(const mobj &Du_ti, @@ -175,9 +175,9 @@ public: const mobj2 &Du_spec, const mobj &Dd_tf, const mobj &Ds_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, robj &result); @@ -186,9 +186,9 @@ public: const mobj2 &Du_spec, const mobj &Dd_tf, const mobj &Ds_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, robj &result); template accelerator_inline static void SigmaToNucleonQ2NonEyeSite(const mobj &Du_ti, @@ -196,9 +196,9 @@ public: const mobj2 &Du_spec, const mobj &Dd_tf, const mobj &Ds_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, robj &result); public: template @@ -209,7 +209,7 @@ public: const Gamma Gamma_H, const Gamma GammaB_sigma, const Gamma GammaB_nucl, - const std::string op, + const std::string op, SpinMatrixField &stn_corr); template static void SigmaToNucleonNonEye(const PropagatorField &qq_ti, @@ -220,7 +220,7 @@ public: const Gamma Gamma_H, const Gamma GammaB_sigma, const Gamma GammaB_nucl, - const std::string op, + const std::string op, SpinMatrixField &stn_corr); }; //This computes a baryon contraction on a lattice site, including the spin-trace of the correlation matrix @@ -229,10 +229,10 @@ template accelerator_inline void BaryonUtils::BaryonSite(const mobj &D1, const mobj &D2, const mobj &D3, - const Gamma GammaA_i, - const Gamma GammaB_i, - const Gamma GammaA_f, - const Gamma GammaB_f, + const Gamma GammaA_i, + const Gamma GammaB_i, + const Gamma GammaA_f, + const Gamma GammaB_f, const int parity, const bool * wick_contraction, robj &result) @@ -349,10 +349,10 @@ template accelerator_inline void BaryonUtils::BaryonSiteMatrix(const mobj &D1, const mobj &D2, const mobj &D3, - const Gamma GammaA_i, - const Gamma GammaB_i, - const Gamma GammaA_f, - const Gamma GammaB_f, + const Gamma GammaA_i, + const Gamma GammaB_i, + const Gamma GammaA_f, + const Gamma GammaB_f, const bool * wick_contraction, robj &result) { @@ -496,10 +496,10 @@ template void BaryonUtils::ContractBaryons(const PropagatorField &q1_left, const PropagatorField &q2_left, const PropagatorField &q3_left, - const Gamma GammaA_left, - const Gamma GammaB_left, - const Gamma GammaA_right, - const Gamma GammaB_right, + const Gamma GammaA_left, + const Gamma GammaB_left, + const Gamma GammaA_right, + const Gamma GammaB_right, const bool* wick_contractions, const int parity, ComplexField &baryon_corr) @@ -548,10 +548,10 @@ template void BaryonUtils::ContractBaryonsMatrix(const PropagatorField &q1_left, const PropagatorField &q2_left, const PropagatorField &q3_left, - const Gamma GammaA_left, - const Gamma GammaB_left, - const Gamma GammaA_right, - const Gamma GammaB_right, + const Gamma GammaA_left, + const Gamma GammaB_left, + const Gamma GammaA_right, + const Gamma GammaB_right, const bool* wick_contractions, SpinMatrixField &baryon_corr) { @@ -612,10 +612,10 @@ template void BaryonUtils::ContractBaryonsSlicedMatrix(const mobj &D1, const mobj &D2, const mobj &D3, - const Gamma GammaA_left, - const Gamma GammaB_left, - const Gamma GammaA_right, - const Gamma GammaB_right, + const Gamma GammaA_left, + const Gamma GammaB_left, + const Gamma GammaA_right, + const Gamma GammaB_right, const bool* wick_contractions, const int nt, robj &result) @@ -646,9 +646,9 @@ void BaryonUtils::BaryonGamma3ptGroup1Site( const mobj2 &Dq2_spec, const mobj2 &Dq3_spec, const mobj &Dq4_tf, - const Gamma GammaJ, - const Gamma GammaBi, - const Gamma GammaBf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, int wick_contraction, robj &result) { @@ -740,9 +740,9 @@ void BaryonUtils::BaryonGamma3ptGroup2Site( const mobj &Dq2_ti, const mobj2 &Dq3_spec, const mobj &Dq4_tf, - const Gamma GammaJ, - const Gamma GammaBi, - const Gamma GammaBf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, int wick_contraction, robj &result) { @@ -831,9 +831,9 @@ void BaryonUtils::BaryonGamma3ptGroup3Site( const mobj2 &Dq2_spec, const mobj &Dq3_ti, const mobj &Dq4_tf, - const Gamma GammaJ, - const Gamma GammaBi, - const Gamma GammaBf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, int wick_contraction, robj &result) { @@ -926,9 +926,9 @@ void BaryonUtils::BaryonGamma3pt( const PropagatorField &q_tf, int group, int wick_contraction, - const Gamma GammaJ, - const Gamma GammaBi, - const Gamma GammaBf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, SpinMatrixField &stn_corr) { assert(Ns==4 && "Baryon code only implemented for N_spin = 4"); @@ -992,9 +992,9 @@ void BaryonUtils::SigmaToNucleonQ1EyeSite(const mobj &Dq_loop, const mobj2 &Du_spec, const mobj &Dd_tf, const mobj &Ds_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, robj &result) { @@ -1054,9 +1054,9 @@ void BaryonUtils::SigmaToNucleonQ1NonEyeSite(const mobj &Du_ti, const mobj2 &Du_spec, const mobj &Dd_tf, const mobj &Ds_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, robj &result) { @@ -1123,9 +1123,9 @@ void BaryonUtils::SigmaToNucleonQ2EyeSite(const mobj &Dq_loop, const mobj2 &Du_spec, const mobj &Dd_tf, const mobj &Ds_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, robj &result) { @@ -1183,9 +1183,9 @@ void BaryonUtils::SigmaToNucleonQ2NonEyeSite(const mobj &Du_ti, const mobj2 &Du_spec, const mobj &Dd_tf, const mobj &Ds_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, robj &result) { @@ -1252,9 +1252,9 @@ void BaryonUtils::SigmaToNucleonEye(const PropagatorField &qq_loop, const mobj &Du_spec, const PropagatorField &qd_tf, const PropagatorField &qs_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, const std::string op, SpinMatrixField &stn_corr) { @@ -1296,9 +1296,9 @@ void BaryonUtils::SigmaToNucleonNonEye(const PropagatorField &qq_ti, const mobj &Du_spec, const PropagatorField &qd_tf, const PropagatorField &qs_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, const std::string op, SpinMatrixField &stn_corr) { From 45fc7ded3a18b8fc905127fb4d987a5bb16654e1 Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Tue, 12 Jan 2021 09:10:37 +0000 Subject: [PATCH 05/23] test for sum --- Grid/qcd/utils/BaryonUtils.h | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index 80a80a76..828b4085 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -950,7 +950,8 @@ void BaryonUtils::BaryonGamma3pt( typedef decltype(coalescedRead(vcorr[0])) spinor; spinor result=Zero(); BaryonGamma3ptGroup1Site(Dq_ti,Dq_spec_p[0],Dq_spec_p[1],Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); - coalescedWrite(vcorr[ss],result); + //coalescedWrite(vcorr[ss],vcorr[ss]+result); //diff by factor 10??? + coalescedWrite(vcorr[ss],vcorr[ss]+result); });//end loop over lattice sites } else if (group == 2) { @@ -1271,6 +1272,9 @@ void BaryonUtils::SigmaToNucleonEye(const PropagatorField &qq_loop, bool doQ1 = (op == "Q1"); bool doQ2 = (op == "Q2"); + + Vector my_Dq_spec{Du_spec}; + mobj * Dq_spec_p = &my_Dq_spec[0]; accelerator_for(ss, grid->oSites(), grid->Nsimd(), { auto Dq_loop = vq_loop(ss); @@ -1279,9 +1283,9 @@ void BaryonUtils::SigmaToNucleonEye(const PropagatorField &qq_loop, typedef decltype(coalescedRead(vcorr[0])) spinor; spinor result=Zero(); if(doQ1){ - SigmaToNucleonQ1EyeSite(Dq_loop,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); + SigmaToNucleonQ1EyeSite(Dq_loop,Dq_spec_p[0],Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); } else if(doQ2){ - SigmaToNucleonQ2EyeSite(Dq_loop,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); + SigmaToNucleonQ2EyeSite(Dq_loop,Dq_spec_p[0],Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); } else { assert(0 && "Weak Operator not correctly specified"); } @@ -1316,6 +1320,9 @@ void BaryonUtils::SigmaToNucleonNonEye(const PropagatorField &qq_ti, bool doQ1 = (op == "Q1"); bool doQ2 = (op == "Q2"); + + Vector my_Dq_spec{Du_spec}; + mobj * Dq_spec_p = &my_Dq_spec[0]; accelerator_for(ss, grid->oSites(), grid->Nsimd(), { auto Dq_ti = vq_ti(ss); @@ -1325,9 +1332,9 @@ void BaryonUtils::SigmaToNucleonNonEye(const PropagatorField &qq_ti, typedef decltype(coalescedRead(vcorr[0])) spinor; spinor result=Zero(); if(doQ1){ - SigmaToNucleonQ1NonEyeSite(Dq_ti,Dq_tf,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); + SigmaToNucleonQ1NonEyeSite(Dq_ti,Dq_tf,Dq_spec_p[0],Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); } else if(doQ2){ - SigmaToNucleonQ2NonEyeSite(Dq_ti,Dq_tf,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); + SigmaToNucleonQ2NonEyeSite(Dq_ti,Dq_tf,Dq_spec_p[0],Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); } else { assert(0 && "Weak Operator not correctly specified"); } From fa12b9a3295f1fdf392b0367de4136e576b5401e Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Wed, 13 Jan 2021 10:04:17 +0000 Subject: [PATCH 06/23] bugfix --- Grid/qcd/utils/BaryonUtils.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index 828b4085..69bf8959 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -961,7 +961,7 @@ void BaryonUtils::BaryonGamma3pt( typedef decltype(coalescedRead(vcorr[0])) spinor; spinor result=Zero(); BaryonGamma3ptGroup2Site(Dq_spec_p[0],Dq_ti,Dq_spec_p[1],Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); - coalescedWrite(vcorr[ss],result); + coalescedWrite(vcorr[ss],vcorr[ss]+result); });//end loop over lattice sites } else if (group == 3) { accelerator_for(ss, grid->oSites(), grid->Nsimd(), { @@ -970,7 +970,7 @@ void BaryonUtils::BaryonGamma3pt( typedef decltype(coalescedRead(vcorr[0])) spinor; spinor result=Zero(); BaryonGamma3ptGroup3Site(Dq_spec_p[0],Dq_spec_p[1],Dq_ti,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); - coalescedWrite(vcorr[ss],result); + coalescedWrite(vcorr[ss],vcorr[ss]+result); });//end loop over lattice sites } From 9b73a937e7dbfca0ff0c49b5003c0416e572bb82 Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Mon, 18 Jan 2021 18:57:05 +0000 Subject: [PATCH 07/23] bugfix --- Grid/qcd/utils/BaryonUtils.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index 69bf8959..edc5c8d5 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -951,7 +951,7 @@ void BaryonUtils::BaryonGamma3pt( spinor result=Zero(); BaryonGamma3ptGroup1Site(Dq_ti,Dq_spec_p[0],Dq_spec_p[1],Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); //coalescedWrite(vcorr[ss],vcorr[ss]+result); //diff by factor 10??? - coalescedWrite(vcorr[ss],vcorr[ss]+result); + coalescedWrite(vcorr[ss],coalescedRead(vcorr[ss])+result); });//end loop over lattice sites } else if (group == 2) { @@ -961,7 +961,7 @@ void BaryonUtils::BaryonGamma3pt( typedef decltype(coalescedRead(vcorr[0])) spinor; spinor result=Zero(); BaryonGamma3ptGroup2Site(Dq_spec_p[0],Dq_ti,Dq_spec_p[1],Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); - coalescedWrite(vcorr[ss],vcorr[ss]+result); + coalescedWrite(vcorr[ss],coalescedRead(vcorr[ss])+result); });//end loop over lattice sites } else if (group == 3) { accelerator_for(ss, grid->oSites(), grid->Nsimd(), { @@ -970,7 +970,7 @@ void BaryonUtils::BaryonGamma3pt( typedef decltype(coalescedRead(vcorr[0])) spinor; spinor result=Zero(); BaryonGamma3ptGroup3Site(Dq_spec_p[0],Dq_spec_p[1],Dq_ti,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); - coalescedWrite(vcorr[ss],vcorr[ss]+result); + coalescedWrite(vcorr[ss],coalescedRead(vcorr[ss])+result); });//end loop over lattice sites } From 8bfa0e74f837c914efaf95929149f9ce5b0a5487 Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Tue, 19 Jan 2021 12:27:57 +0000 Subject: [PATCH 08/23] final version, tested on CPU and GPU --- Grid/qcd/utils/BaryonUtils.h | 1 - 1 file changed, 1 deletion(-) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index edc5c8d5..ca8b66ef 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -950,7 +950,6 @@ void BaryonUtils::BaryonGamma3pt( typedef decltype(coalescedRead(vcorr[0])) spinor; spinor result=Zero(); BaryonGamma3ptGroup1Site(Dq_ti,Dq_spec_p[0],Dq_spec_p[1],Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); - //coalescedWrite(vcorr[ss],vcorr[ss]+result); //diff by factor 10??? coalescedWrite(vcorr[ss],coalescedRead(vcorr[ss])+result); });//end loop over lattice sites From fc6d07897fe4a9f7f1a80bd9a7849fe373648b42 Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Tue, 19 Jan 2021 12:32:48 +0000 Subject: [PATCH 09/23] revert changes --- tests/solver/Test_zMADWF_prec.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/solver/Test_zMADWF_prec.cc b/tests/solver/Test_zMADWF_prec.cc index f18e1d86..d1168764 100644 --- a/tests/solver/Test_zMADWF_prec.cc +++ b/tests/solver/Test_zMADWF_prec.cc @@ -52,7 +52,7 @@ struct TestParams{ bool zmobius_inner; double lambda_max; //upper bound of H_T eigenvalue range required to generate zMobius approximation - TestParams(): load_config(false), config_file("ckpoint_lat.1000"), mass(0.01), + TestParams(): load_config(true), config_file("ckpoint_lat.1000"), mass(0.01), Ls_outer(24), b_plus_c_outer(2.0), resid_outer(1e-8), Ls_inner(12), b_plus_c_inner(1.0), resid_inner(1e-8), zmobius_inner(true), lambda_max(1.42), outer_precon("Standard"), inner_precon("Standard") {} @@ -246,7 +246,7 @@ void run(const TestParams ¶ms){ typename RunParamsInner::SchurSolverType SchurSolver_inner(CG_inner); ZeroGuesser Guess; - MADWF > madwf(D_outer, D_inner, PV_outer, SchurSolver_inner, Guess, params.resid_outer, 10000, &update); + MADWF > madwf(D_outer, D_inner, PV_outer, SchurSolver_inner, Guess, params.resid_outer, 100, &update); LatticeFermionD result_MADWF(FGrid_outer); result_MADWF = Zero(); From df16202865e5ff8277b805336b12a054b316ea08 Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Tue, 19 Jan 2021 19:25:27 +0000 Subject: [PATCH 10/23] weird bug in 2pt function... --- Grid/qcd/utils/BaryonUtils.h | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index ca8b66ef..7393c232 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -513,6 +513,7 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_left, GridBase *grid = q1_left.Grid(); autoView(vbaryon_corr , baryon_corr , AcceleratorWrite); + autoView( vcorr_read , baryon_corr , AcceleratorRead); autoView( v1 , q1_left , AcceleratorRead); autoView( v2 , q2_left , AcceleratorRead); autoView( v3 , q3_left , AcceleratorRead); @@ -533,7 +534,8 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_left, auto D1 = v1(ss); auto D2 = v2(ss); auto D3 = v3(ss); - typedef decltype(coalescedRead(vbaryon_corr[0])) cVec; + //typedef decltype(coalescedRead(vbaryon_corr[0])) cVec; + typedef decltype(coalescedRead(vcorr_read[0])) cVec; cVec result=Zero(); BaryonSite(D1,D2,D3,GammaA_left,GammaB_left,GammaA_right,GammaB_right,parity,wick_contractions,result); coalescedWrite(vbaryon_corr[ss],result); @@ -562,6 +564,7 @@ void BaryonUtils::ContractBaryonsMatrix(const PropagatorField &q1_left, GridBase *grid = q1_left.Grid(); autoView(vbaryon_corr , baryon_corr , AcceleratorWrite); + autoView( vcorr_read , baryon_corr , AcceleratorRead); autoView( v1 , q1_left , AcceleratorRead); autoView( v2 , q2_left , AcceleratorRead); autoView( v3 , q3_left , AcceleratorRead); @@ -570,7 +573,8 @@ void BaryonUtils::ContractBaryonsMatrix(const PropagatorField &q1_left, auto D1 = v1(ss); auto D2 = v2(ss); auto D3 = v3(ss); - typedef decltype(coalescedRead(vbaryon_corr[0])) spinor; + //typedef decltype(coalescedRead(vbaryon_corr[0])) spinor; + typedef decltype(coalescedRead(vcorr_read[0])) spinor; spinor result=Zero(); BaryonSiteMatrix(D1,D2,D3,GammaA_left,GammaB_left,GammaA_right,GammaB_right,wick_contractions,result); coalescedWrite(vbaryon_corr[ss],result); @@ -937,6 +941,7 @@ void BaryonUtils::BaryonGamma3pt( GridBase *grid = q_tf.Grid(); autoView( vcorr , stn_corr , AcceleratorWrite); + autoView( vcorr_read , stn_corr , AcceleratorRead); autoView( vq_ti , q_ti , AcceleratorRead); autoView( vq_tf , q_tf , AcceleratorRead); @@ -947,29 +952,29 @@ void BaryonUtils::BaryonGamma3pt( accelerator_for(ss, grid->oSites(), grid->Nsimd(), { auto Dq_ti = vq_ti(ss); auto Dq_tf = vq_tf(ss); - typedef decltype(coalescedRead(vcorr[0])) spinor; + typedef decltype(coalescedRead(vcorr_read[0])) spinor; spinor result=Zero(); BaryonGamma3ptGroup1Site(Dq_ti,Dq_spec_p[0],Dq_spec_p[1],Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); - coalescedWrite(vcorr[ss],coalescedRead(vcorr[ss])+result); + coalescedWrite(vcorr[ss],coalescedRead(vcorr_read[ss])+result); });//end loop over lattice sites } else if (group == 2) { accelerator_for(ss, grid->oSites(), grid->Nsimd(), { auto Dq_ti = vq_ti(ss); auto Dq_tf = vq_tf(ss); - typedef decltype(coalescedRead(vcorr[0])) spinor; + typedef decltype(coalescedRead(vcorr_read[0])) spinor; spinor result=Zero(); BaryonGamma3ptGroup2Site(Dq_spec_p[0],Dq_ti,Dq_spec_p[1],Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); - coalescedWrite(vcorr[ss],coalescedRead(vcorr[ss])+result); + coalescedWrite(vcorr[ss],coalescedRead(vcorr_read[ss])+result); });//end loop over lattice sites } else if (group == 3) { accelerator_for(ss, grid->oSites(), grid->Nsimd(), { auto Dq_ti = vq_ti(ss); auto Dq_tf = vq_tf(ss); - typedef decltype(coalescedRead(vcorr[0])) spinor; + typedef decltype(coalescedRead(vcorr_read[0])) spinor; spinor result=Zero(); BaryonGamma3ptGroup3Site(Dq_spec_p[0],Dq_spec_p[1],Dq_ti,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); - coalescedWrite(vcorr[ss],coalescedRead(vcorr[ss])+result); + coalescedWrite(vcorr[ss],coalescedRead(vcorr_read[ss])+result); });//end loop over lattice sites } From 81d88d9f4df1c9944c8b33db6dd123c36560757d Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Wed, 27 Jan 2021 21:09:51 +0000 Subject: [PATCH 11/23] fixes --- Grid/qcd/utils/BaryonUtils.h | 63 ++++++++++++++++++------------------ 1 file changed, 32 insertions(+), 31 deletions(-) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index 7393c232..d6b48ba0 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -52,7 +52,7 @@ public: const Gamma GammaA_right, const Gamma GammaB_right, const int parity, - const bool * wick_contractions, + const int wick_contractions, robj &result); template accelerator_inline static void BaryonSiteMatrix(const mobj &D1, @@ -62,12 +62,12 @@ public: const Gamma GammaB_left, const Gamma GammaA_right, const Gamma GammaB_right, - const bool * wick_contractions, + const int wick_contractions, robj &result); public: static void WickContractions(std::string qi, std::string qf, - bool* wick_contractions); + int &wick_contractions); static void ContractBaryons(const PropagatorField &q1_left, const PropagatorField &q2_left, const PropagatorField &q3_left, @@ -75,7 +75,7 @@ public: const Gamma GammaB_left, const Gamma GammaA_right, const Gamma GammaB_right, - const bool* wick_contractions, + const int wick_contractions, const int parity, ComplexField &baryon_corr); static void ContractBaryonsMatrix(const PropagatorField &q1_left, @@ -85,7 +85,7 @@ public: const Gamma GammaB_left, const Gamma GammaA_right, const Gamma GammaB_right, - const bool* wick_contractions, + const int wick_contractions, SpinMatrixField &baryon_corr); template static void ContractBaryonsSliced(const mobj &D1, @@ -95,7 +95,7 @@ public: const Gamma GammaB_left, const Gamma GammaA_right, const Gamma GammaB_right, - const bool* wick_contractions, + const int wick_contractions, const int parity, const int nt, robj &result); @@ -107,7 +107,7 @@ public: const Gamma GammaB_left, const Gamma GammaA_right, const Gamma GammaB_right, - const bool* wick_contractions, + const int wick_contractions, const int nt, robj &result); private: @@ -234,7 +234,7 @@ void BaryonUtils::BaryonSite(const mobj &D1, const Gamma GammaA_f, const Gamma GammaB_f, const int parity, - const bool * wick_contraction, + const int wick_contraction, robj &result) { @@ -268,7 +268,7 @@ void BaryonUtils::BaryonSite(const mobj &D1, ee = Real(eSgn_f * eSgn_i); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s]; //This is the \delta_{456}^{123} part - if (wick_contraction[0]){ + if (wick_contraction & 1){ for (int rho=0; rho::BaryonSite(const mobj &D1, } } //This is the \delta_{456}^{231} part - if (wick_contraction[1]){ + if (wick_contraction & 2){ for (int rho=0; rho::BaryonSite(const mobj &D1, }} } //This is the \delta_{456}^{312} part - if (wick_contraction[2]){ + if (wick_contraction & 4){ for (int rho=0; rho::BaryonSite(const mobj &D1, }} } //This is the \delta_{456}^{132} part - if (wick_contraction[3]){ + if (wick_contraction & 8){ for (int rho=0; rho::BaryonSite(const mobj &D1, } } //This is the \delta_{456}^{321} part - if (wick_contraction[4]){ + if (wick_contraction & 16){ for (int rho=0; rho::BaryonSite(const mobj &D1, }} } //This is the \delta_{456}^{213} part - if (wick_contraction[5]){ + if (wick_contraction & 32){ for (int rho=0; rho::BaryonSiteMatrix(const mobj &D1, const Gamma GammaB_i, const Gamma GammaA_f, const Gamma GammaB_f, - const bool * wick_contraction, + const int wick_contraction, robj &result) { @@ -383,7 +383,7 @@ void BaryonUtils::BaryonSiteMatrix(const mobj &D1, ee = Real(eSgn_f * eSgn_i); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s]; //This is the \delta_{456}^{123} part - if (wick_contraction[0]){ + if (wick_contraction & 1){ for (int rho_i=0; rho_i::BaryonSiteMatrix(const mobj &D1, }} } //This is the \delta_{456}^{231} part - if (wick_contraction[1]){ + if (wick_contraction & 2){ for (int rho_i=0; rho_i::BaryonSiteMatrix(const mobj &D1, }} } //This is the \delta_{456}^{312} part - if (wick_contraction[2]){ + if (wick_contraction & 4){ for (int rho_i=0; rho_i::BaryonSiteMatrix(const mobj &D1, }} } //This is the \delta_{456}^{132} part - if (wick_contraction[3]){ + if (wick_contraction & 8){ for (int rho_i=0; rho_i::BaryonSiteMatrix(const mobj &D1, }} } //This is the \delta_{456}^{321} part - if (wick_contraction[4]){ + if (wick_contraction & 16){ for (int rho_i=0; rho_i::BaryonSiteMatrix(const mobj &D1, }} } //This is the \delta_{456}^{213} part - if (wick_contraction[5]){ + if (wick_contraction & 32){ for (int rho_i=0; rho_i::BaryonSiteMatrix(const mobj &D1, * flavours. * * The array wick_contractions must be of length 6 */ template -void BaryonUtils::WickContractions(std::string qi, std::string qf, bool* wick_contractions) { +void BaryonUtils::WickContractions(std::string qi, std::string qf, int &wick_contractions) { + assert(qi.size() == 3 && qf.size() == 3 && "Only sets of 3 quarks accepted."); const int epsilon[6][3] = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}}; + wick_contractions=0; for (int ie=0; ie < 6 ; ie++) { - wick_contractions[ie] = (qi.size() == 3 && qf.size() == 3 - && qi[0] == qf[epsilon[ie][0]] + wick_contractions += ( ( qi[0] == qf[epsilon[ie][0]] && qi[1] == qf[epsilon[ie][1]] - && qi[2] == qf[epsilon[ie][2]]); + && qi[2] == qf[epsilon[ie][2]]) ? 1 : 0) << ie; } } @@ -500,7 +501,7 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_left, const Gamma GammaB_left, const Gamma GammaA_right, const Gamma GammaB_right, - const bool* wick_contractions, + const int wick_contractions, const int parity, ComplexField &baryon_corr) { @@ -522,9 +523,9 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_left, bytes += grid->oSites() * (432.*sizeof(vComplex) + 126.*sizeof(int) + 36.*sizeof(Real)); for (int ie=0; ie < 6 ; ie++){ if(ie==0 or ie==3){ - bytes += grid->oSites() * (4.*sizeof(int) + 4752.*sizeof(vComplex)) * wick_contractions[ie]; + //bytes += grid->oSites() * (4.*sizeof(int) + 4752.*sizeof(vComplex)) * wick_contractions[ie]; } else{ - bytes += grid->oSites() * (64.*sizeof(int) + 5184.*sizeof(vComplex)) * wick_contractions[ie]; + //bytes += grid->oSites() * (64.*sizeof(int) + 5184.*sizeof(vComplex)) * wick_contractions[ie]; } } Real t=0.; @@ -554,7 +555,7 @@ void BaryonUtils::ContractBaryonsMatrix(const PropagatorField &q1_left, const Gamma GammaB_left, const Gamma GammaA_right, const Gamma GammaB_right, - const bool* wick_contractions, + const int wick_contractions, SpinMatrixField &baryon_corr) { @@ -595,7 +596,7 @@ void BaryonUtils::ContractBaryonsSliced(const mobj &D1, const Gamma GammaB_left, const Gamma GammaA_right, const Gamma GammaB_right, - const bool* wick_contractions, + const int wick_contractions, const int parity, const int nt, robj &result) @@ -620,7 +621,7 @@ void BaryonUtils::ContractBaryonsSlicedMatrix(const mobj &D1, const Gamma GammaB_left, const Gamma GammaA_right, const Gamma GammaB_right, - const bool* wick_contractions, + const int wick_contractions, const int nt, robj &result) { From 712bb406502922fb0ceb733e7a74f0ce2b902e2c Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Tue, 15 Dec 2020 16:33:29 +0000 Subject: [PATCH 12/23] merge develop --- tests/solver/Test_zMADWF_prec.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/solver/Test_zMADWF_prec.cc b/tests/solver/Test_zMADWF_prec.cc index d1168764..f18e1d86 100644 --- a/tests/solver/Test_zMADWF_prec.cc +++ b/tests/solver/Test_zMADWF_prec.cc @@ -52,7 +52,7 @@ struct TestParams{ bool zmobius_inner; double lambda_max; //upper bound of H_T eigenvalue range required to generate zMobius approximation - TestParams(): load_config(true), config_file("ckpoint_lat.1000"), mass(0.01), + TestParams(): load_config(false), config_file("ckpoint_lat.1000"), mass(0.01), Ls_outer(24), b_plus_c_outer(2.0), resid_outer(1e-8), Ls_inner(12), b_plus_c_inner(1.0), resid_inner(1e-8), zmobius_inner(true), lambda_max(1.42), outer_precon("Standard"), inner_precon("Standard") {} @@ -246,7 +246,7 @@ void run(const TestParams ¶ms){ typename RunParamsInner::SchurSolverType SchurSolver_inner(CG_inner); ZeroGuesser Guess; - MADWF > madwf(D_outer, D_inner, PV_outer, SchurSolver_inner, Guess, params.resid_outer, 100, &update); + MADWF > madwf(D_outer, D_inner, PV_outer, SchurSolver_inner, Guess, params.resid_outer, 10000, &update); LatticeFermionD result_MADWF(FGrid_outer); result_MADWF = Zero(); From 7905afa9f5b19b147b59eb809b87d1a27e4fc950 Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Tue, 19 Jan 2021 12:32:48 +0000 Subject: [PATCH 13/23] revert changes --- tests/solver/Test_zMADWF_prec.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/solver/Test_zMADWF_prec.cc b/tests/solver/Test_zMADWF_prec.cc index f18e1d86..d1168764 100644 --- a/tests/solver/Test_zMADWF_prec.cc +++ b/tests/solver/Test_zMADWF_prec.cc @@ -52,7 +52,7 @@ struct TestParams{ bool zmobius_inner; double lambda_max; //upper bound of H_T eigenvalue range required to generate zMobius approximation - TestParams(): load_config(false), config_file("ckpoint_lat.1000"), mass(0.01), + TestParams(): load_config(true), config_file("ckpoint_lat.1000"), mass(0.01), Ls_outer(24), b_plus_c_outer(2.0), resid_outer(1e-8), Ls_inner(12), b_plus_c_inner(1.0), resid_inner(1e-8), zmobius_inner(true), lambda_max(1.42), outer_precon("Standard"), inner_precon("Standard") {} @@ -246,7 +246,7 @@ void run(const TestParams ¶ms){ typename RunParamsInner::SchurSolverType SchurSolver_inner(CG_inner); ZeroGuesser Guess; - MADWF > madwf(D_outer, D_inner, PV_outer, SchurSolver_inner, Guess, params.resid_outer, 10000, &update); + MADWF > madwf(D_outer, D_inner, PV_outer, SchurSolver_inner, Guess, params.resid_outer, 100, &update); LatticeFermionD result_MADWF(FGrid_outer); result_MADWF = Zero(); From 96dd7a8fbd66d438618962ca930e7bd1eef34d08 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Mon, 16 Nov 2020 17:15:34 +0100 Subject: [PATCH 14/23] Flop cout matches DiRAC-ITT-2020 --- benchmarks/Benchmark_ITT.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/Benchmark_ITT.cc b/benchmarks/Benchmark_ITT.cc index 032535b3..5d602ce9 100644 --- a/benchmarks/Benchmark_ITT.cc +++ b/benchmarks/Benchmark_ITT.cc @@ -445,7 +445,7 @@ public: // 1344= 3*(2*8+6)*2*8 + 8*3*2*2 + 3*4*2*8 // 1344 = Nc* (6+(Nc-1)*8)*2*Nd + Nd*Nc*2*2 + Nd*Nc*Ns*2 // double flops=(1344.0*volume)/2; -#if 1 +#if 0 double fps = Nc* (6+(Nc-1)*8)*Ns*Nd + Nd*Nc*Ns + Nd*Nc*Ns*2; #else double fps = Nc* (6+(Nc-1)*8)*Ns*Nd + 2*Nd*Nc*Ns + 2*Nd*Nc*Ns*2; From a673b6a54da17b319d8f6707893a9a3f4005be32 Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Thu, 28 Jan 2021 14:15:09 +0000 Subject: [PATCH 15/23] prettify --- Grid/qcd/utils/BaryonUtils.h | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index d6b48ba0..56c5781d 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -513,11 +513,11 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_left, GridBase *grid = q1_left.Grid(); - autoView(vbaryon_corr , baryon_corr , AcceleratorWrite); - autoView( vcorr_read , baryon_corr , AcceleratorRead); - autoView( v1 , q1_left , AcceleratorRead); - autoView( v2 , q2_left , AcceleratorRead); - autoView( v3 , q3_left , AcceleratorRead); + autoView( vbaryon_corr , baryon_corr , AcceleratorWrite); + autoView( vcorr_read , baryon_corr , AcceleratorRead); + autoView( v1 , q1_left , AcceleratorRead); + autoView( v2 , q2_left , AcceleratorRead); + autoView( v3 , q3_left , AcceleratorRead); Real bytes =0.; bytes += grid->oSites() * (432.*sizeof(vComplex) + 126.*sizeof(int) + 36.*sizeof(Real)); @@ -564,11 +564,11 @@ void BaryonUtils::ContractBaryonsMatrix(const PropagatorField &q1_left, GridBase *grid = q1_left.Grid(); - autoView(vbaryon_corr , baryon_corr , AcceleratorWrite); - autoView( vcorr_read , baryon_corr , AcceleratorRead); - autoView( v1 , q1_left , AcceleratorRead); - autoView( v2 , q2_left , AcceleratorRead); - autoView( v3 , q3_left , AcceleratorRead); + autoView( vbaryon_corr , baryon_corr , AcceleratorWrite); + autoView( vcorr_read , baryon_corr , AcceleratorRead); + autoView( v1 , q1_left , AcceleratorRead); + autoView( v2 , q2_left , AcceleratorRead); + autoView( v3 , q3_left , AcceleratorRead); accelerator_for(ss, grid->oSites(), grid->Nsimd(), { auto D1 = v1(ss); @@ -941,10 +941,10 @@ void BaryonUtils::BaryonGamma3pt( GridBase *grid = q_tf.Grid(); - autoView( vcorr , stn_corr , AcceleratorWrite); + autoView( vcorr , stn_corr , AcceleratorWrite); autoView( vcorr_read , stn_corr , AcceleratorRead); - autoView( vq_ti , q_ti , AcceleratorRead); - autoView( vq_tf , q_tf , AcceleratorRead); + autoView( vq_ti , q_ti , AcceleratorRead); + autoView( vq_tf , q_tf , AcceleratorRead); Vector my_Dq_spec{Dq_spec1,Dq_spec2}; mobj * Dq_spec_p = &my_Dq_spec[0]; From bc496dd8440449a2e9d1bc5dce074192640cd094 Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Thu, 28 Jan 2021 14:29:56 +0000 Subject: [PATCH 16/23] change back benchmark_ITT --- benchmarks/Benchmark_ITT.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/Benchmark_ITT.cc b/benchmarks/Benchmark_ITT.cc index 5d602ce9..032535b3 100644 --- a/benchmarks/Benchmark_ITT.cc +++ b/benchmarks/Benchmark_ITT.cc @@ -445,7 +445,7 @@ public: // 1344= 3*(2*8+6)*2*8 + 8*3*2*2 + 3*4*2*8 // 1344 = Nc* (6+(Nc-1)*8)*2*Nd + Nd*Nc*2*2 + Nd*Nc*Ns*2 // double flops=(1344.0*volume)/2; -#if 0 +#if 1 double fps = Nc* (6+(Nc-1)*8)*Ns*Nd + Nd*Nc*Ns + Nd*Nc*Ns*2; #else double fps = Nc* (6+(Nc-1)*8)*Ns*Nd + 2*Nd*Nc*Ns + 2*Nd*Nc*Ns*2; From a5033329247f6410928965ce564b2989c7bac908 Mon Sep 17 00:00:00 2001 From: Michael Marshall <43034299+mmphys@users.noreply.github.com> Date: Sun, 14 Feb 2021 21:27:54 +0000 Subject: [PATCH 17/23] Seems the intention with AutoConf produced Grid/Config.h was to use sed to translate standard PACKAGE_ #defines into GRID_ however due to missing '' after -i this hasn't been working. Perhaps it is too late to fix this, since we don't know who/what is relying on this downstream? ... but if they are, and AutoConf is being used, then likely these #defines have just been redefined anyway. Seems reasonable to redefine PACKAGE and VERSION as well, as none of these macros are used throughout Grid or Hadrons. --- configure.ac | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/configure.ac b/configure.ac index f077ca93..7d95e4e2 100644 --- a/configure.ac +++ b/configure.ac @@ -7,7 +7,7 @@ AM_INIT_AUTOMAKE([subdir-objects 1.13]) AM_EXTRA_RECURSIVE_TARGETS([tests bench]) AC_CONFIG_MACRO_DIR([m4]) AC_CONFIG_SRCDIR([Grid/Grid.h]) -AC_CONFIG_HEADERS([Grid/Config.h],[sed -i 's|PACKAGE_|GRID_|' Grid/Config.h]) +AC_CONFIG_HEADERS([Grid/Config.h],[[sed -i '' -e 's|PACKAGE_|GRID_|' -e 's|[[:space:]]PACKAGE[[:space:]]| GRID_PACKAGE |' -e 's|[[:space:]]VERSION[[:space:]]| GRID_PACKAGE_VERSION |' Grid/Config.h]]) m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])]) ################ Get git info From 35114c9e629c53546c8e95edd6def7b3e1692c7a Mon Sep 17 00:00:00 2001 From: Michael Marshall <43034299+mmphys@users.noreply.github.com> Date: Wed, 17 Feb 2021 13:24:15 +0000 Subject: [PATCH 18/23] Mac OS (Darwin) sed -i flag for in-place editing differs from posix / gnu --- configure.ac | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/configure.ac b/configure.ac index 7d95e4e2..702ce826 100644 --- a/configure.ac +++ b/configure.ac @@ -7,7 +7,12 @@ AM_INIT_AUTOMAKE([subdir-objects 1.13]) AM_EXTRA_RECURSIVE_TARGETS([tests bench]) AC_CONFIG_MACRO_DIR([m4]) AC_CONFIG_SRCDIR([Grid/Grid.h]) -AC_CONFIG_HEADERS([Grid/Config.h],[[sed -i '' -e 's|PACKAGE_|GRID_|' -e 's|[[:space:]]PACKAGE[[:space:]]| GRID_PACKAGE |' -e 's|[[:space:]]VERSION[[:space:]]| GRID_PACKAGE_VERSION |' Grid/Config.h]]) +AC_CONFIG_HEADERS([Grid/Config.h],[[$SED_INPLACE -e 's|PACKAGE_|GRID_|' -e 's|[[:space:]]PACKAGE[[:space:]]| GRID_PACKAGE |' -e 's|[[:space:]]VERSION[[:space:]]| GRID_PACKAGE_VERSION |' Grid/Config.h]], + [if test x"$host_os" == x"${host_os#darwin}" ; then] + [SED_INPLACE="sed -i"] + [else] + [SED_INPLACE="sed -i .bak"] + [fi]) m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])]) ################ Get git info From 86b58d5aff2adac0b41b4fc7be22e77dfbd583e2 Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Thu, 18 Feb 2021 12:04:32 +0000 Subject: [PATCH 19/23] changed if and accelerator_for - no runtime errors any more --- Grid/qcd/utils/BaryonUtils.h | 129 ++++++++++++++++++----------------- 1 file changed, 68 insertions(+), 61 deletions(-) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index 56c5781d..4ac5f685 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -513,19 +513,18 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_left, GridBase *grid = q1_left.Grid(); - autoView( vbaryon_corr , baryon_corr , AcceleratorWrite); - autoView( vcorr_read , baryon_corr , AcceleratorRead); - autoView( v1 , q1_left , AcceleratorRead); - autoView( v2 , q2_left , AcceleratorRead); - autoView( v3 , q3_left , AcceleratorRead); + autoView(vbaryon_corr , baryon_corr , AcceleratorWrite); + autoView( v1 , q1_left , AcceleratorRead); + autoView( v2 , q2_left , AcceleratorRead); + autoView( v3 , q3_left , AcceleratorRead); Real bytes =0.; bytes += grid->oSites() * (432.*sizeof(vComplex) + 126.*sizeof(int) + 36.*sizeof(Real)); for (int ie=0; ie < 6 ; ie++){ if(ie==0 or ie==3){ - //bytes += grid->oSites() * (4.*sizeof(int) + 4752.*sizeof(vComplex)) * wick_contractions[ie]; + bytes += ( wick_contractions & (1 << ie) ) ? grid->oSites() * (4.*sizeof(int) + 4752.*sizeof(vComplex)) : 0.; } else{ - //bytes += grid->oSites() * (64.*sizeof(int) + 5184.*sizeof(vComplex)) * wick_contractions[ie]; + bytes += ( wick_contractions & (1 << ie) ) ? grid->oSites() * (64.*sizeof(int) + 5184.*sizeof(vComplex)) : 0.; } } Real t=0.; @@ -535,8 +534,7 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_left, auto D1 = v1(ss); auto D2 = v2(ss); auto D3 = v3(ss); - //typedef decltype(coalescedRead(vbaryon_corr[0])) cVec; - typedef decltype(coalescedRead(vcorr_read[0])) cVec; + typedef decltype(coalescedRead(vbaryon_corr[0])) cVec; cVec result=Zero(); BaryonSite(D1,D2,D3,GammaA_left,GammaB_left,GammaA_right,GammaB_right,parity,wick_contractions,result); coalescedWrite(vbaryon_corr[ss],result); @@ -561,21 +559,19 @@ void BaryonUtils::ContractBaryonsMatrix(const PropagatorField &q1_left, assert(Ns==4 && "Baryon code only implemented for N_spin = 4"); assert(Nc==3 && "Baryon code only implemented for N_colour = 3"); - + GridBase *grid = q1_left.Grid(); - autoView( vbaryon_corr , baryon_corr , AcceleratorWrite); - autoView( vcorr_read , baryon_corr , AcceleratorRead); - autoView( v1 , q1_left , AcceleratorRead); - autoView( v2 , q2_left , AcceleratorRead); - autoView( v3 , q3_left , AcceleratorRead); + autoView(vbaryon_corr , baryon_corr , AcceleratorWrite); + autoView( v1 , q1_left , AcceleratorRead); + autoView( v2 , q2_left , AcceleratorRead); + autoView( v3 , q3_left , AcceleratorRead); accelerator_for(ss, grid->oSites(), grid->Nsimd(), { auto D1 = v1(ss); auto D2 = v2(ss); auto D3 = v3(ss); - //typedef decltype(coalescedRead(vbaryon_corr[0])) spinor; - typedef decltype(coalescedRead(vcorr_read[0])) spinor; + typedef decltype(coalescedRead(vbaryon_corr[0])) spinor; spinor result=Zero(); BaryonSiteMatrix(D1,D2,D3,GammaA_left,GammaB_left,GammaA_right,GammaB_right,wick_contractions,result); coalescedWrite(vbaryon_corr[ss],result); @@ -941,10 +937,9 @@ void BaryonUtils::BaryonGamma3pt( GridBase *grid = q_tf.Grid(); - autoView( vcorr , stn_corr , AcceleratorWrite); - autoView( vcorr_read , stn_corr , AcceleratorRead); - autoView( vq_ti , q_ti , AcceleratorRead); - autoView( vq_tf , q_tf , AcceleratorRead); + autoView( vcorr , stn_corr , AcceleratorWrite); + autoView( vq_ti , q_ti , AcceleratorRead); + autoView( vq_tf , q_tf , AcceleratorRead); Vector my_Dq_spec{Dq_spec1,Dq_spec2}; mobj * Dq_spec_p = &my_Dq_spec[0]; @@ -953,29 +948,28 @@ void BaryonUtils::BaryonGamma3pt( accelerator_for(ss, grid->oSites(), grid->Nsimd(), { auto Dq_ti = vq_ti(ss); auto Dq_tf = vq_tf(ss); - typedef decltype(coalescedRead(vcorr_read[0])) spinor; + typedef decltype(coalescedRead(vcorr[0])) spinor; spinor result=Zero(); BaryonGamma3ptGroup1Site(Dq_ti,Dq_spec_p[0],Dq_spec_p[1],Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); - coalescedWrite(vcorr[ss],coalescedRead(vcorr_read[ss])+result); + coalescedWrite(vcorr[ss],coalescedRead(vcorr[ss])+result); });//end loop over lattice sites - } else if (group == 2) { accelerator_for(ss, grid->oSites(), grid->Nsimd(), { auto Dq_ti = vq_ti(ss); auto Dq_tf = vq_tf(ss); - typedef decltype(coalescedRead(vcorr_read[0])) spinor; + typedef decltype(coalescedRead(vcorr[0])) spinor; spinor result=Zero(); BaryonGamma3ptGroup2Site(Dq_spec_p[0],Dq_ti,Dq_spec_p[1],Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); - coalescedWrite(vcorr[ss],coalescedRead(vcorr_read[ss])+result); + coalescedWrite(vcorr[ss],coalescedRead(vcorr[ss])+result); });//end loop over lattice sites } else if (group == 3) { accelerator_for(ss, grid->oSites(), grid->Nsimd(), { auto Dq_ti = vq_ti(ss); auto Dq_tf = vq_tf(ss); - typedef decltype(coalescedRead(vcorr_read[0])) spinor; + typedef decltype(coalescedRead(vcorr[0])) spinor; spinor result=Zero(); BaryonGamma3ptGroup3Site(Dq_spec_p[0],Dq_spec_p[1],Dq_ti,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); - coalescedWrite(vcorr[ss],coalescedRead(vcorr_read[ss])+result); + coalescedWrite(vcorr[ss],coalescedRead(vcorr[ss])+result); });//end loop over lattice sites } @@ -1206,6 +1200,7 @@ void BaryonUtils::SigmaToNucleonQ2NonEyeSite(const mobj &Du_ti, Real ee; + for (int ie_n=0; ie_n < 6 ; ie_n++){ int a_n = (ie_n < 3 ? ie_n : (6-ie_n)%3 ); //epsilon[ie_n][0]; //a int b_n = (ie_n < 3 ? (ie_n+1)%3 : (8-ie_n)%3 ); //epsilon[ie_n][1]; //b @@ -1250,6 +1245,7 @@ void BaryonUtils::SigmaToNucleonQ2NonEyeSite(const mobj &Du_ti, }} } } + } template @@ -1275,27 +1271,32 @@ void BaryonUtils::SigmaToNucleonEye(const PropagatorField &qq_loop, autoView( vd_tf , qd_tf , AcceleratorRead); autoView( vs_ti , qs_ti , AcceleratorRead); - bool doQ1 = (op == "Q1"); - bool doQ2 = (op == "Q2"); - Vector my_Dq_spec{Du_spec}; mobj * Dq_spec_p = &my_Dq_spec[0]; - accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - auto Dq_loop = vq_loop(ss); - auto Dd_tf = vd_tf(ss); - auto Ds_ti = vs_ti(ss); - typedef decltype(coalescedRead(vcorr[0])) spinor; - spinor result=Zero(); - if(doQ1){ + if(op == "Q1"){ + accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + auto Dq_loop = vq_loop(ss); + auto Dd_tf = vd_tf(ss); + auto Ds_ti = vs_ti(ss); + typedef decltype(coalescedRead(vcorr[0])) spinor; + spinor result=Zero(); SigmaToNucleonQ1EyeSite(Dq_loop,Dq_spec_p[0],Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); - } else if(doQ2){ + coalescedWrite(vcorr[ss],result); + });//end loop over lattice sites + } else if(op == "Q2"){ + accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + auto Dq_loop = vq_loop(ss); + auto Dd_tf = vd_tf(ss); + auto Ds_ti = vs_ti(ss); + typedef decltype(coalescedRead(vcorr[0])) spinor; + spinor result=Zero(); SigmaToNucleonQ2EyeSite(Dq_loop,Dq_spec_p[0],Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); - } else { - assert(0 && "Weak Operator not correctly specified"); - } - coalescedWrite(vcorr[ss],result); - });//end loop over lattice sites + coalescedWrite(vcorr[ss],result); + });//end loop over lattice sites + } else { + assert(0 && "Weak Operator not correctly specified"); + } } template @@ -1322,29 +1323,35 @@ void BaryonUtils::SigmaToNucleonNonEye(const PropagatorField &qq_ti, autoView( vq_tf , qq_tf , AcceleratorRead ); autoView( vd_tf , qd_tf , AcceleratorRead ); autoView( vs_ti , qs_ti , AcceleratorRead ); - - bool doQ1 = (op == "Q1"); - bool doQ2 = (op == "Q2"); Vector my_Dq_spec{Du_spec}; mobj * Dq_spec_p = &my_Dq_spec[0]; - accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - auto Dq_ti = vq_ti(ss); - auto Dq_tf = vq_tf(ss); - auto Dd_tf = vd_tf(ss); - auto Ds_ti = vs_ti(ss); - typedef decltype(coalescedRead(vcorr[0])) spinor; - spinor result=Zero(); - if(doQ1){ + if(op == "Q1"){ + accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + auto Dq_ti = vq_ti(ss); + auto Dq_tf = vq_tf(ss); + auto Dd_tf = vd_tf(ss); + auto Ds_ti = vs_ti(ss); + typedef decltype(coalescedRead(vcorr[0])) spinor; + spinor result=Zero(); SigmaToNucleonQ1NonEyeSite(Dq_ti,Dq_tf,Dq_spec_p[0],Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); - } else if(doQ2){ + coalescedWrite(vcorr[ss],result); + });//end loop over lattice sites + } else if(op == "Q2"){ + accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + auto Dq_ti = vq_ti(ss); + auto Dq_tf = vq_tf(ss); + auto Dd_tf = vd_tf(ss); + auto Ds_ti = vs_ti(ss); + typedef decltype(coalescedRead(vcorr[0])) spinor; + spinor result=Zero(); SigmaToNucleonQ2NonEyeSite(Dq_ti,Dq_tf,Dq_spec_p[0],Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); - } else { - assert(0 && "Weak Operator not correctly specified"); - } - coalescedWrite(vcorr[ss],result); - });//end loop over lattice sites + coalescedWrite(vcorr[ss],result); + });//end loop over lattice sites + } else { + assert(0 && "Weak Operator not correctly specified"); + } } NAMESPACE_END(Grid); From 7ae030f5851e8a9d74bedb79a28ef598b4006c26 Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Thu, 18 Feb 2021 13:24:50 +0000 Subject: [PATCH 20/23] changed back A2AUtils warning --- Grid/qcd/utils/A2Autils.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Grid/qcd/utils/A2Autils.h b/Grid/qcd/utils/A2Autils.h index 497927dd..b63d8571 100644 --- a/Grid/qcd/utils/A2Autils.h +++ b/Grid/qcd/utils/A2Autils.h @@ -1047,7 +1047,7 @@ A2Autils::ContractWWVV(std::vector &WWVV, { GridBase *grid = vs[0].Grid(); - //int nd = grid->_ndimension; + int nd = grid->_ndimension; int Nsimd = grid->Nsimd(); int N_t = WW_sd.dimensions()[0]; int N_s = WW_sd.dimensions()[1]; From e3d019bc2f46c8c338bc3306a1954180a8b94dd7 Mon Sep 17 00:00:00 2001 From: Daniel Richtmann Date: Mon, 22 Feb 2021 14:56:52 +0100 Subject: [PATCH 21/23] Enable performance counting in WilsonFermion like in others --- .../fermion/implementation/WilsonFermionImplementation.h | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h b/Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h index 4977ea68..84ac25c1 100644 --- a/Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h +++ b/Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h @@ -397,6 +397,7 @@ void WilsonFermion::DhopDerivEO(GaugeField &mat, const FermionField &U, co template void WilsonFermion::Dhop(const FermionField &in, FermionField &out, int dag) { + DhopCalls+=2; conformable(in.Grid(), _grid); // verifies full grid conformable(in.Grid(), out.Grid()); @@ -408,6 +409,7 @@ void WilsonFermion::Dhop(const FermionField &in, FermionField &out, int da template void WilsonFermion::DhopOE(const FermionField &in, FermionField &out, int dag) { + DhopCalls++; conformable(in.Grid(), _cbgrid); // verifies half grid conformable(in.Grid(), out.Grid()); // drops the cb check @@ -420,6 +422,7 @@ void WilsonFermion::DhopOE(const FermionField &in, FermionField &out, int template void WilsonFermion::DhopEO(const FermionField &in, FermionField &out,int dag) { + DhopCalls++; conformable(in.Grid(), _cbgrid); // verifies half grid conformable(in.Grid(), out.Grid()); // drops the cb check From c073e62e0be2f2e421588ad76f90dce33e7154c0 Mon Sep 17 00:00:00 2001 From: Daniel Richtmann Date: Mon, 22 Feb 2021 15:17:07 +0100 Subject: [PATCH 22/23] Correct misleading ac help string --- configure.ac | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/configure.ac b/configure.ac index 702ce826..fb0c78fc 100644 --- a/configure.ac +++ b/configure.ac @@ -130,7 +130,7 @@ esac ############### fermions AC_ARG_ENABLE([fermion-reps], - [AC_HELP_STRING([--fermion-reps=yes|no], [enable extra fermion representation support])], + [AC_HELP_STRING([--enable-fermion-reps=yes|no], [enable extra fermion representation support])], [ac_FERMION_REPS=${enable_fermion_reps}], [ac_FERMION_REPS=yes]) AM_CONDITIONAL(BUILD_FERMION_REPS, [ test "${ac_FERMION_REPS}X" == "yesX" ]) From d5ab571a894682e36bc24488958f01a18983c7f7 Mon Sep 17 00:00:00 2001 From: Christopher Kelly Date: Tue, 23 Feb 2021 11:49:56 -0500 Subject: [PATCH 23/23] Added the ability to apply a custom "filter" to the conjugate momentum in the Integrator classes, applied both after refresh and after applying the forces Added a conjugate momentum "filter" that applies a phase to each site. With sites set to 1.0 or 0.0 this acts as a mask and enables, for example, the freezing of inactive gauge links in DDHMC Added tests/forces/Test_momentum_filter demonstrating the use of the filter to freeze boundary links --- Grid/qcd/hmc/integrators/Integrator.h | 28 ++++ Grid/qcd/hmc/integrators/MomentumFilter.h | 94 +++++++++++++ tests/forces/Test_momentum_filter.cc | 154 ++++++++++++++++++++++ 3 files changed, 276 insertions(+) create mode 100644 Grid/qcd/hmc/integrators/MomentumFilter.h create mode 100644 tests/forces/Test_momentum_filter.cc diff --git a/Grid/qcd/hmc/integrators/Integrator.h b/Grid/qcd/hmc/integrators/Integrator.h index 70055754..77b7de52 100644 --- a/Grid/qcd/hmc/integrators/Integrator.h +++ b/Grid/qcd/hmc/integrators/Integrator.h @@ -33,6 +33,7 @@ directory #define INTEGRATOR_INCLUDED #include +#include "MomentumFilter.h" NAMESPACE_BEGIN(Grid); @@ -78,8 +79,19 @@ protected: RepresentationPolicy Representations; IntegratorParameters Params; + //Filters allow the user to manipulate the conjugate momentum, for example to freeze links in DDHMC + //It is applied whenever the momentum is updated / refreshed + //The default filter does nothing + MomentumFilterBase const* MomFilter; + const ActionSet as; + //Get a pointer to a shared static instance of the "do-nothing" momentum filter to serve as a default + static MomentumFilterBase const* getDefaultMomFilter(){ + static MomentumFilterNone filter; + return &filter; + } + void update_P(Field& U, int level, double ep) { t_P[level] += ep; @@ -135,6 +147,8 @@ protected: // Force from the other representations as[level].apply(update_P_hireps, Representations, Mom, U, ep); + + MomFilter->applyFilter(Mom); } void update_U(Field& U, double ep) @@ -174,11 +188,23 @@ public: t_P.resize(levels, 0.0); t_U = 0.0; // initialization of smearer delegated outside of Integrator + + //Default the momentum filter to "do-nothing" + MomFilter = getDefaultMomFilter(); }; virtual ~Integrator() {} virtual std::string integrator_name() = 0; + + //Set the momentum filter allowing for manipulation of the conjugate momentum + void setMomentumFilter(const MomentumFilterBase &filter){ + MomFilter = &filter; + } + + //Access the conjugate momentum + const MomentaField & getMomentum() const{ return P; } + void print_parameters() { @@ -249,6 +275,8 @@ public: // Refresh the higher representation actions as[level].apply(refresh_hireps, Representations, pRNG); } + + MomFilter->applyFilter(P); } // to be used by the actionlevel class to iterate diff --git a/Grid/qcd/hmc/integrators/MomentumFilter.h b/Grid/qcd/hmc/integrators/MomentumFilter.h new file mode 100644 index 00000000..2a15d80c --- /dev/null +++ b/Grid/qcd/hmc/integrators/MomentumFilter.h @@ -0,0 +1,94 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/hmc/integrators/MomentumFilter.h + +Copyright (C) 2015 + +Author: Christopher Kelly +Author: Peter Boyle + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution +directory +*************************************************************************************/ +/* END LEGAL */ +//-------------------------------------------------------------------- +#ifndef MOMENTUM_FILTER +#define MOMENTUM_FILTER + +NAMESPACE_BEGIN(Grid); + +//These filter objects allow the user to manipulate the conjugate momentum as part of the update / refresh + +template +struct MomentumFilterBase{ + virtual void applyFilter(MomentaField &P) const; +}; + +//Do nothing +template +struct MomentumFilterNone: public MomentumFilterBase{ + void applyFilter(MomentaField &P) const override{} +}; + +//Multiply each site/direction by a Lorentz vector complex number field +//Can be used to implement a mask, zeroing out sites +template +struct MomentumFilterApplyPhase: public MomentumFilterBase{ + typedef typename MomentaField::vector_type vector_type; //SIMD-vectorized complex type + typedef typename MomentaField::scalar_type scalar_type; //scalar complex type + typedef iVector >, Nd > LorentzScalarType; //complex phase for each site/direction + typedef Lattice LatticeLorentzScalarType; + + LatticeLorentzScalarType phase; + + MomentumFilterApplyPhase(const LatticeLorentzScalarType _phase): phase(_phase){} + + //Default to uniform field of (1,0) + MomentumFilterApplyPhase(GridBase* _grid): phase(_grid){ + LorentzScalarType one; + for(int mu=0;mu + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory + *************************************************************************************/ + /* END LEGAL */ +#include + +using namespace std; +using namespace Grid; + +//Get the mu-directected links on the upper boundary and the bulk remainder +template +void getLinksBoundaryBulk(Field &bound, Field &bulk, Field &from, const Coordinate &latt_size){ + bound = Zero(); bulk = Zero(); + for(int mu=0;mu seeds({1,2,3,4}); + + GridParallelRNG pRNG(&Grid); + pRNG.SeedFixedIntegers(seeds); + + typedef PeriodicGimplR Gimpl; + typedef WilsonGaugeAction GaugeAction; + typedef NoHirep Representation; //fundamental + typedef NoSmearing Smearing; + typedef MinimumNorm2 Omelyan; + typedef Gimpl::Field Field; + typedef MomentumFilterApplyPhase Filter; + Filter filter(&Grid); + + //Setup a filter that disables link update on links passing through the global lattice boundary + typedef Filter::LatticeLorentzScalarType MaskType; + typedef Filter::LorentzScalarType MaskSiteType; + + MaskSiteType zero, one; + for(int mu=0;mu::HotConfiguration(pRNG,U); + + //Get the original links on the bulk and boundary for later use + Field Ubnd_orig(&Grid), Ubulk_orig(&Grid); + getLinksBoundaryBulk(Ubnd_orig, Ubulk_orig, U, latt_size); + + ActionSet actions(1); + double beta=6; + GaugeAction gauge_action(beta); + actions[0].push_back(&gauge_action); + + Smearing smear; + IntegratorParameters params(1,1.); //1 MD step + Omelyan integrator(&Grid, params, actions, smear); + + integrator.setMomentumFilter(filter); + + integrator.refresh(U, pRNG); //doesn't actually change the gauge field + + //Check the momentum is zero on the boundary + const auto &P = integrator.getMomentum(); + Field Pbnd(&Grid), Pbulk(&Grid); + getLinksBoundaryBulk(Pbnd, Pbulk, const_cast(P), latt_size); + + RealD Pbnd_nrm = norm2(Pbnd); //expect zero + std::cout << GridLogMessage << "After refresh, norm2 of mu-directed conjugate momentum on boundary is: " << Pbnd_nrm << " (expect 0)" << std::endl; + RealD Pbulk_nrm = norm2(Pbulk); //expect non-zero + std::cout << GridLogMessage << "After refresh, norm2 of bulk conjugate momentum is: " << Pbulk_nrm << " (expect non-zero)" << std::endl; + + //Evolve the gauge field + integrator.integrate(U); + + //Check momentum is still zero on boundary + getLinksBoundaryBulk(Pbnd, Pbulk, const_cast(P), latt_size); + + Pbnd_nrm = norm2(Pbnd); //expect zero + std::cout << GridLogMessage << "After integrate, norm2 of mu-directed conjugate momentum on boundary is: " << Pbnd_nrm << " (expect 0)" << std::endl; + Pbulk_nrm = norm2(Pbulk); //expect non-zero + std::cout << GridLogMessage << "After integrate, norm2 of bulk conjugate momentum is: " << Pbulk_nrm << " (expect non-zero)" << std::endl; + + //Get the new bulk and bound links + Field Ubnd_new(&Grid), Ubulk_new(&Grid); + getLinksBoundaryBulk(Ubnd_new, Ubulk_new, U, latt_size); + + Field Ubnd_diff = Ubnd_new - Ubnd_orig; + Field Ubulk_diff = Ubulk_new - Ubulk_orig; + + RealD Ubnd_change = norm2( Ubnd_diff ); + RealD Ubulk_change = norm2( Ubulk_diff ); + std::cout << GridLogMessage << "After integrate, norm2 of change in mu-directed boundary links is : " << Ubnd_change << " (expect 0)" << std::endl; + std::cout << GridLogMessage << "After integrate, norm2 of change in bulk links is : " << Ubulk_change << " (expect non-zero)" << std::endl; + + Grid_finalize(); +}