diff --git a/lib/qcd/action/ActionParams.h b/lib/qcd/action/ActionParams.h index dcbdfce8..91e94741 100644 --- a/lib/qcd/action/ActionParams.h +++ b/lib/qcd/action/ActionParams.h @@ -45,6 +45,10 @@ namespace QCD { WilsonImplParams() : overlapCommsCompute(false) {}; }; + struct StaggeredImplParams { + StaggeredImplParams() {}; + }; + struct OneFlavourRationalParams { RealD lo; RealD hi; diff --git a/lib/qcd/action/Actions.h b/lib/qcd/action/Actions.h index ba6e577d..1c4acd86 100644 --- a/lib/qcd/action/Actions.h +++ b/lib/qcd/action/Actions.h @@ -53,6 +53,7 @@ Author: paboyle #include #include #include //used by all wilson type fermions +#include //used by all wilson type fermions //////////////////////////////////////////// // Gauge Actions @@ -108,6 +109,10 @@ typedef SymanzikGaugeAction ConjugateSymanzikGaugeAction //////////////////////////////////////////////////////////////////////////////////////////////////// +#define FermOpStaggeredTemplateInstantiate(A) \ + template class A; \ + template class A; + #define FermOp4dVecTemplateInstantiate(A) \ template class A; \ template class A; \ @@ -147,6 +152,8 @@ typedef SymanzikGaugeAction ConjugateSymanzikGaugeAction //#include +#include + #include // Cayley types #include #include diff --git a/lib/qcd/action/fermion/FermionOperatorImpl.h b/lib/qcd/action/fermion/FermionOperatorImpl.h index 0800dea6..0b162f77 100644 --- a/lib/qcd/action/fermion/FermionOperatorImpl.h +++ b/lib/qcd/action/fermion/FermionOperatorImpl.h @@ -343,7 +343,7 @@ class GparityWilsonImpl : public ConjugateGaugeImpl + class StaggeredImpl : public PeriodicGaugeImpl > { + + public: + + typedef RealD _Coeff_t ; + static const int Dimension = Representation::Dimension; + typedef PeriodicGaugeImpl > Gimpl; + + //Necessary? + constexpr bool is_fundamental() const{return Dimension == Nc ? 1 : 0;} + + const bool LsVectorised=false; + typedef _Coeff_t Coeff_t; + + INHERIT_GIMPL_TYPES(Gimpl); + + template using iImplSpinor = iScalar > >; + template using iImplHalfSpinor = iVector >, Ngp>; + template using iImplDoubledGaugeField = iVector >, Nds>; + + typedef iImplSpinor SiteSpinor; + typedef iImplHalfSpinor SiteHalfSpinor; + typedef iImplDoubledGaugeField SiteDoubledGaugeField; + + typedef Lattice FermionField; + typedef Lattice DoubledGaugeField; + + typedef SimpleCompressor Compressor; + typedef StaggeredImplParams ImplParams; + typedef CartesianStencil StencilImpl; + + ImplParams Params; + + StaggeredImpl(const ImplParams &p = ImplParams()) : Params(p){}; + + inline void multLink(SiteSpinor &phi, + const SiteDoubledGaugeField &U, + const SiteSpinor &chi, + int mu){ + mult(&phi(), &U(mu), &chi()); + } + inline void multLinkAdd(SiteSpinor &phi, + const SiteDoubledGaugeField &U, + const SiteSpinor &chi, + int mu){ + mac(&phi(), &U(mu), &chi()); + } + + template + inline void loadLinkElement(Simd ®, ref &memory) { + reg = memory; + } + + inline void DoubleStore(GridBase *GaugeGrid, + DoubledGaugeField &Uds, + DoubledGaugeField &UUUds, // for Naik term + const GaugeField &Umu) { + conformable(Uds._grid, GaugeGrid); + conformable(Umu._grid, GaugeGrid); + GaugeLinkField U(GaugeGrid); + for (int mu = 0; mu < Nd; mu++) { + U = PeekIndex(Umu, mu); + PokeIndex(Uds, U, mu); + PokeIndex(UUUds, U, mu); + std::cout << GridLogMessage << " NOT created the treble links for staggered yet" <(Uds, U, mu + 4); + PokeIndex(UUUds, U, mu+4); + } + } + + inline void InsertForce4D(GaugeField &mat, FermionField &Btilde, FermionField &A,int mu){ + GaugeLinkField link(mat._grid); + link = TraceIndex(outerProduct(Btilde,A)); + PokeIndex(mat,link,mu); + } + + inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField Ã,int mu){ + assert (0); + // Must never hit + } + }; + + + typedef WilsonImpl WilsonImplR; // Real.. whichever prec typedef WilsonImpl WilsonImplF; // Float typedef WilsonImpl WilsonImplD; // Double @@ -527,6 +618,10 @@ PARALLEL_FOR_LOOP typedef GparityWilsonImpl GparityWilsonImplF; // Float typedef GparityWilsonImpl GparityWilsonImplD; // Double + typedef StaggeredImpl StaggeredImplR; // Real.. whichever prec + typedef StaggeredImpl StaggeredImplF; // Float + typedef StaggeredImpl StaggeredImplD; // Double + }} #endif diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc b/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc new file mode 100644 index 00000000..73cc272a --- /dev/null +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc @@ -0,0 +1,309 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/action/fermion/ImprovedStaggeredFermion.cc + +Copyright (C) 2015 + +Author: Azusa Yamaguchi, 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 */ +#include + +namespace Grid { +namespace QCD { + +const std::vector +ImprovedStaggeredFermionStatic::directions({0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3}); +const std::vector +ImprovedStaggeredFermionStatic::displacements({1, 1, 1, 1, -1, -1, -1, -1, 3, 3, 3, 3, -3, -3, -3, -3}); + +///////////////////////////////// +// Constructor and gauge import +///////////////////////////////// + +template +ImprovedStaggeredFermion::ImprovedStaggeredFermion(GaugeField &_Umu, GridCartesian &Fgrid, + GridRedBlackCartesian &Hgrid, RealD _mass, + const ImplParams &p) + : Kernels(p), + _grid(&Fgrid), + _cbgrid(&Hgrid), + Stencil(&Fgrid, npoint, Even, directions, displacements), + StencilEven(&Hgrid, npoint, Even, directions, displacements), // source is Even + StencilOdd(&Hgrid, npoint, Odd, directions, displacements), // source is Odd + mass(_mass), + Lebesgue(_grid), + LebesgueEvenOdd(_cbgrid), + Umu(&Fgrid), + UmuEven(&Hgrid), + UmuOdd(&Hgrid), + UUUmu(&Fgrid), + UUUmuEven(&Hgrid), + UUUmuOdd(&Hgrid) { + // Allocate the required comms buffer + ImportGauge(_Umu); +} + +template +void ImprovedStaggeredFermion::ImportGauge(const GaugeField &_Umu) { + GaugeField HUmu(_Umu._grid); + HUmu = _Umu * (-0.5); + Impl::DoubleStore(GaugeGrid(), Umu, UUUmu, HUmu); + pickCheckerboard(Even, UmuEven, Umu); + pickCheckerboard(Odd, UmuOdd, Umu); + pickCheckerboard(Even, UUUmuEven, UUUmu); + pickCheckerboard(Odd, UUUmuOdd, UUUmu); +} + +///////////////////////////// +// Implement the interface +///////////////////////////// + +template +RealD ImprovedStaggeredFermion::M(const FermionField &in, FermionField &out) { + out.checkerboard = in.checkerboard; + Dhop(in, out, DaggerNo); + return axpy_norm(out, mass, in, out); +} + +template +RealD ImprovedStaggeredFermion::Mdag(const FermionField &in, FermionField &out) { + out.checkerboard = in.checkerboard; + Dhop(in, out, DaggerYes); + return axpy_norm(out, mass, in, out); +} + +template +void ImprovedStaggeredFermion::Meooe(const FermionField &in, FermionField &out) { + if (in.checkerboard == Odd) { + DhopEO(in, out, DaggerNo); + } else { + DhopOE(in, out, DaggerNo); + } +} +template +void ImprovedStaggeredFermion::MeooeDag(const FermionField &in, FermionField &out) { + if (in.checkerboard == Odd) { + DhopEO(in, out, DaggerYes); + } else { + DhopOE(in, out, DaggerYes); + } +} + +template +void ImprovedStaggeredFermion::Mooee(const FermionField &in, FermionField &out) { + out.checkerboard = in.checkerboard; + typename FermionField::scalar_type scal(mass); + out = scal * in; +} + +template +void ImprovedStaggeredFermion::MooeeDag(const FermionField &in, FermionField &out) { + out.checkerboard = in.checkerboard; + Mooee(in, out); +} + +template +void ImprovedStaggeredFermion::MooeeInv(const FermionField &in, FermionField &out) { + out.checkerboard = in.checkerboard; + out = (1.0 / (mass)) * in; +} + +template +void ImprovedStaggeredFermion::MooeeInvDag(const FermionField &in, + FermionField &out) { + out.checkerboard = in.checkerboard; + MooeeInv(in, out); +} + +/////////////////////////////////// +// Internal +/////////////////////////////////// + +template +void ImprovedStaggeredFermion::DerivInternal(StencilImpl &st, DoubledGaugeField &U, DoubledGaugeField &UUU, + GaugeField & mat, + const FermionField &A, const FermionField &B, int dag) { + assert((dag == DaggerNo) || (dag == DaggerYes)); + + Compressor compressor; + + FermionField Btilde(B._grid); + FermionField Atilde(B._grid); + Atilde = A; + + st.HaloExchange(B, compressor); + + for (int mu = 0; mu < Nd; mu++) { + + //////////////////////// + // Call the single hop + //////////////////////// + PARALLEL_FOR_LOOP + for (int sss = 0; sss < B._grid->oSites(); sss++) { + Kernels::DhopDir(st, U, UUU, st.CommBuf(), sss, sss, B, Btilde, mu,1); + } + + // Force in three link terms + // + // Impl::InsertForce4D(mat, Btilde, Atilde, mu); + // + // dU_ac(x)/dt = i p_ab U_bc(x) + // + // => dS_f/dt = dS_f/dU_ac(x) . dU_ac(x)/dt = i p_ab U_bc(x) dS_f/dU_ac(x) + // + // One link: form fragments S_f = A U B + // + // write Btilde = U(x) B(x+mu) + // + // mat+= TraceIndex(outerProduct(Btilde,A)); + // + // Three link: form fragments S_f = A UUU B + // + // mat+= outer ( A, UUUB) <-- Best take DhopDeriv with one linke or identity matrix + // mat+= outer ( AU, UUB) <-- and then use covariant cshift? + // mat+= outer ( AUU, UB) <-- Returned from call to DhopDir + + assert(0);// need to figure out the force interface with a blasted three link term. + + } +} + +template +void ImprovedStaggeredFermion::DhopDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) { + + conformable(U._grid, _grid); + conformable(U._grid, V._grid); + conformable(U._grid, mat._grid); + + mat.checkerboard = U.checkerboard; + + DerivInternal(Stencil, Umu, UUUmu, mat, U, V, dag); +} + +template +void ImprovedStaggeredFermion::DhopDerivOE(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) { + + conformable(U._grid, _cbgrid); + conformable(U._grid, V._grid); + conformable(U._grid, mat._grid); + + assert(V.checkerboard == Even); + assert(U.checkerboard == Odd); + mat.checkerboard = Odd; + + DerivInternal(StencilEven, UmuOdd, UUUmuOdd, mat, U, V, dag); +} + +template +void ImprovedStaggeredFermion::DhopDerivEO(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) { + + conformable(U._grid, _cbgrid); + conformable(U._grid, V._grid); + conformable(U._grid, mat._grid); + + assert(V.checkerboard == Odd); + assert(U.checkerboard == Even); + mat.checkerboard = Even; + + DerivInternal(StencilOdd, UmuEven, UUUmuEven, mat, U, V, dag); +} + +template +void ImprovedStaggeredFermion::Dhop(const FermionField &in, FermionField &out, int dag) { + conformable(in._grid, _grid); // verifies full grid + conformable(in._grid, out._grid); + + out.checkerboard = in.checkerboard; + + DhopInternal(Stencil, Lebesgue, Umu, UUUmu, in, out, dag); +} + +template +void ImprovedStaggeredFermion::DhopOE(const FermionField &in, FermionField &out, int dag) { + conformable(in._grid, _cbgrid); // verifies half grid + conformable(in._grid, out._grid); // drops the cb check + + assert(in.checkerboard == Even); + out.checkerboard = Odd; + + DhopInternal(StencilEven, LebesgueEvenOdd, UmuOdd, UUUmuOdd, in, out, dag); +} + +template +void ImprovedStaggeredFermion::DhopEO(const FermionField &in, FermionField &out, int dag) { + conformable(in._grid, _cbgrid); // verifies half grid + conformable(in._grid, out._grid); // drops the cb check + + assert(in.checkerboard == Odd); + out.checkerboard = Even; + + DhopInternal(StencilOdd, LebesgueEvenOdd, UmuEven, UUUmuEven, in, out, dag); +} + +template +void ImprovedStaggeredFermion::Mdir(const FermionField &in, FermionField &out, int dir, int disp) { + DhopDir(in, out, dir, disp); +} + +template +void ImprovedStaggeredFermion::DhopDir(const FermionField &in, FermionField &out, int dir, int disp) { + + Compressor compressor; + Stencil.HaloExchange(in, compressor); + +PARALLEL_FOR_LOOP + for (int sss = 0; sss < in._grid->oSites(); sss++) { + Kernels::DhopDir(Stencil, Umu, UUUmu, Stencil.CommBuf(), sss, sss, in, out, dir, disp); + } +}; + +template +void ImprovedStaggeredFermion::DhopInternal(StencilImpl &st, LebesgueOrder &lo, + DoubledGaugeField &U, + DoubledGaugeField &UUU, + const FermionField &in, + FermionField &out, int dag) { + assert((dag == DaggerNo) || (dag == DaggerYes)); + + Compressor compressor; + st.HaloExchange(in, compressor); + + if (dag == DaggerYes) { + PARALLEL_FOR_LOOP + for (int sss = 0; sss < in._grid->oSites(); sss++) { + Kernels::DhopSiteDag(st, lo, U, UUU, st.CommBuf(), sss, sss, in, out); + } + } else { + PARALLEL_FOR_LOOP + for (int sss = 0; sss < in._grid->oSites(); sss++) { + Kernels::DhopSite(st, lo, U, UUU, st.CommBuf(), sss, sss, in, out); + } + } +}; + +FermOpStaggeredTemplateInstantiate(ImprovedStaggeredFermion); + + //AdjointFermOpTemplateInstantiate(ImprovedStaggeredFermion); + //TwoIndexFermOpTemplateInstantiate(ImprovedStaggeredFermion); + +}} diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion.h b/lib/qcd/action/fermion/ImprovedStaggeredFermion.h new file mode 100644 index 00000000..bf05240e --- /dev/null +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion.h @@ -0,0 +1,152 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/action/fermion/ImprovedStaggered.h + +Copyright (C) 2015 + +Author: Azusa Yamaguchi, 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 GRID_QCD_IMPR_STAG_FERMION_H +#define GRID_QCD_IMPR_STAG_FERMION_H + +namespace Grid { + +namespace QCD { + +class ImprovedStaggeredFermionStatic { + public: + static const std::vector directions; + static const std::vector displacements; + static const int npoint = 16; +}; + +template +class ImprovedStaggeredFermion : public StaggeredKernels, public ImprovedStaggeredFermionStatic { + public: + INHERIT_IMPL_TYPES(Impl); + typedef StaggeredKernels Kernels; + + /////////////////////////////////////////////////////////////// + // Implement the abstract base + /////////////////////////////////////////////////////////////// + GridBase *GaugeGrid(void) { return _grid; } + GridBase *GaugeRedBlackGrid(void) { return _cbgrid; } + GridBase *FermionGrid(void) { return _grid; } + GridBase *FermionRedBlackGrid(void) { return _cbgrid; } + + ////////////////////////////////////////////////////////////////// + // override multiply; cut number routines if pass dagger argument + // and also make interface more uniformly consistent + ////////////////////////////////////////////////////////////////// + RealD M(const FermionField &in, FermionField &out); + RealD Mdag(const FermionField &in, FermionField &out); + + ///////////////////////////////////////////////////////// + // half checkerboard operations + ///////////////////////////////////////////////////////// + void Meooe(const FermionField &in, FermionField &out); + void MeooeDag(const FermionField &in, FermionField &out); + + // allow override for twisted mass and clover + virtual void Mooee(const FermionField &in, FermionField &out); + virtual void MooeeDag(const FermionField &in, FermionField &out); + virtual void MooeeInv(const FermionField &in, FermionField &out); + virtual void MooeeInvDag(const FermionField &in, FermionField &out); + + //////////////////////// + // Derivative interface + //////////////////////// + // Interface calls an internal routine + void DhopDeriv (GaugeField &mat, const FermionField &U, const FermionField &V, int dag); + void DhopDerivOE(GaugeField &mat, const FermionField &U, const FermionField &V, int dag); + void DhopDerivEO(GaugeField &mat, const FermionField &U, const FermionField &V, int dag); + + /////////////////////////////////////////////////////////////// + // non-hermitian hopping term; half cb or both + /////////////////////////////////////////////////////////////// + void Dhop (const FermionField &in, FermionField &out, int dag); + void DhopOE(const FermionField &in, FermionField &out, int dag); + void DhopEO(const FermionField &in, FermionField &out, int dag); + + /////////////////////////////////////////////////////////////// + // Multigrid assistance; force term uses too + /////////////////////////////////////////////////////////////// + void Mdir(const FermionField &in, FermionField &out, int dir, int disp); + void DhopDir(const FermionField &in, FermionField &out, int dir, int disp); + + /////////////////////////////////////////////////////////////// + // Extra methods added by derived + /////////////////////////////////////////////////////////////// + void DerivInternal(StencilImpl &st, + DoubledGaugeField &U,DoubledGaugeField &UUU, + GaugeField &mat, + const FermionField &A, const FermionField &B, int dag); + + void DhopInternal(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,DoubledGaugeField &UUU, + const FermionField &in, FermionField &out, int dag); + + // Constructor + ImprovedStaggeredFermion(GaugeField &_Umu, GridCartesian &Fgrid, + GridRedBlackCartesian &Hgrid, RealD _mass, + const ImplParams &p = ImplParams()); + + // DoubleStore impl dependent + void ImportGauge(const GaugeField &_Umu); + + /////////////////////////////////////////////////////////////// + // Data members require to support the functionality + /////////////////////////////////////////////////////////////// + + // protected: + public: + // any other parameters of action ??? + + RealD mass; + + GridBase *_grid; + GridBase *_cbgrid; + + // Defines the stencils for even and odd + StencilImpl Stencil; + StencilImpl StencilEven; + StencilImpl StencilOdd; + + // Copy of the gauge field , with even and odd subsets + DoubledGaugeField Umu; + DoubledGaugeField UmuEven; + DoubledGaugeField UmuOdd; + + DoubledGaugeField UUUmu; + DoubledGaugeField UUUmuEven; + DoubledGaugeField UUUmuOdd; + + LebesgueOrder Lebesgue; + LebesgueOrder LebesgueEvenOdd; +}; + +typedef ImprovedStaggeredFermion ImprovedStaggeredFermionF; +typedef ImprovedStaggeredFermion ImprovedStaggeredFermionD; + +} +} +#endif diff --git a/lib/qcd/action/fermion/StaggeredKernels.cc b/lib/qcd/action/fermion/StaggeredKernels.cc new file mode 100644 index 00000000..8df7f3e4 --- /dev/null +++ b/lib/qcd/action/fermion/StaggeredKernels.cc @@ -0,0 +1,223 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/action/fermion/WilsonKernels.cc + +Copyright (C) 2015 + +Author: Azusa Yamaguchi, 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 */ +#include +namespace Grid { +namespace QCD { + +template +StaggeredKernels::StaggeredKernels(const ImplParams &p) : Base(p){}; + +//////////////////////////////////////////// +// Generic implementation; move to different file? +//////////////////////////////////////////// + +template +void StaggeredKernels::DhopSiteDepth(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, + SiteSpinor *buf, int sF, + int sU, const FermionField &in, SiteSpinor &out,int threeLink) { + const SiteSpinor *chi_p; + SiteSpinor chi; + SiteSpinor Uchi; + StencilEntry *SE; + int ptype; + int skew = 0; + if (threeLink) skew=8; + /////////////////////////// + // Xp + /////////////////////////// + + SE = st.GetEntry(ptype, Xp+skew, sF); + if (SE->_is_local) { + if (SE->_permute) { + chi_p = χ + permute(chi, in._odata[SE->_offset], ptype); + } else { + chi_p = &in._odata[SE->_offset]; + } + } else { + chi_p = &buf[SE->_offset]; + } + Impl::multLink(Uchi, U._odata[sU], *chi_p, Xp); + + /////////////////////////// + // Yp + /////////////////////////// + SE = st.GetEntry(ptype, Yp+skew, sF); + if (SE->_is_local) { + if (SE->_permute) { + chi_p = χ + permute(chi, in._odata[SE->_offset], ptype); + } else { + chi_p = &in._odata[SE->_offset]; + } + } else { + chi_p = &buf[SE->_offset]; + } + Impl::multLinkAdd(Uchi, U._odata[sU], *chi_p, Yp); + + /////////////////////////// + // Zp + /////////////////////////// + SE = st.GetEntry(ptype, Zp+skew, sF); + if (SE->_is_local) { + if (SE->_permute) { + chi_p = χ + permute(chi, in._odata[SE->_offset], ptype); + } else { + chi_p = &in._odata[SE->_offset]; + } + } else { + chi_p = &buf[SE->_offset]; + } + Impl::multLinkAdd(Uchi, U._odata[sU], *chi_p, Zp); + + /////////////////////////// + // Tp + /////////////////////////// + SE = st.GetEntry(ptype, Tp+skew, sF); + if (SE->_is_local) { + if (SE->_permute) { + chi_p = χ + permute(chi, in._odata[SE->_offset], ptype); + } else { + chi_p = &in._odata[SE->_offset]; + } + } else { + chi_p = &buf[SE->_offset]; + } + Impl::multLinkAdd(Uchi, U._odata[sU], *chi_p, Tp); + + /////////////////////////// + // Xm + /////////////////////////// + SE = st.GetEntry(ptype, Xm+skew, sF); + if (SE->_is_local) { + if (SE->_permute) { + chi_p = χ + permute(chi, in._odata[SE->_offset], ptype); + } else { + chi_p = &in._odata[SE->_offset]; + } + } else { + chi_p = &buf[SE->_offset]; + } + Impl::multLinkAdd(Uchi, U._odata[sU], *chi_p, Xm); + + /////////////////////////// + // Ym + /////////////////////////// + SE = st.GetEntry(ptype, Ym+skew, sF); + if (SE->_is_local) { + if (SE->_permute) { + chi_p = χ + permute(chi, in._odata[SE->_offset], ptype); + } else { + chi_p = &in._odata[SE->_offset]; + } + } else { + chi_p = &buf[SE->_offset]; + } + Impl::multLinkAdd(Uchi, U._odata[sU], *chi_p, Ym); + + /////////////////////////// + // Zm + /////////////////////////// + SE = st.GetEntry(ptype, Zm+skew, sF); + if (SE->_is_local) { + if (SE->_permute) { + chi_p = χ + permute(chi, in._odata[SE->_offset], ptype); + } else { + chi_p = &in._odata[SE->_offset]; + } + } else { + chi_p = &buf[SE->_offset]; + } + Impl::multLinkAdd(Uchi, U._odata[sU], *chi_p, Zm); + + /////////////////////////// + // Tm + /////////////////////////// + SE = st.GetEntry(ptype, Tm+skew, sF); + if (SE->_is_local) { + if (SE->_permute) { + chi_p = χ + permute(chi, in._odata[SE->_offset], ptype); + } else { + chi_p = &in._odata[SE->_offset]; + } + } else { + chi_p = &buf[SE->_offset]; + } + Impl::multLinkAdd(Uchi, U._odata[sU], *chi_p, Tm); + + vstream(out, Uchi); +}; + +// Need controls to do interior, exterior, or both +template +void StaggeredKernels::DhopSiteDag(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, DoubledGaugeField &UUU, + SiteSpinor *buf, int sF, + int sU, const FermionField &in, FermionField &out) { + SiteSpinor naik; + SiteSpinor naive; + int oneLink =0; + int threeLink=1; + DhopSiteDepth(st,lo,U,buf,sF,sU,in,naive,oneLink); + DhopSiteDepth(st,lo,UUU,buf,sF,sU,in,naik,threeLink); + out._odata[sF] =naive+naik; +}; +template +void StaggeredKernels::DhopSite(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, DoubledGaugeField &UUU, + SiteSpinor *buf, int sF, + int sU, const FermionField &in, FermionField &out) { + SiteSpinor naik; + SiteSpinor naive; + int oneLink =0; + int threeLink=1; + DhopSiteDepth(st,lo,U,buf,sF,sU,in,naive,oneLink); + DhopSiteDepth(st,lo,UUU,buf,sF,sU,in,naik,threeLink); + out._odata[sF] =-naive-naik; +}; + +template +void StaggeredKernels::DhopDir( StencilImpl &st, DoubledGaugeField &U, DoubledGaugeField &UUU, SiteSpinor *buf, int sF, + int sU, const FermionField &in, FermionField &out, int dir, int disp) +{ + // Disp should be either +1,-1,+3,-3 + // What about "dag" ? + // Because we work out pU . dS/dU + // U + assert(0); +} + +FermOpStaggeredTemplateInstantiate(StaggeredKernels); + +}} + diff --git a/lib/qcd/action/fermion/StaggeredKernels.h b/lib/qcd/action/fermion/StaggeredKernels.h new file mode 100644 index 00000000..f51a4b37 --- /dev/null +++ b/lib/qcd/action/fermion/StaggeredKernels.h @@ -0,0 +1,70 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/action/fermion/StaggeredKernels.h + +Copyright (C) 2015 + +Author: Azusa Yamaguchi, 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 GRID_QCD_STAGGERED_KERNELS_H +#define GRID_QCD_STAGGERED_KERNELS_H + +namespace Grid { +namespace QCD { + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // Helper routines that implement Staggered stencil for a single site. + //////////////////////////////////////////////////////////////////////////////////////////////////////////////// +class StaggeredKernelsStatic { + public: +}; + +template class StaggeredKernels : public FermionOperator , public StaggeredKernelsStatic { + public: + + INHERIT_IMPL_TYPES(Impl); + typedef FermionOperator Base; + +public: + + void DhopDir(StencilImpl &st, DoubledGaugeField &U, DoubledGaugeField &UUU, SiteSpinor * buf, + int sF, int sU, const FermionField &in, FermionField &out, int dir,int disp); + + void DhopSiteDepth(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, SiteSpinor * buf, + int sF, int sU, const FermionField &in, SiteSpinor &out,int threeLink); + + void DhopSite(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, DoubledGaugeField &UUU, SiteSpinor * buf, + int sF, int sU, const FermionField &in, FermionField &out); + + void DhopSiteDag(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, DoubledGaugeField &UUU, SiteSpinor * buf, + int sF, int sU, const FermionField &in, FermionField &out); + +public: + + StaggeredKernels(const ImplParams &p = ImplParams()); + +}; + +}} + +#endif diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc index 4bc28bc7..845250fc 100644 --- a/lib/qcd/action/fermion/WilsonFermion.cc +++ b/lib/qcd/action/fermion/WilsonFermion.cc @@ -34,10 +34,9 @@ directory namespace Grid { namespace QCD { -const std::vector WilsonFermionStatic::directions({0, 1, 2, 3, 0, 1, 2, - 3}); -const std::vector WilsonFermionStatic::displacements({1, 1, 1, 1, -1, -1, - -1, -1}); +const std::vector WilsonFermionStatic::directions({0, 1, 2, 3, 0, 1, 2, 3}); +const std::vector WilsonFermionStatic::displacements({1, 1, 1, 1, -1, -1, -1, -1}); + int WilsonFermionStatic::HandOptDslash; ///////////////////////////////// @@ -166,8 +165,7 @@ void WilsonFermion::DerivInternal(StencilImpl &st, DoubledGaugeField &U, //////////////////////// PARALLEL_FOR_LOOP for (int sss = 0; sss < B._grid->oSites(); sss++) { - Kernels::DiracOptDhopDir(st, U, st.CommBuf(), sss, sss, B, Btilde, mu, - gamma); + Kernels::DiracOptDhopDir(st, U, st.CommBuf(), sss, sss, B, Btilde, mu, gamma); } ////////////////////////////////////////////////// @@ -277,8 +275,7 @@ void WilsonFermion::DhopDirDisp(const FermionField &in, FermionField &out, PARALLEL_FOR_LOOP for (int sss = 0; sss < in._grid->oSites(); sss++) { - Kernels::DiracOptDhopDir(Stencil, Umu, Stencil.CommBuf(), sss, sss, in, out, - dirdisp, gamma); + Kernels::DiracOptDhopDir(Stencil, Umu, Stencil.CommBuf(), sss, sss, in, out, dirdisp, gamma); } }; @@ -295,14 +292,12 @@ void WilsonFermion::DhopInternal(StencilImpl &st, LebesgueOrder &lo, if (dag == DaggerYes) { PARALLEL_FOR_LOOP for (int sss = 0; sss < in._grid->oSites(); sss++) { - Kernels::DiracOptDhopSiteDag(st, lo, U, st.CommBuf(), sss, sss, 1, 1, in, - out); + Kernels::DiracOptDhopSiteDag(st, lo, U, st.CommBuf(), sss, sss, 1, 1, in, out); } } else { PARALLEL_FOR_LOOP for (int sss = 0; sss < in._grid->oSites(); sss++) { - Kernels::DiracOptDhopSite(st, lo, U, st.CommBuf(), sss, sss, 1, 1, in, - out); + Kernels::DiracOptDhopSite(st, lo, U, st.CommBuf(), sss, sss, 1, 1, in, out); } } }; diff --git a/lib/qcd/action/pseudofermion/TwoFlavour.h b/lib/qcd/action/pseudofermion/TwoFlavour.h index 6b65a95d..ddc17d42 100644 --- a/lib/qcd/action/pseudofermion/TwoFlavour.h +++ b/lib/qcd/action/pseudofermion/TwoFlavour.h @@ -63,8 +63,7 @@ class TwoFlavourPseudoFermionAction : public Action { Phi(Op.FermionGrid()){}; ////////////////////////////////////////////////////////////////////////////////////// - // Push the gauge field in to the dops. Assume any BC's and smearing already - // applied + // Push the gauge field in to the dops. Assume any BC's and smearing already applied ////////////////////////////////////////////////////////////////////////////////////// virtual void refresh(const GaugeField &U, GridParallelRNG &pRNG) { // P(phi) = e^{- phi^dag (MdagM)^-1 phi} @@ -107,8 +106,7 @@ class TwoFlavourPseudoFermionAction : public Action { MdagMOp.Op(X, Y); RealD action = norm2(Y); - std::cout << GridLogMessage << "Pseudofermion action " << action - << std::endl; + std::cout << GridLogMessage << "Pseudofermion action " << action << std::endl; return action; }; @@ -119,6 +117,7 @@ class TwoFlavourPseudoFermionAction : public Action { // // = - Ydag dM X - Xdag dMdag Y // + // ////////////////////////////////////////////////////// virtual void deriv(const GaugeField &U, GaugeField &dSdU) { FermOp.ImportGauge(U); @@ -133,8 +132,7 @@ class TwoFlavourPseudoFermionAction : public Action { DerivativeSolver(MdagMOp, Phi, X); // X = (MdagM)^-1 phi MdagMOp.Op(X, Y); // Y = M X = (Mdag)^-1 phi - // Our conventions really make this UdSdU; We do not differentiate wrt Udag - // here. + // Our conventions really make this UdSdU; We do not differentiate wrt Udag here. // So must take dSdU - adj(dSdU) and left multiply by mom to get dS/dt. FermOp.MDeriv(tmp, Y, X, DaggerNo);