diff --git a/lib/qcd/action/fermion/MobiusEOFAFermiondense.cc b/lib/qcd/action/fermion/MobiusEOFAFermiondense.cc index d66b8cd9..bbf06585 100644 --- a/lib/qcd/action/fermion/MobiusEOFAFermiondense.cc +++ b/lib/qcd/action/fermion/MobiusEOFAFermiondense.cc @@ -28,157 +28,156 @@ 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 MobiusEOFAFermion::MooeeInv(const FermionField& psi, FermionField& chi) - { - this->MooeeInternal(psi, chi, DaggerNo, InverseYes); +/* + * Dense matrix versions of routines + */ +template +void MobiusEOFAFermion::MooeeInv(const FermionField& psi, FermionField& chi) +{ + this->MooeeInternal(psi, chi, DaggerNo, InverseYes); +} + +template +void MobiusEOFAFermion::MooeeInv_shift(const FermionField& psi, FermionField& chi) +{ + this->MooeeInternal(psi, chi, DaggerNo, InverseYes); +} + +template +void MobiusEOFAFermion::MooeeInvDag(const FermionField& psi, FermionField& chi) +{ + this->MooeeInternal(psi, chi, DaggerYes, InverseYes); +} + +template +void MobiusEOFAFermion::MooeeInvDag_shift(const FermionField& psi, FermionField& chi) +{ + this->MooeeInternal(psi, chi, DaggerYes, InverseYes); +} + +template +void MobiusEOFAFermion::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; + + int pm = this->pm; + RealD shift = this->shift; + RealD alpha = this->alpha; + RealD k = this->k; + RealD mq1 = this->mq1; + + 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]; } - template - void MobiusEOFAFermion::MooeeInv_shift(const FermionField& psi, FermionField& chi) - { - this->MooeeInternal(psi, chi, DaggerNo, InverseYes); + for(int s=0; scee[s]; } - template - void MobiusEOFAFermion::MooeeInvDag(const FermionField& psi, FermionField& chi) - { - this->MooeeInternal(psi, chi, DaggerYes, InverseYes); + for(int s=0; scee[s+1]; } + Pplus (0,Ls-1) = mq1*this->cee[0]; + Pminus(Ls-1,0) = mq1*this->cee[Ls-1]; - template - void MobiusEOFAFermion::MooeeInvDag_shift(const FermionField& psi, FermionField& chi) - { - this->MooeeInternal(psi, chi, DaggerYes, InverseYes); - } - - template - void MobiusEOFAFermion::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; - - int pm = this->pm; - RealD shift = this->shift; - RealD alpha = this->alpha; - RealD k = this->k; - RealD mq1 = this->mq1; - - 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]; + if(shift != 0.0){ + Coeff_t N = 2.0 * ( std::pow(alpha+1.0,Ls) + mq1*std::pow(alpha-1.0,Ls) ); + for(int s=0; scee[s]; - } + Eigen::MatrixXd PplusMat ; + Eigen::MatrixXd PminusMat; - for(int s=0; scee[s+1]; - } - Pplus (0,Ls-1) = mq1*this->cee[0]; - Pminus(Ls-1,0) = mq1*this->cee[Ls-1]; + if(inv){ + PplusMat = Pplus.inverse(); + PminusMat = Pminus.inverse(); + } else { + PplusMat = Pplus; + PminusMat = Pminus; + } - if(shift != 0.0){ - Coeff_t N = 2.0 * ( std::pow(alpha+1.0,Ls) + mq1*std::pow(alpha-1.0,Ls) ); - for(int s=0; s::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); - template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); - template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); - template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); - template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); - template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); +template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); +template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); +template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); +template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); +template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); +template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); - INSTANTIATE_DPERP_MOBIUS_EOFA(GparityWilsonImplFH); - INSTANTIATE_DPERP_MOBIUS_EOFA(GparityWilsonImplDF); - INSTANTIATE_DPERP_MOBIUS_EOFA(WilsonImplFH); - INSTANTIATE_DPERP_MOBIUS_EOFA(WilsonImplDF); - INSTANTIATE_DPERP_MOBIUS_EOFA(ZWilsonImplFH); - INSTANTIATE_DPERP_MOBIUS_EOFA(ZWilsonImplDF); +INSTANTIATE_DPERP_MOBIUS_EOFA(GparityWilsonImplFH); +INSTANTIATE_DPERP_MOBIUS_EOFA(GparityWilsonImplDF); +INSTANTIATE_DPERP_MOBIUS_EOFA(WilsonImplFH); +INSTANTIATE_DPERP_MOBIUS_EOFA(WilsonImplDF); +INSTANTIATE_DPERP_MOBIUS_EOFA(ZWilsonImplFH); +INSTANTIATE_DPERP_MOBIUS_EOFA(ZWilsonImplDF); - template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); - template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); - template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); - template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); - template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); - template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); +template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); +template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); +template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); +template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); +template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); +template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); - #endif +#endif -}} +NAMESPACE_END(Grid);