diff --git a/HMC/Mobius2f_DDHMC.cc b/HMC/Mobius2f_DDHMC.cc index 7815c25a..96e19bde 100644 --- a/HMC/Mobius2f_DDHMC.cc +++ b/HMC/Mobius2f_DDHMC.cc @@ -62,6 +62,7 @@ int main(int argc, char **argv) { // Typedefs to simplify notation typedef WilsonImplR FermionImplPolicy; typedef MobiusFermionR FermionAction; + typedef typename FermionAction::Impl_t Fimpl; typedef DirichletFermionOperator DirichletFermion; typedef MobiusEOFAFermionR FermionEOFAAction; @@ -81,8 +82,8 @@ int main(int argc, char **argv) { MD.trajL = 1.0; HMCparameters HMCparams; - HMCparams.StartTrajectory = 114; - HMCparams.Trajectories = 1000; + HMCparams.StartTrajectory = 26; + HMCparams.Trajectories = 1; HMCparams.NoMetropolisUntil= 0; // "[HotStart, ColdStart, TepidStart, CheckpointStart]\n"; // HMCparams.StartingType =std::string("ColdStart"); @@ -106,9 +107,9 @@ int main(int argc, char **argv) { TheHMC.Resources.SetRNGSeeds(RNGpar); // Momentum Dirichlet - Coordinate Block({16,16,16,4}); + Coordinate Block({16,16,16,16}); - // TheHMC.Resources.SetMomentumFilter(new DirichletFilter(Block)); + // TheHMC.Resources.SetMomentumFilter(new DirichletFilter(Block)); TheHMC.Resources.SetMomentumFilter(new DDHMCFilter(Block,1)); // Construct observables // here there is too much indirection @@ -125,7 +126,7 @@ int main(int argc, char **argv) { RealD b = 1.0; RealD c = 0.0; - std::vector hasenbusch({ 0.1, 0.3, 0.6 }); + std::vector hasenbusch({ 0.04, 0.2 }); auto GridPtr = TheHMC.Resources.GetCartesian(); auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); @@ -149,9 +150,11 @@ int main(int argc, char **argv) { // These lines are unecessary if BC are all periodic std::vector boundary = {1,1,1,-1}; FermionAction::ImplParams Params(boundary); + FermionAction::ImplParams DirichletParams(boundary); + DirichletParams.locally_periodic=true; double ActionStoppingCondition = 1e-10; - double DerivativeStoppingCondition = 1e-8; + double DerivativeStoppingCondition = 1e-10; double MaxCGIterations = 30000; //////////////////////////////////// @@ -190,29 +193,29 @@ int main(int argc, char **argv) { std::vector DDenominators; std::vector DirichletNumerators; std::vector DirichletDenominators; - - std::vector *> Quotients; - std::vector *> BoundaryNumerators; - std::vector *> BoundaryDenominators; - std::vector *> BoundaryQuotients; - for(int h=0;h *> BoundaryNumerators; + std::vector *> BoundaryDenominators; - std::cout << GridLogMessage << " 2f quotient Action "<< light_num[h] << " / " << light_den[h]<< std::endl; + std::vector *> Quotients; + std::vector *> BoundaryQuotients; + std::vector *> BoundaryNums; + std::vector *> BoundaryDens; + for(int h=0;h + BoundaryNumerators.push_back (new SchurFactoredFermionOperator (*PeriNumerators[h], *DirichletNumerators[h], ActionCG,Block)); - BoundaryDenominators.push_back (new SchurFactoredFermionOperator + + BoundaryDenominators.push_back (new SchurFactoredFermionOperator (*PeriDenominators[h], *DirichletDenominators[h], ActionCG,Block)); @@ -220,24 +223,30 @@ int main(int argc, char **argv) { //////////////////////////////////////////////////////////////////////////// // Standard CG for 2f force //////////////////////////////////////////////////////////////////////////// + std::cout << GridLogMessage << " 2f quotient Action "<< light_num[h] << " / " << light_den[h]<< std::endl; + Quotients.push_back (new - TwoFlavourEvenOddRatioPseudoFermionAction + TwoFlavourEvenOddRatioPseudoFermionAction (*DirichletNumerators[h], *DirichletDenominators[h], ActionCG, ActionCG)); + std::cout << GridLogMessage << " 2f Boundary quotient Action "<< light_num[h] << " / " << light_den[h]<< std::endl; + BoundaryQuotients.push_back(new - DomainDecomposedBoundaryTwoFlavourRatioPseudoFermion + DomainDecomposedBoundaryTwoFlavourRatioPseudoFermion (*BoundaryNumerators[h], *BoundaryDenominators[h], ActionCG,ActionCG)); - + } for(int h=0;h + +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 +#include +#include +#include +#include + +#include +#include +#include + + + +NAMESPACE_BEGIN(Grid); + +#define MIXED_PRECISION + /* + * Need a plan for gauge field update for mixed precision in HMC (2x speed up) + * -- Store the single prec action operator. + * -- Clone the gauge field from the operator function argument. + * -- Build the mixed precision operator dynamically from the passed operator and single prec clone. + */ + 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.) + { }; + + void operator()(LinearOperatorBase &LinOpU, const FieldD &src, FieldD &psi) + { + + SchurOperatorD * SchurOpU = static_cast(&LinOpU); + + // Assumption made in code to extract gauge field + // We could avoid storing LinopD reference alltogether ? + assert(&(SchurOpU->_Mat)==&(LinOpD._Mat)); + + //////////////////////////////////////////////////////////////////////////////////// + // Must snarf a single precision copy of the gauge field in Linop_d argument + //////////////////////////////////////////////////////////////////////////////////// + typedef typename FermionOperatorF::GaugeField GaugeFieldF; + typedef typename FermionOperatorD::GaugeField GaugeFieldD; + typedef typename FermionOperatorF::GaugeLinkField GaugeLinkFieldF; + typedef typename FermionOperatorD::GaugeLinkField GaugeLinkFieldD; + + GridBase * GridPtrF = SinglePrecGrid4; + GridBase * GridPtrD = FermOpD.GaugeGrid(); + + //////////////////////////////////////////////////////////////////////////////////// + // Moving this to a Clone method of fermion operator would allow to duplicate the + // physics parameters and decrease gauge field copies + //////////////////////////////////////////////////////////////////////////////////// + auto &Umu_d = FermOpD.GetDoubledGaugeField(); + auto &Umu_f = FermOpF.GetDoubledGaugeField(); + auto &Umu_fe= FermOpF.GetDoubledGaugeFieldE(); + auto &Umu_fo= FermOpF.GetDoubledGaugeFieldO(); + precisionChange(Umu_f,Umu_d); + pickCheckerboard(Even,Umu_fe,Umu_f); + pickCheckerboard(Odd ,Umu_fo,Umu_f); + + //////////////////////////////////////////////////////////////////////////////////// + // Could test to make sure that LinOpF and LinOpD agree to single prec? + //////////////////////////////////////////////////////////////////////////////////// + /* + FieldD srcD(FermOpD.FermionRedBlackGrid()); + FieldD tmpD(FermOpD.FermionRedBlackGrid()); + FieldF tmpF(FermOpF.FermionRedBlackGrid()); + FieldF srcF(FermOpF.FermionRedBlackGrid()); + srcD = 1.0; + precisionChange(srcF,srcD); + std::cout << GridLogMessage << " Prec Src "< MPCG(Tolerance,MaxInnerIterations,MaxOuterIterations,SinglePrecGrid5,LinOpF,LinOpD); + std::cout << GridLogMessage << "Calling mixed precision Conjugate Gradient" < DirichletFermion; + typedef DirichletFermionOperator DirichletFermionF; + + typedef MobiusEOFAFermionR FermionEOFAAction; + 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 = 12; + MD.trajL = 1.0; + + HMCparameters HMCparams; + HMCparams.StartTrajectory = 26; + HMCparams.Trajectories = 1000; + 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_EOFA4D_lat"; + CPparams.rng_prefix = "ckpoint_EOFA4D_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); + + // Momentum Dirichlet + Coordinate Block({16,16,16,16}); + + // TheHMC.Resources.SetMomentumFilter(new DirichletFilter(Block)); + TheHMC.Resources.SetMomentumFilter(new DDHMCFilter(Block,1)); + // Construct observables + // here there is too much indirection + typedef PlaquetteMod PlaqObs; + TheHMC.Resources.AddObservable(); + ////////////////////////////////////////////// + + const int Ls = 16; + Real beta = 2.13; + Real light_mass = 0.01; + Real strange_mass = 0.04; + Real pv_mass = 1.0; + RealD M5 = 1.8; + RealD b = 1.0; + RealD c = 0.0; + + std::vector hasenbusch({ 0.016, 0.04, 0.4 }); + + auto GridPtr = TheHMC.Resources.GetCartesian(); + auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); + auto FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,GridPtr); + auto FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,GridPtr); + + Coordinate latt = GridDefaultLatt(); + Coordinate mpi = GridDefaultMpi(); + Coordinate simdF = GridDefaultSimd(Nd,vComplexF::Nsimd()); + Coordinate simdD = GridDefaultSimd(Nd,vComplexD::Nsimd()); + auto GridPtrF = SpaceTimeGrid::makeFourDimGrid(latt,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); + + // These lines are unecessary if BC are all periodic + std::vector boundary = {1,1,1,-1}; + FermionAction::ImplParams Params(boundary); + FermionAction::ImplParams DirichletParams(boundary); + DirichletParams.locally_periodic=true; + + double ActionStoppingCondition = 1e-10; + double DerivativeStoppingCondition = 1e-8; + double MaxCGIterations = 30000; + + //////////////////////////////////// + // Collect actions + //////////////////////////////////// + ActionLevel Level1(1); + ActionLevel Level2(8); + + ConjugateGradient ActionCG(ActionStoppingCondition,MaxCGIterations); + ConjugateGradient DerivativeCG(DerivativeStoppingCondition,MaxCGIterations); + + //////////////////////////////////// + // up down action + //////////////////////////////////// + std::vector light_den; + std::vector light_num; + + int n_hasenbusch = hasenbusch.size(); + light_den.push_back(light_mass); + for(int h=0;h LinearOperatorF; + typedef SchurDiagMooeeOperator LinearOperatorD; + // typedef SchurDiagMooeeDagOperator LinearOperatorDagF; + // typedef SchurDiagMooeeDagOperator LinearOperatorDagD; + typedef MixedPrecisionConjugateGradientOperatorFunction MxPCG; + // typedef MixedPrecisionConjugateGradientOperatorFunction MxDagPCG; + + std::vector PeriNumerators; + std::vector PeriDenominators; + + std::vector DNumerators; + std::vector DDenominators; + std::vector DDenominatorsF; + std::vector DirichletNumerators; + std::vector DirichletDenominators; + std::vector DirichletDenominatorsF; + + std::vector *> Quotients; + std::vector *> BoundaryNumerators; + std::vector *> BoundaryDenominators; + std::vector *> BoundaryQuotients; + + std::vector ActionMPCG; + std::vector MPCG; + std::vector DenominatorsF; + std::vector LinOpD; + std::vector LinOpF; + + for(int h=0;h + (*PeriNumerators[h], + *DirichletNumerators[h], + ActionCG,Block)); + BoundaryDenominators.push_back (new SchurFactoredFermionOperator + (*PeriDenominators[h], + *DirichletDenominators[h], + ActionCG,Block)); + + // Dirichlet Schur even odd MpsDagMpc operators on local domains + LinOpD.push_back(new LinearOperatorD(*DirichletDenominators[h])); + LinOpF.push_back(new LinearOperatorF(*DirichletDenominatorsF[h])); + + const int MX_inner = 1000; + + MPCG.push_back(new MxPCG(DerivativeStoppingCondition, + MX_inner, + MaxCGIterations, + GridPtrF, + FrbGridF, + *DirichletDenominatorsF[h],*DirichletDenominators[h], + *LinOpF[h], *LinOpD[h]) ); + + ActionMPCG.push_back(new MxPCG(ActionStoppingCondition, + MX_inner, + MaxCGIterations, + GridPtrF, + FrbGridF, + *DirichletDenominatorsF[h],*DirichletDenominators[h], + *LinOpF[h], *LinOpD[h]) ); + + //////////////////////////////////////////////////////////////////////////// + // Standard CG for 2f force + //////////////////////////////////////////////////////////////////////////// + // Quotients.push_back (new TwoFlavourEvenOddRatioPseudoFermionAction(*Numerators[h],*Denominators[h],*MPCG[h],*ActionMPCG[h],ActionCG)); + Quotients.push_back (new + TwoFlavourEvenOddRatioPseudoFermionAction + (*DirichletNumerators[h], + *DirichletDenominators[h], + *ActionMPCG[h], + *ActionMPCG[h], + ActionCG)); + + BoundaryQuotients.push_back(new + DomainDecomposedBoundaryTwoFlavourRatioPseudoFermion + (*BoundaryNumerators[h], + *BoundaryDenominators[h], + ActionCG,ActionCG)); + + } + + for(int h=0;h