From 3544965f5497a18c9e2c84f5ae9874e8b80adcce Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Thu, 7 Jul 2022 17:49:20 +0100 Subject: [PATCH 1/8] Stream doesn't work --- Grid/threads/Accelerator.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Grid/threads/Accelerator.h b/Grid/threads/Accelerator.h index d88c08b4..3be3bbe7 100644 --- a/Grid/threads/Accelerator.h +++ b/Grid/threads/Accelerator.h @@ -190,7 +190,7 @@ void Lambda6Apply(uint64_t num1, uint64_t num2, uint64_t num3, #define accelerator_barrier(dummy) \ { \ - cudaStreamSynchronize(cpuStream); \ + cudaDeviceSynchronize(); \ cudaError err = cudaGetLastError(); \ if ( cudaSuccess != err ) { \ printf("accelerator_barrier(): Cuda error %s \n", \ @@ -362,11 +362,11 @@ accelerator_inline int acceleratorSIMTlane(int Nsimd) { dim3 hip_blocks ((num1+nt-1)/nt,num2,1); \ if(hip_threads.x * hip_threads.y * hip_threads.z <= 64){ \ hipLaunchKernelGGL(LambdaApply64,hip_blocks,hip_threads, \ - 0,cpuStream, \ + 0,0/*cpuStream*/, \ num1,num2,nsimd, lambda); \ } else { \ hipLaunchKernelGGL(LambdaApply,hip_blocks,hip_threads, \ - 0,cpuStream, \ + 0,0/*cpuStream*/, \ num1,num2,nsimd, lambda); \ } \ } From c0f84824025c17a0bb97ece42833f671e3ce6770 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Thu, 7 Jul 2022 17:49:36 +0100 Subject: [PATCH 2/8] Remove SSC marks --- benchmarks/Benchmark_gparity.cc | 4 ---- 1 file changed, 4 deletions(-) diff --git a/benchmarks/Benchmark_gparity.cc b/benchmarks/Benchmark_gparity.cc index ce84ecbc..c10e7849 100644 --- a/benchmarks/Benchmark_gparity.cc +++ b/benchmarks/Benchmark_gparity.cc @@ -98,9 +98,7 @@ int main (int argc, char ** argv) std::cout<Barrier(); @@ -141,9 +139,7 @@ int main (int argc, char ** argv) std::cout<Barrier(); From b0fe664e9d74b634ca880c56095c8f90a34ca6aa Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sun, 10 Jul 2022 21:31:25 +0100 Subject: [PATCH 3/8] Better force log info --- Grid/qcd/action/ActionBase.h | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/Grid/qcd/action/ActionBase.h b/Grid/qcd/action/ActionBase.h index fa69b4e5..b8c81f99 100644 --- a/Grid/qcd/action/ActionBase.h +++ b/Grid/qcd/action/ActionBase.h @@ -42,6 +42,8 @@ public: bool is_smeared = false; RealD deriv_norm_sum; RealD deriv_max_sum; + RealD Fdt_norm_sum; + RealD Fdt_max_sum; int deriv_num; RealD deriv_us; RealD S_us; @@ -51,12 +53,17 @@ public: deriv_num=0; deriv_norm_sum = deriv_max_sum=0.0; } - void deriv_log(RealD nrm, RealD max) { deriv_max_sum+=max; deriv_norm_sum+=nrm; deriv_num++;} - RealD deriv_max_average(void) { return deriv_max_sum/deriv_num; }; - RealD deriv_norm_average(void) { return deriv_norm_sum/deriv_num; }; + void deriv_log(RealD nrm, RealD max,RealD Fdt_nrm,RealD Fdt_max) { + deriv_max_sum+=max; deriv_norm_sum+=nrm; + Fdt_max_sum+=Fdt_max; Fdt_norm_sum+=Fdt_nrm; deriv_num++; + } + RealD deriv_max_average(void) { return deriv_max_sum/deriv_num; }; + RealD deriv_norm_average(void) { return deriv_norm_sum/deriv_num; }; + RealD Fdt_max_average(void) { return Fdt_max_sum/deriv_num; }; + RealD Fdt_norm_average(void) { return Fdt_norm_sum/deriv_num; }; RealD deriv_timer(void) { return deriv_us; }; - RealD S_timer(void) { return deriv_us; }; - RealD refresh_timer(void) { return deriv_us; }; + RealD S_timer(void) { return S_us; }; + RealD refresh_timer(void) { return refresh_us; }; void deriv_timer_start(void) { deriv_us-=usecond(); } void deriv_timer_stop(void) { deriv_us+=usecond(); } void refresh_timer_start(void) { refresh_us-=usecond(); } From 1f907d330da8496d369877d7c93bd1379178acdc Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sun, 10 Jul 2022 21:31:48 +0100 Subject: [PATCH 4/8] Different default params for dirichlet --- Grid/qcd/action/ActionParams.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Grid/qcd/action/ActionParams.h b/Grid/qcd/action/ActionParams.h index b203f27a..274ff318 100644 --- a/Grid/qcd/action/ActionParams.h +++ b/Grid/qcd/action/ActionParams.h @@ -39,7 +39,7 @@ struct GparityWilsonImplParams { Coordinate twists; //mu=Nd-1 is assumed to be the time direction and a twist value of 1 indicates antiperiodic BCs Coordinate dirichlet; // Blocksize of dirichlet BCs - GparityWilsonImplParams() : twists(Nd, 0), dirichlet(Nd, 0) {}; + GparityWilsonImplParams() : twists(Nd, 0) { dirichlet.resize(0); }; }; struct WilsonImplParams { @@ -48,13 +48,13 @@ struct WilsonImplParams { AcceleratorVector twist_n_2pi_L; AcceleratorVector boundary_phases; WilsonImplParams() { - dirichlet.resize(Nd,0); + dirichlet.resize(0); boundary_phases.resize(Nd, 1.0); twist_n_2pi_L.resize(Nd, 0.0); }; WilsonImplParams(const AcceleratorVector phi) : boundary_phases(phi), overlapCommsCompute(false) { twist_n_2pi_L.resize(Nd, 0.0); - dirichlet.resize(Nd,0); + dirichlet.resize(0); } }; @@ -62,7 +62,7 @@ struct StaggeredImplParams { Coordinate dirichlet; // Blocksize of dirichlet BCs StaggeredImplParams() { - dirichlet.resize(Nd,0); + dirichlet.resize(0); }; }; From 58182fe3455a0d78de6074e6c5e38bdd88e67751 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sun, 10 Jul 2022 21:32:58 +0100 Subject: [PATCH 5/8] Different approach to default dirichlet params --- Grid/stencil/Stencil.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Grid/stencil/Stencil.h b/Grid/stencil/Stencil.h index eb73ba5f..9d005289 100644 --- a/Grid/stencil/Stencil.h +++ b/Grid/stencil/Stencil.h @@ -711,7 +711,9 @@ public: this->_comms_recv.resize(npoints); this->same_node.resize(npoints); - if ( p.dirichlet.size() ) DirichletBlock(p.dirichlet); // comms send/recv set up + if ( p.dirichlet.size() ==0 ) p.dirichlet.resize(grid->Nd(),0); + + DirichletBlock(p.dirichlet); // comms send/recv set up _unified_buffer_size=0; surface_list.resize(0); From 177b1a7ec630bc02f661e2e10299b73b47e9281a Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sun, 10 Jul 2022 21:34:10 +0100 Subject: [PATCH 6/8] 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 Date: Sun, 10 Jul 2022 21:35:18 +0100 Subject: [PATCH 7/8] MixedPrec support --- .../OneFlavourEvenOddRationalRatio.h | 30 +++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h b/Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h index 1b36ae0f..078bf845 100644 --- a/Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h +++ b/Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h @@ -67,6 +67,36 @@ NAMESPACE_BEGIN(Grid); virtual std::string action_name(){return "OneFlavourEvenOddRatioRationalPseudoFermionAction";} }; + template + class OneFlavourEvenOddRatioRationalMixedPrecPseudoFermionAction : public GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction { + public: + typedef OneFlavourRationalParams Params; + private: + static RationalActionParams transcribe(const Params &in){ + RationalActionParams out; + out.inv_pow = 2; + out.lo = in.lo; + out.hi = in.hi; + out.MaxIter = in.MaxIter; + out.action_tolerance = out.md_tolerance = in.tolerance; + out.action_degree = out.md_degree = in.degree; + out.precision = in.precision; + out.BoundsCheckFreq = in.BoundsCheckFreq; + return out; + } + + public: + OneFlavourEvenOddRatioRationalMixedPrecPseudoFermionAction(FermionOperator &_NumOp, + FermionOperator &_DenOp, + FermionOperator &_NumOpF, + FermionOperator &_DenOpF, + const Params & p, Integer ReliableUpdateFreq + ) : + GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction(_NumOp, _DenOp,_NumOpF, _DenOpF, transcribe(p),ReliableUpdateFreq){} + + virtual std::string action_name(){return "OneFlavourEvenOddRatioRationalPseudoFermionAction";} + }; + NAMESPACE_END(Grid); #endif From fab50c57d9d14be6cc0f32de93160c94eb563fec Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Mon, 11 Jul 2022 18:42:27 +0100 Subject: [PATCH 8/8] More loggin --- Grid/qcd/hmc/integrators/Integrator.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Grid/qcd/hmc/integrators/Integrator.h b/Grid/qcd/hmc/integrators/Integrator.h index 5bdef14d..33a77f32 100644 --- a/Grid/qcd/hmc/integrators/Integrator.h +++ b/Grid/qcd/hmc/integrators/Integrator.h @@ -153,7 +153,7 @@ protected: Real force_max = std::sqrt(maxLocalNorm2(force)); Real impulse_max = force_max * ep * HMC_MOMENTUM_DENOMINATOR; - as[level].actions.at(a)->deriv_log(force_abs,force_max); + as[level].actions.at(a)->deriv_log(force_abs,force_max,impulse_abs,impulse_max); std::cout << GridLogIntegrator<< "["<deriv_max_average() <<" norm " << as[level].actions.at(actionID)->deriv_norm_average() + <<" Fdt max " << as[level].actions.at(actionID)->Fdt_max_average() + <<" norm " << as[level].actions.at(actionID)->Fdt_norm_average() <<" calls " << as[level].actions.at(actionID)->deriv_num << std::endl; }