From d185f2eaa7c55811db517caa8c783cfee04364e5 Mon Sep 17 00:00:00 2001 From: Christopher Kelly Date: Wed, 23 Dec 2020 16:26:10 -0500 Subject: [PATCH] OneFlavourEvenOddRatioRationalPseudoFermionAction now derives from GeneralEvenOddRatioRationalPseudoFermionAction, simply performs transcription of parameters --- .../GeneralEvenOddRationalRatio.h | 2 +- .../OneFlavourEvenOddRationalRatio.h | 251 ++---------------- Grid/qcd/action/pseudofermion/PseudoFermion.h | 2 +- 3 files changed, 20 insertions(+), 235 deletions(-) diff --git a/Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatio.h b/Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatio.h index af64fba7..a92b6f7d 100644 --- a/Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatio.h +++ b/Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatio.h @@ -83,7 +83,7 @@ NAMESPACE_BEGIN(Grid); GeneralEvenOddRatioRationalPseudoFermionAction(FermionOperator &_NumOp, FermionOperator &_DenOp, - Params & p + const Params & p ) : NumOp(_NumOp), DenOp(_DenOp), diff --git a/Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h b/Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h index e5f0b602..1b36ae0f 100644 --- a/Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h +++ b/Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h @@ -40,246 +40,31 @@ NAMESPACE_BEGIN(Grid); // Here N/D \sim R_{-1/2} ~ (M^dagM)^{-1/2} template - class OneFlavourEvenOddRatioRationalPseudoFermionAction : public Action { + class OneFlavourEvenOddRatioRationalPseudoFermionAction : public GeneralEvenOddRatioRationalPseudoFermionAction { public: - - INHERIT_IMPL_TYPES(Impl); - typedef OneFlavourRationalParams Params; - Params param; - - MultiShiftFunction PowerHalf ; - MultiShiftFunction PowerNegHalf; - MultiShiftFunction PowerQuarter; - MultiShiftFunction PowerNegQuarter; - private: - - FermionOperator & NumOp;// the basic operator - FermionOperator & DenOp;// the basic operator - FermionField PhiEven; // the pseudo fermion field for this trajectory - FermionField PhiOdd; // the pseudo fermion field for this trajectory + static RationalActionParams transcribe(const Params &in){ + RationalActionParams out; + out.inv_pow = 2; + out.lo = in.lo; + out.hi = in.hi; + out.MaxIter = in.MaxIter; + out.action_tolerance = out.md_tolerance = in.tolerance; + out.action_degree = out.md_degree = in.degree; + out.precision = in.precision; + out.BoundsCheckFreq = in.BoundsCheckFreq; + return out; + } public: - OneFlavourEvenOddRatioRationalPseudoFermionAction(FermionOperator &_NumOp, - FermionOperator &_DenOp, - Params & p - ) : - NumOp(_NumOp), - DenOp(_DenOp), - PhiOdd (_NumOp.FermionRedBlackGrid()), - PhiEven(_NumOp.FermionRedBlackGrid()), - param(p) - { - AlgRemez remez(param.lo,param.hi,param.precision); + FermionOperator &_DenOp, + const Params & p + ) : + GeneralEvenOddRatioRationalPseudoFermionAction(_NumOp, _DenOp, transcribe(p)){} - // MdagM^(+- 1/2) - std::cout< sig^2 = 0.5. - // - // So eta should be of width sig = 1/sqrt(2). - - RealD scale = std::sqrt(0.5); - - FermionField eta(NumOp.FermionGrid()); - FermionField etaOdd (NumOp.FermionRedBlackGrid()); - FermionField etaEven(NumOp.FermionRedBlackGrid()); - FermionField tmp(NumOp.FermionRedBlackGrid()); - - gaussian(pRNG,eta); eta=eta*scale; - - pickCheckerboard(Even,etaEven,eta); - pickCheckerboard(Odd,etaOdd,eta); - - NumOp.ImportGauge(U); - DenOp.ImportGauge(U); - - - // MdagM^1/4 eta - SchurDifferentiableOperator MdagM(DenOp); - ConjugateGradientMultiShift msCG_M(param.MaxIter,PowerQuarter); - msCG_M(MdagM,etaOdd,tmp); - - // VdagV^-1/4 MdagM^1/4 eta - SchurDifferentiableOperator VdagV(NumOp); - ConjugateGradientMultiShift msCG_V(param.MaxIter,PowerNegQuarter); - msCG_V(VdagV,tmp,PhiOdd); - - assert(NumOp.ConstEE() == 1); - assert(DenOp.ConstEE() == 1); - PhiEven = Zero(); - - }; - - ////////////////////////////////////////////////////// - // S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi - ////////////////////////////////////////////////////// - virtual RealD S(const GaugeField &U) { - - NumOp.ImportGauge(U); - DenOp.ImportGauge(U); - - FermionField X(NumOp.FermionRedBlackGrid()); - FermionField Y(NumOp.FermionRedBlackGrid()); - - // VdagV^1/4 Phi - SchurDifferentiableOperator VdagV(NumOp); - ConjugateGradientMultiShift msCG_V(param.MaxIter,PowerQuarter); - msCG_V(VdagV,PhiOdd,X); - - // MdagM^-1/4 VdagV^1/4 Phi - SchurDifferentiableOperator MdagM(DenOp); - ConjugateGradientMultiShift msCG_M(param.MaxIter,PowerNegQuarter); - msCG_M(MdagM,X,Y); - - // Randomly apply rational bounds checks. - if ( (rand()%param.BoundsCheckFreq)==0 ) { - FermionField gauss(NumOp.FermionRedBlackGrid()); - gauss = PhiOdd; - HighBoundCheck(MdagM,gauss,param.hi); - InverseSqrtBoundsCheck(param.MaxIter,param.tolerance*100,MdagM,gauss,PowerNegHalf); - } - - // Phidag VdagV^1/4 MdagM^-1/4 MdagM^-1/4 VdagV^1/4 Phi - RealD action = norm2(Y); - - return action; - }; - - // S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi - // - // Here, M is some 5D operator and V is the Pauli-Villars field - // N and D makeup the rat. poly of the M term and P and & makeup the rat.poly of the denom term - // - // Need - // dS_f/dU = chi^dag d[P/Q] N/D P/Q chi - // + chi^dag P/Q d[N/D] P/Q chi - // + chi^dag P/Q N/D d[P/Q] chi - // - // P/Q is expressed as partial fraction expansion: - // - // a0 + \sum_k ak/(V^dagV + bk) - // - // d[P/Q] is then - // - // \sum_k -ak [V^dagV+bk]^{-1} [ dV^dag V + V^dag dV ] [V^dag V + bk]^{-1} - // - // and similar for N/D. - // - // Need - // MpvPhi_k = [Vdag V + bk]^{-1} chi - // MpvPhi = {a0 + \sum_k ak [Vdag V + bk]^{-1} }chi - // - // MfMpvPhi_k = [MdagM+bk]^{-1} MpvPhi - // MfMpvPhi = {a0 + \sum_k ak [Mdag M + bk]^{-1} } MpvPhi - // - // MpvMfMpvPhi_k = [Vdag V + bk]^{-1} MfMpvchi - // - - virtual void deriv(const GaugeField &U,GaugeField & dSdU) { - - const int n_f = PowerNegHalf.poles.size(); - const int n_pv = PowerQuarter.poles.size(); - - std::vector MpvPhi_k (n_pv,NumOp.FermionRedBlackGrid()); - std::vector MpvMfMpvPhi_k(n_pv,NumOp.FermionRedBlackGrid()); - std::vector MfMpvPhi_k (n_f ,NumOp.FermionRedBlackGrid()); - - FermionField MpvPhi(NumOp.FermionRedBlackGrid()); - FermionField MfMpvPhi(NumOp.FermionRedBlackGrid()); - FermionField MpvMfMpvPhi(NumOp.FermionRedBlackGrid()); - FermionField Y(NumOp.FermionRedBlackGrid()); - - GaugeField tmp(NumOp.GaugeGrid()); - - NumOp.ImportGauge(U); - DenOp.ImportGauge(U); - - SchurDifferentiableOperator VdagV(NumOp); - SchurDifferentiableOperator MdagM(DenOp); - - ConjugateGradientMultiShift msCG_V(param.MaxIter,PowerQuarter); - ConjugateGradientMultiShift msCG_M(param.MaxIter,PowerNegHalf); - - msCG_V(VdagV,PhiOdd,MpvPhi_k,MpvPhi); - msCG_M(MdagM,MpvPhi,MfMpvPhi_k,MfMpvPhi); - msCG_V(VdagV,MfMpvPhi,MpvMfMpvPhi_k,MpvMfMpvPhi); - - RealD ak; - - dSdU = Zero(); - - // With these building blocks - // - // dS/dU = - // \sum_k -ak MfMpvPhi_k^dag [ dM^dag M + M^dag dM ] MfMpvPhi_k (1) - // + \sum_k -ak MpvMfMpvPhi_k^\dag [ dV^dag V + V^dag dV ] MpvPhi_k (2) - // -ak MpvPhi_k^dag [ dV^dag V + V^dag dV ] MpvMfMpvPhi_k (3) - - //(1) - for(int k=0;k #include #include -#include #include +#include #include #endif