From 3e139b52d3b84ef9b318d315ffdff2ce31104488 Mon Sep 17 00:00:00 2001 From: paboyle Date: Sun, 14 Jan 2018 22:47:24 +0000 Subject: [PATCH] Namespace --- .../OneFlavourEvenOddRationalRatio.h | 416 +++++++++--------- 1 file changed, 207 insertions(+), 209 deletions(-) diff --git a/lib/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h b/lib/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h index e1758e03..9aa27762 100644 --- a/lib/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h +++ b/lib/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h @@ -1,4 +1,4 @@ - /************************************************************************************* +/************************************************************************************* Grid physics library, www.github.com/paboyle/Grid @@ -23,259 +23,257 @@ 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_EVEN_ODD_RATIONAL_RATIO_H #define QCD_PSEUDOFERMION_ONE_FLAVOUR_EVEN_ODD_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 OneFlavourEvenOddRatioRationalPseudoFermionAction : 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 PhiEven; // the pseudo fermion field for this trajectory + FermionField PhiOdd; // the pseudo fermion field for this trajectory + +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); + + // MdagM^(+- 1/2) + std::cout< - class OneFlavourEvenOddRatioRationalPseudoFermionAction : 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 eta(NumOp.FermionGrid()); + FermionField etaOdd (NumOp.FermionRedBlackGrid()); + FermionField etaEven(NumOp.FermionRedBlackGrid()); + FermionField tmp(NumOp.FermionRedBlackGrid()); - MultiShiftFunction PowerHalf ; - MultiShiftFunction PowerNegHalf; - MultiShiftFunction PowerQuarter; - MultiShiftFunction PowerNegQuarter; + gaussian(pRNG,eta); eta=eta*scale; - 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 + pickCheckerboard(Even,etaEven,eta); + pickCheckerboard(Odd,etaOdd,eta); - 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); - - // 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); + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); - // MdagM^1/4 eta - SchurDifferentiableOperator MdagM(DenOp); - ConjugateGradientMultiShift msCG_M(param.MaxIter,PowerQuarter); - msCG_M(MdagM,etaOdd,tmp); + // 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); + // 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; + 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) { + ////////////////////////////////////////////////////// + // 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.FermionRedBlackGrid()); - FermionField Y(NumOp.FermionRedBlackGrid()); + 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); + // 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); + // MdagM^-1/4 VdagV^1/4 Phi + SchurDifferentiableOperator 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.FermionRedBlackGrid()); - std::vector MpvMfMpvPhi_k(n_pv,NumOp.FermionRedBlackGrid()); - std::vector MfMpvPhi_k (n_f ,NumOp.FermionRedBlackGrid()); + 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()); + FermionField MpvPhi(NumOp.FermionRedBlackGrid()); + FermionField MfMpvPhi(NumOp.FermionRedBlackGrid()); + FermionField MpvMfMpvPhi(NumOp.FermionRedBlackGrid()); + FermionField Y(NumOp.FermionRedBlackGrid()); - GaugeField tmp(NumOp.GaugeGrid()); + GaugeField tmp(NumOp.GaugeGrid()); - NumOp.ImportGauge(U); - DenOp.ImportGauge(U); + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); - SchurDifferentiableOperator VdagV(NumOp); - SchurDifferentiableOperator MdagM(DenOp); + SchurDifferentiableOperator VdagV(NumOp); + SchurDifferentiableOperator MdagM(DenOp); - 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,PhiOdd,MpvPhi_k,MpvPhi); - msCG_M(MdagM,MpvPhi,MfMpvPhi_k,MfMpvPhi); - msCG_V(VdagV,MfMpvPhi,MpvMfMpvPhi_k,MpvMfMpvPhi); + msCG_V(VdagV,PhiOdd,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