#ifndef Hadrons_MContraction_A2AMesonField_hpp_ #define Hadrons_MContraction_A2AMesonField_hpp_ #include #include #include #include #include BEGIN_HADRONS_NAMESPACE /****************************************************************************** * A2AMesonField * ******************************************************************************/ BEGIN_MODULE_NAMESPACE(MContraction) typedef std::pair GammaPair; class A2AMesonFieldPar : Serializable { public: GRID_SERIALIZABLE_CLASS_MEMBERS(A2AMesonFieldPar, int, cacheBlock, int, schurBlock, int, N, int, Nl, std::string, A2A, std::string, output); }; template class TA2AMesonField : public Module { public: FERM_TYPE_ALIASES(FImpl, ); SOLVER_TYPE_ALIASES(FImpl, ); typedef A2AModesSchurDiagTwo A2ABase; public: // constructor TA2AMesonField(const std::string name); // destructor virtual ~TA2AMesonField(void){}; // dependency relation virtual std::vector getInput(void); virtual std::vector getOutput(void); // setup virtual void setup(void); // execution virtual void execute(void); // Arithmetic help. Move to Grid?? virtual void MesonField(Eigen::Tensor &mat, const std::vector &lhs, const std::vector &rhs, std::vector gammas, const std::vector &mom, int orthogdim) ; }; MODULE_REGISTER(A2AMesonField, ARG(TA2AMesonField), MContraction); MODULE_REGISTER(ZA2AMesonField, ARG(TA2AMesonField), MContraction); /****************************************************************************** * TA2AMesonField implementation * ******************************************************************************/ // constructor ///////////////////////////////////////////////////////////////// template TA2AMesonField::TA2AMesonField(const std::string name) : Module(name) { } // dependencies/products /////////////////////////////////////////////////////// template std::vector TA2AMesonField::getInput(void) { std::vector in = {par().A2A + "_class"}; in.push_back(par().A2A + "_w_high_4d"); in.push_back(par().A2A + "_v_high_4d"); return in; } template std::vector TA2AMesonField::getOutput(void) { std::vector out = {}; return out; } // setup /////////////////////////////////////////////////////////////////////// template void TA2AMesonField::setup(void) { auto &a2a = envGet(A2ABase, par().A2A + "_class"); int nt = env().getDim(Tp); int Nl = par().Nl; int N = par().N; int Ls_ = env().getObjectLs(par().A2A + "_class"); // Four D fields envTmp(std::vector, "w", 1, par().schurBlock, FermionField(env().getGrid(1))); envTmp(std::vector, "v", 1, par().schurBlock, FermionField(env().getGrid(1))); // 5D tmp envTmpLat(FermionField, "tmp_5d", Ls_); } ////////////////////////////////////////////////////////////////////////////////// // Cache blocked arithmetic routine // Could move to Grid ??? ////////////////////////////////////////////////////////////////////////////////// template void TA2AMesonField::MesonField(Eigen::Tensor &mat, const std::vector &lhs, const std::vector &rhs, std::vector gammas, const std::vector &mom, int orthogdim) { typedef typename FImpl::SiteSpinor vobj; typedef typename vobj::scalar_object sobj; typedef typename vobj::scalar_type scalar_type; typedef typename vobj::vector_type vector_type; typedef iSpinMatrix SpinMatrix_v; typedef iSpinMatrix SpinMatrix_s; int Lblock = lhs.size(); int Rblock = rhs.size(); GridBase *grid = lhs[0]._grid; const int Nd = grid->_ndimension; const int Nsimd = grid->Nsimd(); int Nt = grid->GlobalDimensions()[orthogdim]; int Ngamma = gammas.size(); int Nmom = mom.size(); int fd=grid->_fdimensions[orthogdim]; int ld=grid->_ldimensions[orthogdim]; int rd=grid->_rdimensions[orthogdim]; // will locally sum vectors first // sum across these down to scalars // splitting the SIMD int MFrvol = rd*Lblock*Rblock*Nmom; int MFlvol = ld*Lblock*Rblock*Nmom; Vector lvSum(MFrvol); parallel_for (int r = 0; r < MFrvol; r++){ lvSum[r] = zero; } Vector lsSum(MFlvol); parallel_for (int r = 0; r < MFlvol; r++){ lsSum[r]=scalar_type(0.0); } int e1= grid->_slice_nblock[orthogdim]; int e2= grid->_slice_block [orthogdim]; int stride=grid->_slice_stride[orthogdim]; std::cout << GridLogMessage << " Entering first parallel loop "<_ostride[orthogdim]; // base offset for start of plane for(int n=0;n > phase(Nmom); for(int m=0;m icoor(Nd); std::vector extracted(Nsimd); for(int i=0;iiCoorFromIindex(icoor,idx); int ldx = rt+icoor[orthogdim]*rd; int ij_ldx = m+Nmom*i+Nmom*Lblock*j+Nmom*Lblock*Rblock*ldx; lsSum[ij_ldx]=lsSum[ij_ldx]+extracted[idx]; } }}} } assert(mat.dimension(0) == Nt); assert(mat.dimension(1) == Nmom); assert(mat.dimension(2) == Ngamma); assert(mat.dimension(3) == Lblock); assert(mat.dimension(4) == Rblock); mat.setZero(); parallel_for(int t=0;t_processor_coor[orthogdim]){ for(int i=0;iGlobalSumVector(&mat(0,0,0,0,0),Nmom*Rblock*Lblock*Nt*Ngamma); return; } // execution /////////////////////////////////////////////////////////////////// template void TA2AMesonField::execute(void) { LOG(Message) << "Computing A2A meson field" << std::endl; auto &a2a = envGet(A2ABase, par().A2A + "_class"); // 2+6+4+4 = 16 gammas // Ordering defined here std::vector gammas ( { Gamma::Algebra::Identity, Gamma::Algebra::Gamma5, Gamma::Algebra::GammaX, Gamma::Algebra::GammaY, Gamma::Algebra::GammaZ, Gamma::Algebra::GammaT, Gamma::Algebra::GammaXGamma5, Gamma::Algebra::GammaYGamma5, Gamma::Algebra::GammaZGamma5, Gamma::Algebra::GammaTGamma5, Gamma::Algebra::SigmaXY, Gamma::Algebra::SigmaXZ, Gamma::Algebra::SigmaXT, Gamma::Algebra::SigmaYZ, Gamma::Algebra::SigmaYT, Gamma::Algebra::SigmaZT }); /////////////////////////////////////////////// // Square assumption for now Nl = Nr = N /////////////////////////////////////////////// int nt = env().getDim(Tp); int N = par().N; int Nl = par().Nl; int ngamma = gammas.size(); /////////////////////////////////////////////// // Momentum setup /////////////////////////////////////////////// std::vector phases(1,env().getGrid(1)); int nmom = phases.size(); phases[0] = Complex(1.0); Eigen::Tensor mesonField (nmom,ngamma,nt,N,N); LOG(Message) << "N = Nh+Nl for A2A MesonField is " << N << std::endl; envGetTmp(std::vector, w); envGetTmp(std::vector, v); envGetTmp(FermionField, tmp_5d); LOG(Message) << "Finding v and w vectors for N = " << N << std::endl; int schurBlock = par().schurBlock; int cacheBlock = par().cacheBlock; for(int i_base=0;i_base mesonFieldBlocked(nmom,ngamma,nt,N_i,N_j); MesonField(mesonFieldBlocked, w, v, gammas, phases,Tp); /////////////////////////////////////////////////////////////// // Copy out to full meson field tensor /////////////////////////////////////////////////////////////// for(int ii =0;ii< N_i;ii++) { for(int jj =0;jj< N_j;jj++) { for(int m =0;m< nmom;m++) { for(int g =0;g< ngamma;g++) { for(int t =0;t< nt;t++) { mesonField(m,g,t,i_base+ii,j_base+jj) = mesonFieldBlocked(m,g,t,ii,jj); }}}}} LOG(Message) << "Contracted MesonFields " <