From 177b1a7ec630bc02f661e2e10299b73b47e9281a Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sun, 10 Jul 2022 21:34:10 +0100 Subject: [PATCH] Mixed prec --- HMC/Mobius2p1f_DD_RHMC_96I_mixed.cc | 13 +- HMC/Mobius2p1f_DD_RHMC_96I_mixedmshift.cc | 470 ++++++++++++++++++++++ 2 files changed, 480 insertions(+), 3 deletions(-) create mode 100644 HMC/Mobius2p1f_DD_RHMC_96I_mixedmshift.cc diff --git a/HMC/Mobius2p1f_DD_RHMC_96I_mixed.cc b/HMC/Mobius2p1f_DD_RHMC_96I_mixed.cc index a9b5dc7e..13bda00b 100644 --- a/HMC/Mobius2p1f_DD_RHMC_96I_mixed.cc +++ b/HMC/Mobius2p1f_DD_RHMC_96I_mixed.cc @@ -128,8 +128,14 @@ template MPCG(Tolerance,MaxInnerIterations,MaxOuterIterations,SinglePrecGrid5,LinOpF,LinOpD); +#if 1 + RealD delta=1.e-4; + std::cout << GridLogMessage << "Calling reliable update Conjugate Gradient" < MPCG(Tolerance,MaxInnerIterations*MaxOuterIterations,delta,SinglePrecGrid5,LinOpF,LinOpD); +#else std::cout << GridLogMessage << "Calling mixed precision Conjugate Gradient" < MPCG(Tolerance,MaxInnerIterations,MaxOuterIterations,SinglePrecGrid5,LinOpF,LinOpD); +#endif MPCG(src,psi); } }; @@ -161,7 +167,7 @@ int main(int argc, char **argv) { // MD.name = std::string("Force Gradient"); typedef GenericHMCRunner HMCWrapper; MD.name = std::string("MinimumNorm2"); - MD.MDsteps = 4; + MD.MDsteps = 6; MD.trajL = 1.0; HMCparameters HMCparams; @@ -204,7 +210,8 @@ int main(int argc, char **argv) { Real light_mass = 7.8e-4; Real strange_mass = 0.02132; Real pv_mass = 1.0; - std::vector hasenbusch({ light_mass, 3.8e-3, 0.0145, 0.045, 0.108, 0.25, 0.51 , pv_mass }); + // std::vector hasenbusch({ light_mass, 3.8e-3, 0.0145, 0.045, 0.108, 0.25, 0.51 , pv_mass }); + std::vector hasenbusch({ light_mass, 5e-3, 0.0145, 0.045, 0.108, 0.25, 0.51 , pv_mass }); // FIXME: // Same in MC and MD diff --git a/HMC/Mobius2p1f_DD_RHMC_96I_mixedmshift.cc b/HMC/Mobius2p1f_DD_RHMC_96I_mixedmshift.cc new file mode 100644 index 00000000..869f41f8 --- /dev/null +++ b/HMC/Mobius2p1f_DD_RHMC_96I_mixedmshift.cc @@ -0,0 +1,470 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/Test_hmc_EODWFRatio.cc + +Copyright (C) 2015-2016 + +Author: Peter Boyle +Author: Guido Cossu + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +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 */ +#include + +NAMESPACE_BEGIN(Grid); + +template + class MixedPrecisionConjugateGradientOperatorFunction : public OperatorFunction { + public: + typedef typename FermionOperatorD::FermionField FieldD; + typedef typename FermionOperatorF::FermionField FieldF; + + using OperatorFunction::operator(); + + RealD Tolerance; + RealD InnerTolerance; //Initial tolerance for inner CG. Defaults to Tolerance but can be changed + Integer MaxInnerIterations; + Integer MaxOuterIterations; + GridBase* SinglePrecGrid4; //Grid for single-precision fields + GridBase* SinglePrecGrid5; //Grid for single-precision fields + RealD OuterLoopNormMult; //Stop the outer loop and move to a final double prec solve when the residual is OuterLoopNormMult * Tolerance + + FermionOperatorF &FermOpF; + FermionOperatorD &FermOpD;; + SchurOperatorF &LinOpF; + SchurOperatorD &LinOpD; + + Integer TotalInnerIterations; //Number of inner CG iterations + Integer TotalOuterIterations; //Number of restarts + Integer TotalFinalStepIterations; //Number of CG iterations in final patch-up step + + MixedPrecisionConjugateGradientOperatorFunction(RealD tol, + Integer maxinnerit, + Integer maxouterit, + GridBase* _sp_grid4, + GridBase* _sp_grid5, + FermionOperatorF &_FermOpF, + FermionOperatorD &_FermOpD, + SchurOperatorF &_LinOpF, + SchurOperatorD &_LinOpD): + LinOpF(_LinOpF), + LinOpD(_LinOpD), + FermOpF(_FermOpF), + FermOpD(_FermOpD), + Tolerance(tol), + InnerTolerance(tol), + MaxInnerIterations(maxinnerit), + MaxOuterIterations(maxouterit), + SinglePrecGrid4(_sp_grid4), + SinglePrecGrid5(_sp_grid5), + OuterLoopNormMult(100.) + { + /* Debugging instances of objects; references are stored + std::cout << GridLogMessage << " Mixed precision CG wrapper LinOpF " < &LinOpU, const FieldD &src, FieldD &psi) { + + std::cout << GridLogMessage << " Mixed precision CG wrapper operator() "<(&LinOpU); + + // std::cout << GridLogMessage << " Mixed precision CG wrapper operator() FermOpU " <_Mat)<_Mat)==&(LinOpD._Mat)); + + //////////////////////////////////////////////////////////////////////////////////// + // Must snarf a single precision copy of the gauge field in Linop_d argument + //////////////////////////////////////////////////////////////////////////////////// + typedef typename FermionOperatorF::GaugeField GaugeFieldF; + typedef typename FermionOperatorF::GaugeLinkField GaugeLinkFieldF; + typedef typename FermionOperatorD::GaugeField GaugeFieldD; + typedef typename FermionOperatorD::GaugeLinkField GaugeLinkFieldD; + + GridBase * GridPtrF = SinglePrecGrid4; + GridBase * GridPtrD = FermOpD.Umu.Grid(); + GaugeFieldF U_f (GridPtrF); + GaugeLinkFieldF Umu_f(GridPtrF); + // std::cout << " Dim gauge field "<Nd()<Nd()<(FermOpD.Umu, mu); + precisionChange(Umu_f,Umu_d); + PokeIndex(FermOpF.Umu, Umu_f, mu); + } + pickCheckerboard(Even,FermOpF.UmuEven,FermOpF.Umu); + pickCheckerboard(Odd ,FermOpF.UmuOdd ,FermOpF.Umu); + + //////////////////////////////////////////////////////////////////////////////////// + // Make a mixed precision conjugate gradient + //////////////////////////////////////////////////////////////////////////////////// +#if 1 + RealD delta=1.e-4; + std::cout << GridLogMessage << "Calling reliable update Conjugate Gradient" < MPCG(Tolerance,MaxInnerIterations*MaxOuterIterations,delta,SinglePrecGrid5,LinOpF,LinOpD); +#else + std::cout << GridLogMessage << "Calling mixed precision Conjugate Gradient" < MPCG(Tolerance,MaxInnerIterations,MaxOuterIterations,SinglePrecGrid5,LinOpF,LinOpD); +#endif + MPCG(src,psi); + } + }; + +NAMESPACE_END(Grid); + + +int main(int argc, char **argv) { + using namespace Grid; + + Grid_init(&argc, &argv); + int threads = GridThread::GetThreads(); + + // Typedefs to simplify notation + typedef WilsonImplR FermionImplPolicy; + typedef WilsonImplF FermionImplPolicyF; + + typedef MobiusFermionR FermionAction; + typedef MobiusFermionF FermionActionF; + typedef typename FermionAction::FermionField FermionField; + typedef typename FermionActionF::FermionField FermionFieldF; + + typedef Grid::XmlReader Serialiser; + + //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + IntegratorParameters MD; + // typedef GenericHMCRunner HMCWrapper; + // MD.name = std::string("Leap Frog"); + // typedef GenericHMCRunner HMCWrapper; + // MD.name = std::string("Force Gradient"); + typedef GenericHMCRunner HMCWrapper; + MD.name = std::string("MinimumNorm2"); + MD.MDsteps = 6; + MD.trajL = 1.0; + + HMCparameters HMCparams; + HMCparams.StartTrajectory = 1077; + HMCparams.Trajectories = 1; + HMCparams.NoMetropolisUntil= 0; + // "[HotStart, ColdStart, TepidStart, CheckpointStart]\n"; + // HMCparams.StartingType =std::string("ColdStart"); + HMCparams.StartingType =std::string("CheckpointStart"); + HMCparams.MD = MD; + HMCWrapper TheHMC(HMCparams); + + // Grid from the command line arguments --grid and --mpi + TheHMC.Resources.AddFourDimGrid("gauge"); // use default simd lanes decomposition + + CheckpointerParameters CPparams; + CPparams.config_prefix = "ckpoint_DDHMC_lat"; + CPparams.rng_prefix = "ckpoint_DDHMC_rng"; + CPparams.saveInterval = 1; + CPparams.format = "IEEE64BIG"; + TheHMC.Resources.LoadNerscCheckpointer(CPparams); + + RNGModuleParameters RNGpar; + RNGpar.serial_seeds = "1 2 3 4 5"; + RNGpar.parallel_seeds = "6 7 8 9 10"; + TheHMC.Resources.SetRNGSeeds(RNGpar); + + // Construct observables + // here there is too much indirection + typedef PlaquetteMod PlaqObs; + TheHMC.Resources.AddObservable(); + ////////////////////////////////////////////// + + const int Ls = 12; + RealD M5 = 1.8; + RealD b = 1.5; + RealD c = 0.5; + Real beta = 2.31; + // Real light_mass = 5.4e-4; + Real light_mass = 7.8e-4; + Real strange_mass = 0.02132; + Real pv_mass = 1.0; + // std::vector hasenbusch({ light_mass, 3.8e-3, 0.0145, 0.045, 0.108, 0.25, 0.51 , pv_mass }); + std::vector hasenbusch({ light_mass, 5e-3, 0.0145, 0.045, 0.108, 0.25, 0.51 , pv_mass }); + + // FIXME: + // Same in MC and MD + // Need to mix precision too + OneFlavourRationalParams SFRp; // Strange + SFRp.lo = 4.0e-3; + SFRp.hi = 90.0; + SFRp.MaxIter = 60000; + SFRp.tolerance= 1.0e-8; + SFRp.mdtolerance= 1.0e-6; + SFRp.degree = 12; + SFRp.precision= 50; + SFRp.BoundsCheckFreq=0; + + OneFlavourRationalParams OFRp; // Up/down + OFRp.lo = 2.0e-5; + OFRp.hi = 90.0; + OFRp.MaxIter = 60000; + OFRp.tolerance= 1.0e-8; + OFRp.mdtolerance= 1.0e-6; + // OFRp.degree = 20; converges + // OFRp.degree = 16; + OFRp.degree = 12; + OFRp.precision= 80; + OFRp.BoundsCheckFreq=0; + + auto GridPtr = TheHMC.Resources.GetCartesian(); + auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); + + typedef SchurDiagMooeeOperator LinearOperatorF; + typedef SchurDiagMooeeOperator LinearOperatorD; + typedef MixedPrecisionConjugateGradientOperatorFunction MxPCG; + + //////////////////////////////////////////////////////////////// + // Domain decomposed + //////////////////////////////////////////////////////////////// + Coordinate latt4 = GridPtr->GlobalDimensions(); + Coordinate mpi = GridPtr->ProcessorGrid(); + Coordinate shm; + + GlobalSharedMemory::GetShmDims(mpi,shm); + + Coordinate CommDim(Nd); + for(int d=0;d1 ? 1 : 0; + + Coordinate NonDirichlet(Nd+1,0); + Coordinate Dirichlet(Nd+1,0); + Dirichlet[1] = CommDim[0]*latt4[0]/mpi[0] * shm[0]; + Dirichlet[2] = CommDim[1]*latt4[1]/mpi[1] * shm[1]; + Dirichlet[3] = CommDim[2]*latt4[2]/mpi[2] * shm[2]; + Dirichlet[4] = CommDim[3]*latt4[3]/mpi[3] * shm[3]; + + Coordinate Block4(Nd); + Block4[0] = Dirichlet[1]; + Block4[1] = Dirichlet[2]; + Block4[2] = Dirichlet[3]; + Block4[3] = Dirichlet[4]; + + int Width=3; + TheHMC.Resources.SetMomentumFilter(new DDHMCFilter(Block4,Width)); + + ////////////////////////// + // Fermion Grids + ////////////////////////// + auto FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,GridPtr); + auto FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,GridPtr); + + Coordinate simdF = GridDefaultSimd(Nd,vComplexF::Nsimd()); + auto GridPtrF = SpaceTimeGrid::makeFourDimGrid(latt4,simdF,mpi); + auto GridRBPtrF = SpaceTimeGrid::makeFourDimRedBlackGrid(GridPtrF); + auto FGridF = SpaceTimeGrid::makeFiveDimGrid(Ls,GridPtrF); + auto FrbGridF = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,GridPtrF); + + IwasakiGaugeActionR GaugeAction(beta); + + // temporarily need a gauge field + LatticeGaugeField U(GridPtr); + LatticeGaugeFieldF UF(GridPtrF); + + std::cout << GridLogMessage << " Running the HMC "<< std::endl; + TheHMC.ReadCommandLine(argc,argv); // params on CML or from param file + TheHMC.initializeGaugeFieldAndRNGs(U); + + + // These lines are unecessary if BC are all periodic + std::vector boundary = {1,1,1,-1}; + FermionAction::ImplParams Params(boundary); + Params.dirichlet=NonDirichlet; + FermionAction::ImplParams ParamsDir(boundary); + ParamsDir.dirichlet=Dirichlet; + + // double StoppingCondition = 1e-14; + // double MDStoppingCondition = 1e-9; + double StoppingCondition = 1e-10; + double MDStoppingCondition = 1e-7; + double MDStoppingConditionLoose = 1e-6; + double MaxCGIterations = 300000; + ConjugateGradient CG(StoppingCondition,MaxCGIterations); + ConjugateGradient MDCG(MDStoppingCondition,MaxCGIterations); + + //////////////////////////////////// + // Collect actions + //////////////////////////////////// + ActionLevel Level1(1); + ActionLevel Level2(4); + ActionLevel Level3(8); + + //////////////////////////////////// + // Strange action + //////////////////////////////////// + FermionAction StrangeOp (U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,strange_mass,M5,b,c, Params); + FermionAction StrangePauliVillarsOp(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,pv_mass, M5,b,c, Params); + + FermionAction StrangeOpDir (U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,strange_mass,M5,b,c, ParamsDir); + FermionAction StrangePauliVillarsOpDir(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,pv_mass, M5,b,c, ParamsDir); + + OneFlavourEvenOddRatioRationalPseudoFermionAction StrangePseudoFermionBdy(StrangeOpDir,StrangeOp,SFRp); + OneFlavourEvenOddRatioRationalPseudoFermionAction StrangePseudoFermionLocal(StrangePauliVillarsOpDir,StrangeOpDir,SFRp); + OneFlavourEvenOddRatioRationalPseudoFermionAction StrangePseudoFermionPVBdy(StrangePauliVillarsOp,StrangePauliVillarsOpDir,SFRp); + Level1.push_back(&StrangePseudoFermionBdy); + Level2.push_back(&StrangePseudoFermionLocal); + Level1.push_back(&StrangePseudoFermionPVBdy); + + //////////////////////////////////// + // up down action + //////////////////////////////////// + std::vector light_den; + std::vector light_num; + std::vector dirichlet_den; + std::vector dirichlet_num; + + int n_hasenbusch = hasenbusch.size(); + light_den.push_back(light_mass); dirichlet_den.push_back(0); + for(int h=0;h Numerators; + std::vector NumeratorsF; + std::vector Denominators; + std::vector DenominatorsF; + std::vector *> Quotients; + +#define MIXED_PRECISION +#ifdef MIXED_PRECISION + std::vector *> Bdys; +#else + std::vector *> Bdys; +#endif + std::vector ActionMPCG; + std::vector MPCG; + + typedef SchurDiagMooeeOperator LinearOperatorF; + typedef SchurDiagMooeeOperator LinearOperatorD; + std::vector LinOpD; + std::vector LinOpF; + + for(int h=0;h(*Numerators[h],*Denominators[h],MDCG,CG)); + Quotients.push_back (new TwoFlavourEvenOddRatioPseudoFermionAction(*Numerators[h],*Denominators[h],*MPCG[h],*ActionMPCG[h],CG)); + } else { +#ifdef MIXED_PRECISION + Bdys.push_back( new OneFlavourEvenOddRatioRationalMixedPrecPseudoFermionAction( + *Numerators[h],*Denominators[h], + *NumeratorsF[h],*DenominatorsF[h], + OFRp, 500) ); +#else + Bdys.push_back( new OneFlavourEvenOddRatioRationalPseudoFermionAction(*Numerators[h],*Denominators[h],OFRp)); + Bdys.push_back( new OneFlavourEvenOddRatioRationalPseudoFermionAction(*Numerators[h],*Denominators[h],OFRp)); +#endif + } + } + + int nquo=Quotients.size(); + Level1.push_back(Bdys[0]); + Level1.push_back(Bdys[1]); + for(int h=0;h