diff --git a/lib/qcd/action/fermion/AbstractEOFAFermion.h b/lib/qcd/action/fermion/AbstractEOFAFermion.h index abe06b8c..15faa401 100644 --- a/lib/qcd/action/fermion/AbstractEOFAFermion.h +++ b/lib/qcd/action/fermion/AbstractEOFAFermion.h @@ -35,66 +35,66 @@ See the full license in the file "LICENSE" in the top level distribution directo 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); + // 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; + 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 + RealD alpha; // Mobius scale + RealD k; // EOFA normalization constant - virtual void Instantiatable(void) = 0; + 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; + // 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); - }; + // Implement derivatives in base class: + // 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; + // 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) ); - } - }; + // 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(this->alpha+1.0,Ls) + _mq2*std::pow(this->alpha-1.0,Ls) ) / + ( std::pow(this->alpha+1.0,Ls) + _mq3*std::pow(this->alpha-1.0,Ls) ); + }; + }; }} #endif diff --git a/lib/qcd/action/fermion/DomainWallEOFAFermion.h b/lib/qcd/action/fermion/DomainWallEOFAFermion.h index 179736ba..d48e3b8f 100644 --- a/lib/qcd/action/fermion/DomainWallEOFAFermion.h +++ b/lib/qcd/action/fermion/DomainWallEOFAFermion.h @@ -35,75 +35,93 @@ See the full license in the file "LICENSE" in the top level distribution directo 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); + /*template struct switcheroo { + static inline int iscomplex() { return 0; } + template + static inline vec mult(vec a, vec b){ return real_mult(a,b); } + }; - public: - // Modified (0,Ls-1) and (Ls-1,0) elements of Mooee for red-black preconditioned Shamir EOFA - Coeff_t dm; - Coeff_t dp; + template<> struct switcheroo { + static inline int iscomplex() { return 1; } + template + static inline vec mult(vec a, vec b){ return a*b; } + }; - virtual void Instantiatable(void) {}; + template<> struct switcheroo { + static inline int iscomplex() { return 1; } + template + static inline vec mult(vec a, vec b) { return a*b; } + };*/ - // 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); + template + class DomainWallEOFAFermion : public AbstractEOFAFermion + { + public: + INHERIT_IMPL_TYPES(Impl); - // override multiply - virtual RealD M (const FermionField& in, FermionField& out); - virtual RealD Mdag (const FermionField& in, FermionField& out); + public: + // Modified (0,Ls-1) and (Ls-1,0) elements of Mooee + // for red-black preconditioned Shamir EOFA + Coeff_t dm; + Coeff_t dp; - // 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 Instantiatable(void) {}; - virtual void M5D (const FermionField &psi, FermionField &chi); - virtual void M5Ddag (const FermionField &psi, FermionField &chi); + // EOFA-specific operations + 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); - ///////////////////////////////////////////////////// - // 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); + // override multiply + virtual RealD M (const FermionField& in, FermionField& out); + virtual RealD Mdag (const FermionField& in, FermionField& out); - virtual void RefreshShiftCoefficients(RealD new_shift); + // 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); - // 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()); + virtual void M5D (const FermionField& psi, FermionField& chi); + virtual void M5Ddag (const FermionField& psi, FermionField& chi); - protected: - virtual void SetCoefficientsInternal(RealD zolo_hi, std::vector &gamma, RealD b, RealD c); - }; + ///////////////////////////////////////////////////// + // 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: + 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::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); diff --git a/tests/debug/Test_reweight_dwf_eofa_gparity.cc b/tests/debug/Test_reweight_dwf_eofa_gparity.cc new file mode 100644 index 00000000..b77ec33e --- /dev/null +++ b/tests/debug/Test_reweight_dwf_eofa_gparity.cc @@ -0,0 +1,209 @@ +/************************************************************************************* + +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; + +typedef typename GparityDomainWallFermionR::FermionField FermionField; + +// parameters for test +const std::vector grid_dim = { 8, 8, 8, 8 }; +const int Ls = 8; +const int Nhits = 10; +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 + GparityDomainWallFermionR::ImplParams params; + GparityDomainWallFermionR Ddwf_f(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, M5, params); + GparityDomainWallFermionR Ddwf_b(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mb, M5, params); + 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; + GparityDomainWallEOFAFermionR Deofa_L(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, mf, mb, shift_L, pm, M5, params); + GparityDomainWallEOFAFermionR Deofa_R(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mb, mf, mb, shift_R, pm, M5, params); + 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(); +}