From 4f0631615fb8351d3e091d935e044a00a2bf75ff Mon Sep 17 00:00:00 2001 From: Vera Guelpers Date: Tue, 30 Apr 2019 12:04:59 +0100 Subject: [PATCH] A2A Lepton-Meson Field contraction --- Grid/qcd/utils/A2Autils.h | 208 +++++++++++++ Hadrons/Modules.hpp | 1 + .../Modules/MContraction/A2ALeptonField.cc | 34 +++ .../Modules/MContraction/A2ALeptonField.hpp | 283 ++++++++++++++++++ Hadrons/modules.inc | 2 + 5 files changed, 528 insertions(+) create mode 100644 Hadrons/Modules/MContraction/A2ALeptonField.cc create mode 100644 Hadrons/Modules/MContraction/A2ALeptonField.hpp diff --git a/Grid/qcd/utils/A2Autils.h b/Grid/qcd/utils/A2Autils.h index 89b4d4bd..6750cbff 100644 --- a/Grid/qcd/utils/A2Autils.h +++ b/Grid/qcd/utils/A2Autils.h @@ -68,6 +68,16 @@ public: const std::vector &emB1, int orthogdim, double *t_kernel = nullptr, double *t_gsum = nullptr); + template + static void LeptonField(TensorType &mat, + const FermionField *lhs_wi, + const FermionField *rhs_vj, + const std::vector &lepton02, + const std::vector &lepton03, + const std::vector &lepton12, + const std::vector &lepton13, + int orthogdim, double *t_kernel = nullptr, double *t_gsum = nullptr); + static void ContractWWVV(std::vector &WWVV, const Eigen::Tensor &WW_sd, const FermionField *vs, @@ -809,6 +819,204 @@ void A2Autils::AslashField(TensorType &mat, if (t_gsum) *t_gsum += usecond(); } +// "Lepton" field w_i(x)^dag * L_mu * g^L_mu * v_j(x) +// +// With V-A current: +// +// g^L_mu = g_mu * (1 - g5) +// +// then in spin space +// +// ( 0 0 L_02 L_03 ) +// L_mu * g^L_mu = ( 0 0 L_12 L_13 ) +// ( 0 0 0 0 ) +// ( 0 0 0 0 ) +// +// +// With +// +// L_02 = 2*L3 + 2*i*L2 +// L_03 = -2*L1 + 2*i*L0 +// L_12 = 2*L1 + 2*i*L0 +// L_13 = 2*L3 - 2*i*L2 +// +// (where mu = (x,y,z,t) ) +// +template +template +void A2Autils::LeptonField(TensorType &mat, + const FermionField *lhs_wi, + const FermionField *rhs_vj, + const std::vector &lepton02, + const std::vector &lepton03, + const std::vector &lepton12, + const std::vector &lepton13, + int orthogdim, double *t_kernel, double *t_gsum) +{ + typedef typename FermionField::vector_object 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; + typedef iSinglet Singlet_v; + typedef iSinglet Singlet_s; + + int Lblock = mat.dimension(3); + int Rblock = mat.dimension(4); + + GridBase *grid = lhs_wi[0]._grid; + + const int Nd = grid->_ndimension; + const int Nsimd = grid->Nsimd(); + + int Nt = grid->GlobalDimensions()[orthogdim]; + + int Nlep = lepton02.size(); + assert(lepton03.size() == Nlep); + assert(lepton12.size() == Nlep); + assert(lepton13.size() == Nlep); + + 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*Nlep; + int MFlvol = ld*Lblock*Rblock*Nlep; + + 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]; + + // Nested parallelism would be ok + // Wasting cores here. Test case r + if (t_kernel) *t_kernel = -usecond(); + parallel_for(int r=0;r_ostride[orthogdim]; // base offset for start of plane + + for(int n=0;n icoor(Nd); + std::vector extracted(Nsimd); + + for(int i=0;i(lvSum[ij_rdx],extracted); + for(int idx=0;idxiCoorFromIindex(icoor,idx); + + int ldx = rt+icoor[orthogdim]*rd; + int ij_ldx = m+Nlep*i+Nlep*Lblock*j+Nlep*Lblock*Rblock*ldx; + + lsSum[ij_ldx]=lsSum[ij_ldx]+extracted[idx]; + } + } + } + if (t_kernel) *t_kernel += usecond(); + + // ld loop and local only?? + int pd = grid->_processors[orthogdim]; + int pc = grid->_processor_coor[orthogdim]; + parallel_for_nest2(int lt=0;ltGlobalSumVector(&mat(0,0,0,0,0),Nlep*Nt*Lblock*Rblock); + if (t_gsum) *t_gsum += usecond(); +} + + //////////////////////////////////////////// // Schematic thoughts about more generalised four quark insertion // diff --git a/Hadrons/Modules.hpp b/Hadrons/Modules.hpp index 898cb532..b0aae2a8 100644 --- a/Hadrons/Modules.hpp +++ b/Hadrons/Modules.hpp @@ -3,6 +3,7 @@ #include #include #include +#include #include #include #include diff --git a/Hadrons/Modules/MContraction/A2ALeptonField.cc b/Hadrons/Modules/MContraction/A2ALeptonField.cc new file mode 100644 index 00000000..3cc3d86a --- /dev/null +++ b/Hadrons/Modules/MContraction/A2ALeptonField.cc @@ -0,0 +1,34 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MContraction/A2ALeptonField.cc + +Copyright (C) 2015-2019 + +Author: Vera Guelpers + +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 Grid; +using namespace Hadrons; +using namespace MContraction; + +template class Grid::Hadrons::MContraction::TA2ALeptonField; diff --git a/Hadrons/Modules/MContraction/A2ALeptonField.hpp b/Hadrons/Modules/MContraction/A2ALeptonField.hpp new file mode 100644 index 00000000..7e041cf5 --- /dev/null +++ b/Hadrons/Modules/MContraction/A2ALeptonField.hpp @@ -0,0 +1,283 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MContraction/A2ALeptonField.cc + +Copyright (C) 2015-2019 + +Author: Vera Guelpers + +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 Hadrons_MContraction_A2ALeptonField_hpp_ +#define Hadrons_MContraction_A2ALeptonField_hpp_ + +#include +#include +#include +#include + +#ifndef MLF_IO_TYPE +#define MLF_IO_TYPE ComplexF +#endif + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * A2ALeptonField * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MContraction) + +class A2ALeptonFieldPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(A2ALeptonFieldPar, + int, cacheBlock, + int, block, + std::string, left, + std::string, right, + std::string, output, + std::vector, lepton); +}; + +class A2ALeptonFieldMetadata: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(A2ALeptonFieldMetadata, + std::string, leptonName, + int, sidx1, + int, sidx2); +}; + +template +class LeptonFieldKernel: public A2AKernel +{ +public: + typedef typename FImpl::FermionField FermionField; +public: + LeptonFieldKernel(const std::vector &lepton02, + const std::vector &lepton03, + const std::vector &lepton12, + const std::vector &lepton13, + GridBase *grid) + : lepton02_(lepton02), lepton03_(lepton03), lepton12_(lepton12), lepton13_(lepton13), grid_(grid) + { + vol_ = 1.; + for (auto &d: grid_->GlobalDimensions()) + { + vol_ *= d; + } + } + + virtual ~LeptonFieldKernel(void) = default; + virtual void operator()(A2AMatrixSet &m, const FermionField *left, + const FermionField *right, + const unsigned int orthogDim, double &t) + { + A2Autils::LeptonField(m, left, right, lepton02_, lepton03_, lepton12_, lepton13_, orthogDim, &t); + } + + virtual double flops(const unsigned int blockSizei, const unsigned int blockSizej) + { + return 0.; + } + + virtual double bytes(const unsigned int blockSizei, const unsigned int blockSizej) + { + return 0.; + } +private: + const std::vector &lepton02_, &lepton03_, &lepton12_, &lepton13_; + GridBase *grid_; + double vol_; +}; + +template +class TA2ALeptonField: public Module +{ +public: + FERM_TYPE_ALIASES(FImpl,); + typedef A2AMatrixBlockComputation Computation; + typedef LeptonFieldKernel Kernel; +public: + // constructor + TA2ALeptonField(const std::string name); + // destructor + virtual ~TA2ALeptonField(void) {}; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_TMP(A2ALeptonField, ARG(TA2ALeptonField), MContraction); + +/****************************************************************************** + * TA2ALeptonField implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TA2ALeptonField::TA2ALeptonField(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TA2ALeptonField::getInput(void) +{ + std::vector in = par().lepton; + in.push_back(par().left); + in.push_back(par().right); + + return in; +} + +template +std::vector TA2ALeptonField::getOutput(void) +{ + std::vector out = {}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TA2ALeptonField::setup(void) +{ + envTmp(Computation, "computation", 1, envGetGrid(FermionField), + env().getNd() - 1, par().lepton.size()*8, 1, par().block, + par().cacheBlock, this); + envTmp(std::vector, "Lmu", 1, + 4, envGetGrid(LatticeComplex)); + envTmp(std::vector, "L02", 1, + par().lepton.size()*8, envGetGrid(LatticeComplex)); + envTmp(std::vector, "L03", 1, + par().lepton.size()*8, envGetGrid(LatticeComplex)); + envTmp(std::vector, "L12", 1, + par().lepton.size()*8, envGetGrid(LatticeComplex)); + envTmp(std::vector, "L13", 1, + par().lepton.size()*8, envGetGrid(LatticeComplex)); + envTmpLat(PropagatorField, "prop_buf"); + +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TA2ALeptonField::execute(void) +{ + auto &left = envGet(std::vector, par().left); + auto &right = envGet(std::vector, par().right); + + int nt = env().getDim().back(); + int N_i = left.size(); + int N_j = right.size(); + int block = par().block; + int cacheBlock = par().cacheBlock; + + LOG(Message) << "Computing all-to-all lepton insertion fields" << std::endl; + LOG(Message) << "Left: '" << par().left << "' Right: '" << par().right << "'" << std::endl; + LOG(Message) << "leptons: " << std::endl; + for (auto &name: par().lepton) + { + LOG(Message) << " " << name << std::endl; + } + LOG(Message) << "lepton insertion field size: " << nt << "*" << N_i << "*" << N_j + << " (filesize " << sizeString(nt*N_i*N_j*sizeof(MLF_IO_TYPE)) + << "/spin index/lepton)" << std::endl; + + envGetTmp(std::vector, L02); + envGetTmp(std::vector, L03); + envGetTmp(std::vector, L12); + envGetTmp(std::vector, L13); + + unsigned int ii = 0; + envGetTmp(PropagatorField, prop_buf); + envGetTmp(std::vector, Lmu); + for (unsigned int i = 0; i < par().lepton.size(); ++i) + { + auto &lepton = envGet(PropagatorField, par().lepton[i]); + // need only half of the spin indices (s1s2) of the lepton field + // the other ones are zero because of the left-handed current + for (unsigned int s1 = 0; s1 < 2; ++s1) + for (unsigned int s2 = 0; s2 < 4; ++s2) + { + for (unsigned int mu = 0; mu < 4; ++mu) + { + prop_buf = GammaL(Gamma::gmu[mu]) * lepton; + //lepon is unit matix in color space, so just pick one diagonal enty + Lmu[mu] = peekColour(peekSpin(prop_buf,s1,s2),0,0); + } + // build the required combinations of lepton fields for the a2a kernel + // [g^L_mu * L_mu]_{ab} = L_{ab} (with mu = (x,y,z,t)) + L02[ii] = 2.0*Lmu[3] + 2.0*timesI(Lmu[2]); + L03[ii] = 2.0*timesI(Lmu[0]) - 2.0*Lmu[1]; + L12[ii] = 2.0*timesI(Lmu[0]) + 2.0*Lmu[1]; + L13[ii] = 2.0*Lmu[3] - 2.0*timesI(Lmu[2]); + ii++; + } + } + + auto ionameFn = [this](const unsigned int index, const unsigned int dummy) + { + //index = 8*l + 4*sindex1 + sindex2 + unsigned int sindex = index % 8; + unsigned int sindex2 = sindex % 4; + unsigned int sindex1 = (sindex - sindex2)/4; + unsigned int l = (index - sindex)/8; + return par().lepton[l] + "_" + std::to_string(sindex1) + std::to_string(sindex2); + }; + + auto filenameFn = [this, &ionameFn](const unsigned int index, const unsigned int dummy) + { + return par().output + "." + std::to_string(vm().getTrajectory()) + + "/" + ionameFn(index, dummy) + ".h5"; + }; + + auto metadataFn = [this](const unsigned int index, const unsigned int dummy) + { + A2ALeptonFieldMetadata md; + + unsigned int sindex = index % 8; + unsigned int sindex2 = sindex % 4; + unsigned int sindex1 = (sindex - sindex2)/4; + unsigned int l = (index - sindex)/8; + md.leptonName = par().lepton[l]; + md.sidx1 = sindex1; + md.sidx2 = sindex2; + + return md; + }; + + // executing computation + Kernel kernel(L02, L03, L12, L13, envGetGrid(FermionField)); + envGetTmp(Computation, computation); + computation.execute(left, right, kernel, ionameFn, filenameFn, metadataFn); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MContraction_A2ALeptonField_hpp_ diff --git a/Hadrons/modules.inc b/Hadrons/modules.inc index 3b9e6563..b428b5eb 100644 --- a/Hadrons/modules.inc +++ b/Hadrons/modules.inc @@ -7,6 +7,7 @@ modules_cc =\ Modules/MContraction/WeakNonEye3pt.cc \ Modules/MContraction/A2AAslashField.cc \ Modules/MContraction/A2AMesonField.cc \ + Modules/MContraction/A2ALeptonField.cc \ Modules/MContraction/DiscLoop.cc \ Modules/MContraction/Gamma3pt.cc \ Modules/MFermion/FreeProp.cc \ @@ -72,6 +73,7 @@ modules_hpp =\ Modules/MContraction/A2AAslashField.hpp \ Modules/MContraction/A2ALoop.hpp \ Modules/MContraction/A2AMesonField.hpp \ + Modules/MContraction/A2ALeptonField.hpp \ Modules/MContraction/Meson.hpp \ Modules/MContraction/DiscLoop.hpp \ Modules/MContraction/Gamma3pt.hpp \