diff --git a/lib/qcd/action/fermion/DomainWallEOFAFermionssp.cc b/lib/qcd/action/fermion/DomainWallEOFAFermionssp.cc index 80a4bf09..f49eb2aa 100644 --- a/lib/qcd/action/fermion/DomainWallEOFAFermionssp.cc +++ b/lib/qcd/action/fermion/DomainWallEOFAFermionssp.cc @@ -28,141 +28,140 @@ with this program; if not, write to the Free Software Foundation, Inc., See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ -/* END LEGAL */ + /* END LEGAL */ #include #include -namespace Grid { -namespace QCD { +NAMESPACE_BEGIN(Grid); - // FIXME -- make a version of these routines with site loop outermost for cache reuse. - // Pminus fowards - // Pplus backwards - template - void DomainWallEOFAFermion::M5D(const FermionField& psi, const FermionField& phi, - FermionField& chi, std::vector& lower, std::vector& diag, std::vector& upper) - { - Coeff_t one(1.0); - int Ls = this->Ls; - for(int s=0; s +void DomainWallEOFAFermion::M5D(const FermionField& psi, const FermionField& phi, + FermionField& chi, std::vector& lower, std::vector& diag, std::vector& upper) +{ + Coeff_t one(1.0); + int Ls = this->Ls; + for(int s=0; s - void DomainWallEOFAFermion::M5Ddag(const FermionField& psi, const FermionField& phi, - FermionField& chi, std::vector& lower, std::vector& diag, std::vector& upper) - { - Coeff_t one(1.0); - int Ls = this->Ls; - for(int s=0; s +void DomainWallEOFAFermion::M5Ddag(const FermionField& psi, const FermionField& phi, + FermionField& chi, std::vector& lower, std::vector& diag, std::vector& upper) +{ + Coeff_t one(1.0); + int Ls = this->Ls; + for(int s=0; s - void DomainWallEOFAFermion::MooeeInv(const FermionField& psi, FermionField& chi) - { - Coeff_t one(1.0); - Coeff_t czero(0.0); - chi.checkerboard = psi.checkerboard; - int Ls = this->Ls; +template +void DomainWallEOFAFermion::MooeeInv(const FermionField& psi, FermionField& chi) +{ + Coeff_t one(1.0); + Coeff_t czero(0.0); + chi.checkerboard = psi.checkerboard; + int Ls = this->Ls; - FermionField tmp(psi._grid); + FermionField tmp(psi._grid); - // Apply (L^{\prime})^{-1} - axpby_ssp(chi, one, psi, czero, psi, 0, 0); // chi[0]=psi[0] - for(int s=1; slee[s-1], chi, s, s-1);// recursion Psi[s] -lee P_+ chi[s-1] - } + // Apply (L^{\prime})^{-1} + axpby_ssp(chi, one, psi, czero, psi, 0, 0); // chi[0]=psi[0] + for(int s=1; slee[s-1], chi, s, s-1);// recursion Psi[s] -lee P_+ chi[s-1] + } - // L_m^{-1} - for(int s=0; sleem[s], chi, Ls-1, s); - } + // L_m^{-1} + for(int s=0; sleem[s], chi, Ls-1, s); + } - // U_m^{-1} D^{-1} - for(int s=0; sdee[s], chi, -this->ueem[s]/this->dee[Ls], chi, s, Ls-1); - } - axpby_ssp_pminus(tmp, czero, chi, one/this->dee[Ls-1], chi, Ls-1, Ls-1); - axpby_ssp_pplus(chi, one, tmp, one/this->dee[Ls], chi, Ls-1, Ls-1); + // U_m^{-1} D^{-1} + for(int s=0; sdee[s], chi, -this->ueem[s]/this->dee[Ls], chi, s, Ls-1); + } + axpby_ssp_pminus(tmp, czero, chi, one/this->dee[Ls-1], chi, Ls-1, Ls-1); + axpby_ssp_pplus(chi, one, tmp, one/this->dee[Ls], chi, Ls-1, Ls-1); - // Apply U^{-1} - for(int s=Ls-2; s>=0; s--){ - axpby_ssp_pminus(chi, one, chi, -this->uee[s], chi, s, s+1); // chi[Ls] - } - } + // Apply U^{-1} + for(int s=Ls-2; s>=0; s--){ + axpby_ssp_pminus(chi, one, chi, -this->uee[s], chi, s, s+1); // chi[Ls] + } +} - template - void DomainWallEOFAFermion::MooeeInvDag(const FermionField& psi, FermionField& chi) - { - Coeff_t one(1.0); - Coeff_t czero(0.0); - chi.checkerboard = psi.checkerboard; - int Ls = this->Ls; +template +void DomainWallEOFAFermion::MooeeInvDag(const FermionField& psi, FermionField& chi) +{ + Coeff_t one(1.0); + Coeff_t czero(0.0); + chi.checkerboard = psi.checkerboard; + int Ls = this->Ls; - FermionField tmp(psi._grid); + FermionField tmp(psi._grid); - // Apply (U^{\prime})^{-dagger} - axpby_ssp(chi, one, psi, czero, psi, 0, 0); // chi[0]=psi[0] - for(int s=1; suee[s-1]), chi, s, s-1); - } + // Apply (U^{\prime})^{-dagger} + axpby_ssp(chi, one, psi, czero, psi, 0, 0); // chi[0]=psi[0] + for(int s=1; suee[s-1]), chi, s, s-1); + } - // U_m^{-\dagger} - for(int s=0; sueem[s]), chi, Ls-1, s); - } + // U_m^{-\dagger} + for(int s=0; sueem[s]), chi, Ls-1, s); + } - // L_m^{-\dagger} D^{-dagger} - for(int s=0; sdee[s]), chi, -conjugate(this->leem[s]/this->dee[Ls-1]), chi, s, Ls-1); - } - axpby_ssp_pminus(tmp, czero, chi, one/conjugate(this->dee[Ls-1]), chi, Ls-1, Ls-1); - axpby_ssp_pplus(chi, one, tmp, one/conjugate(this->dee[Ls]), chi, Ls-1, Ls-1); + // L_m^{-\dagger} D^{-dagger} + for(int s=0; sdee[s]), chi, -conjugate(this->leem[s]/this->dee[Ls-1]), chi, s, Ls-1); + } + axpby_ssp_pminus(tmp, czero, chi, one/conjugate(this->dee[Ls-1]), chi, Ls-1, Ls-1); + axpby_ssp_pplus(chi, one, tmp, one/conjugate(this->dee[Ls]), chi, Ls-1, Ls-1); - // Apply L^{-dagger} - for(int s=Ls-2; s>=0; s--){ - axpby_ssp_pplus(chi, one, chi, -conjugate(this->lee[s]), chi, s, s+1); // chi[Ls] - } - } + // Apply L^{-dagger} + for(int s=Ls-2; s>=0; s--){ + axpby_ssp_pplus(chi, one, chi, -conjugate(this->lee[s]), chi, s, s+1); // chi[Ls] + } +} - #ifdef DOMAIN_WALL_EOFA_DPERP_LINALG +#ifdef DOMAIN_WALL_EOFA_DPERP_LINALG - INSTANTIATE_DPERP_DWF_EOFA(WilsonImplF); - INSTANTIATE_DPERP_DWF_EOFA(WilsonImplD); - INSTANTIATE_DPERP_DWF_EOFA(GparityWilsonImplF); - INSTANTIATE_DPERP_DWF_EOFA(GparityWilsonImplD); - INSTANTIATE_DPERP_DWF_EOFA(ZWilsonImplF); - INSTANTIATE_DPERP_DWF_EOFA(ZWilsonImplD); +INSTANTIATE_DPERP_DWF_EOFA(WilsonImplF); +INSTANTIATE_DPERP_DWF_EOFA(WilsonImplD); +INSTANTIATE_DPERP_DWF_EOFA(GparityWilsonImplF); +INSTANTIATE_DPERP_DWF_EOFA(GparityWilsonImplD); +INSTANTIATE_DPERP_DWF_EOFA(ZWilsonImplF); +INSTANTIATE_DPERP_DWF_EOFA(ZWilsonImplD); - INSTANTIATE_DPERP_DWF_EOFA(WilsonImplFH); - INSTANTIATE_DPERP_DWF_EOFA(WilsonImplDF); - INSTANTIATE_DPERP_DWF_EOFA(GparityWilsonImplFH); - INSTANTIATE_DPERP_DWF_EOFA(GparityWilsonImplDF); - INSTANTIATE_DPERP_DWF_EOFA(ZWilsonImplFH); - INSTANTIATE_DPERP_DWF_EOFA(ZWilsonImplDF); +INSTANTIATE_DPERP_DWF_EOFA(WilsonImplFH); +INSTANTIATE_DPERP_DWF_EOFA(WilsonImplDF); +INSTANTIATE_DPERP_DWF_EOFA(GparityWilsonImplFH); +INSTANTIATE_DPERP_DWF_EOFA(GparityWilsonImplDF); +INSTANTIATE_DPERP_DWF_EOFA(ZWilsonImplFH); +INSTANTIATE_DPERP_DWF_EOFA(ZWilsonImplDF); - #endif +#endif -}} +NAMESPACE_END(Grid);