portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-09 21:50:45 +01:00
Code Releases Activity

Namespace

Browse Source
This commit is contained in:
paboyle 2018-01-14 22:42:31 +00:00
parent 8d52e0a349
commit f5e74033f9
1 changed files with 108 additions and 108 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

216
lib/qcd/action/pseudofermion/TwoFlavourRatio.h
View File

@ -1,4 +1,4 @@
/************************************************************************************* /*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
@ -25,149 +25,149 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory 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 #ifndef QCD_PSEUDOFERMION_TWO_FLAVOUR_RATIO_H
#define QCD_PSEUDOFERMION_TWO_FLAVOUR_RATIO_H #define QCD_PSEUDOFERMION_TWO_FLAVOUR_RATIO_H
namespace Grid{ NAMESPACE_BEGIN(Grid);
namespace QCD{
/////////////////////////////////////// ///////////////////////////////////////
// Two flavour ratio // Two flavour ratio
/////////////////////////////////////// ///////////////////////////////////////
template<class Impl> template<class Impl>
class TwoFlavourRatioPseudoFermionAction : public Action<typename Impl::GaugeField> { class TwoFlavourRatioPseudoFermionAction : public Action<typename Impl::GaugeField> {
public: public:
INHERIT_IMPL_TYPES(Impl); INHERIT_IMPL_TYPES(Impl);
private: private:
FermionOperator<Impl> & NumOp;// the basic operator FermionOperator<Impl> & NumOp;// the basic operator
FermionOperator<Impl> & DenOp;// the basic operator FermionOperator<Impl> & DenOp;// the basic operator
OperatorFunction<FermionField> &DerivativeSolver; OperatorFunction<FermionField> &DerivativeSolver;
OperatorFunction<FermionField> &ActionSolver; OperatorFunction<FermionField> &ActionSolver;
FermionField Phi; // the pseudo fermion field for this trajectory FermionField Phi; // the pseudo fermion field for this trajectory
public: public:
TwoFlavourRatioPseudoFermionAction(FermionOperator<Impl> &_NumOp, TwoFlavourRatioPseudoFermionAction(FermionOperator<Impl> &_NumOp,
FermionOperator<Impl> &_DenOp, FermionOperator<Impl> &_DenOp,
OperatorFunction<FermionField> & DS, OperatorFunction<FermionField> & DS,
OperatorFunction<FermionField> & AS OperatorFunction<FermionField> & AS
) : NumOp(_NumOp), DenOp(_DenOp), DerivativeSolver(DS), ActionSolver(AS), Phi(_NumOp.FermionGrid()) {}; ) : 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(){ virtual std::string LogParameters(){
std::stringstream sstream; std::stringstream sstream;
sstream << GridLogMessage << "["<<action_name()<<"] has no parameters" << std::endl; sstream << GridLogMessage << "["<<action_name()<<"] has no parameters" << std::endl;
return sstream.str(); return sstream.str();
} }
virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) { virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) {
// P(phi) = e^{- phi^dag V (MdagM)^-1 Vdag phi} // P(phi) = e^{- phi^dag V (MdagM)^-1 Vdag phi}
// //
// NumOp == V // NumOp == V
// DenOp == M // DenOp == M
// //
// Take phi = Vdag^{-1} Mdag eta ; eta = Mdag^{-1} Vdag Phi // Take phi = Vdag^{-1} Mdag eta ; eta = Mdag^{-1} Vdag Phi
// //
// P(eta) = e^{- eta^dag eta} // P(eta) = e^{- eta^dag eta}
// //
// e^{x^2/2 sig^2} => sig^2 = 0.5. // 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.... // So eta should be of width sig = 1/sqrt(2) and must multiply by 0.707....
// //
RealD scale = std::sqrt(0.5); RealD scale = std::sqrt(0.5);
FermionField eta(NumOp.FermionGrid()); FermionField eta(NumOp.FermionGrid());
FermionField tmp(NumOp.FermionGrid()); FermionField tmp(NumOp.FermionGrid());
gaussian(pRNG,eta); gaussian(pRNG,eta);
NumOp.ImportGauge(U); NumOp.ImportGauge(U);
DenOp.ImportGauge(U); DenOp.ImportGauge(U);
// Note: this hard codes normal equations type solvers; alternate implementation needed for // Note: this hard codes normal equations type solvers; alternate implementation needed for
// non-herm style solvers. // non-herm style solvers.
MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(NumOp); MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(NumOp);
DenOp.Mdag(eta,Phi); // Mdag eta DenOp.Mdag(eta,Phi); // Mdag eta
tmp = zero; tmp = zero;
ActionSolver(MdagMOp,Phi,tmp); // (VdagV)^-1 Mdag eta = V^-1 Vdag^-1 Mdag eta ActionSolver(MdagMOp,Phi,tmp); // (VdagV)^-1 Mdag eta = V^-1 Vdag^-1 Mdag eta
NumOp.M(tmp,Phi); // 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 // S = phi^dag V (Mdag M)^-1 Vdag phi
////////////////////////////////////////////////////// //////////////////////////////////////////////////////
virtual RealD S(const GaugeField &U) { virtual RealD S(const GaugeField &U) {
NumOp.ImportGauge(U); NumOp.ImportGauge(U);
DenOp.ImportGauge(U); DenOp.ImportGauge(U);
FermionField X(NumOp.FermionGrid()); FermionField X(NumOp.FermionGrid());
FermionField Y(NumOp.FermionGrid()); FermionField Y(NumOp.FermionGrid());
MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(DenOp); MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(DenOp);
NumOp.Mdag(Phi,Y); // Y= Vdag phi NumOp.Mdag(Phi,Y); // Y= Vdag phi
X=zero; X=zero;
ActionSolver(MdagMOp,Y,X); // X= (MdagM)^-1 Vdag phi ActionSolver(MdagMOp,Y,X); // X= (MdagM)^-1 Vdag phi
DenOp.M(X,Y); // Y= Mdag^-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 // 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 [ Mdag dM + dMdag M ] (Mdag M)^-1 V^dag phi
// + phi^dag V (Mdag M)^-1 dV^dag phi // + phi^dag V (Mdag M)^-1 dV^dag phi
////////////////////////////////////////////////////// //////////////////////////////////////////////////////
virtual void deriv(const GaugeField &U,GaugeField & dSdU) { virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
NumOp.ImportGauge(U); NumOp.ImportGauge(U);
DenOp.ImportGauge(U); DenOp.ImportGauge(U);
MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(DenOp); MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(DenOp);
FermionField X(NumOp.FermionGrid()); FermionField X(NumOp.FermionGrid());
FermionField Y(NumOp.FermionGrid()); FermionField Y(NumOp.FermionGrid());
GaugeField force(NumOp.GaugeGrid()); GaugeField force(NumOp.GaugeGrid());
//Y=Vdag phi //Y=Vdag phi
//X = (Mdag M)^-1 V^dag phi //X = (Mdag M)^-1 V^dag phi
//Y = (Mdag)^-1 V^dag phi //Y = (Mdag)^-1 V^dag phi
NumOp.Mdag(Phi,Y); // Y= Vdag phi NumOp.Mdag(Phi,Y); // Y= Vdag phi
X=zero; X=zero;
DerivativeSolver(MdagMOp,Y,X); // X= (MdagM)^-1 Vdag phi DerivativeSolver(MdagMOp,Y,X); // X= (MdagM)^-1 Vdag phi
DenOp.M(X,Y); // Y= Mdag^-1 Vdag phi DenOp.M(X,Y); // Y= Mdag^-1 Vdag phi
// phi^dag V (Mdag M)^-1 dV^dag phi // phi^dag V (Mdag M)^-1 dV^dag phi
NumOp.MDeriv(force , X, Phi, DaggerYes ); dSdU=force; NumOp.MDeriv(force , X, Phi, DaggerYes ); dSdU=force;
// phi^dag dV (Mdag M)^-1 V^dag phi // phi^dag dV (Mdag M)^-1 V^dag phi
NumOp.MDeriv(force , Phi, X ,DaggerNo ); dSdU=dSdU+force; 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 Mdag dM (Mdag M)^-1 V^dag phi
// - phi^dag V (Mdag M)^-1 dMdag M (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,Y,X,DaggerNo); dSdU=dSdU-force;
DenOp.MDeriv(force,X,Y,DaggerYes); dSdU=dSdU-force; DenOp.MDeriv(force,X,Y,DaggerYes); dSdU=dSdU-force;
dSdU *= -1.0; dSdU *= -1.0;
//dSdU = - Ta(dSdU); //dSdU = - Ta(dSdU);
};
};
NAMESPACE_END(Grid):
};
};
}
}
#endif #endif