diff --git a/HMC/Mobius2p1fIDSDRGparityEOFA_40ID.cc b/HMC/Mobius2p1fIDSDRGparityEOFA_40ID.cc index 63bc2db4..927a073b 100644 --- a/HMC/Mobius2p1fIDSDRGparityEOFA_40ID.cc +++ b/HMC/Mobius2p1fIDSDRGparityEOFA_40ID.cc @@ -82,7 +82,7 @@ struct EOFAparameters: Serializable { rat_params.lo = 0.1; rat_params.hi = 25.0; - rat_params.MaxIter = 10000; + rat_params.MaxIter = 50000; rat_params.tolerance= 1.0e-9; rat_params.degree = 14; rat_params.precision= 50; @@ -180,7 +180,7 @@ void computeEigenvalues(std::string param_file, FunctionHermOp Cheb_wrap(Cheb, hermop); std::cout << "IRL: alpha=" << params.alpha << " beta=" << params.beta << " mu=" << params.mu << " ord=" << params.ord << std::endl; - ImplicitlyRestartedLanczos IRL(Cheb_wrap, hermop_wrap, params.n_stop, params.n_want, params.n_use, params.tolerance, 10000); + ImplicitlyRestartedLanczos IRL(Cheb_wrap, hermop_wrap, params.n_stop, params.n_want, params.n_use, params.tolerance, 50000); std::vector eval(params.n_use); std::vector evec(params.n_use, rbGrid); @@ -228,19 +228,19 @@ void checkRHMC(GridCartesian* Grid, GridRedBlackCartesian* rbGrid, const Lattice if(action_or_md == 0 || action_or_md == 2){ std::cout << "Starting: Checking quality of RHMC action approx for " << quark_descr << " quark numerator and power -1/" << inv_pow << std::endl; - InversePowerBoundsCheck(inv_pow, 10000, 1e16, MdagM,gauss_o, rhmc.ApproxNegPowerAction); //use large tolerance to prevent exit on fail; we are trying to tune here! + InversePowerBoundsCheck(inv_pow, 50000, 1e16, MdagM,gauss_o, rhmc.ApproxNegPowerAction); //use large tolerance to prevent exit on fail; we are trying to tune here! std::cout << "Finished: Checking quality of RHMC action approx for " << quark_descr << " quark numerator and power -1/" << inv_pow << std::endl; std::cout << "Starting: Checking quality of RHMC action approx for " << quark_descr << " quark numerator and power -1/" << 2*inv_pow << std::endl; - InversePowerBoundsCheck(2*inv_pow, 10000, 1e16, MdagM,gauss_o, rhmc.ApproxNegHalfPowerAction); + InversePowerBoundsCheck(2*inv_pow, 50000, 1e16, MdagM,gauss_o, rhmc.ApproxNegHalfPowerAction); std::cout << "Finished: Checking quality of RHMC action approx for " << quark_descr << " quark numerator and power -1/" << 2*inv_pow << std::endl; std::cout << "Starting: Checking quality of RHMC action approx for " << quark_descr << " quark denominator and power -1/" << inv_pow << std::endl; - InversePowerBoundsCheck(inv_pow, 10000, 1e16, VdagV,gauss_o, rhmc.ApproxNegPowerAction); + InversePowerBoundsCheck(inv_pow, 50000, 1e16, VdagV,gauss_o, rhmc.ApproxNegPowerAction); std::cout << "Finished: Checking quality of RHMC action approx for " << quark_descr << " quark denominator and power -1/" << inv_pow << std::endl; std::cout << "Starting: Checking quality of RHMC action approx for " << quark_descr << " quark denominator and power -1/" << 2*inv_pow << std::endl; - InversePowerBoundsCheck(2*inv_pow, 10000, 1e16, VdagV,gauss_o, rhmc.ApproxNegHalfPowerAction); + InversePowerBoundsCheck(2*inv_pow, 50000, 1e16, VdagV,gauss_o, rhmc.ApproxNegHalfPowerAction); std::cout << "Finished: Checking quality of RHMC action approx for " << quark_descr << " quark denominator and power -1/" << 2*inv_pow << std::endl; } @@ -248,19 +248,19 @@ void checkRHMC(GridCartesian* Grid, GridRedBlackCartesian* rbGrid, const Lattice if(action_or_md == 1 || action_or_md == 2){ std::cout << "Starting: Checking quality of RHMC MD approx for " << quark_descr << " quark numerator and power -1/" << inv_pow << std::endl; - InversePowerBoundsCheck(inv_pow, 10000, 1e16, MdagM,gauss_o, rhmc.ApproxNegPowerMD); + InversePowerBoundsCheck(inv_pow, 50000, 1e16, MdagM,gauss_o, rhmc.ApproxNegPowerMD); std::cout << "Finished: Checking quality of RHMC MD approx for " << quark_descr << " quark numerator and power -1/" << inv_pow << std::endl; std::cout << "Starting: Checking quality of RHMC MD approx for " << quark_descr << " quark numerator and power -1/" << 2*inv_pow << std::endl; - InversePowerBoundsCheck(2*inv_pow, 10000, 1e16, MdagM,gauss_o, rhmc.ApproxNegHalfPowerMD); + InversePowerBoundsCheck(2*inv_pow, 50000, 1e16, MdagM,gauss_o, rhmc.ApproxNegHalfPowerMD); std::cout << "Finished: Checking quality of RHMC MD approx for " << quark_descr << " quark numerator and power -1/" << 2*inv_pow << std::endl; std::cout << "Starting: Checking quality of RHMC MD approx for " << quark_descr << " quark denominator and power -1/" << inv_pow << std::endl; - InversePowerBoundsCheck(inv_pow, 10000, 1e16, VdagV,gauss_o, rhmc.ApproxNegPowerMD); + InversePowerBoundsCheck(inv_pow, 50000, 1e16, VdagV,gauss_o, rhmc.ApproxNegPowerMD); std::cout << "Finished: Checking quality of RHMC MD approx for " << quark_descr << " quark denominator and power -1/" << inv_pow << std::endl; std::cout << "Starting: Checking quality of RHMC MD approx for " << quark_descr << " quark denominator and power -1/" << 2*inv_pow << std::endl; - InversePowerBoundsCheck(2*inv_pow, 10000, 1e16, VdagV,gauss_o, rhmc.ApproxNegHalfPowerMD); + InversePowerBoundsCheck(2*inv_pow, 50000, 1e16, VdagV,gauss_o, rhmc.ApproxNegHalfPowerMD); std::cout << "Finished: Checking quality of RHMC MD approx for " << quark_descr << " quark denominator and power -1/" << 2*inv_pow << std::endl; } } @@ -554,11 +554,16 @@ int main(int argc, char **argv) { //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: IntegratorParameters MD; typedef ConjugateHMCRunnerD HMCWrapper; //NB: This is the "Omelyan integrator" - typedef HMCWrapper::ImplPolicy GaugeImplPolicy; MD.name = std::string("MinimumNorm2"); + + // typedef ConjugateHMCRunnerD HMCWrapper; + // MD.name = std::string("ForceGradient"); + MD.MDsteps = user_params.Steps; MD.trajL = 1.0; + typedef HMCWrapper::ImplPolicy GaugeImplPolicy; + HMCparameters HMCparams; HMCparams.StartTrajectory = user_params.StartTrajectory; HMCparams.Trajectories = user_params.Trajectories; @@ -660,11 +665,11 @@ int main(int argc, char **argv) { typedef MixedPrecisionConjugateGradientOperatorFunction EOFA_mxCG; typedef MixedPrecisionReliableUpdateConjugateGradientOperatorFunction EOFA_relupCG; - std::vector eofa_light_masses = { light_mass , 0.01 }; - std::vector eofa_pv_masses = { 0.01 , 1.0 }; - int n_light_hsb = 2; + std::vector eofa_light_masses = { light_mass , 0.01, 0.1 }; + std::vector eofa_pv_masses = { 0.01 , 0.1, 1.0 }; + int n_light_hsb = 3; - EOFAmixPrecPFaction* EOFA_pfactions[2]; + EOFAmixPrecPFaction* EOFA_pfactions[n_light_hsb]; for(int i=0;iInnerTolerance = user_params.eofa_l.action_mixcg_inner_tolerance; - EOFA_mxCG* ActionMCG_R = new EOFA_mxCG(user_params.eofa_l.action_tolerance, 10000, 1000, UGridF, FrbGridF, *RopF, *RopD, *linopR_F, *linopR_D); + EOFA_mxCG* ActionMCG_R = new EOFA_mxCG(user_params.eofa_l.action_tolerance, 50000, 1000, UGridF, FrbGridF, *RopF, *RopD, *linopR_F, *linopR_D); ActionMCG_R->InnerTolerance = user_params.eofa_l.action_mixcg_inner_tolerance; - EOFA_mxCG* DerivMCG_L = new EOFA_mxCG(user_params.eofa_l.md_tolerance, 10000, 1000, UGridF, FrbGridF, *LopF, *LopD, *linopL_F, *linopL_D); + EOFA_mxCG* DerivMCG_L = new EOFA_mxCG(user_params.eofa_l.md_tolerance, 50000, 1000, UGridF, FrbGridF, *LopF, *LopD, *linopL_F, *linopL_D); DerivMCG_L->InnerTolerance = user_params.eofa_l.md_mixcg_inner_tolerance; - EOFA_mxCG* DerivMCG_R = new EOFA_mxCG(user_params.eofa_l.md_tolerance, 10000, 1000, UGridF, FrbGridF, *RopF, *RopD, *linopR_F, *linopR_D); + EOFA_mxCG* DerivMCG_R = new EOFA_mxCG(user_params.eofa_l.md_tolerance, 50000, 1000, UGridF, FrbGridF, *RopF, *RopD, *linopR_F, *linopR_D); DerivMCG_R->InnerTolerance = user_params.eofa_l.md_mixcg_inner_tolerance; std::cout << GridLogMessage << "Set EOFA action solver action tolerance outer=" << ActionMCG_L->Tolerance << " inner=" << ActionMCG_L->InnerTolerance << std::endl; @@ -728,7 +733,7 @@ int main(int argc, char **argv) { RationalActionParams rat_act_params_s; rat_act_params_s.inv_pow = 4; // (M^dag M)^{1/4} rat_act_params_s.precision= 60; - rat_act_params_s.MaxIter = 10000; + rat_act_params_s.MaxIter = 50000; user_params.rat_quo_s.Export(rat_act_params_s); std::cout << GridLogMessage << " Heavy quark bounds check every " << rat_act_params_s.BoundsCheckFreq << " trajectories (avg)" << std::endl; @@ -752,7 +757,7 @@ int main(int argc, char **argv) { RationalActionParams rat_act_params_DSDR; rat_act_params_DSDR.inv_pow = 2; // (M^dag M)^{1/2} rat_act_params_DSDR.precision= 60; - rat_act_params_DSDR.MaxIter = 10000; + rat_act_params_DSDR.MaxIter = 50000; user_params.rat_quo_DSDR.Export(rat_act_params_DSDR); std::cout << GridLogMessage << "DSDR quark bounds check every " << rat_act_params_DSDR.BoundsCheckFreq << " trajectories (avg)" << std::endl; diff --git a/HMC/Mobius2p1fIDSDRGparityEOFA_48ID.cc b/HMC/Mobius2p1fIDSDRGparityEOFA_48ID.cc index eeb5b4b7..09f7a325 100644 --- a/HMC/Mobius2p1fIDSDRGparityEOFA_48ID.cc +++ b/HMC/Mobius2p1fIDSDRGparityEOFA_48ID.cc @@ -80,7 +80,7 @@ struct EOFAparameters: Serializable { md_tolerance = 1e-8; md_mixcg_inner_tolerance = 1e-8; - rat_params.lo = 0.1; + rat_params.lo = 1.0; rat_params.hi = 25.0; rat_params.MaxIter = 10000; rat_params.tolerance= 1.0e-9; @@ -98,7 +98,7 @@ struct EvolParameters: Serializable { bool, MetropolisTest, std::string, StartingType, std::vector, GparityDirs, - EOFAparameters, eofa_l, + std::vector, eofa_l, RatQuoParameters, rat_quo_s, RatQuoParameters, rat_quo_DSDR); @@ -578,8 +578,8 @@ int main(int argc, char **argv) { // Collect actions //////////////////////////////////// ActionLevel Level1(1); //light quark + strange quark - ActionLevel Level2(1); //DSDR - ActionLevel Level3(8); //gauge (8 increments per step) + ActionLevel Level2(2); //DSDR + ActionLevel Level3(4); //gauge ///////////////////////////////////////////////////////////// @@ -591,11 +591,12 @@ int main(int argc, char **argv) { typedef ExactOneFlavourRatioMixedPrecHeatbathPseudoFermionAction EOFAmixPrecPFaction; typedef MixedPrecisionConjugateGradientOperatorFunction EOFA_mxCG; - std::vector eofa_light_masses = { light_mass , 0.1 }; - std::vector eofa_pv_masses = { 0.1 , 1.0 }; - int n_light_hsb = 2; - - EOFAmixPrecPFaction* EOFA_pfactions[2]; + std::vector eofa_light_masses = { light_mass , 0.004, 0.016, 0.064, 0.256 }; + std::vector eofa_pv_masses = { 0.004 , 0.016, 0.064, 0.256, 1.0 }; + int n_light_hsb = 5; + assert(user_params.eofa_l.size() == n_light_hsb); + + EOFAmixPrecPFaction* EOFA_pfactions[n_light_hsb]; for(int i=0;iInnerTolerance = user_params.eofa_l.action_mixcg_inner_tolerance; + EOFA_mxCG* ActionMCG_L = new EOFA_mxCG(user_params.eofa_l[i].action_tolerance, 10000, 1000, UGridF, FrbGridF, *LopF, *LopD, *linopL_F, *linopL_D); + ActionMCG_L->InnerTolerance = user_params.eofa_l[i].action_mixcg_inner_tolerance; - EOFA_mxCG* ActionMCG_R = new EOFA_mxCG(user_params.eofa_l.action_tolerance, 10000, 1000, UGridF, FrbGridF, *RopF, *RopD, *linopR_F, *linopR_D); - ActionMCG_R->InnerTolerance = user_params.eofa_l.action_mixcg_inner_tolerance; + EOFA_mxCG* ActionMCG_R = new EOFA_mxCG(user_params.eofa_l[i].action_tolerance, 10000, 1000, UGridF, FrbGridF, *RopF, *RopD, *linopR_F, *linopR_D); + ActionMCG_R->InnerTolerance = user_params.eofa_l[i].action_mixcg_inner_tolerance; - EOFA_mxCG* DerivMCG_L = new EOFA_mxCG(user_params.eofa_l.md_tolerance, 10000, 1000, UGridF, FrbGridF, *LopF, *LopD, *linopL_F, *linopL_D); - DerivMCG_L->InnerTolerance = user_params.eofa_l.md_mixcg_inner_tolerance; + EOFA_mxCG* DerivMCG_L = new EOFA_mxCG(user_params.eofa_l[i].md_tolerance, 10000, 1000, UGridF, FrbGridF, *LopF, *LopD, *linopL_F, *linopL_D); + DerivMCG_L->InnerTolerance = user_params.eofa_l[i].md_mixcg_inner_tolerance; - EOFA_mxCG* DerivMCG_R = new EOFA_mxCG(user_params.eofa_l.md_tolerance, 10000, 1000, UGridF, FrbGridF, *RopF, *RopD, *linopR_F, *linopR_D); - DerivMCG_R->InnerTolerance = user_params.eofa_l.md_mixcg_inner_tolerance; + EOFA_mxCG* DerivMCG_R = new EOFA_mxCG(user_params.eofa_l[i].md_tolerance, 10000, 1000, UGridF, FrbGridF, *RopF, *RopD, *linopR_F, *linopR_D); + DerivMCG_R->InnerTolerance = user_params.eofa_l[i].md_mixcg_inner_tolerance; std::cout << GridLogMessage << "Set EOFA action solver action tolerance outer=" << ActionMCG_L->Tolerance << " inner=" << ActionMCG_L->InnerTolerance << std::endl; std::cout << GridLogMessage << "Set EOFA MD solver tolerance outer=" << DerivMCG_L->Tolerance << " inner=" << DerivMCG_L->InnerTolerance << std::endl; @@ -632,7 +633,7 @@ int main(int argc, char **argv) { *ActionMCG_L, *ActionMCG_R, *ActionMCG_L, *ActionMCG_R, *DerivMCG_L, *DerivMCG_R, - user_params.eofa_l.rat_params, true); + user_params.eofa_l[i].rat_params, true); EOFA_pfactions[i] = EOFA; Level1.push_back(EOFA); } @@ -729,8 +730,8 @@ int main(int argc, char **argv) { else if(sarg == "--check_eofa"){ assert(i < argc-1); check_eofa = true; - eofa_which_hsb = std::stoi(argv[i+1]); - assert(eofa_which_hsb >= 0 && eofa_which_hsb < n_light_hsb); + eofa_which_hsb = std::stoi(argv[i+1]); //-1 indicates all hasenbusch + assert(eofa_which_hsb == -1 || (eofa_which_hsb >= 0 && eofa_which_hsb < n_light_hsb) ); } else if(sarg == "--upper_bound_eofa"){ assert(i < argc-1); @@ -753,7 +754,19 @@ int main(int argc, char **argv) { //std::cout << GridLogMessage << "EOFA MD solver tolerance outer=" << DerivMCG_L.Tolerance << " inner=" << DerivMCG_L.InnerTolerance << std::endl; - if(check_eofa) checkEOFA(*EOFA_pfactions[eofa_which_hsb], FGridD, TheHMC.Resources.GetParallelRNG(), Ud); + if(check_eofa){ + if(eofa_which_hsb >= 0){ + std::cout << GridLogMessage << "Starting checking EOFA Hasenbusch " << eofa_which_hsb << std::endl; + checkEOFA(*EOFA_pfactions[eofa_which_hsb], FGridD, TheHMC.Resources.GetParallelRNG(), Ud); + std::cout << GridLogMessage << "Finished checking EOFA Hasenbusch " << eofa_which_hsb << std::endl; + }else{ + for(int i=0;i(lanc_params_s, FGridD, FrbGridD, Ud, Numerator_sD, TheHMC.Resources.GetParallelRNG());