diff --git a/lib/qcd/action/fermion/DomainWallEOFAFermiondense.cc b/lib/qcd/action/fermion/DomainWallEOFAFermiondense.cc index c27074d9..5c269c86 100644 --- a/lib/qcd/action/fermion/DomainWallEOFAFermiondense.cc +++ b/lib/qcd/action/fermion/DomainWallEOFAFermiondense.cc @@ -28,132 +28,131 @@ 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 #include -namespace Grid { -namespace QCD { +NAMESPACE_BEGIN(Grid); - /* - * Dense matrix versions of routines - */ - template - void DomainWallEOFAFermion::MooeeInvDag(const FermionField& psi, FermionField& chi) - { - this->MooeeInternal(psi, chi, DaggerYes, InverseYes); +/* + * Dense matrix versions of routines + */ +template +void DomainWallEOFAFermion::MooeeInvDag(const FermionField& psi, FermionField& chi) +{ + this->MooeeInternal(psi, chi, DaggerYes, InverseYes); +} + +template +void DomainWallEOFAFermion::MooeeInv(const FermionField& psi, FermionField& chi) +{ + this->MooeeInternal(psi, chi, DaggerNo, InverseYes); +} + +template +void DomainWallEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv) +{ + int Ls = this->Ls; + int LLs = psi._grid->_rdimensions[0]; + int vol = psi._grid->oSites()/LLs; + + chi.checkerboard = psi.checkerboard; + + assert(Ls==LLs); + + Eigen::MatrixXd Pplus = Eigen::MatrixXd::Zero(Ls,Ls); + Eigen::MatrixXd Pminus = Eigen::MatrixXd::Zero(Ls,Ls); + + for(int s=0;sbee[s]; + Pminus(s,s) = this->bee[s]; + } + + for(int s=0; scee[s]; + } + + for(int s=0; scee[s+1]; + } + + Pplus (0,Ls-1) = this->dp; + Pminus(Ls-1,0) = this->dm; + + Eigen::MatrixXd PplusMat ; + Eigen::MatrixXd PminusMat; + + if(inv) { + PplusMat = Pplus.inverse(); + PminusMat = Pminus.inverse(); + } else { + PplusMat = Pplus; + PminusMat = Pminus; + } + + if(dag){ + PplusMat.adjointInPlace(); + PminusMat.adjointInPlace(); + } + + // For the non-vectorised s-direction this is simple + + for(auto site=0; site - void DomainWallEOFAFermion::MooeeInv(const FermionField& psi, FermionField& chi) - { - this->MooeeInternal(psi, chi, DaggerNo, InverseYes); - } +#ifdef DOMAIN_WALL_EOFA_DPERP_DENSE - template - void DomainWallEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv) - { - int Ls = this->Ls; - int LLs = psi._grid->_rdimensions[0]; - int vol = psi._grid->oSites()/LLs; +INSTANTIATE_DPERP_DWF_EOFA(GparityWilsonImplF); +INSTANTIATE_DPERP_DWF_EOFA(GparityWilsonImplD); +INSTANTIATE_DPERP_DWF_EOFA(WilsonImplF); +INSTANTIATE_DPERP_DWF_EOFA(WilsonImplD); +INSTANTIATE_DPERP_DWF_EOFA(ZWilsonImplF); +INSTANTIATE_DPERP_DWF_EOFA(ZWilsonImplD); - chi.checkerboard = psi.checkerboard; +template void DomainWallEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); +template void DomainWallEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); +template void DomainWallEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); +template void DomainWallEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); +template void DomainWallEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); +template void DomainWallEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); - assert(Ls==LLs); +INSTANTIATE_DPERP_DWF_EOFA(GparityWilsonImplFH); +INSTANTIATE_DPERP_DWF_EOFA(GparityWilsonImplDF); +INSTANTIATE_DPERP_DWF_EOFA(WilsonImplFH); +INSTANTIATE_DPERP_DWF_EOFA(WilsonImplDF); +INSTANTIATE_DPERP_DWF_EOFA(ZWilsonImplFH); +INSTANTIATE_DPERP_DWF_EOFA(ZWilsonImplDF); - Eigen::MatrixXd Pplus = Eigen::MatrixXd::Zero(Ls,Ls); - Eigen::MatrixXd Pminus = Eigen::MatrixXd::Zero(Ls,Ls); +template void DomainWallEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); +template void DomainWallEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); +template void DomainWallEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); +template void DomainWallEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); +template void DomainWallEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); +template void DomainWallEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); - for(int s=0;sbee[s]; - Pminus(s,s) = this->bee[s]; - } +#endif - for(int s=0; scee[s]; - } - - for(int s=0; scee[s+1]; - } - - Pplus (0,Ls-1) = this->dp; - Pminus(Ls-1,0) = this->dm; - - Eigen::MatrixXd PplusMat ; - Eigen::MatrixXd PminusMat; - - if(inv) { - PplusMat = Pplus.inverse(); - PminusMat = Pminus.inverse(); - } else { - PplusMat = Pplus; - PminusMat = Pminus; - } - - if(dag){ - PplusMat.adjointInPlace(); - PminusMat.adjointInPlace(); - } - - // For the non-vectorised s-direction this is simple - - for(auto site=0; site::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); - template void DomainWallEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); - template void DomainWallEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); - template void DomainWallEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); - template void DomainWallEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); - template void DomainWallEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); - - INSTANTIATE_DPERP_DWF_EOFA(GparityWilsonImplFH); - INSTANTIATE_DPERP_DWF_EOFA(GparityWilsonImplDF); - INSTANTIATE_DPERP_DWF_EOFA(WilsonImplFH); - INSTANTIATE_DPERP_DWF_EOFA(WilsonImplDF); - INSTANTIATE_DPERP_DWF_EOFA(ZWilsonImplFH); - INSTANTIATE_DPERP_DWF_EOFA(ZWilsonImplDF); - - template void DomainWallEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); - template void DomainWallEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); - template void DomainWallEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); - template void DomainWallEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); - template void DomainWallEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); - template void DomainWallEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); - - #endif - -}} +NAMESPACE_END(Grid);