From 202a7fe900461bf4f65947abc637b6c679839c15 Mon Sep 17 00:00:00 2001 From: David Murphy Date: Tue, 15 Aug 2017 13:36:08 -0400 Subject: [PATCH] Re-import DWF and abstract base EOFA fermion classes and tests --- lib/qcd/action/fermion/AbstractEOFAFermion.h | 100 +++ .../action/fermion/DomainWallEOFAFermion.cc | 440 +++++++++++++ .../action/fermion/DomainWallEOFAFermion.h | 115 ++++ .../fermion/DomainWallEOFAFermioncache.cc | 248 +++++++ .../fermion/DomainWallEOFAFermiondense.cc | 159 +++++ .../fermion/DomainWallEOFAFermionssp.cc | 168 +++++ .../fermion/DomainWallEOFAFermionvec.cc | 605 ++++++++++++++++++ lib/qcd/action/fermion/Fermion.h | 26 +- lib/qcd/modules/FermionOperatorModules.h | 49 +- tests/core/Test_dwf_eofa_even_odd.cc | 239 +++++++ tests/debug/Test_reweight_dwf_eofa.cc | 206 ++++++ 11 files changed, 2346 insertions(+), 9 deletions(-) create mode 100644 lib/qcd/action/fermion/AbstractEOFAFermion.h create mode 100644 lib/qcd/action/fermion/DomainWallEOFAFermion.cc create mode 100644 lib/qcd/action/fermion/DomainWallEOFAFermion.h create mode 100644 lib/qcd/action/fermion/DomainWallEOFAFermioncache.cc create mode 100644 lib/qcd/action/fermion/DomainWallEOFAFermiondense.cc create mode 100644 lib/qcd/action/fermion/DomainWallEOFAFermionssp.cc create mode 100644 lib/qcd/action/fermion/DomainWallEOFAFermionvec.cc create mode 100644 tests/core/Test_dwf_eofa_even_odd.cc create mode 100644 tests/debug/Test_reweight_dwf_eofa.cc diff --git a/lib/qcd/action/fermion/AbstractEOFAFermion.h b/lib/qcd/action/fermion/AbstractEOFAFermion.h new file mode 100644 index 00000000..abe06b8c --- /dev/null +++ b/lib/qcd/action/fermion/AbstractEOFAFermion.h @@ -0,0 +1,100 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/action/fermion/AbstractEOFAFermion.h + +Copyright (C) 2017 + +Author: Peter Boyle +Author: Peter Boyle +Author: David Murphy + +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_ABSTRACT_EOFA_FERMION_H +#define GRID_QCD_ABSTRACT_EOFA_FERMION_H + +#include + +namespace Grid { +namespace QCD { + + // DJM: Abstract base class for EOFA fermion types. + // Defines layout of additional EOFA-specific parameters and operators. + // Use to construct EOFA pseudofermion actions that are agnostic to Shamir / Mobius / etc., + // and ensure that no one can construct EOFA pseudofermion action with non-EOFA fermion type. + template + class AbstractEOFAFermion : public CayleyFermion5D { + public: + INHERIT_IMPL_TYPES(Impl); + + public: + // Fermion operator: D(mq1) + shift*\gamma_{5}*R_{5}*\Delta_{\pm}(mq2,mq3)*P_{\pm} + RealD mq1; + RealD mq2; + RealD mq3; + RealD shift; + int pm; + + RealD alpha; // Mobius scale + RealD k; // EOFA normalization constant + + virtual void Instantiatable(void) = 0; + + // EOFA-specific operations + // Force user to implement in derived classes + virtual void Omega (const FermionField& in, FermionField& out, int sign, int dag) = 0; + virtual void Dtilde (const FermionField& in, FermionField& out) = 0; + virtual void DtildeInv(const FermionField& in, FermionField& out) = 0; + + // Implement derivatives in base clcass: for EOFA both DWF and Mobius just need d(Dw)/dU + virtual void MDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag){ + this->DhopDeriv(mat, U, V, dag); + }; + virtual void MoeDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag){ + this->DhopDerivOE(mat, U, V, dag); + }; + virtual void MeoDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag){ + this->DhopDerivEO(mat, U, V, dag); + }; + + // Recompute 5D coefficients for different value of shift constant (needed for heatbath loop over poles) + virtual void RefreshShiftCoefficients(RealD new_shift) = 0; + + // Constructors + AbstractEOFAFermion(GaugeField& _Umu, + GridCartesian& FiveDimGrid, + GridRedBlackCartesian& FiveDimRedBlackGrid, + GridCartesian& FourDimGrid, + GridRedBlackCartesian& FourDimRedBlackGrid, + RealD _mq1, RealD _mq2, RealD _mq3, RealD _shift, int _pm, + RealD _M5, RealD _b, RealD _c, const ImplParams &p= ImplParams()) : + CayleyFermion5D(_Umu, FiveDimGrid, FiveDimRedBlackGrid, FourDimGrid, FourDimRedBlackGrid, _mq1, _M5, p), + mq1(_mq1), mq2(_mq2), mq3(_mq3), shift(_shift), pm(_pm) + { + int Ls = this->Ls; + this->alpha = _b + _c; + this->k = this->alpha * (_mq3 - _mq2) * std::pow(this->alpha+1.0,2*Ls) / + ( std::pow(alpha+1.0,Ls) + _mq2*std::pow(alpha-1.0,Ls) ) / + ( std::pow(alpha+1.0,Ls) + _mq3*std::pow(alpha-1.0,Ls) ); + } + }; +}} + +#endif diff --git a/lib/qcd/action/fermion/DomainWallEOFAFermion.cc b/lib/qcd/action/fermion/DomainWallEOFAFermion.cc new file mode 100644 index 00000000..dc4f6504 --- /dev/null +++ b/lib/qcd/action/fermion/DomainWallEOFAFermion.cc @@ -0,0 +1,440 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/action/fermion/DomainWallEOFAFermion.cc + +Copyright (C) 2017 + +Author: Peter Boyle +Author: Peter Boyle +Author: Peter Boyle +Author: paboyle +Author: David Murphy + +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 +#include +#include + +namespace Grid { +namespace QCD { + + template + DomainWallEOFAFermion::DomainWallEOFAFermion( + GaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _mq1, RealD _mq2, RealD _mq3, + RealD _shift, int _pm, RealD _M5, const ImplParams &p) : + AbstractEOFAFermion(_Umu, FiveDimGrid, FiveDimRedBlackGrid, + FourDimGrid, FourDimRedBlackGrid, _mq1, _mq2, _mq3, + _shift, _pm, _M5, 1.0, 0.0, p) + { + RealD eps = 1.0; + Approx::zolotarev_data *zdata = Approx::higham(eps,this->Ls); + assert(zdata->n == this->Ls); + + std::cout << GridLogMessage << "DomainWallEOFAFermion with Ls=" << this->Ls << std::endl; + this->SetCoefficientsTanh(zdata, 1.0, 0.0); + + Approx::zolotarev_free(zdata); + } + + /*************************************************************** + /* Additional EOFA operators only called outside the inverter. + /* Since speed is not essential, simple axpby-style + /* implementations should be fine. + /***************************************************************/ + template + void DomainWallEOFAFermion::Omega(const FermionField& psi, FermionField& Din, int sign, int dag) + { + int Ls = this->Ls; + + Din = zero; + if((sign == 1) && (dag == 0)){ axpby_ssp(Din, 0.0, psi, 1.0, psi, Ls-1, 0); } + else if((sign == -1) && (dag == 0)){ axpby_ssp(Din, 0.0, psi, 1.0, psi, 0, 0); } + else if((sign == 1 ) && (dag == 1)){ axpby_ssp(Din, 0.0, psi, 1.0, psi, 0, Ls-1); } + else if((sign == -1) && (dag == 1)){ axpby_ssp(Din, 0.0, psi, 1.0, psi, 0, 0); } + } + + // This is just the identity for DWF + template + void DomainWallEOFAFermion::Dtilde(const FermionField& psi, FermionField& chi){ chi = psi; } + + // This is just the identity for DWF + template + void DomainWallEOFAFermion::DtildeInv(const FermionField& psi, FermionField& chi){ chi = psi; } + + /*****************************************************************************************************/ + + template + RealD DomainWallEOFAFermion::M(const FermionField& psi, FermionField& chi) + { + int Ls = this->Ls; + + FermionField Din(psi._grid); + + this->Meooe5D(psi, Din); + this->DW(Din, chi, DaggerNo); + axpby(chi, 1.0, 1.0, chi, psi); + this->M5D(psi, chi); + return(norm2(chi)); + } + + template + RealD DomainWallEOFAFermion::Mdag(const FermionField& psi, FermionField& chi) + { + int Ls = this->Ls; + + FermionField Din(psi._grid); + + this->DW(psi, Din, DaggerYes); + this->MeooeDag5D(Din, chi); + this->M5Ddag(psi, chi); + axpby(chi, 1.0, 1.0, chi, psi); + return(norm2(chi)); + } + + /******************************************************************** + /* Performance critical fermion operators called inside the inverter + /********************************************************************/ + + template + void DomainWallEOFAFermion::M5D(const FermionField& psi, FermionField& chi) + { + int Ls = this->Ls; + int pm = this->pm; + RealD shift = this->shift; + RealD mq1 = this->mq1; + RealD mq2 = this->mq2; + RealD mq3 = this->mq3; + + // coefficients for shift operator ( = shift*\gamma_{5}*R_{5}*\Delta_{\pm}(mq2,mq3)*P_{\pm} ) + Coeff_t shiftp(0.0), shiftm(0.0); + if(shift != 0.0){ + if(pm == 1){ shiftp = shift*(mq3-mq2); } + else{ shiftm = -shift*(mq3-mq2); } + } + + std::vector diag(Ls,1.0); + std::vector upper(Ls,-1.0); upper[Ls-1] = mq1 + shiftm; + std::vector lower(Ls,-1.0); lower[0] = mq1 + shiftp; + + #if(0) + std::cout << GridLogMessage << "DomainWallEOFAFermion::M5D(FF&,FF&):" << std::endl; + for(int i=0; i::iscomplex()) { + sp[l] = PplusMat (l*istride+s1*ostride,s2); + sm[l] = PminusMat(l*istride+s1*ostride,s2); + } else { + // if real + scalar_type tmp; + tmp = PplusMat (l*istride+s1*ostride,s2); + sp[l] = scalar_type(tmp.real(),tmp.real()); + tmp = PminusMat(l*istride+s1*ostride,s2); + sm[l] = scalar_type(tmp.real(),tmp.real()); + } + } + Matp[LLs*s2+s1] = Vp; + Matm[LLs*s2+s1] = Vm; + }} + } + + FermOpTemplateInstantiate(DomainWallEOFAFermion); + GparityFermOpTemplateInstantiate(DomainWallEOFAFermion); + +}} diff --git a/lib/qcd/action/fermion/DomainWallEOFAFermion.h b/lib/qcd/action/fermion/DomainWallEOFAFermion.h new file mode 100644 index 00000000..179736ba --- /dev/null +++ b/lib/qcd/action/fermion/DomainWallEOFAFermion.h @@ -0,0 +1,115 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/action/fermion/DomainWallEOFAFermion.h + +Copyright (C) 2017 + +Author: Peter Boyle +Author: Peter Boyle +Author: David Murphy + +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_DOMAIN_WALL_EOFA_FERMION_H +#define GRID_QCD_DOMAIN_WALL_EOFA_FERMION_H + +#include + +namespace Grid { +namespace QCD { + + // DJM: EOFA with (Shamir) domain wall fermions. + // We overload and re-implement only the routines which have a different operator + // structure than the CayleyFermion5D base class. + template + class DomainWallEOFAFermion : public AbstractEOFAFermion + { + public: + INHERIT_IMPL_TYPES(Impl); + + public: + // Modified (0,Ls-1) and (Ls-1,0) elements of Mooee for red-black preconditioned Shamir EOFA + Coeff_t dm; + Coeff_t dp; + + virtual void Instantiatable(void) {}; + + // EOFA specific operators + virtual void Omega (const FermionField& in, FermionField &out, int sign, int dag); + virtual void Dtilde (const FermionField& in, FermionField &out); + virtual void DtildeInv (const FermionField& in, FermionField &out); + + // override multiply + virtual RealD M (const FermionField& in, FermionField& out); + virtual RealD Mdag (const FermionField& in, FermionField& out); + + // half checkerboard operations + 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); + + virtual void M5D (const FermionField &psi, FermionField &chi); + virtual void M5Ddag (const FermionField &psi, FermionField &chi); + + ///////////////////////////////////////////////////// + // Instantiate different versions depending on Impl + ///////////////////////////////////////////////////// + void M5D (const FermionField &psi, const FermionField &phi, FermionField &chi, + std::vector &lower, std::vector &diag, std::vector &upper); + void M5Ddag (const FermionField &psi, const FermionField &phi, FermionField &chi, + std::vector &lower, std::vector &diag, std::vector &upper); + void MooeeInternal (const FermionField &in, FermionField &out, int dag, int inv); + void MooeeInternalCompute(int dag, int inv, Vector >& Matp, Vector >& Matm); + void MooeeInternalAsm (const FermionField &in, FermionField &out, int LLs, int site, + Vector >& Matp, Vector >& Matm); + void MooeeInternalZAsm (const FermionField &in, FermionField &out, int LLs, int site, + Vector >& Matp, Vector >& Matm); + + virtual void RefreshShiftCoefficients(RealD new_shift); + + // Constructors + DomainWallEOFAFermion(GaugeField& _Umu, + GridCartesian& FiveDimGrid, + GridRedBlackCartesian& FiveDimRedBlackGrid, + GridCartesian& FourDimGrid, + GridRedBlackCartesian& FourDimRedBlackGrid, + RealD _mq1, RealD _mq2, RealD _mq3, RealD _shift, + int pm, RealD _M5, const ImplParams& p=ImplParams()); + + protected: + virtual void SetCoefficientsInternal(RealD zolo_hi, std::vector &gamma, RealD b, RealD c); + }; +}} + +#define INSTANTIATE_DPERP_DWF_EOFA(A)\ +template void DomainWallEOFAFermion::M5D(const FermionField& psi, const FermionField& phi, FermionField& chi,\ + std::vector& lower, std::vector& diag, std::vector& upper); \ +template void DomainWallEOFAFermion::M5Ddag(const FermionField& psi, const FermionField& phi, FermionField& chi,\ + std::vector& lower, std::vector& diag, std::vector& upper); \ +template void DomainWallEOFAFermion::MooeeInv(const FermionField& psi, FermionField& chi); \ +template void DomainWallEOFAFermion::MooeeInvDag(const FermionField& psi, FermionField& chi); + +#undef DOMAIN_WALL_EOFA_DPERP_DENSE +#define DOMAIN_WALL_EOFA_DPERP_CACHE +#undef DOMAIN_WALL_EOFA_DPERP_LINALG +#define DOMAIN_WALL_EOFA_DPERP_VEC + +#endif diff --git a/lib/qcd/action/fermion/DomainWallEOFAFermioncache.cc b/lib/qcd/action/fermion/DomainWallEOFAFermioncache.cc new file mode 100644 index 00000000..0b214d31 --- /dev/null +++ b/lib/qcd/action/fermion/DomainWallEOFAFermioncache.cc @@ -0,0 +1,248 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/action/fermion/DomainWallEOFAFermioncache.cc + +Copyright (C) 2017 + +Author: Peter Boyle +Author: Peter Boyle +Author: Peter Boyle +Author: paboyle +Author: David Murphy + +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 +#include + +namespace Grid { +namespace QCD { + + // 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) + { + int Ls = this->Ls; + GridBase* grid = psi._grid; + + assert(phi.checkerboard == psi.checkerboard); + chi.checkerboard = psi.checkerboard; + // Flops = 6.0*(Nc*Ns) *Ls*vol + this->M5Dcalls++; + this->M5Dtime -= usecond(); + + parallel_for(int ss=0; ssoSites(); ss+=Ls){ // adds Ls + for(int s=0; sM5Dtime += usecond(); + } + + template + void DomainWallEOFAFermion::M5Ddag(const FermionField& psi, const FermionField& phi, + FermionField& chi, std::vector& lower, std::vector& diag, std::vector& upper) + { + int Ls = this->Ls; + GridBase* grid = psi._grid; + assert(phi.checkerboard == psi.checkerboard); + chi.checkerboard=psi.checkerboard; + + // Flops = 6.0*(Nc*Ns) *Ls*vol + this->M5Dcalls++; + this->M5Dtime -= usecond(); + + parallel_for(int ss=0; ssoSites(); ss+=Ls){ // adds Ls + auto tmp = psi._odata[0]; + for(int s=0; sM5Dtime += usecond(); + } + + template + void DomainWallEOFAFermion::MooeeInv(const FermionField& psi, FermionField& chi) + { + GridBase* grid = psi._grid; + int Ls = this->Ls; + + chi.checkerboard = psi.checkerboard; + + this->MooeeInvCalls++; + this->MooeeInvTime -= usecond(); + + parallel_for(int ss=0; ssoSites(); ss+=Ls){ // adds Ls + + auto tmp1 = psi._odata[0]; + auto tmp2 = psi._odata[0]; + + // flops = 12*2*Ls + 12*2*Ls + 3*12*Ls + 12*2*Ls = 12*Ls * (9) = 108*Ls flops + // Apply (L^{\prime})^{-1} + chi[ss] = psi[ss]; // chi[0]=psi[0] + for(int s=1; slee[s-1]*tmp1; + } + + // L_m^{-1} + for(int s=0; sleem[s]*tmp1; + } + + // U_m^{-1} D^{-1} + for(int s=0; sdee[s])*chi[ss+s] - (this->ueem[s]/this->dee[Ls])*tmp1; + } + spProj5m(tmp2, chi[ss+Ls-1]); + chi[ss+Ls-1] = (1.0/this->dee[Ls])*tmp1 + (1.0/this->dee[Ls-1])*tmp2; + + // Apply U^{-1} + for(int s=Ls-2; s>=0; s--){ + spProj5m(tmp1, chi[ss+s+1]); + chi[ss+s] = chi[ss+s] - this->uee[s]*tmp1; + } + } + + this->MooeeInvTime += usecond(); + } + + template + void DomainWallEOFAFermion::MooeeInvDag(const FermionField& psi, FermionField& chi) + { + GridBase* grid = psi._grid; + int Ls = this->Ls; + + assert(psi.checkerboard == psi.checkerboard); + chi.checkerboard = psi.checkerboard; + + std::vector ueec(Ls); + std::vector deec(Ls+1); + std::vector leec(Ls); + std::vector ueemc(Ls); + std::vector leemc(Ls); + + for(int s=0; suee[s]); + deec[s] = conjugate(this->dee[s]); + leec[s] = conjugate(this->lee[s]); + ueemc[s] = conjugate(this->ueem[s]); + leemc[s] = conjugate(this->leem[s]); + } + deec[Ls] = conjugate(this->dee[Ls]); + + this->MooeeInvCalls++; + this->MooeeInvTime -= usecond(); + + parallel_for(int ss=0; ssoSites(); ss+=Ls){ // adds Ls + + auto tmp1 = psi._odata[0]; + auto tmp2 = psi._odata[0]; + + // Apply (U^{\prime})^{-dagger} + chi[ss] = psi[ss]; + for(int s=1; s=0; s--){ + spProj5p(tmp1, chi[ss+s+1]); + chi[ss+s] = chi[ss+s] - leec[s]*tmp1; + } + } + + this->MooeeInvTime += usecond(); + } + + #ifdef DOMAIN_WALL_EOFA_DPERP_CACHE + + 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); + + #endif + +}} diff --git a/lib/qcd/action/fermion/DomainWallEOFAFermiondense.cc b/lib/qcd/action/fermion/DomainWallEOFAFermiondense.cc new file mode 100644 index 00000000..c27074d9 --- /dev/null +++ b/lib/qcd/action/fermion/DomainWallEOFAFermiondense.cc @@ -0,0 +1,159 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/action/fermion/DomainWallEOFAFermiondense.cc + +Copyright (C) 2017 + +Author: Peter Boyle +Author: Peter Boyle +Author: Peter Boyle +Author: paboyle +Author: David Murphy + +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 +#include +#include + +namespace Grid { +namespace QCD { + + /* + * 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::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 + +}} diff --git a/lib/qcd/action/fermion/DomainWallEOFAFermionssp.cc b/lib/qcd/action/fermion/DomainWallEOFAFermionssp.cc new file mode 100644 index 00000000..80a4bf09 --- /dev/null +++ b/lib/qcd/action/fermion/DomainWallEOFAFermionssp.cc @@ -0,0 +1,168 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/action/fermion/DomainWallEOFAFermionssp.cc + +Copyright (C) 2017 + +Author: Peter Boyle +Author: Peter Boyle +Author: Peter Boyle +Author: paboyle +Author: David Murphy + +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 +#include + +namespace Grid { +namespace QCD { + + // 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::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; + + 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] + } + + // 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); + + // 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; + + 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); + } + + // 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); + + // 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 + + 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); + + #endif + +}} diff --git a/lib/qcd/action/fermion/DomainWallEOFAFermionvec.cc b/lib/qcd/action/fermion/DomainWallEOFAFermionvec.cc new file mode 100644 index 00000000..81ce448c --- /dev/null +++ b/lib/qcd/action/fermion/DomainWallEOFAFermionvec.cc @@ -0,0 +1,605 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/action/fermion/DomainWallEOFAFermionvec.cc + +Copyright (C) 2017 + +Author: Peter Boyle +Author: Peter Boyle +Author: Peter Boyle +Author: paboyle +Author: David Murphy + +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 +#include + +namespace Grid { +namespace QCD { + + /* + * 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::M5D(const FermionField& psi, const FermionField& phi, + FermionField& chi, std::vector& lower, std::vector& diag, std::vector& upper) + { + GridBase* grid = psi._grid; + int Ls = this->Ls; + int LLs = grid->_rdimensions[0]; + const int nsimd = Simd::Nsimd(); + + Vector > u(LLs); + Vector > l(LLs); + Vector > d(LLs); + + assert(Ls/LLs == nsimd); + assert(phi.checkerboard == psi.checkerboard); + + chi.checkerboard = psi.checkerboard; + + // just directly address via type pun + typedef typename Simd::scalar_type scalar_type; + scalar_type* u_p = (scalar_type*) &u[0]; + scalar_type* l_p = (scalar_type*) &l[0]; + scalar_type* d_p = (scalar_type*) &d[0]; + + for(int o=0;oM5Dcalls++; + this->M5Dtime -= usecond(); + + assert(Nc == 3); + + parallel_for(int ss=0; ssoSites(); ss+=LLs){ // adds LLs + + #if 0 + + alignas(64) SiteHalfSpinor hp; + alignas(64) SiteHalfSpinor hm; + alignas(64) SiteSpinor fp; + alignas(64) SiteSpinor fm; + + for(int v=0; v= v){ rotate(hm, hm, nsimd-1); } + + hp = 0.5*hp; + hm = 0.5*hm; + + spRecon5m(fp, hp); + spRecon5p(fm, hm); + + chi[ss+v] = d[v]*phi[ss+v]; + chi[ss+v] = chi[ss+v] + u[v]*fp; + chi[ss+v] = chi[ss+v] + l[v]*fm; + + } + + #else + + for(int v=0; v(hp_00.v); + hp_01.v = Optimization::Rotate::tRotate<2>(hp_01.v); + hp_02.v = Optimization::Rotate::tRotate<2>(hp_02.v); + hp_10.v = Optimization::Rotate::tRotate<2>(hp_10.v); + hp_11.v = Optimization::Rotate::tRotate<2>(hp_11.v); + hp_12.v = Optimization::Rotate::tRotate<2>(hp_12.v); + } + + if(vm >= v){ + hm_00.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_00.v); + hm_01.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_01.v); + hm_02.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_02.v); + hm_10.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_10.v); + hm_11.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_11.v); + hm_12.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_12.v); + } + + // Can force these to real arithmetic and save 2x. + Simd p_00 = switcheroo::mult(d[v]()()(), phi[ss+v]()(0)(0)) + switcheroo::mult(l[v]()()(), hm_00); + Simd p_01 = switcheroo::mult(d[v]()()(), phi[ss+v]()(0)(1)) + switcheroo::mult(l[v]()()(), hm_01); + Simd p_02 = switcheroo::mult(d[v]()()(), phi[ss+v]()(0)(2)) + switcheroo::mult(l[v]()()(), hm_02); + Simd p_10 = switcheroo::mult(d[v]()()(), phi[ss+v]()(1)(0)) + switcheroo::mult(l[v]()()(), hm_10); + Simd p_11 = switcheroo::mult(d[v]()()(), phi[ss+v]()(1)(1)) + switcheroo::mult(l[v]()()(), hm_11); + Simd p_12 = switcheroo::mult(d[v]()()(), phi[ss+v]()(1)(2)) + switcheroo::mult(l[v]()()(), hm_12); + Simd p_20 = switcheroo::mult(d[v]()()(), phi[ss+v]()(2)(0)) + switcheroo::mult(u[v]()()(), hp_00); + Simd p_21 = switcheroo::mult(d[v]()()(), phi[ss+v]()(2)(1)) + switcheroo::mult(u[v]()()(), hp_01); + Simd p_22 = switcheroo::mult(d[v]()()(), phi[ss+v]()(2)(2)) + switcheroo::mult(u[v]()()(), hp_02); + Simd p_30 = switcheroo::mult(d[v]()()(), phi[ss+v]()(3)(0)) + switcheroo::mult(u[v]()()(), hp_10); + Simd p_31 = switcheroo::mult(d[v]()()(), phi[ss+v]()(3)(1)) + switcheroo::mult(u[v]()()(), hp_11); + Simd p_32 = switcheroo::mult(d[v]()()(), phi[ss+v]()(3)(2)) + switcheroo::mult(u[v]()()(), hp_12); + + vstream(chi[ss+v]()(0)(0), p_00); + vstream(chi[ss+v]()(0)(1), p_01); + vstream(chi[ss+v]()(0)(2), p_02); + vstream(chi[ss+v]()(1)(0), p_10); + vstream(chi[ss+v]()(1)(1), p_11); + vstream(chi[ss+v]()(1)(2), p_12); + vstream(chi[ss+v]()(2)(0), p_20); + vstream(chi[ss+v]()(2)(1), p_21); + vstream(chi[ss+v]()(2)(2), p_22); + vstream(chi[ss+v]()(3)(0), p_30); + vstream(chi[ss+v]()(3)(1), p_31); + vstream(chi[ss+v]()(3)(2), p_32); + } + + #endif + } + + this->M5Dtime += usecond(); + } + + template + void DomainWallEOFAFermion::M5Ddag(const FermionField& psi, const FermionField& phi, + FermionField& chi, std::vector& lower, std::vector& diag, std::vector& upper) + { + GridBase* grid = psi._grid; + int Ls = this->Ls; + int LLs = grid->_rdimensions[0]; + int nsimd = Simd::Nsimd(); + + Vector > u(LLs); + Vector > l(LLs); + Vector > d(LLs); + + assert(Ls/LLs == nsimd); + assert(phi.checkerboard == psi.checkerboard); + + chi.checkerboard = psi.checkerboard; + + // just directly address via type pun + typedef typename Simd::scalar_type scalar_type; + scalar_type* u_p = (scalar_type*) &u[0]; + scalar_type* l_p = (scalar_type*) &l[0]; + scalar_type* d_p = (scalar_type*) &d[0]; + + for(int o=0; oM5Dcalls++; + this->M5Dtime -= usecond(); + + parallel_for(int ss=0; ssoSites(); ss+=LLs){ // adds LLs + + #if 0 + + alignas(64) SiteHalfSpinor hp; + alignas(64) SiteHalfSpinor hm; + alignas(64) SiteSpinor fp; + alignas(64) SiteSpinor fm; + + for(int v=0; v= v){ rotate(hm, hm, nsimd-1); } + + hp = hp*0.5; + hm = hm*0.5; + spRecon5p(fp, hp); + spRecon5m(fm, hm); + + chi[ss+v] = d[v]*phi[ss+v]+u[v]*fp; + chi[ss+v] = chi[ss+v] +l[v]*fm; + } + + #else + + for(int v=0; v(hp_00.v); + hp_01.v = Optimization::Rotate::tRotate<2>(hp_01.v); + hp_02.v = Optimization::Rotate::tRotate<2>(hp_02.v); + hp_10.v = Optimization::Rotate::tRotate<2>(hp_10.v); + hp_11.v = Optimization::Rotate::tRotate<2>(hp_11.v); + hp_12.v = Optimization::Rotate::tRotate<2>(hp_12.v); + } + + if(vm >= v){ + hm_00.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_00.v); + hm_01.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_01.v); + hm_02.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_02.v); + hm_10.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_10.v); + hm_11.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_11.v); + hm_12.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_12.v); + } + + Simd p_00 = switcheroo::mult(d[v]()()(), phi[ss+v]()(0)(0)) + switcheroo::mult(u[v]()()(), hp_00); + Simd p_01 = switcheroo::mult(d[v]()()(), phi[ss+v]()(0)(1)) + switcheroo::mult(u[v]()()(), hp_01); + Simd p_02 = switcheroo::mult(d[v]()()(), phi[ss+v]()(0)(2)) + switcheroo::mult(u[v]()()(), hp_02); + Simd p_10 = switcheroo::mult(d[v]()()(), phi[ss+v]()(1)(0)) + switcheroo::mult(u[v]()()(), hp_10); + Simd p_11 = switcheroo::mult(d[v]()()(), phi[ss+v]()(1)(1)) + switcheroo::mult(u[v]()()(), hp_11); + Simd p_12 = switcheroo::mult(d[v]()()(), phi[ss+v]()(1)(2)) + switcheroo::mult(u[v]()()(), hp_12); + Simd p_20 = switcheroo::mult(d[v]()()(), phi[ss+v]()(2)(0)) + switcheroo::mult(l[v]()()(), hm_00); + Simd p_21 = switcheroo::mult(d[v]()()(), phi[ss+v]()(2)(1)) + switcheroo::mult(l[v]()()(), hm_01); + Simd p_22 = switcheroo::mult(d[v]()()(), phi[ss+v]()(2)(2)) + switcheroo::mult(l[v]()()(), hm_02); + Simd p_30 = switcheroo::mult(d[v]()()(), phi[ss+v]()(3)(0)) + switcheroo::mult(l[v]()()(), hm_10); + Simd p_31 = switcheroo::mult(d[v]()()(), phi[ss+v]()(3)(1)) + switcheroo::mult(l[v]()()(), hm_11); + Simd p_32 = switcheroo::mult(d[v]()()(), phi[ss+v]()(3)(2)) + switcheroo::mult(l[v]()()(), hm_12); + + vstream(chi[ss+v]()(0)(0), p_00); + vstream(chi[ss+v]()(0)(1), p_01); + vstream(chi[ss+v]()(0)(2), p_02); + vstream(chi[ss+v]()(1)(0), p_10); + vstream(chi[ss+v]()(1)(1), p_11); + vstream(chi[ss+v]()(1)(2), p_12); + vstream(chi[ss+v]()(2)(0), p_20); + vstream(chi[ss+v]()(2)(1), p_21); + vstream(chi[ss+v]()(2)(2), p_22); + vstream(chi[ss+v]()(3)(0), p_30); + vstream(chi[ss+v]()(3)(1), p_31); + vstream(chi[ss+v]()(3)(2), p_32); + } + #endif + + } + + this->M5Dtime += usecond(); + } + + #ifdef AVX512 + #include + #include + #include + #endif + + template + void DomainWallEOFAFermion::MooeeInternalAsm(const FermionField& psi, FermionField& chi, + int LLs, int site, Vector >& Matp, Vector >& Matm) + { + #ifndef AVX512 + { + SiteHalfSpinor BcastP; + SiteHalfSpinor BcastM; + SiteHalfSpinor SiteChiP; + SiteHalfSpinor SiteChiM; + + // Ls*Ls * 2 * 12 * vol flops + for(int s1=0; s1); + for(int s1=0; s1 + void DomainWallEOFAFermion::MooeeInternalZAsm(const FermionField& psi, FermionField& chi, + int LLs, int site, Vector >& Matp, Vector >& Matm) + { + std::cout << "Error: zMobius not implemented for EOFA" << std::endl; + exit(-1); + }; + + 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; + + Vector > Matp; + Vector > Matm; + Vector > *_Matp; + Vector > *_Matm; + + // MooeeInternalCompute(dag,inv,Matp,Matm); + if(inv && dag){ + _Matp = &this->MatpInvDag; + _Matm = &this->MatmInvDag; + } + + if(inv && (!dag)){ + _Matp = &this->MatpInv; + _Matm = &this->MatmInv; + } + + if(!inv){ + MooeeInternalCompute(dag, inv, Matp, Matm); + _Matp = &Matp; + _Matm = &Matm; + } + + assert(_Matp->size() == Ls*LLs); + + this->MooeeInvCalls++; + this->MooeeInvTime -= usecond(); + + if(switcheroo::iscomplex()){ + parallel_for(auto site=0; siteMooeeInvTime += usecond(); + } + + #ifdef DOMAIN_WALL_EOFA_DPERP_VEC + + INSTANTIATE_DPERP_DWF_EOFA(DomainWallVec5dImplD); + INSTANTIATE_DPERP_DWF_EOFA(DomainWallVec5dImplF); + INSTANTIATE_DPERP_DWF_EOFA(ZDomainWallVec5dImplD); + INSTANTIATE_DPERP_DWF_EOFA(ZDomainWallVec5dImplF); + + INSTANTIATE_DPERP_DWF_EOFA(DomainWallVec5dImplDF); + INSTANTIATE_DPERP_DWF_EOFA(DomainWallVec5dImplFH); + INSTANTIATE_DPERP_DWF_EOFA(ZDomainWallVec5dImplDF); + INSTANTIATE_DPERP_DWF_EOFA(ZDomainWallVec5dImplFH); + + 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); + 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 + +}} diff --git a/lib/qcd/action/fermion/Fermion.h b/lib/qcd/action/fermion/Fermion.h index 0f803f44..ac2a94b2 100644 --- a/lib/qcd/action/fermion/Fermion.h +++ b/lib/qcd/action/fermion/Fermion.h @@ -55,7 +55,7 @@ Author: Peter Boyle #include #include // Cayley types #include -#include +#include #include #include #include @@ -113,6 +113,14 @@ typedef DomainWallFermion DomainWallFermionRL; typedef DomainWallFermion DomainWallFermionFH; typedef DomainWallFermion DomainWallFermionDF; +typedef DomainWallEOFAFermion DomainWallEOFAFermionR; +typedef DomainWallEOFAFermion DomainWallEOFAFermionF; +typedef DomainWallEOFAFermion DomainWallEOFAFermionD; + +typedef DomainWallEOFAFermion DomainWallEOFAFermionRL; +typedef DomainWallEOFAFermion DomainWallEOFAFermionFH; +typedef DomainWallEOFAFermion DomainWallEOFAFermionDF; + typedef MobiusFermion MobiusFermionR; typedef MobiusFermion MobiusFermionF; typedef MobiusFermion MobiusFermionD; @@ -138,6 +146,14 @@ typedef DomainWallFermion DomainWallFermionVec5dRL; typedef DomainWallFermion DomainWallFermionVec5dFH; typedef DomainWallFermion DomainWallFermionVec5dDF; +typedef DomainWallEOFAFermion DomainWallEOFAFermionVec5dR; +typedef DomainWallEOFAFermion DomainWallEOFAFermionVec5dF; +typedef DomainWallEOFAFermion DomainWallEOFAFermionVec5dD; + +typedef DomainWallEOFAFermion DomainWallEOFAFermionVec5dRL; +typedef DomainWallEOFAFermion DomainWallEOFAFermionVec5dFH; +typedef DomainWallEOFAFermion DomainWallEOFAFermionVec5dDF; + typedef MobiusFermion MobiusFermionVec5dR; typedef MobiusFermion MobiusFermionVec5dF; typedef MobiusFermion MobiusFermionVec5dD; @@ -206,6 +222,14 @@ typedef DomainWallFermion GparityDomainWallFermionRL; typedef DomainWallFermion GparityDomainWallFermionFH; typedef DomainWallFermion GparityDomainWallFermionDF; +typedef DomainWallEOFAFermion GparityDomainWallEOFAFermionR; +typedef DomainWallEOFAFermion GparityDomainWallEOFAFermionF; +typedef DomainWallEOFAFermion GparityDomainWallEOFAFermionD; + +typedef DomainWallEOFAFermion GparityDomainWallEOFAFermionRL; +typedef DomainWallEOFAFermion GparityDomainWallEOFAFermionFH; +typedef DomainWallEOFAFermion GparityDomainWallEOFAFermionDF; + typedef WilsonTMFermion GparityWilsonTMFermionR; typedef WilsonTMFermion GparityWilsonTMFermionF; typedef WilsonTMFermion GparityWilsonTMFermionD; diff --git a/lib/qcd/modules/FermionOperatorModules.h b/lib/qcd/modules/FermionOperatorModules.h index c66842c6..fc9d96a7 100644 --- a/lib/qcd/modules/FermionOperatorModules.h +++ b/lib/qcd/modules/FermionOperatorModules.h @@ -72,7 +72,7 @@ protected: } virtual unsigned int Ls(){ - return 0; + return 0; } virtual void print_parameters(){ @@ -97,7 +97,7 @@ class HMC_FermionOperatorModuleFactory : public Factory < FermionOperatorModuleBase > , Reader > { public: // use SINGLETON FUNCTOR MACRO HERE - typedef Reader TheReader; + typedef Reader TheReader; HMC_FermionOperatorModuleFactory(const HMC_FermionOperatorModuleFactory& e) = delete; void operator=(const HMC_FermionOperatorModuleFactory& e) = delete; @@ -122,7 +122,7 @@ namespace QCD{ // Modules class WilsonFermionParameters : Serializable { public: - GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonFermionParameters, + GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonFermionParameters, RealD, mass); }; @@ -144,7 +144,7 @@ class WilsonFermionModule: public FermionOperatorModuleGridRefs[0]; auto GridMod5d = this->GridRefs[1]; typename FermionImpl::GaugeField U(GridMod->get_full()); - this->FOPtr.reset(new MobiusFermion( U, *(GridMod->get_full()), *(GridMod->get_rb()), + this->FOPtr.reset(new MobiusFermion( U, *(GridMod->get_full()), *(GridMod->get_rb()), *(GridMod5d->get_full()), *(GridMod5d->get_rb()), this->Par_.mass, this->Par_.M5, this->Par_.b, this->Par_.c)); } @@ -175,7 +175,7 @@ class MobiusFermionModule: public FermionOperatorModuleGridRefs[0]; auto GridMod5d = this->GridRefs[1]; typename FermionImpl::GaugeField U(GridMod->get_full()); - this->FOPtr.reset(new DomainWallFermion( U, *(GridMod->get_full()), *(GridMod->get_rb()), + this->FOPtr.reset(new DomainWallFermion( U, *(GridMod->get_full()), *(GridMod->get_rb()), *(GridMod5d->get_full()), *(GridMod5d->get_rb()), this->Par_.mass, this->Par_.M5)); } }; +class DomainWallEOFAFermionParameters : Serializable { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(DomainWallEOFAFermionParameters, + RealD, mq1, + RealD, mq2, + RealD, mq3, + RealD, shift, + int, pm, + RealD, M5, + unsigned int, Ls); +}; + +template +class DomainWallEOFAFermionModule: public FermionOperatorModule { + typedef FermionOperatorModule FermBase; + using FermBase::FermBase; // for constructors + + virtual unsigned int Ls(){ + return this->Par_.Ls; + } + + // acquire resource + virtual void initialize(){ + auto GridMod = this->GridRefs[0]; + auto GridMod5d = this->GridRefs[1]; + typename FermionImpl::GaugeField U(GridMod->get_full()); + this->FOPtr.reset(new DomainWallEOFAFermion( U, *(GridMod->get_full()), *(GridMod->get_rb()), + *(GridMod5d->get_full()), *(GridMod5d->get_rb()), + this->Par_.mq1, this->Par_.mq2, this->Par_.mq3, + this->Par_.shift, this->Par_.pm, this->Par_.M5)); + } +}; + } // QCD } // Grid -#endif //FERMIONOPERATOR_MODULES_H \ No newline at end of file +#endif //FERMIONOPERATOR_MODULES_H diff --git a/tests/core/Test_dwf_eofa_even_odd.cc b/tests/core/Test_dwf_eofa_even_odd.cc new file mode 100644 index 00000000..5fe0f653 --- /dev/null +++ b/tests/core/Test_dwf_eofa_even_odd.cc @@ -0,0 +1,239 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/core/Test_dwf_eofa_even_odd.cc + +Copyright (C) 2017 + +Author: Peter Boyle +Author: paboyle +Author: David Murphy + +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 std; +using namespace Grid; +using namespace Grid::QCD; + +template +struct scal { + d internal; +}; + +Gamma::Algebra Gmu [] = { + Gamma::Algebra::GammaX, + Gamma::Algebra::GammaY, + Gamma::Algebra::GammaZ, + Gamma::Algebra::GammaT +}; + +int main (int argc, char ** argv) +{ + Grid_init(&argc, &argv); + + int threads = GridThread::GetThreads(); + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + const int Ls = 8; + // GridCartesian* UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()), GridDefaultMpi()); + GridCartesian* UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()), GridDefaultMpi()); + GridCartesian* FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls, UGrid); + GridRedBlackCartesian* UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridRedBlackCartesian* FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, UGrid); + + std::vector seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + + LatticeFermion src (FGrid); random(RNG5, src); + LatticeFermion phi (FGrid); random(RNG5, phi); + LatticeFermion chi (FGrid); random(RNG5, chi); + LatticeFermion result(FGrid); result = zero; + LatticeFermion ref (FGrid); ref = zero; + LatticeFermion tmp (FGrid); tmp = zero; + LatticeFermion err (FGrid); err = zero; + LatticeGaugeField Umu (UGrid); SU3::HotConfiguration(RNG4, Umu); + std::vector U(4,UGrid); + + // Only one non-zero (y) + Umu = zero; + for(int nn=0; nn0){ U[nn] = zero; } + PokeIndex(Umu, U[nn], nn); + } + + RealD mq1 = 0.1; + RealD mq2 = 0.5; + RealD mq3 = 1.0; + RealD shift = 0.1234; + RealD M5 = 1.8; + int pm = 1; + DomainWallEOFAFermionR Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mq1, mq2, mq3, shift, pm, M5); + + LatticeFermion src_e (FrbGrid); + LatticeFermion src_o (FrbGrid); + LatticeFermion r_e (FrbGrid); + LatticeFermion r_o (FrbGrid); + LatticeFermion r_eo (FGrid); + LatticeFermion r_eeoo(FGrid); + + std::cout << GridLogMessage << "==========================================================" << std::endl; + std::cout << GridLogMessage << "= Testing that Meo + Moe + Moo + Mee = Munprec " << std::endl; + std::cout << GridLogMessage << "==========================================================" << std::endl; + + pickCheckerboard(Even, src_e, src); + pickCheckerboard(Odd, src_o, src); + + Ddwf.Meooe(src_e, r_o); std::cout << GridLogMessage << "Applied Meo" << std::endl; + Ddwf.Meooe(src_o, r_e); std::cout << GridLogMessage << "Applied Moe" << std::endl; + setCheckerboard(r_eo, r_o); + setCheckerboard(r_eo, r_e); + + Ddwf.Mooee(src_e, r_e); std::cout << GridLogMessage << "Applied Mee" << std::endl; + Ddwf.Mooee(src_o, r_o); std::cout << GridLogMessage << "Applied Moo" << std::endl; + setCheckerboard(r_eeoo, r_e); + setCheckerboard(r_eeoo, r_o); + + r_eo = r_eo + r_eeoo; + Ddwf.M(src, ref); + + // std::cout << GridLogMessage << r_eo << std::endl; + // std::cout << GridLogMessage << ref << std::endl; + + err = ref - r_eo; + std::cout << GridLogMessage << "EO norm diff " << norm2(err) << " " << norm2(ref) << " " << norm2(r_eo) << std::endl; + + LatticeComplex cerr(FGrid); + cerr = localInnerProduct(err,err); + // std::cout << GridLogMessage << cerr << std::endl; + + std::cout << GridLogMessage << "==============================================================" << std::endl; + std::cout << GridLogMessage << "= Test Ddagger is the dagger of D by requiring " << std::endl; + std::cout << GridLogMessage << "= < phi | Deo | chi > * = < chi | Deo^dag| phi> " << std::endl; + std::cout << GridLogMessage << "==============================================================" << std::endl; + + LatticeFermion chi_e (FrbGrid); + LatticeFermion chi_o (FrbGrid); + + LatticeFermion dchi_e(FrbGrid); + LatticeFermion dchi_o(FrbGrid); + + LatticeFermion phi_e (FrbGrid); + LatticeFermion phi_o (FrbGrid); + + LatticeFermion dphi_e(FrbGrid); + LatticeFermion dphi_o(FrbGrid); + + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd , chi_o, chi); + pickCheckerboard(Even, phi_e, phi); + pickCheckerboard(Odd , phi_o, phi); + + Ddwf.Meooe (chi_e, dchi_o); + Ddwf.Meooe (chi_o, dchi_e); + Ddwf.MeooeDag(phi_e, dphi_o); + Ddwf.MeooeDag(phi_o, dphi_e); + + ComplexD pDce = innerProduct(phi_e, dchi_e); + ComplexD pDco = innerProduct(phi_o, dchi_o); + ComplexD cDpe = innerProduct(chi_e, dphi_e); + ComplexD cDpo = innerProduct(chi_o, dphi_o); + + std::cout << GridLogMessage << "e " << pDce << " " << cDpe << std::endl; + std::cout << GridLogMessage << "o " << pDco << " " << cDpo << std::endl; + + std::cout << GridLogMessage << "pDce - conj(cDpo) " << pDce-conj(cDpo) << std::endl; + std::cout << GridLogMessage << "pDco - conj(cDpe) " << pDco-conj(cDpe) << std::endl; + + std::cout << GridLogMessage << "==============================================================" << std::endl; + std::cout << GridLogMessage << "= Test MeeInv Mee = 1 " << std::endl; + std::cout << GridLogMessage << "==============================================================" << std::endl; + + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd , chi_o, chi); + + Ddwf.Mooee (chi_e, src_e); + Ddwf.MooeeInv(src_e, phi_e); + + Ddwf.Mooee (chi_o, src_o); + Ddwf.MooeeInv(src_o, phi_o); + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = phi - chi; + std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + + std::cout << GridLogMessage << "==============================================================" << std::endl; + std::cout << GridLogMessage << "= Test MeeInvDag MeeDag = 1 " << std::endl; + std::cout << GridLogMessage << "==============================================================" << std::endl; + + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd , chi_o, chi); + + Ddwf.MooeeDag (chi_e, src_e); + Ddwf.MooeeInvDag(src_e, phi_e); + + Ddwf.MooeeDag (chi_o, src_o); + Ddwf.MooeeInvDag(src_o, phi_o); + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = phi - chi; + std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + + std::cout << GridLogMessage << "==============================================================" << std::endl; + std::cout << GridLogMessage << "= Test MpcDagMpc is Hermitian " << std::endl; + std::cout << GridLogMessage << "==============================================================" << std::endl; + + random(RNG5, phi); + random(RNG5, chi); + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd , chi_o, chi); + pickCheckerboard(Even, phi_e, phi); + pickCheckerboard(Odd , phi_o, phi); + RealD t1,t2; + + SchurDiagMooeeOperator HermOpEO(Ddwf); + HermOpEO.MpcDagMpc(chi_e, dchi_e, t1, t2); + HermOpEO.MpcDagMpc(chi_o, dchi_o, t1, t2); + + HermOpEO.MpcDagMpc(phi_e, dphi_e, t1, t2); + HermOpEO.MpcDagMpc(phi_o, dphi_o, t1, t2); + + pDce = innerProduct(phi_e, dchi_e); + pDco = innerProduct(phi_o, dchi_o); + cDpe = innerProduct(chi_e, dphi_e); + cDpo = innerProduct(chi_o, dphi_o); + + std::cout << GridLogMessage << "e " << pDce << " " << cDpe << std::endl; + std::cout << GridLogMessage << "o " << pDco << " " << cDpo << std::endl; + + std::cout << GridLogMessage << "pDce - conj(cDpo) " << pDco-conj(cDpo) << std::endl; + std::cout << GridLogMessage << "pDco - conj(cDpe) " << pDce-conj(cDpe) << std::endl; + + Grid_finalize(); +} diff --git a/tests/debug/Test_reweight_dwf_eofa.cc b/tests/debug/Test_reweight_dwf_eofa.cc new file mode 100644 index 00000000..98a17e2f --- /dev/null +++ b/tests/debug/Test_reweight_dwf_eofa.cc @@ -0,0 +1,206 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/debug/Test_reweight_dwf_eofa.cc + +Copyright (C) 2017 + +Author: Peter Boyle +Author: paboyle +Author: David Murphy + +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 std; +using namespace Grid; +using namespace Grid::QCD; + +// parameters for test +const std::vector grid_dim = { 8, 8, 8, 8 }; +const int Ls = 8; +const int Nhits = 25; +const int max_iter = 5000; +const RealD mf = 0.1; +const RealD mb = 0.11; +const RealD M5 = 1.8; +const RealD stop_tol = 1.0e-12; + +RealD mean(const std::vector& data) +{ + int N = data.size(); + RealD mean(0.0); + for(int i=0; i& data, int sample) +{ + int N = data.size(); + RealD mean(0.0); + for(int i=0; i& jacks, RealD mean) +{ + int N = jacks.size(); + RealD std(0.0); + for(int i=0; i jack_stats(const std::vector& data) +{ + int N = data.size(); + std::vector jack_samples(N); + std::vector jack_stats(2); + + jack_stats[0] = mean(data); + for(int i=0; i seeds4({1, 2, 3, 4}); + std::vector seeds5({5, 6, 7, 8}); + GridParallelRNG RNG5(FGrid); + RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); + RNG4.SeedFixedIntegers(seeds4); + + // Random gauge field + LatticeGaugeField Umu(UGrid); + SU3::HotConfiguration(RNG4, Umu); + + // Initialize RHMC fermion operators + DomainWallFermionR Ddwf_f(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, M5); + DomainWallFermionR Ddwf_b(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mb, M5); + SchurDiagMooeeOperator MdagM(Ddwf_f); + SchurDiagMooeeOperator VdagV(Ddwf_b); + + // Degree 12 rational approximations to x^(1/4) and x^(-1/4) + double lo = 0.0001; + double hi = 95.0; + int precision = 64; + int degree = 12; + AlgRemez remez(lo, hi, precision); + std::cout << GridLogMessage << "Generating degree " << degree << " for x^(1/4)" << std::endl; + remez.generateApprox(degree, 1, 4); + MultiShiftFunction PowerQuarter(remez, stop_tol, false); + MultiShiftFunction PowerNegQuarter(remez, stop_tol, true); + + // Stochastically estimate reweighting factor via RHMC + RealD scale = std::sqrt(0.5); + std::vector rw_rhmc(Nhits); + ConjugateGradientMultiShift msCG_V(max_iter, PowerQuarter); + ConjugateGradientMultiShift msCG_M(max_iter, PowerNegQuarter); + std::cout.precision(12); + + for(int hit=0; hit tmp(2, Ddwf_f.FermionRedBlackGrid()); + gaussian(RNG5, Phi); + Phi = Phi*scale; + + pickCheckerboard(Odd, PhiOdd, Phi); + + // evaluate -log(rw) + msCG_V(VdagV, PhiOdd, tmp[0]); + msCG_M(MdagM, tmp[0], tmp[1]); + rw_rhmc[hit] = norm2(tmp[1]) - norm2(PhiOdd); + std::cout << std::endl << "==================================================" << std::endl; + std::cout << " --- RHMC: Hit " << hit << ": rw = " << rw_rhmc[hit]; + std::cout << std::endl << "==================================================" << std::endl << std::endl; + + } + + // Initialize EOFA fermion operators + RealD shift_L = 0.0; + RealD shift_R = -1.0; + int pm = 1; + DomainWallEOFAFermionR Deofa_L(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, mf, mb, shift_L, pm, M5); + DomainWallEOFAFermionR Deofa_R(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mb, mf, mb, shift_R, pm, M5); + MdagMLinearOperator LdagL(Deofa_L); + MdagMLinearOperator RdagR(Deofa_R); + + // Stochastically estimate reweighting factor via EOFA + RealD k = Deofa_L.k; + std::vector rw_eofa(Nhits); + ConjugateGradient CG(stop_tol, max_iter); + SchurRedBlackDiagMooeeSolve SchurSolver(CG); + + for(int hit=0; hit tmp(2, Deofa_L.FermionGrid()); + gaussian(RNG5, Phi); + Phi = Phi*scale; + + // evaluate -log(rw) + // LH term + for(int s=0; s rhmc_result = jack_stats(rw_rhmc); + std::vector eofa_result = jack_stats(rw_eofa); + std::cout << std::endl << "RHMC: rw = " << rhmc_result[0] << " +/- " << rhmc_result[1] << std::endl; + std::cout << std::endl << "EOFA: rw = " << eofa_result[0] << " +/- " << eofa_result[1] << std::endl; + + Grid_finalize(); +}