diff --git a/lib/qcd/action/pseudofermion/TwoFlavourRatio.h b/lib/qcd/action/pseudofermion/TwoFlavourRatio.h index bcbf9364..7897bb12 100644 --- a/lib/qcd/action/pseudofermion/TwoFlavourRatio.h +++ b/lib/qcd/action/pseudofermion/TwoFlavourRatio.h @@ -1,4 +1,4 @@ - /************************************************************************************* +/************************************************************************************* Grid physics library, www.github.com/paboyle/Grid @@ -25,149 +25,149 @@ Author: paboyle 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_TWO_FLAVOUR_RATIO_H #define QCD_PSEUDOFERMION_TWO_FLAVOUR_RATIO_H -namespace Grid{ - namespace QCD{ +NAMESPACE_BEGIN(Grid); - /////////////////////////////////////// - // Two flavour ratio - /////////////////////////////////////// - template - class TwoFlavourRatioPseudoFermionAction : public Action { - public: - INHERIT_IMPL_TYPES(Impl); +/////////////////////////////////////// +// Two flavour ratio +/////////////////////////////////////// +template +class TwoFlavourRatioPseudoFermionAction : public Action { +public: + INHERIT_IMPL_TYPES(Impl); - private: - FermionOperator & NumOp;// the basic operator - FermionOperator & DenOp;// the basic operator +private: + FermionOperator & NumOp;// the basic operator + FermionOperator & DenOp;// the basic operator - OperatorFunction &DerivativeSolver; - OperatorFunction &ActionSolver; + OperatorFunction &DerivativeSolver; + OperatorFunction &ActionSolver; - FermionField Phi; // the pseudo fermion field for this trajectory + FermionField Phi; // the pseudo fermion field for this trajectory - public: - TwoFlavourRatioPseudoFermionAction(FermionOperator &_NumOp, - FermionOperator &_DenOp, - OperatorFunction & DS, - OperatorFunction & AS - ) : NumOp(_NumOp), DenOp(_DenOp), DerivativeSolver(DS), ActionSolver(AS), Phi(_NumOp.FermionGrid()) {}; +public: + TwoFlavourRatioPseudoFermionAction(FermionOperator &_NumOp, + FermionOperator &_DenOp, + OperatorFunction & DS, + OperatorFunction & AS + ) : NumOp(_NumOp), DenOp(_DenOp), DerivativeSolver(DS), ActionSolver(AS), Phi(_NumOp.FermionGrid()) {}; - virtual std::string action_name(){return "TwoFlavourRatioPseudoFermionAction";} + virtual std::string action_name(){return "TwoFlavourRatioPseudoFermionAction";} - virtual std::string LogParameters(){ - std::stringstream sstream; - sstream << GridLogMessage << "["< sig^2 = 0.5. - // - // So eta should be of width sig = 1/sqrt(2) and must multiply by 0.707.... - // - RealD scale = std::sqrt(0.5); + // P(phi) = e^{- phi^dag V (MdagM)^-1 Vdag phi} + // + // NumOp == V + // DenOp == M + // + // Take phi = Vdag^{-1} Mdag eta ; eta = Mdag^{-1} Vdag Phi + // + // 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) and must multiply by 0.707.... + // + RealD scale = std::sqrt(0.5); - FermionField eta(NumOp.FermionGrid()); - FermionField tmp(NumOp.FermionGrid()); + FermionField eta(NumOp.FermionGrid()); + FermionField tmp(NumOp.FermionGrid()); - gaussian(pRNG,eta); + gaussian(pRNG,eta); - NumOp.ImportGauge(U); - DenOp.ImportGauge(U); + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); - // Note: this hard codes normal equations type solvers; alternate implementation needed for - // non-herm style solvers. - MdagMLinearOperator ,FermionField> MdagMOp(NumOp); + // Note: this hard codes normal equations type solvers; alternate implementation needed for + // non-herm style solvers. + MdagMLinearOperator ,FermionField> MdagMOp(NumOp); - DenOp.Mdag(eta,Phi); // Mdag eta - tmp = zero; - ActionSolver(MdagMOp,Phi,tmp); // (VdagV)^-1 Mdag eta = V^-1 Vdag^-1 Mdag eta - NumOp.M(tmp,Phi); // Vdag^-1 Mdag eta + DenOp.Mdag(eta,Phi); // Mdag eta + tmp = zero; + ActionSolver(MdagMOp,Phi,tmp); // (VdagV)^-1 Mdag eta = V^-1 Vdag^-1 Mdag eta + NumOp.M(tmp,Phi); // Vdag^-1 Mdag eta - Phi=Phi*scale; + Phi=Phi*scale; - }; + }; - ////////////////////////////////////////////////////// - // S = phi^dag V (Mdag M)^-1 Vdag phi - ////////////////////////////////////////////////////// - virtual RealD S(const GaugeField &U) { + ////////////////////////////////////////////////////// + // S = phi^dag V (Mdag M)^-1 Vdag phi + ////////////////////////////////////////////////////// + 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()); - MdagMLinearOperator ,FermionField> MdagMOp(DenOp); + MdagMLinearOperator ,FermionField> MdagMOp(DenOp); - NumOp.Mdag(Phi,Y); // Y= Vdag phi - X=zero; - ActionSolver(MdagMOp,Y,X); // X= (MdagM)^-1 Vdag phi - DenOp.M(X,Y); // Y= Mdag^-1 Vdag phi + NumOp.Mdag(Phi,Y); // Y= Vdag phi + X=zero; + ActionSolver(MdagMOp,Y,X); // X= (MdagM)^-1 Vdag phi + DenOp.M(X,Y); // Y= Mdag^-1 Vdag phi - RealD action = norm2(Y); + RealD action = norm2(Y); - return action; - }; + return action; + }; - ////////////////////////////////////////////////////// - // dS/du = phi^dag dV (Mdag M)^-1 V^dag phi - // - phi^dag V (Mdag M)^-1 [ Mdag dM + dMdag M ] (Mdag M)^-1 V^dag phi - // + phi^dag V (Mdag M)^-1 dV^dag phi - ////////////////////////////////////////////////////// - virtual void deriv(const GaugeField &U,GaugeField & dSdU) { + ////////////////////////////////////////////////////// + // dS/du = phi^dag dV (Mdag M)^-1 V^dag phi + // - phi^dag V (Mdag M)^-1 [ Mdag dM + dMdag M ] (Mdag M)^-1 V^dag phi + // + phi^dag V (Mdag M)^-1 dV^dag phi + ////////////////////////////////////////////////////// + virtual void deriv(const GaugeField &U,GaugeField & dSdU) { - NumOp.ImportGauge(U); - DenOp.ImportGauge(U); + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); - MdagMLinearOperator ,FermionField> MdagMOp(DenOp); + MdagMLinearOperator ,FermionField> MdagMOp(DenOp); - FermionField X(NumOp.FermionGrid()); - FermionField Y(NumOp.FermionGrid()); + FermionField X(NumOp.FermionGrid()); + FermionField Y(NumOp.FermionGrid()); - GaugeField force(NumOp.GaugeGrid()); + GaugeField force(NumOp.GaugeGrid()); - //Y=Vdag phi - //X = (Mdag M)^-1 V^dag phi - //Y = (Mdag)^-1 V^dag phi - NumOp.Mdag(Phi,Y); // Y= Vdag phi - X=zero; - DerivativeSolver(MdagMOp,Y,X); // X= (MdagM)^-1 Vdag phi - DenOp.M(X,Y); // Y= Mdag^-1 Vdag phi + //Y=Vdag phi + //X = (Mdag M)^-1 V^dag phi + //Y = (Mdag)^-1 V^dag phi + NumOp.Mdag(Phi,Y); // Y= Vdag phi + X=zero; + DerivativeSolver(MdagMOp,Y,X); // X= (MdagM)^-1 Vdag phi + DenOp.M(X,Y); // Y= Mdag^-1 Vdag phi - // phi^dag V (Mdag M)^-1 dV^dag phi - NumOp.MDeriv(force , X, Phi, DaggerYes ); dSdU=force; + // phi^dag V (Mdag M)^-1 dV^dag phi + NumOp.MDeriv(force , X, Phi, DaggerYes ); dSdU=force; - // phi^dag dV (Mdag M)^-1 V^dag phi - NumOp.MDeriv(force , Phi, X ,DaggerNo ); dSdU=dSdU+force; + // phi^dag dV (Mdag M)^-1 V^dag phi + NumOp.MDeriv(force , Phi, X ,DaggerNo ); dSdU=dSdU+force; - // - phi^dag V (Mdag M)^-1 Mdag dM (Mdag M)^-1 V^dag phi - // - phi^dag V (Mdag M)^-1 dMdag M (Mdag M)^-1 V^dag phi - DenOp.MDeriv(force,Y,X,DaggerNo); dSdU=dSdU-force; - DenOp.MDeriv(force,X,Y,DaggerYes); dSdU=dSdU-force; + // - phi^dag V (Mdag M)^-1 Mdag dM (Mdag M)^-1 V^dag phi + // - phi^dag V (Mdag M)^-1 dMdag M (Mdag M)^-1 V^dag phi + DenOp.MDeriv(force,Y,X,DaggerNo); dSdU=dSdU-force; + DenOp.MDeriv(force,X,Y,DaggerYes); dSdU=dSdU-force; - dSdU *= -1.0; - //dSdU = - Ta(dSdU); + dSdU *= -1.0; + //dSdU = - Ta(dSdU); + + }; +}; + +NAMESPACE_END(Grid): - }; - }; - } -} #endif