diff --git a/lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h b/lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h index 0bbc0ae6..2c66be78 100644 --- a/lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h +++ b/lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h @@ -26,164 +26,163 @@ with this program; if not, write to the Free Software Foundation, Inc., See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ -/* END LEGAL */ + /* END LEGAL */ #ifndef QCD_PSEUDOFERMION_TWO_FLAVOUR_EVEN_ODD_H #define QCD_PSEUDOFERMION_TWO_FLAVOUR_EVEN_ODD_H -namespace Grid { - namespace QCD { - - //////////////////////////////////////////////////////////////////////// - // Two flavour pseudofermion action for any EO prec dop - //////////////////////////////////////////////////////////////////////// - template - class TwoFlavourEvenOddPseudoFermionAction - : public Action { - public: - INHERIT_IMPL_TYPES(Impl); - - private: - FermionOperator &FermOp; // the basic operator - - OperatorFunction &DerivativeSolver; - OperatorFunction &ActionSolver; - - FermionField PhiOdd; // the pseudo fermion field for this trajectory - FermionField PhiEven; // the pseudo fermion field for this trajectory - - public: - ///////////////////////////////////////////////// - // Pass in required objects. - ///////////////////////////////////////////////// - TwoFlavourEvenOddPseudoFermionAction(FermionOperator &Op, - OperatorFunction &DS, - OperatorFunction &AS) - : FermOp(Op), - DerivativeSolver(DS), - ActionSolver(AS), - PhiEven(Op.FermionRedBlackGrid()), - PhiOdd(Op.FermionRedBlackGrid()) - {}; +NAMESPACE_BEGIN(Grid); - virtual std::string action_name(){return "TwoFlavourEvenOddPseudoFermionAction";} + +//////////////////////////////////////////////////////////////////////// +// Two flavour pseudofermion action for any EO prec dop +//////////////////////////////////////////////////////////////////////// +template +class TwoFlavourEvenOddPseudoFermionAction + : public Action { +public: + INHERIT_IMPL_TYPES(Impl); + +private: + FermionOperator &FermOp; // the basic operator + + OperatorFunction &DerivativeSolver; + OperatorFunction &ActionSolver; + + FermionField PhiOdd; // the pseudo fermion field for this trajectory + FermionField PhiEven; // the pseudo fermion field for this trajectory + +public: + ///////////////////////////////////////////////// + // Pass in required objects. + ///////////////////////////////////////////////// + TwoFlavourEvenOddPseudoFermionAction(FermionOperator &Op, + OperatorFunction &DS, + OperatorFunction &AS) + : FermOp(Op), + DerivativeSolver(DS), + ActionSolver(AS), + PhiEven(Op.FermionRedBlackGrid()), + PhiOdd(Op.FermionRedBlackGrid()) + {}; + + virtual std::string action_name(){return "TwoFlavourEvenOddPseudoFermionAction";} - virtual std::string LogParameters(){ - std::stringstream sstream; - sstream << GridLogMessage << "["< sig^2 = 0.5. + // P(phi) = e^{- phi^dag (MpcdagMpc)^-1 phi} + // Phi = McpDag eta + // P(eta) = e^{- eta^dag eta} + // + // e^{x^2/2 sig^2} => sig^2 = 0.5. - RealD scale = std::sqrt(0.5); + RealD scale = std::sqrt(0.5); - FermionField eta (FermOp.FermionGrid()); - FermionField etaOdd (FermOp.FermionRedBlackGrid()); - FermionField etaEven(FermOp.FermionRedBlackGrid()); + FermionField eta (FermOp.FermionGrid()); + FermionField etaOdd (FermOp.FermionRedBlackGrid()); + FermionField etaEven(FermOp.FermionRedBlackGrid()); - gaussian(pRNG,eta); - pickCheckerboard(Even,etaEven,eta); - pickCheckerboard(Odd,etaOdd,eta); + gaussian(pRNG,eta); + pickCheckerboard(Even,etaEven,eta); + pickCheckerboard(Odd,etaOdd,eta); - FermOp.ImportGauge(U); - SchurDifferentiableOperator PCop(FermOp); + FermOp.ImportGauge(U); + SchurDifferentiableOperator PCop(FermOp); - PCop.MpcDag(etaOdd,PhiOdd); + PCop.MpcDag(etaOdd,PhiOdd); - FermOp.MooeeDag(etaEven,PhiEven); + FermOp.MooeeDag(etaEven,PhiEven); - PhiOdd =PhiOdd*scale; - PhiEven=PhiEven*scale; + PhiOdd =PhiOdd*scale; + PhiEven=PhiEven*scale; - }; + }; - ////////////////////////////////////////////////////// - // S = phi^dag (Mdag M)^-1 phi (odd) - // + phi^dag (Mdag M)^-1 phi (even) - ////////////////////////////////////////////////////// - virtual RealD S(const GaugeField &U) { + ////////////////////////////////////////////////////// + // S = phi^dag (Mdag M)^-1 phi (odd) + // + phi^dag (Mdag M)^-1 phi (even) + ////////////////////////////////////////////////////// + virtual RealD S(const GaugeField &U) { - FermOp.ImportGauge(U); + FermOp.ImportGauge(U); - FermionField X(FermOp.FermionRedBlackGrid()); - FermionField Y(FermOp.FermionRedBlackGrid()); + FermionField X(FermOp.FermionRedBlackGrid()); + FermionField Y(FermOp.FermionRedBlackGrid()); - SchurDifferentiableOperator PCop(FermOp); + SchurDifferentiableOperator PCop(FermOp); - X=zero; - ActionSolver(PCop,PhiOdd,X); - PCop.Op(X,Y); - RealD action = norm2(Y); + X=zero; + ActionSolver(PCop,PhiOdd,X); + PCop.Op(X,Y); + RealD action = norm2(Y); - // The EE factorised block; normally can replace with zero if det is constant (gauge field indept) - // Only really clover term that creates this. - FermOp.MooeeInvDag(PhiEven,Y); - action = action + norm2(Y); + // The EE factorised block; normally can replace with zero if det is constant (gauge field indept) + // Only really clover term that creates this. + FermOp.MooeeInvDag(PhiEven,Y); + action = action + norm2(Y); - std::cout << GridLogMessage << "Pseudofermion EO action "< Mpc(FermOp); + SchurDifferentiableOperator Mpc(FermOp); - // Our conventions really make this UdSdU; We do not differentiate wrt Udag here. - // So must take dSdU - adj(dSdU) and left multiply by mom to get dS/dt. + // Our conventions really make this UdSdU; We do not differentiate wrt Udag here. + // So must take dSdU - adj(dSdU) and left multiply by mom to get dS/dt. - X=zero; - DerivativeSolver(Mpc,PhiOdd,X); - Mpc.Mpc(X,Y); - Mpc.MpcDeriv(tmp , Y, X ); dSdU=tmp; - Mpc.MpcDagDeriv(tmp , X, Y); dSdU=dSdU+tmp; + X=zero; + DerivativeSolver(Mpc,PhiOdd,X); + Mpc.Mpc(X,Y); + Mpc.MpcDeriv(tmp , Y, X ); dSdU=tmp; + Mpc.MpcDagDeriv(tmp , X, Y); dSdU=dSdU+tmp; - // Treat the EE case. (MdagM)^-1 = Minv Minvdag - // Deriv defaults to zero. - // FermOp.MooeeInvDag(PhiOdd,Y); - // FermOp.MooeeInv(Y,X); - // FermOp.MeeDeriv(tmp , Y, X,DaggerNo ); dSdU=tmp; - // FermOp.MeeDeriv(tmp , X, Y,DaggerYes); dSdU=dSdU+tmp; + // Treat the EE case. (MdagM)^-1 = Minv Minvdag + // Deriv defaults to zero. + // FermOp.MooeeInvDag(PhiOdd,Y); + // FermOp.MooeeInv(Y,X); + // FermOp.MeeDeriv(tmp , Y, X,DaggerNo ); dSdU=tmp; + // FermOp.MeeDeriv(tmp , X, Y,DaggerYes); dSdU=dSdU+tmp; - assert(FermOp.ConstEE() == 1); + assert(FermOp.ConstEE() == 1); - /* - FermOp.MooeeInvDag(PhiOdd,Y); - FermOp.MooeeInv(Y,X); - FermOp.MeeDeriv(tmp , Y, X,DaggerNo ); dSdU=tmp; - FermOp.MeeDeriv(tmp , X, Y,DaggerYes); dSdU=dSdU+tmp; - */ + /* + FermOp.MooeeInvDag(PhiOdd,Y); + FermOp.MooeeInv(Y,X); + FermOp.MeeDeriv(tmp , Y, X,DaggerNo ); dSdU=tmp; + FermOp.MeeDeriv(tmp , X, Y,DaggerYes); dSdU=dSdU+tmp; + */ - //dSdU = Ta(dSdU); + //dSdU = Ta(dSdU); - }; + }; - }; +}; - } -} +NAMESPACE_END(Grid); #endif