diff --git a/lib/qcd/action/pseudofermion/OneFlavourRational.h b/lib/qcd/action/pseudofermion/OneFlavourRational.h index 078aebf1..ca324e94 100644 --- a/lib/qcd/action/pseudofermion/OneFlavourRational.h +++ b/lib/qcd/action/pseudofermion/OneFlavourRational.h @@ -1,4 +1,4 @@ - /************************************************************************************* +/************************************************************************************* Grid physics library, www.github.com/paboyle/Grid @@ -23,191 +23,190 @@ 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_H #define QCD_PSEUDOFERMION_ONE_FLAVOUR_RATIONAL_H -namespace Grid{ - namespace QCD{ +NAMESPACE_BEGIN(Grid); - /////////////////////////////////////// - // One flavour rational - /////////////////////////////////////// +/////////////////////////////////////// +// One flavour rational +/////////////////////////////////////// - // S_f = chi^dag * N(M^dag*M)/D(M^dag*M) * chi - // - // Here, M is some operator - // N and D makeup the rat. poly - // +// S_f = chi^dag * N(M^dag*M)/D(M^dag*M) * chi +// +// Here, M is some operator +// N and D makeup the rat. poly +// - template - class OneFlavourRationalPseudoFermionAction : public Action { - public: - INHERIT_IMPL_TYPES(Impl); +template +class OneFlavourRationalPseudoFermionAction : public Action { +public: + INHERIT_IMPL_TYPES(Impl); - typedef OneFlavourRationalParams Params; - Params param; + typedef OneFlavourRationalParams Params; + Params param; - MultiShiftFunction PowerHalf ; - MultiShiftFunction PowerNegHalf; - MultiShiftFunction PowerQuarter; - MultiShiftFunction PowerNegQuarter; + MultiShiftFunction PowerHalf ; + MultiShiftFunction PowerNegHalf; + MultiShiftFunction PowerQuarter; + MultiShiftFunction PowerNegQuarter; - private: +private: - FermionOperator & FermOp;// the basic operator + FermionOperator & FermOp;// the basic operator - // NOT using "Nroots"; IroIro is -- perhaps later, but this wasn't good for us historically - // and hasenbusch works better + // NOT using "Nroots"; IroIro is -- perhaps later, but this wasn't good for us historically + // and hasenbusch works better - FermionField Phi; // the pseudo fermion field for this trajectory + FermionField Phi; // the pseudo fermion field for this trajectory - public: +public: - OneFlavourRationalPseudoFermionAction(FermionOperator &Op, - Params & p - ) : FermOp(Op), Phi(Op.FermionGrid()), param(p) - { - AlgRemez remez(param.lo,param.hi,param.precision); + OneFlavourRationalPseudoFermionAction(FermionOperator &Op, + Params & p + ) : FermOp(Op), Phi(Op.FermionGrid()), 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). + // P(phi) = e^{- phi^dag (MdagM)^-1/2 phi} + // = e^{- phi^dag (MdagM)^-1/4 (MdagM)^-1/4 phi} + // Phi = 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). - RealD scale = std::sqrt(0.5); + RealD scale = std::sqrt(0.5); - FermionField eta(FermOp.FermionGrid()); + FermionField eta(FermOp.FermionGrid()); - gaussian(pRNG,eta); + gaussian(pRNG,eta); - FermOp.ImportGauge(U); + FermOp.ImportGauge(U); - // mutishift CG - MdagMLinearOperator ,FermionField> MdagMOp(FermOp); - ConjugateGradientMultiShift msCG(param.MaxIter,PowerQuarter); - msCG(MdagMOp,eta,Phi); + // mutishift CG + MdagMLinearOperator ,FermionField> MdagMOp(FermOp); + ConjugateGradientMultiShift msCG(param.MaxIter,PowerQuarter); + msCG(MdagMOp,eta,Phi); - Phi=Phi*scale; + Phi=Phi*scale; - }; + }; - ////////////////////////////////////////////////////// - // S = phi^dag (Mdag M)^-1/2 phi - ////////////////////////////////////////////////////// - virtual RealD S(const GaugeField &U) { + ////////////////////////////////////////////////////// + // S = phi^dag (Mdag M)^-1/2 phi + ////////////////////////////////////////////////////// + virtual RealD S(const GaugeField &U) { - FermOp.ImportGauge(U); + FermOp.ImportGauge(U); - FermionField Y(FermOp.FermionGrid()); + FermionField Y(FermOp.FermionGrid()); - MdagMLinearOperator ,FermionField> MdagMOp(FermOp); + MdagMLinearOperator ,FermionField> MdagMOp(FermOp); - ConjugateGradientMultiShift msCG(param.MaxIter,PowerNegQuarter); + ConjugateGradientMultiShift msCG(param.MaxIter,PowerNegQuarter); - msCG(MdagMOp,Phi,Y); + msCG(MdagMOp,Phi,Y); - RealD action = norm2(Y); - std::cout << GridLogMessage << "Pseudofermion action FIXME -- is -1/4 solve or -1/2 solve faster??? "< MPhi_k (Npole,FermOp.FermionGrid()); + std::vector MPhi_k (Npole,FermOp.FermionGrid()); - FermionField X(FermOp.FermionGrid()); - FermionField Y(FermOp.FermionGrid()); + FermionField X(FermOp.FermionGrid()); + FermionField Y(FermOp.FermionGrid()); - GaugeField tmp(FermOp.GaugeGrid()); + GaugeField tmp(FermOp.GaugeGrid()); - FermOp.ImportGauge(U); + FermOp.ImportGauge(U); - MdagMLinearOperator ,FermionField> MdagMOp(FermOp); + MdagMLinearOperator ,FermionField> MdagMOp(FermOp); - ConjugateGradientMultiShift msCG(param.MaxIter,PowerNegHalf); + ConjugateGradientMultiShift msCG(param.MaxIter,PowerNegHalf); - msCG(MdagMOp,Phi,MPhi_k); + msCG(MdagMOp,Phi,MPhi_k); - dSdU = zero; - for(int k=0;k