From e3c18ce87259232eb1663c22bb4a843457a6b4cf Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sat, 15 May 2021 09:06:58 -0400 Subject: [PATCH] Adding boundary det --- HMC/Mobius2f_DDHMC.cc | 217 +++++++----------------------------------- 1 file changed, 35 insertions(+), 182 deletions(-) diff --git a/HMC/Mobius2f_DDHMC.cc b/HMC/Mobius2f_DDHMC.cc index 7bd7b2fb..ebd99f26 100644 --- a/HMC/Mobius2f_DDHMC.cc +++ b/HMC/Mobius2f_DDHMC.cc @@ -30,12 +30,15 @@ directory #include #include #include +#include + +#include +#include +#include + //#include //#include -#ifdef GRID_DEFAULT_PRECISION_DOUBLE -#define MIXED_PRECISION -#endif NAMESPACE_BEGIN(Grid); @@ -46,126 +49,6 @@ NAMESPACE_BEGIN(Grid); * -- 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.) - { - /* 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); - - //////////////////////////////////////////////////////////////////////////////////// - // Could test to make sure that LinOpF and LinOpD agree to single prec? - //////////////////////////////////////////////////////////////////////////////////// - /* - GridBase *Fgrid = psi._grid; - FieldD tmp2(Fgrid); - FieldD tmp1(Fgrid); - LinOpU.Op(src,tmp1); - LinOpD.Op(src,tmp2); - std::cout << " Double gauge field "<< norm2(FermOpD.Umu)< MPCG(Tolerance,MaxInnerIterations,MaxOuterIterations,SinglePrecGrid5,LinOpF,LinOpD); - std::cout << GridLogMessage << "Calling mixed precision Conjugate Gradient" < DirichletFermion; - typedef DirichletFermionOperator DirichletFermionF; typedef MobiusEOFAFermionR FermionEOFAAction; - typedef MobiusEOFAFermionF FermionEOFAActionF; typedef typename FermionAction::FermionField FermionField; - typedef typename FermionActionF::FermionField FermionFieldF; typedef Grid::XmlReader Serialiser; @@ -204,10 +83,10 @@ int main(int argc, char **argv) { HMCparameters HMCparams; HMCparams.StartTrajectory = 0; HMCparams.Trajectories = 1000; - HMCparams.NoMetropolisUntil= 0; + HMCparams.NoMetropolisUntil= 10; // "[HotStart, ColdStart, TepidStart, CheckpointStart]\n"; - // HMCparams.StartingType =std::string("ColdStart"); HMCparams.StartingType =std::string("ColdStart"); + // HMCparams.StartingType =std::string("CheckpointStart"); HMCparams.MD = MD; HMCWrapper TheHMC(HMCparams); @@ -227,7 +106,7 @@ int main(int argc, char **argv) { TheHMC.Resources.SetRNGSeeds(RNGpar); // Momentum Dirichlet - Coordinate Block({4,4,4,4}); + Coordinate Block({8,8,8,4}); TheHMC.Resources.SetMomentumFilter(new DDHMCFilter(Block,0)); // Construct observables @@ -265,13 +144,10 @@ int main(int argc, char **argv) { // 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); - FermionActionF::ImplParams ParamsF(boundary); - double ActionStoppingCondition = 1e-10; double DerivativeStoppingCondition = 1e-8; @@ -283,46 +159,8 @@ int main(int argc, char **argv) { ActionLevel Level1(1); ActionLevel Level2(8); - //////////////////////////////////// - // Strange action - //////////////////////////////////// - typedef SchurDiagMooeeOperator LinearOperatorF; - typedef SchurDiagMooeeOperator LinearOperatorD; - - typedef SchurDiagMooeeDagOperator LinearOperatorDagF; - typedef SchurDiagMooeeDagOperator LinearOperatorDagD; - - typedef SchurDiagMooeeOperator LinearOperatorEOFAF; - typedef SchurDiagMooeeOperator LinearOperatorEOFAD; - - typedef MixedPrecisionConjugateGradientOperatorFunction MxPCG; - typedef MixedPrecisionConjugateGradientOperatorFunction MxDagPCG; - typedef MixedPrecisionConjugateGradientOperatorFunction MxPCG_EOFA; - - // DJM: setup for EOFA ratio (Mobius) - OneFlavourRationalParams OFRp; - OFRp.lo = 0.1; - OFRp.hi = 25.0; - OFRp.MaxIter = 10000; - OFRp.tolerance= 1.0e-9; - OFRp.degree = 14; - OFRp.precision= 50; - - MobiusEOFAFermionR Strange_Op_L (U , *FGrid , *FrbGrid , *GridPtr , *GridRBPtr , strange_mass, strange_mass, pv_mass, 0.0, -1, M5, b, c); - MobiusEOFAFermionF Strange_Op_LF(UF, *FGridF, *FrbGridF, *GridPtrF, *GridRBPtrF, strange_mass, strange_mass, pv_mass, 0.0, -1, M5, b, c); - MobiusEOFAFermionR Strange_Op_R (U , *FGrid , *FrbGrid , *GridPtr , *GridRBPtr , pv_mass, strange_mass, pv_mass, -1.0, 1, M5, b, c); - MobiusEOFAFermionF Strange_Op_RF(UF, *FGridF, *FrbGridF, *GridPtrF, *GridRBPtrF, pv_mass, strange_mass, pv_mass, -1.0, 1, M5, b, c); - ConjugateGradient ActionCG(ActionStoppingCondition,MaxCGIterations); ConjugateGradient DerivativeCG(DerivativeStoppingCondition,MaxCGIterations); - ExactOneFlavourRatioPseudoFermionAction - EOFA(Strange_Op_L, Strange_Op_R, - ActionCG, - ActionCG, ActionCG, - DerivativeCG, DerivativeCG, - OFRp, true); - - // Level1.push_back(&EOFA); //////////////////////////////////// // up down action @@ -346,30 +184,38 @@ int main(int argc, char **argv) { ////////////////////////////////////////////////////////////// std::vector Numerators; std::vector Denominators; + std::vector PeriNumerators; + std::vector PeriDenominators; std::vector DirichletNumerators; std::vector DirichletDenominators; std::vector *> Quotients; - // std::vector *> Quotients; - std::vector ActionMPCG; - std::vector MPCG; - std::vector MPCGdag; - std::vector DenominatorsF; - std::vector LinOpD; - std::vector LinOpF; - std::vector LinOpDagD; - std::vector LinOpDagF; - + std::vector *> BoundaryNumerators; + std::vector *> BoundaryDenominators; + std::vector *> BoundaryQuotients; for(int h=0;h + (*PeriNumerators[h], + *DirichletNumerators[h], + ActionCG,Block)); + BoundaryDenominators.push_back (new SchurFactoredFermionOperator + (*PeriDenominators[h], + *DirichletDenominators[h], + ActionCG,Block)); + //////////////////////////////////////////////////////////////////////////// // Standard CG for 2f force //////////////////////////////////////////////////////////////////////////// @@ -380,10 +226,17 @@ int main(int argc, char **argv) { ActionCG, ActionCG)); + BoundaryQuotients.push_back(new + DomainDecomposedBoundaryTwoFlavourRatioPseudoFermion + (*BoundaryNumerators[h], + *BoundaryDenominators[h], + ActionCG,ActionCG)); + } for(int h=0;h