portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-04 11:15:55 +01:00
Code Releases Activity

General RHMC pseudofermion action now allows for different rational approximations to be used in the MD and action evaluation

Browse Source
This commit is contained in:
Christopher Kelly 2020-12-23 11:19:26 -05:00
parent 220ad5e3ee
commit 75c6c6b173
2 changed files with 88 additions and 41 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
Grid/qcd/action/ActionParams.h
View File

@ -97,28 +97,34 @@ struct StaggeredImplParams {
struct RationalActionParams : Serializable { struct RationalActionParams : Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(RationalActionParams, GRID_SERIALIZABLE_CLASS_MEMBERS(RationalActionParams,
int, inv_pow, int, inv_pow,
RealD, lo, RealD, lo, //low eigenvalue bound of rational approx
RealD, hi, RealD, hi, //high eigenvalue bound of rational approx
int, MaxIter, int, MaxIter, //maximum iterations in msCG
RealD, tolerance, RealD, action_tolerance, //msCG tolerance in action evaluation
int, degree, int, action_degree, //rational approx tolerance in action evaluation
int, precision, RealD, md_tolerance, //msCG tolerance in MD integration
int, BoundsCheckFreq); int, md_degree, //rational approx tolerance in MD integration
int, precision, //precision of floating point arithmetic
int, BoundsCheckFreq); //frequency the approximation is tested (with Metropolis degree/tolerance); 0 disables the check
// constructor // constructor
RationalActionParams(int _inv_pow = 2, RationalActionParams(int _inv_pow = 2,
RealD _lo = 0.0, RealD _lo = 0.0,
RealD _hi = 1.0, RealD _hi = 1.0,
int _maxit = 1000, int _maxit = 1000,
RealD tol = 1.0e-8, RealD _action_tolerance = 1.0e-8,
int _degree = 10, int _action_degree = 10,
RealD _md_tolerance = 1.0e-8,
int _md_degree = 10,
int _precision = 64, int _precision = 64,
int _BoundsCheckFreq=20) int _BoundsCheckFreq=20)
: inv_pow(_inv_pow), : inv_pow(_inv_pow),
lo(_lo), lo(_lo),
hi(_hi), hi(_hi),
MaxIter(_maxit), MaxIter(_maxit),
tolerance(tol), action_tolerance(_action_tolerance),
degree(_degree), action_degree(_action_degree),
md_tolerance(_md_tolerance),
md_degree(_md_degree),
precision(_precision), precision(_precision),
BoundsCheckFreq(_BoundsCheckFreq){}; BoundsCheckFreq(_BoundsCheckFreq){};
}; };

101
Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatio.h
View File

@ -60,10 +60,17 @@ NAMESPACE_BEGIN(Grid);
typedef RationalActionParams Params; typedef RationalActionParams Params;
Params param; Params param;
MultiShiftFunction ApproxPower ; //rational approx for X^{1/inv_pow} //For action evaluation
MultiShiftFunction ApproxNegPower; //rational approx for X^{-1/inv_pow} MultiShiftFunction ApproxPowerAction ; //rational approx for X^{1/inv_pow}
MultiShiftFunction ApproxHalfPower; //rational approx for X^{1/(2*inv_pow)} MultiShiftFunction ApproxNegPowerAction; //rational approx for X^{-1/inv_pow}
MultiShiftFunction ApproxNegHalfPower; //rational approx for X^{-1/(2*inv_pow)} MultiShiftFunction ApproxHalfPowerAction; //rational approx for X^{1/(2*inv_pow)}
MultiShiftFunction ApproxNegHalfPowerAction; //rational approx for X^{-1/(2*inv_pow)}
//For the MD integration
MultiShiftFunction ApproxPowerMD ; //rational approx for X^{1/inv_pow}
MultiShiftFunction ApproxNegPowerMD; //rational approx for X^{-1/inv_pow}
MultiShiftFunction ApproxHalfPowerMD; //rational approx for X^{1/(2*inv_pow)}
MultiShiftFunction ApproxNegHalfPowerMD; //rational approx for X^{-1/(2*inv_pow)}
private: private:
@ -90,17 +97,49 @@ NAMESPACE_BEGIN(Grid);
int inv_pow = param.inv_pow; int inv_pow = param.inv_pow;
int _2_inv_pow = 2*inv_pow; int _2_inv_pow = 2*inv_pow;
//Generate approximations for action eval
// MdagM^(+- 1/inv_pow) // MdagM^(+- 1/inv_pow)
std::cout<<GridLogMessage << "Generating degree "<<param.degree<<" for x^(1/" << inv_pow << ")"<<std::endl; std::cout<<GridLogMessage << "Generating degree "<<param.action_degree<<" and tolerance " << param.action_tolerance << " for x^(1/" << inv_pow << ")"<<std::endl;
remez.generateApprox(param.degree,1,inv_pow); remez.generateApprox(param.action_degree,1,inv_pow);
ApproxPower.Init(remez,param.tolerance,false); ApproxPowerAction.Init(remez,param.action_tolerance,false);
ApproxNegPower.Init(remez,param.tolerance,true); ApproxNegPowerAction.Init(remez,param.action_tolerance,true);
// VdagV^(+- 1/(2*inv_pow)) // VdagV^(+- 1/(2*inv_pow))
std::cout<<GridLogMessage << "Generating degree "<<param.degree<<" for x^(1/" << _2_inv_pow << ")"<<std::endl; std::cout<<GridLogMessage << "Generating degree "<<param.action_degree<<" and tolerance " << param.action_tolerance <<" for x^(1/" << _2_inv_pow << ")"<<std::endl;
remez.generateApprox(param.degree,1,_2_inv_pow); remez.generateApprox(param.action_degree,1,_2_inv_pow);
ApproxHalfPower.Init(remez,param.tolerance,false); ApproxHalfPowerAction.Init(remez,param.action_tolerance,false);
ApproxNegHalfPower.Init(remez,param.tolerance,true); ApproxNegHalfPowerAction.Init(remez,param.action_tolerance,true);
//Generate approximations for MD
if(param.md_degree != param.action_degree ||
param.md_tolerance < param.action_tolerance //no point in finding less precise polynomial if the degree is the same
){
// MdagM^(+- 1/inv_pow)
std::cout<<GridLogMessage << "Generating degree "<<param.md_degree<<" and tolerance " << param.md_tolerance <<" for x^(1/" << inv_pow << ")"<<std::endl;
remez.generateApprox(param.md_degree,1,inv_pow);
ApproxPowerMD.Init(remez,param.md_tolerance,false);
ApproxNegPowerMD.Init(remez,param.md_tolerance,true);
// VdagV^(+- 1/(2*inv_pow))
std::cout<<GridLogMessage << "Generating degree "<<param.md_degree<<" and tolerance " << param.md_tolerance <<" for x^(1/" << _2_inv_pow << ")"<<std::endl;
remez.generateApprox(param.md_degree,1,_2_inv_pow);
ApproxHalfPowerMD.Init(remez,param.md_tolerance,false);
ApproxNegHalfPowerMD.Init(remez,param.md_tolerance,true);
}else{
std::cout<<GridLogMessage << "Using same rational approximations for MD as for action evaluation" << std::endl;
ApproxPowerMD = ApproxPowerAction;
ApproxNegPowerMD = ApproxNegPowerAction;
for(int i=0;i<ApproxPowerMD.tolerances.size();i++)
ApproxNegPowerMD.tolerances[i] = ApproxPowerMD.tolerances[i] = param.md_tolerance; //used for multishift
ApproxHalfPowerMD = ApproxHalfPowerAction;
ApproxNegHalfPowerMD = ApproxNegHalfPowerAction;
for(int i=0;i<ApproxPowerMD.tolerances.size();i++)
ApproxNegHalfPowerMD.tolerances[i] = ApproxHalfPowerMD.tolerances[i] = param.md_tolerance;
}
std::cout<<GridLogMessage << action_name() << " initialize: complete" << std::endl; std::cout<<GridLogMessage << action_name() << " initialize: complete" << std::endl;
}; };
@ -108,13 +147,15 @@ NAMESPACE_BEGIN(Grid);
virtual std::string LogParameters(){ virtual std::string LogParameters(){
std::stringstream sstream; std::stringstream sstream;
sstream << GridLogMessage << "["<<action_name()<<"] Power : 1/" << param.inv_pow << std::endl; sstream << GridLogMessage << "["<<action_name()<<"] Power : 1/" << param.inv_pow << std::endl;
sstream << GridLogMessage << "["<<action_name()<<"] Low :" << param.lo << std::endl; sstream << GridLogMessage << "["<<action_name()<<"] Low :" << param.lo << std::endl;
sstream << GridLogMessage << "["<<action_name()<<"] High :" << param.hi << std::endl; sstream << GridLogMessage << "["<<action_name()<<"] High :" << param.hi << std::endl;
sstream << GridLogMessage << "["<<action_name()<<"] Max iterations :" << param.MaxIter << std::endl; sstream << GridLogMessage << "["<<action_name()<<"] Max iterations :" << param.MaxIter << std::endl;
sstream << GridLogMessage << "["<<action_name()<<"] Tolerance :" << param.tolerance << std::endl; sstream << GridLogMessage << "["<<action_name()<<"] Tolerance (Action) :" << param.action_tolerance << std::endl;
sstream << GridLogMessage << "["<<action_name()<<"] Degree :" << param.degree << std::endl; sstream << GridLogMessage << "["<<action_name()<<"] Degree (Action) :" << param.action_degree << std::endl;
sstream << GridLogMessage << "["<<action_name()<<"] Precision :" << param.precision << std::endl; sstream << GridLogMessage << "["<<action_name()<<"] Tolerance (MD) :" << param.md_tolerance << std::endl;
sstream << GridLogMessage << "["<<action_name()<<"] Degree (MD) :" << param.md_degree << std::endl;
sstream << GridLogMessage << "["<<action_name()<<"] Precision :" << param.precision << std::endl;
return sstream.str(); return sstream.str();
} }
@ -154,13 +195,13 @@ NAMESPACE_BEGIN(Grid);
// MdagM^1/(2*inv_pow) eta // MdagM^1/(2*inv_pow) eta
std::cout<<GridLogMessage << action_name() << " refresh: doing (M^dag M)^{1/" << 2*param.inv_pow << "} eta" << std::endl; std::cout<<GridLogMessage << action_name() << " refresh: doing (M^dag M)^{1/" << 2*param.inv_pow << "} eta" << std::endl;
SchurDifferentiableOperator<Impl> MdagM(DenOp); SchurDifferentiableOperator<Impl> MdagM(DenOp);
ConjugateGradientMultiShift<FermionField> msCG_M(param.MaxIter,ApproxHalfPower); ConjugateGradientMultiShift<FermionField> msCG_M(param.MaxIter,ApproxHalfPowerAction);
msCG_M(MdagM,etaOdd,tmp); msCG_M(MdagM,etaOdd,tmp);
// VdagV^-1/(2*inv_pow) MdagM^1/(2*inv_pow) eta // VdagV^-1/(2*inv_pow) MdagM^1/(2*inv_pow) eta
std::cout<<GridLogMessage << action_name() << " refresh: doing (V^dag V)^{-1/" << 2*param.inv_pow << "} ( (M^dag M)^{1/" << 2*param.inv_pow << "} eta)" << std::endl; std::cout<<GridLogMessage << action_name() << " refresh: doing (V^dag V)^{-1/" << 2*param.inv_pow << "} ( (M^dag M)^{1/" << 2*param.inv_pow << "} eta)" << std::endl;
SchurDifferentiableOperator<Impl> VdagV(NumOp); SchurDifferentiableOperator<Impl> VdagV(NumOp);
ConjugateGradientMultiShift<FermionField> msCG_V(param.MaxIter,ApproxNegHalfPower); ConjugateGradientMultiShift<FermionField> msCG_V(param.MaxIter,ApproxNegHalfPowerAction);
msCG_V(VdagV,tmp,PhiOdd); msCG_V(VdagV,tmp,PhiOdd);
assert(NumOp.ConstEE() == 1); assert(NumOp.ConstEE() == 1);
@ -183,13 +224,13 @@ NAMESPACE_BEGIN(Grid);
// VdagV^1/(2*inv_pow) Phi // VdagV^1/(2*inv_pow) Phi
std::cout<<GridLogMessage << action_name() << " compute action: doing (V^dag V)^{1/" << 2*param.inv_pow << "} Phi" << std::endl; std::cout<<GridLogMessage << action_name() << " compute action: doing (V^dag V)^{1/" << 2*param.inv_pow << "} Phi" << std::endl;
SchurDifferentiableOperator<Impl> VdagV(NumOp); SchurDifferentiableOperator<Impl> VdagV(NumOp);
ConjugateGradientMultiShift<FermionField> msCG_V(param.MaxIter,ApproxHalfPower); ConjugateGradientMultiShift<FermionField> msCG_V(param.MaxIter,ApproxHalfPowerAction);
msCG_V(VdagV,PhiOdd,X); msCG_V(VdagV,PhiOdd,X);
// MdagM^-1/(2*inv_pow) VdagV^1/(2*inv_pow) Phi // MdagM^-1/(2*inv_pow) VdagV^1/(2*inv_pow) Phi
std::cout<<GridLogMessage << action_name() << " compute action: doing (M^dag M)^{-1/" << 2*param.inv_pow << "} ( (V^dag V)^{1/" << 2*param.inv_pow << "} Phi)" << std::endl; std::cout<<GridLogMessage << action_name() << " compute action: doing (M^dag M)^{-1/" << 2*param.inv_pow << "} ( (V^dag V)^{1/" << 2*param.inv_pow << "} Phi)" << std::endl;
SchurDifferentiableOperator<Impl> MdagM(DenOp); SchurDifferentiableOperator<Impl> MdagM(DenOp);
ConjugateGradientMultiShift<FermionField> msCG_M(param.MaxIter,ApproxNegHalfPower); ConjugateGradientMultiShift<FermionField> msCG_M(param.MaxIter,ApproxNegHalfPowerAction);
msCG_M(MdagM,X,Y); msCG_M(MdagM,X,Y);
// Randomly apply rational bounds checks. // Randomly apply rational bounds checks.
@ -198,7 +239,7 @@ NAMESPACE_BEGIN(Grid);
FermionField gauss(NumOp.FermionRedBlackGrid()); FermionField gauss(NumOp.FermionRedBlackGrid());
gauss = PhiOdd; gauss = PhiOdd;
HighBoundCheck(MdagM,gauss,param.hi); HighBoundCheck(MdagM,gauss,param.hi);
InversePowerBoundsCheck(param.inv_pow,param.MaxIter,param.tolerance*100,MdagM,gauss,ApproxNegPower); InversePowerBoundsCheck(param.inv_pow,param.MaxIter,param.action_tolerance*100,MdagM,gauss,ApproxNegPowerAction);
} }
// Phidag VdagV^1/(2*inv_pow) MdagM^-1/(2*inv_pow) MdagM^-1/(2*inv_pow) VdagV^1/(2*inv_pow) Phi // Phidag VdagV^1/(2*inv_pow) MdagM^-1/(2*inv_pow) MdagM^-1/(2*inv_pow) VdagV^1/(2*inv_pow) Phi
@ -240,8 +281,8 @@ NAMESPACE_BEGIN(Grid);
virtual void deriv(const GaugeField &U,GaugeField & dSdU) { virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
std::cout<<GridLogMessage << action_name() << " deriv: starting" << std::endl; std::cout<<GridLogMessage << action_name() << " deriv: starting" << std::endl;
const int n_f = ApproxNegPower.poles.size(); const int n_f = ApproxNegPowerMD.poles.size();
const int n_pv = ApproxHalfPower.poles.size(); const int n_pv = ApproxHalfPowerMD.poles.size();
std::vector<FermionField> MpvPhi_k (n_pv,NumOp.FermionRedBlackGrid()); std::vector<FermionField> MpvPhi_k (n_pv,NumOp.FermionRedBlackGrid());
std::vector<FermionField> MpvMfMpvPhi_k(n_pv,NumOp.FermionRedBlackGrid()); std::vector<FermionField> MpvMfMpvPhi_k(n_pv,NumOp.FermionRedBlackGrid());
@ -260,8 +301,8 @@ NAMESPACE_BEGIN(Grid);
SchurDifferentiableOperator<Impl> VdagV(NumOp); SchurDifferentiableOperator<Impl> VdagV(NumOp);
SchurDifferentiableOperator<Impl> MdagM(DenOp); SchurDifferentiableOperator<Impl> MdagM(DenOp);
ConjugateGradientMultiShift<FermionField> msCG_V(param.MaxIter,ApproxHalfPower); ConjugateGradientMultiShift<FermionField> msCG_V(param.MaxIter,ApproxHalfPowerMD);
ConjugateGradientMultiShift<FermionField> msCG_M(param.MaxIter,ApproxNegPower); ConjugateGradientMultiShift<FermionField> msCG_M(param.MaxIter,ApproxNegPowerMD);
std::cout<<GridLogMessage << action_name() << " deriv: doing (V^dag V)^{1/" << 2*param.inv_pow << "} Phi" << std::endl; std::cout<<GridLogMessage << action_name() << " deriv: doing (V^dag V)^{1/" << 2*param.inv_pow << "} Phi" << std::endl;
msCG_V(VdagV,PhiOdd,MpvPhi_k,MpvPhi); msCG_V(VdagV,PhiOdd,MpvPhi_k,MpvPhi);
@ -284,7 +325,7 @@ NAMESPACE_BEGIN(Grid);
//(1) //(1)
std::cout<<GridLogMessage << action_name() << " deriv: doing dS/dU part (1)" << std::endl; std::cout<<GridLogMessage << action_name() << " deriv: doing dS/dU part (1)" << std::endl;
for(int k=0;k<n_f;k++){ for(int k=0;k<n_f;k++){
ak = ApproxNegPower.residues[k]; ak = ApproxNegPowerMD.residues[k];
MdagM.Mpc(MfMpvPhi_k[k],Y); MdagM.Mpc(MfMpvPhi_k[k],Y);
MdagM.MpcDagDeriv(tmp , MfMpvPhi_k[k], Y ); dSdU=dSdU+ak*tmp; MdagM.MpcDagDeriv(tmp , MfMpvPhi_k[k], Y ); dSdU=dSdU+ak*tmp;
MdagM.MpcDeriv(tmp , Y, MfMpvPhi_k[k] ); dSdU=dSdU+ak*tmp; MdagM.MpcDeriv(tmp , Y, MfMpvPhi_k[k] ); dSdU=dSdU+ak*tmp;
@ -295,7 +336,7 @@ NAMESPACE_BEGIN(Grid);
std::cout<<GridLogMessage << action_name() << " deriv: doing dS/dU part (2)+(3)" << std::endl; std::cout<<GridLogMessage << action_name() << " deriv: doing dS/dU part (2)+(3)" << std::endl;
for(int k=0;k<n_pv;k++){ for(int k=0;k<n_pv;k++){
ak = ApproxHalfPower.residues[k]; ak = ApproxHalfPowerMD.residues[k];
VdagV.Mpc(MpvPhi_k[k],Y); VdagV.Mpc(MpvPhi_k[k],Y);
VdagV.MpcDagDeriv(tmp,MpvMfMpvPhi_k[k],Y); dSdU=dSdU+ak*tmp; VdagV.MpcDagDeriv(tmp,MpvMfMpvPhi_k[k],Y); dSdU=dSdU+ak*tmp;