From fe44fc50d9d52d5b7e0253448036985fbd4f6466 Mon Sep 17 00:00:00 2001 From: paboyle Date: Sun, 14 Jan 2018 22:45:29 +0000 Subject: [PATCH] Namespace --- .../pseudofermion/OneFlavourRationalRatio.h | 392 +++++++++--------- 1 file changed, 195 insertions(+), 197 deletions(-) diff --git a/lib/qcd/action/pseudofermion/OneFlavourRationalRatio.h b/lib/qcd/action/pseudofermion/OneFlavourRationalRatio.h index 00cfa2e3..7b7eeb19 100644 --- a/lib/qcd/action/pseudofermion/OneFlavourRationalRatio.h +++ b/lib/qcd/action/pseudofermion/OneFlavourRationalRatio.h @@ -1,4 +1,4 @@ - /************************************************************************************* +/************************************************************************************* Grid physics library, www.github.com/paboyle/Grid @@ -23,245 +23,243 @@ Author: Peter Boyle 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 */ +*************************************************************************************/ +/* END LEGAL */ #ifndef QCD_PSEUDOFERMION_ONE_FLAVOUR_RATIONAL_RATIO_H #define QCD_PSEUDOFERMION_ONE_FLAVOUR_RATIONAL_RATIO_H -namespace Grid{ - namespace QCD{ +NAMESPACE_BEGIN(Grid); + +/////////////////////////////////////// +// One flavour rational +/////////////////////////////////////// - /////////////////////////////////////// - // One flavour rational - /////////////////////////////////////// +// 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 P/Q \sim R_{1/4} ~ (V^dagV)^{1/4} +// Here N/D \sim R_{-1/2} ~ (M^dagM)^{-1/2} + +template +class OneFlavourRatioRationalPseudoFermionAction : public Action { +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 Phi; // the pseudo fermion field for this trajectory + +public: + + OneFlavourRatioRationalPseudoFermionAction(FermionOperator &_NumOp, + FermionOperator &_DenOp, + Params & p + ) : NumOp(_NumOp), DenOp(_DenOp), Phi(_NumOp.FermionGrid()), param(p) + { + AlgRemez remez(param.lo,param.hi,param.precision); + + // MdagM^(+- 1/2) + std::cout< - class OneFlavourRatioRationalPseudoFermionAction : public Action { - public: + // P(phi) = e^{- phi^dag (VdagV)^1/4 (MdagM)^-1/2 (VdagV)^1/4 phi} + // = e^{- phi^dag (VdagV)^1/4 (MdagM)^-1/4 (MdagM)^-1/4 (VdagV)^1/4 phi} + // + // Phi = (VdagV)^-1/4 Mdag^{1/4} eta + // + // P(eta) = e^{- eta^dag eta} + // + // e^{x^2/2 sig^2} => sig^2 = 0.5. + // + // So eta should be of width sig = 1/sqrt(2). - INHERIT_IMPL_TYPES(Impl); + RealD scale = std::sqrt(0.5); - typedef OneFlavourRationalParams Params; - Params param; + FermionField tmp(NumOp.FermionGrid()); + FermionField eta(NumOp.FermionGrid()); - MultiShiftFunction PowerHalf ; - MultiShiftFunction PowerNegHalf; - MultiShiftFunction PowerQuarter; - MultiShiftFunction PowerNegQuarter; + gaussian(pRNG,eta); - private: - - FermionOperator & NumOp;// the basic operator - FermionOperator & DenOp;// the basic operator - FermionField Phi; // the pseudo fermion field for this trajectory + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); - public: + // MdagM^1/4 eta + MdagMLinearOperator ,FermionField> MdagM(DenOp); + ConjugateGradientMultiShift msCG_M(param.MaxIter,PowerQuarter); + msCG_M(MdagM,eta,tmp); - OneFlavourRatioRationalPseudoFermionAction(FermionOperator &_NumOp, - FermionOperator &_DenOp, - Params & p - ) : NumOp(_NumOp), DenOp(_DenOp), Phi(_NumOp.FermionGrid()), param(p) - { - AlgRemez remez(param.lo,param.hi,param.precision); + // VdagV^-1/4 MdagM^1/4 eta + MdagMLinearOperator ,FermionField> VdagV(NumOp); + ConjugateGradientMultiShift msCG_V(param.MaxIter,PowerNegQuarter); + msCG_V(VdagV,tmp,Phi); - // 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 tmp(NumOp.FermionGrid()); - FermionField eta(NumOp.FermionGrid()); - - gaussian(pRNG,eta); - - NumOp.ImportGauge(U); - DenOp.ImportGauge(U); - - // MdagM^1/4 eta - MdagMLinearOperator ,FermionField> MdagM(DenOp); - ConjugateGradientMultiShift msCG_M(param.MaxIter,PowerQuarter); - msCG_M(MdagM,eta,tmp); - - // VdagV^-1/4 MdagM^1/4 eta - MdagMLinearOperator ,FermionField> VdagV(NumOp); - ConjugateGradientMultiShift msCG_V(param.MaxIter,PowerNegQuarter); - msCG_V(VdagV,tmp,Phi); - - Phi=Phi*scale; + Phi=Phi*scale; - }; + }; - ////////////////////////////////////////////////////// - // 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) { + ////////////////////////////////////////////////////// + // 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); + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); - FermionField X(NumOp.FermionGrid()); - FermionField Y(NumOp.FermionGrid()); + FermionField X(NumOp.FermionGrid()); + FermionField Y(NumOp.FermionGrid()); - // VdagV^1/4 Phi - MdagMLinearOperator ,FermionField> VdagV(NumOp); - ConjugateGradientMultiShift msCG_V(param.MaxIter,PowerQuarter); - msCG_V(VdagV,Phi,X); + // VdagV^1/4 Phi + MdagMLinearOperator ,FermionField> VdagV(NumOp); + ConjugateGradientMultiShift msCG_V(param.MaxIter,PowerQuarter); + msCG_V(VdagV,Phi,X); - // MdagM^-1/4 VdagV^1/4 Phi - MdagMLinearOperator ,FermionField> MdagM(DenOp); - ConjugateGradientMultiShift msCG_M(param.MaxIter,PowerNegQuarter); - msCG_M(MdagM,X,Y); + // MdagM^-1/4 VdagV^1/4 Phi + MdagMLinearOperator ,FermionField> MdagM(DenOp); + ConjugateGradientMultiShift msCG_M(param.MaxIter,PowerNegQuarter); + msCG_M(MdagM,X,Y); - // Phidag VdagV^1/4 MdagM^-1/4 MdagM^-1/4 VdagV^1/4 Phi - RealD action = norm2(Y); + // Phidag VdagV^1/4 MdagM^-1/4 MdagM^-1/4 VdagV^1/4 Phi + RealD action = norm2(Y); - return action; - }; + 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 - // + // 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) { + virtual void deriv(const GaugeField &U,GaugeField & dSdU) { - const int n_f = PowerNegHalf.poles.size(); - const int n_pv = PowerQuarter.poles.size(); + const int n_f = PowerNegHalf.poles.size(); + const int n_pv = PowerQuarter.poles.size(); - std::vector MpvPhi_k (n_pv,NumOp.FermionGrid()); - std::vector MpvMfMpvPhi_k(n_pv,NumOp.FermionGrid()); - std::vector MfMpvPhi_k (n_f,NumOp.FermionGrid()); + std::vector MpvPhi_k (n_pv,NumOp.FermionGrid()); + std::vector MpvMfMpvPhi_k(n_pv,NumOp.FermionGrid()); + std::vector MfMpvPhi_k (n_f,NumOp.FermionGrid()); - FermionField MpvPhi(NumOp.FermionGrid()); - FermionField MfMpvPhi(NumOp.FermionGrid()); - FermionField MpvMfMpvPhi(NumOp.FermionGrid()); - FermionField Y(NumOp.FermionGrid()); + FermionField MpvPhi(NumOp.FermionGrid()); + FermionField MfMpvPhi(NumOp.FermionGrid()); + FermionField MpvMfMpvPhi(NumOp.FermionGrid()); + FermionField Y(NumOp.FermionGrid()); - GaugeField tmp(NumOp.GaugeGrid()); + GaugeField tmp(NumOp.GaugeGrid()); - NumOp.ImportGauge(U); - DenOp.ImportGauge(U); + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); - MdagMLinearOperator ,FermionField> MdagM(DenOp); - MdagMLinearOperator ,FermionField> VdagV(NumOp); + MdagMLinearOperator ,FermionField> MdagM(DenOp); + MdagMLinearOperator ,FermionField> VdagV(NumOp); - ConjugateGradientMultiShift msCG_V(param.MaxIter,PowerQuarter); - ConjugateGradientMultiShift msCG_M(param.MaxIter,PowerNegHalf); + ConjugateGradientMultiShift msCG_V(param.MaxIter,PowerQuarter); + ConjugateGradientMultiShift msCG_M(param.MaxIter,PowerNegHalf); - msCG_V(VdagV,Phi,MpvPhi_k,MpvPhi); - msCG_M(MdagM,MpvPhi,MfMpvPhi_k,MfMpvPhi); - msCG_V(VdagV,MfMpvPhi,MpvMfMpvPhi_k,MpvMfMpvPhi); + msCG_V(VdagV,Phi,MpvPhi_k,MpvPhi); + msCG_M(MdagM,MpvPhi,MfMpvPhi_k,MfMpvPhi); + msCG_V(VdagV,MfMpvPhi,MpvMfMpvPhi_k,MpvMfMpvPhi); - RealD ak; + RealD ak; - dSdU = zero; + 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) + // 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