From fd933420c65b927da2c12f6c355e26f7d610db1f Mon Sep 17 00:00:00 2001 From: Christopher Kelly Date: Wed, 22 Jun 2022 10:27:48 -0400 Subject: [PATCH] Imported changes from feature/gparity_HMC branch: Added a bounds-check function for the RHMC with arbitrary power Added a pseudofermion action for the rational ratio with an arbitrary power and a mixed-precision variant of the same. The existing one-flavor rational ratio class now uses the general class under the hood To support testing of the two-flavor even-odd ratio pseudofermion, separated the functionality of generating the random field and performing the heatbath step, and added a method to obtain the pseudofermion field Added a new HMC runner start type: CheckpointStartReseed, which reseeds the RNG from scratch, allowing for the creation of new evolution streams from an existing checkpoint. Added log output of seeds used when the RNG is seeded. EOFA changes: To support mixed-precision inversion, generalized the class to maintain a separate solver for the L and R operators in the heatbath (separate solvers are already implemented for the other stages) To support mixed-precision, the action of setting the operator shift coefficients is now maintained in a virtual function. A derived class for mixed-precision solvers ensures the coefficients are applied to both the double and single-prec operators The ||^2 of the random source is now stored by the heatbath and compared to the initial action when it is computed. These should be equal but may differ if the rational bounds are not chosen correctly, hence serving as a useful and free test Fixed calculation of M_eofa (previously incomplete and #if'd out) Added functionality to compute M_eofa^-1 to complement the calculation of M_eofa (both are equally expensive!) To support testing, separated the functionality of generating the random field and performing the heatbath step, and added a method to obtain the pseudofermion field Added a test program which computes the G-parity force using the 1 and 2 flavor implementations and compares the result. Test supports DWF, EOFA and DSDR actions, chosen by a command line option. The Mobius EOFA force test now also checks the rational approximation used for the heatbath Added a test program for the mixed precision EOFA compared to the double-prec implementation, G-parity HMC test now applied GPBC in the y direction and not the t direction (GPBC in t are no longer supported) and checkpoints after every configuration Added a test program which computes the two-flavor G-parity action (via RHMC) with both the 1 and 2 flavor implementations and checks they agree Added a test program to check the implementation of M_eofa^{-1} --- Grid/qcd/action/ActionParams.h | 52 +- Grid/qcd/action/pseudofermion/Bounds.h | 60 ++- .../pseudofermion/ExactOneFlavourRatio.h | 229 +++++++-- .../GeneralEvenOddRationalRatio.h | 372 +++++++++++++++ .../GeneralEvenOddRationalRatioMixedPrec.h | 93 ++++ .../OneFlavourEvenOddRationalRatio.h | 262 +--------- Grid/qcd/action/pseudofermion/PseudoFermion.h | 2 + .../pseudofermion/TwoFlavourEvenOddRatio.h | 28 +- Grid/qcd/hmc/GenericHMCrunner.h | 12 +- Grid/qcd/hmc/HMCModules.h | 2 + Grid/qcd/hmc/integrators/Integrator.h | 6 +- tests/forces/Test_gpdwf_force_1f_2f.cc | 446 ++++++++++++++++++ tests/forces/Test_mobius_force_eofa.cc | 42 ++ .../forces/Test_mobius_gparity_eofa_mixed.cc | 233 +++++++++ tests/hmc/Test_action_dwf_gparity2fvs1f.cc | 257 ++++++++++ tests/hmc/Test_hmc_GparityIwasakiGauge.cc | 4 +- tests/solver/Test_eofa_inv.cc | 125 +++++ 17 files changed, 1920 insertions(+), 305 deletions(-) create mode 100644 Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatio.h create mode 100644 Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatioMixedPrec.h create mode 100644 tests/forces/Test_gpdwf_force_1f_2f.cc create mode 100644 tests/forces/Test_mobius_gparity_eofa_mixed.cc create mode 100644 tests/hmc/Test_action_dwf_gparity2fvs1f.cc create mode 100644 tests/solver/Test_eofa_inv.cc diff --git a/Grid/qcd/action/ActionParams.h b/Grid/qcd/action/ActionParams.h index 6a3f053a..13f2e594 100644 --- a/Grid/qcd/action/ActionParams.h +++ b/Grid/qcd/action/ActionParams.h @@ -37,6 +37,7 @@ NAMESPACE_BEGIN(Grid); // These can move into a params header and be given MacroMagic serialisation struct GparityWilsonImplParams { Coordinate twists; + //mu=Nd-1 is assumed to be the time direction and a twist value of 1 indicates antiperiodic BCs GparityWilsonImplParams() : twists(Nd, 0) {}; }; @@ -66,7 +67,8 @@ struct StaggeredImplParams { RealD, mdtolerance, int, degree, int, precision, - int, BoundsCheckFreq); + int, BoundsCheckFreq, + RealD, BoundsCheckTol); // MaxIter and tolerance, vectors?? @@ -78,7 +80,8 @@ struct StaggeredImplParams { int _degree = 10, int _precision = 64, int _BoundsCheckFreq=20, - RealD mdtol = 1.0e-6) + RealD mdtol = 1.0e-6, + double _BoundsCheckTol=1e-6) : lo(_lo), hi(_hi), MaxIter(_maxit), @@ -86,9 +89,52 @@ struct StaggeredImplParams { mdtolerance(mdtol), degree(_degree), precision(_precision), - BoundsCheckFreq(_BoundsCheckFreq){}; + BoundsCheckFreq(_BoundsCheckFreq), + BoundsCheckTol(_BoundsCheckTol){}; }; + /*Action parameters for the generalized rational action + The approximation is for (M^dag M)^{1/inv_pow} + where inv_pow is the denominator of the fractional power. + Default inv_pow=2 for square root, making this equivalent to + the OneFlavourRational action + */ + struct RationalActionParams : Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(RationalActionParams, + int, inv_pow, + RealD, lo, //low eigenvalue bound of rational approx + RealD, hi, //high eigenvalue bound of rational approx + int, MaxIter, //maximum iterations in msCG + RealD, action_tolerance, //msCG tolerance in action evaluation + int, action_degree, //rational approx tolerance in action evaluation + RealD, md_tolerance, //msCG tolerance in MD integration + int, md_degree, //rational approx tolerance in MD integration + int, precision, //precision of floating point arithmetic + int, BoundsCheckFreq); //frequency the approximation is tested (with Metropolis degree/tolerance); 0 disables the check + // constructor + RationalActionParams(int _inv_pow = 2, + RealD _lo = 0.0, + RealD _hi = 1.0, + int _maxit = 1000, + RealD _action_tolerance = 1.0e-8, + int _action_degree = 10, + RealD _md_tolerance = 1.0e-8, + int _md_degree = 10, + int _precision = 64, + int _BoundsCheckFreq=20) + : inv_pow(_inv_pow), + lo(_lo), + hi(_hi), + MaxIter(_maxit), + action_tolerance(_action_tolerance), + action_degree(_action_degree), + md_tolerance(_md_tolerance), + md_degree(_md_degree), + precision(_precision), + BoundsCheckFreq(_BoundsCheckFreq){}; + }; + + NAMESPACE_END(Grid); #endif diff --git a/Grid/qcd/action/pseudofermion/Bounds.h b/Grid/qcd/action/pseudofermion/Bounds.h index b9621f24..8864b1d7 100644 --- a/Grid/qcd/action/pseudofermion/Bounds.h +++ b/Grid/qcd/action/pseudofermion/Bounds.h @@ -65,13 +65,65 @@ NAMESPACE_BEGIN(Grid); X=X-Y; RealD Nd = norm2(X); std::cout << "************************* "< void InversePowerBoundsCheck(int inv_pow, + int MaxIter,double tol, + LinearOperatorBase &HermOp, + Field &GaussNoise, + MultiShiftFunction &ApproxNegPow) + { + GridBase *FermionGrid = GaussNoise.Grid(); + + Field X(FermionGrid); + Field Y(FermionGrid); + Field Z(FermionGrid); + + Field tmp1(FermionGrid), tmp2(FermionGrid); + + X=GaussNoise; + RealD Nx = norm2(X); + + ConjugateGradientMultiShift msCG(MaxIter,ApproxNegPow); + + tmp1 = X; + + Field* in = &tmp1; + Field* out = &tmp2; + for(int i=0;i class ExactOneFlavourRatioPseudoFermionAction : public Action { @@ -57,37 +61,60 @@ NAMESPACE_BEGIN(Grid); bool use_heatbath_forecasting; AbstractEOFAFermion& Lop; // the basic LH operator AbstractEOFAFermion& Rop; // the basic RH operator - SchurRedBlackDiagMooeeSolve SolverHB; + SchurRedBlackDiagMooeeSolve SolverHBL; + SchurRedBlackDiagMooeeSolve SolverHBR; SchurRedBlackDiagMooeeSolve SolverL; SchurRedBlackDiagMooeeSolve SolverR; SchurRedBlackDiagMooeeSolve DerivativeSolverL; SchurRedBlackDiagMooeeSolve DerivativeSolverR; FermionField Phi; // the pseudofermion field for this trajectory + RealD norm2_eta; //|eta|^2 where eta is the random gaussian field used to generate the pseudofermion field + bool initial_action; //true for the first call to S after refresh, for which the identity S = |eta|^2 holds provided the rational approx is good public: + //Used in the heatbath, refresh the shift coefficients of the L (LorR=0) or R (LorR=1) operator + virtual void heatbathRefreshShiftCoefficients(int LorR, RealD to){ + AbstractEOFAFermion&op = LorR == 0 ? Lop : Rop; + op.RefreshShiftCoefficients(to); + } + + + //Use the same solver for L,R in all cases ExactOneFlavourRatioPseudoFermionAction(AbstractEOFAFermion& _Lop, AbstractEOFAFermion& _Rop, OperatorFunction& CG, Params& p, bool use_fc=false) - : ExactOneFlavourRatioPseudoFermionAction(_Lop,_Rop,CG,CG,CG,CG,CG,p,use_fc) {}; - + : ExactOneFlavourRatioPseudoFermionAction(_Lop,_Rop,CG,CG,CG,CG,CG,CG,p,use_fc) {}; + + //Use the same solver for L,R in the heatbath but different solvers elsewhere ExactOneFlavourRatioPseudoFermionAction(AbstractEOFAFermion& _Lop, AbstractEOFAFermion& _Rop, - OperatorFunction& HeatbathCG, + OperatorFunction& HeatbathCG, + OperatorFunction& ActionCGL, OperatorFunction& ActionCGR, + OperatorFunction& DerivCGL , OperatorFunction& DerivCGR, + Params& p, + bool use_fc=false) + : ExactOneFlavourRatioPseudoFermionAction(_Lop,_Rop,HeatbathCG,HeatbathCG, ActionCGL, ActionCGR, DerivCGL,DerivCGR,p,use_fc) {}; + + //Use different solvers for L,R in all cases + ExactOneFlavourRatioPseudoFermionAction(AbstractEOFAFermion& _Lop, + AbstractEOFAFermion& _Rop, + OperatorFunction& HeatbathCGL, OperatorFunction& HeatbathCGR, OperatorFunction& ActionCGL, OperatorFunction& ActionCGR, OperatorFunction& DerivCGL , OperatorFunction& DerivCGR, Params& p, bool use_fc=false) : Lop(_Lop), Rop(_Rop), - SolverHB(HeatbathCG,false,true), + SolverHBL(HeatbathCGL,false,true), SolverHBR(HeatbathCGR,false,true), SolverL(ActionCGL, false, true), SolverR(ActionCGR, false, true), DerivativeSolverL(DerivCGL, false, true), DerivativeSolverR(DerivCGR, false, true), Phi(_Lop.FermionGrid()), param(p), - use_heatbath_forecasting(use_fc) + use_heatbath_forecasting(use_fc), + initial_action(false) { AlgRemez remez(param.lo, param.hi, param.precision); @@ -97,6 +124,8 @@ NAMESPACE_BEGIN(Grid); PowerNegHalf.Init(remez, param.tolerance, true); }; + const FermionField &getPhi() const{ return Phi; } + virtual std::string action_name() { return "ExactOneFlavourRatioPseudoFermionAction"; } virtual std::string LogParameters() { @@ -117,6 +146,19 @@ NAMESPACE_BEGIN(Grid); else{ for(int s=0; s sig^2 = 0.5. + // + RealD scale = std::sqrt(0.5); + + FermionField eta (Lop.FermionGrid()); + gaussian(pRNG,eta); eta = eta * scale; + + refresh(U,eta); + } + // EOFA heatbath: see Eqn. (29) of arXiv:1706.05843 // We generate a Gaussian noise vector \eta, and then compute // \Phi = M_{\rm EOFA}^{-1/2} * \eta @@ -124,12 +166,10 @@ NAMESPACE_BEGIN(Grid); // // As a check of rational require \Phi^dag M_{EOFA} \Phi == eta^dag M^-1/2^dag M M^-1/2 eta = eta^dag eta // - virtual void refresh(const GaugeField& U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) - { + void refresh(const GaugeField &U, const FermionField &eta) { Lop.ImportGauge(U); Rop.ImportGauge(U); - FermionField eta (Lop.FermionGrid()); FermionField CG_src (Lop.FermionGrid()); FermionField CG_soln (Lop.FermionGrid()); FermionField Forecast_src(Lop.FermionGrid()); @@ -140,11 +180,6 @@ NAMESPACE_BEGIN(Grid); if(use_heatbath_forecasting){ prev_solns.reserve(param.degree); } ChronoForecast, FermionField> Forecast; - // Seed with Gaussian noise vector (var = 0.5) - RealD scale = std::sqrt(0.5); - gaussian(pRNG,eta); - eta = eta * scale; - // \Phi = ( \alpha_{0} + \sum_{k=1}^{N_{p}} \alpha_{l} * \gamma_{l} ) * \eta RealD N(PowerNegHalf.norm); for(int k=0; k tmp(2, Lop.FermionGrid()); - mPhi = phi; + out = in; // LH term: S = S - k <\Phi| P_{-} \Omega_{-}^{\dagger} H(mf)^{-1} \Omega_{-} P_{-} |\Phi> - spProj(Phi, spProj_Phi, -1, Lop.Ls); - Lop.Omega(spProj_Phi, tmp[0], -1, 0); + spProj(in, spProj_in, -1, Lop.Ls); + Lop.Omega(spProj_in, tmp[0], -1, 0); G5R5(tmp[1], tmp[0]); tmp[0] = Zero(); SolverL(Lop, tmp[1], tmp[0]); Lop.Dtilde(tmp[0], tmp[1]); // We actually solved Cayley preconditioned system: transform back Lop.Omega(tmp[1], tmp[0], -1, 1); - mPhi = mPhi - Lop.k * innerProduct(spProj_Phi, tmp[0]).real(); + spProj(tmp[0], tmp[1], -1, Lop.Ls); + + out = out - Lop.k * tmp[1]; // RH term: S = S + k <\Phi| P_{+} \Omega_{+}^{\dagger} ( H(mb) - // - \Delta_{+}(mf,mb) P_{+} )^{-1} \Omega_{-} P_{-} |\Phi> - spProj(Phi, spProj_Phi, 1, Rop.Ls); - Rop.Omega(spProj_Phi, tmp[0], 1, 0); + // - \Delta_{+}(mf,mb) P_{+} )^{-1} \Omega_{+} P_{+} |\Phi> + spProj(in, spProj_in, 1, Rop.Ls); + Rop.Omega(spProj_in, tmp[0], 1, 0); G5R5(tmp[1], tmp[0]); tmp[0] = Zero(); SolverR(Rop, tmp[1], tmp[0]); Rop.Dtilde(tmp[0], tmp[1]); Rop.Omega(tmp[1], tmp[0], 1, 1); - action += Rop.k * innerProduct(spProj_Phi, tmp[0]).real(); -#endif + spProj(tmp[0], tmp[1], 1, Rop.Ls); + + out = out + Rop.k * tmp[1]; } + //Due to the structure of EOFA, it is no more expensive to compute the inverse of Meofa + //To ensure correctness we can simply reuse the heatbath code but use the rational approx + //f(x) = 1/x which corresponds to alpha_0=0, alpha_1=1, beta_1=0 => gamma_1=1 + void MeofaInv(const GaugeField &U, const FermionField &in, FermionField &out) { + Lop.ImportGauge(U); + Rop.ImportGauge(U); + + FermionField CG_src (Lop.FermionGrid()); + FermionField CG_soln (Lop.FermionGrid()); + std::vector tmp(2, Lop.FermionGrid()); + + // \Phi = ( \alpha_{0} + \sum_{k=1}^{N_{p}} \alpha_{l} * \gamma_{l} ) * \eta + // = 1 * \eta + out = in; + + // LH terms: + // \Phi = \Phi + k \sum_{k=1}^{N_{p}} P_{-} \Omega_{-}^{\dagger} ( H(mf) + // - \gamma_{l} \Delta_{-}(mf,mb) P_{-} )^{-1} \Omega_{-} P_{-} \eta + spProj(in, tmp[0], -1, Lop.Ls); + Lop.Omega(tmp[0], tmp[1], -1, 0); + G5R5(CG_src, tmp[1]); + { + heatbathRefreshShiftCoefficients(0, -1.); //-gamma_1 = -1. + + CG_soln = Zero(); // Just use zero as the initial guess + SolverHBL(Lop, CG_src, CG_soln); + + Lop.Dtilde(CG_soln, tmp[0]); // We actually solved Cayley preconditioned system: transform back + tmp[1] = Lop.k * tmp[0]; + } + Lop.Omega(tmp[1], tmp[0], -1, 1); + spProj(tmp[0], tmp[1], -1, Lop.Ls); + out = out + tmp[1]; + + // RH terms: + // \Phi = \Phi - k \sum_{k=1}^{N_{p}} P_{+} \Omega_{+}^{\dagger} ( H(mb) + // - \beta_l\gamma_{l} \Delta_{+}(mf,mb) P_{+} )^{-1} \Omega_{+} P_{+} \eta + spProj(in, tmp[0], 1, Rop.Ls); + Rop.Omega(tmp[0], tmp[1], 1, 0); + G5R5(CG_src, tmp[1]); + { + heatbathRefreshShiftCoefficients(1, 0.); //-gamma_1 * beta_1 = 0 + + CG_soln = Zero(); + SolverHBR(Rop, CG_src, CG_soln); + + Rop.Dtilde(CG_soln, tmp[0]); // We actually solved Cayley preconditioned system: transform back + tmp[1] = - Rop.k * tmp[0]; + } + Rop.Omega(tmp[1], tmp[0], 1, 1); + spProj(tmp[0], tmp[1], 1, Rop.Ls); + out = out + tmp[1]; + + // Reset shift coefficients for energy and force evals + heatbathRefreshShiftCoefficients(0, 0.0); + heatbathRefreshShiftCoefficients(1, -1.0); + }; + + + + // EOFA action: see Eqn. (10) of arXiv:1706.05843 virtual RealD S(const GaugeField& U) { @@ -271,7 +374,7 @@ NAMESPACE_BEGIN(Grid); action -= Lop.k * innerProduct(spProj_Phi, tmp[0]).real(); // RH term: S = S + k <\Phi| P_{+} \Omega_{+}^{\dagger} ( H(mb) - // - \Delta_{+}(mf,mb) P_{+} )^{-1} \Omega_{-} P_{-} |\Phi> + // - \Delta_{+}(mf,mb) P_{+} )^{-1} \Omega_{+} P_{+} |\Phi> spProj(Phi, spProj_Phi, 1, Rop.Ls); Rop.Omega(spProj_Phi, tmp[0], 1, 0); G5R5(tmp[1], tmp[0]); @@ -281,6 +384,26 @@ NAMESPACE_BEGIN(Grid); Rop.Omega(tmp[1], tmp[0], 1, 1); action += Rop.k * innerProduct(spProj_Phi, tmp[0]).real(); + if(initial_action){ + //For the first call to S after refresh, S = |eta|^2. We can use this to ensure the rational approx is good + RealD diff = action - norm2_eta; + + //S_init = eta^dag M^{-1/2} M M^{-1/2} eta + //S_init - eta^dag eta = eta^dag ( M^{-1/2} M M^{-1/2} - 1 ) eta + + //If approximate solution + //S_init - eta^dag eta = eta^dag ( [M^{-1/2}+\delta M^{-1/2}] M [M^{-1/2}+\delta M^{-1/2}] - 1 ) eta + // \approx eta^dag ( \delta M^{-1/2} M^{1/2} + M^{1/2}\delta M^{-1/2} ) eta + // We divide out |eta|^2 to remove source scaling but the tolerance on this check should still be somewhat higher than the actual approx tolerance + RealD test = fabs(diff)/norm2_eta; //test the quality of the rational approx + + std::cout << GridLogMessage << action_name() << " initial action " << action << " expect " << norm2_eta << "; diff " << diff << std::endl; + std::cout << GridLogMessage << action_name() << "[ eta^dag ( M^{-1/2} M M^{-1/2} - 1 ) eta ]/|eta^2| = " << test << " expect 0 (tol " << param.BoundsCheckTol << ")" << std::endl; + + assert( ( test < param.BoundsCheckTol ) && " Initial action check failed" ); + initial_action = false; + } + return action; }; @@ -329,6 +452,40 @@ NAMESPACE_BEGIN(Grid); }; }; + template + class ExactOneFlavourRatioMixedPrecHeatbathPseudoFermionAction : public ExactOneFlavourRatioPseudoFermionAction{ + public: + INHERIT_IMPL_TYPES(ImplD); + typedef OneFlavourRationalParams Params; + + private: + AbstractEOFAFermion& LopF; // the basic LH operator + AbstractEOFAFermion& RopF; // the basic RH operator + + public: + + virtual std::string action_name() { return "ExactOneFlavourRatioMixedPrecHeatbathPseudoFermionAction"; } + + //Used in the heatbath, refresh the shift coefficients of the L (LorR=0) or R (LorR=1) operator + virtual void heatbathRefreshShiftCoefficients(int LorR, RealD to){ + AbstractEOFAFermion &op = LorR == 0 ? LopF : RopF; + op.RefreshShiftCoefficients(to); + this->ExactOneFlavourRatioPseudoFermionAction::heatbathRefreshShiftCoefficients(LorR,to); + } + + ExactOneFlavourRatioMixedPrecHeatbathPseudoFermionAction(AbstractEOFAFermion& _LopF, + AbstractEOFAFermion& _RopF, + AbstractEOFAFermion& _LopD, + AbstractEOFAFermion& _RopD, + OperatorFunction& HeatbathCGL, OperatorFunction& HeatbathCGR, + OperatorFunction& ActionCGL, OperatorFunction& ActionCGR, + OperatorFunction& DerivCGL , OperatorFunction& DerivCGR, + Params& p, + bool use_fc=false) : + LopF(_LopF), RopF(_RopF), ExactOneFlavourRatioPseudoFermionAction(_LopD, _RopD, HeatbathCGL, HeatbathCGR, ActionCGL, ActionCGR, DerivCGL, DerivCGR, p, use_fc){} + }; + + NAMESPACE_END(Grid); #endif diff --git a/Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatio.h b/Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatio.h new file mode 100644 index 00000000..2b08cf49 --- /dev/null +++ b/Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatio.h @@ -0,0 +1,372 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/qcd/action/pseudofermion/GeneralEvenOddRationalRatio.h + + Copyright (C) 2015 + + Author: Christopher Kelly + Author: Peter Boyle + + 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 */ +#ifndef QCD_PSEUDOFERMION_GENERAL_EVEN_ODD_RATIONAL_RATIO_H +#define QCD_PSEUDOFERMION_GENERAL_EVEN_ODD_RATIONAL_RATIO_H + +NAMESPACE_BEGIN(Grid); + + ///////////////////////////////////////////////////////// + // Generic rational approximation for ratios of operators + ///////////////////////////////////////////////////////// + + /* S_f = -log( det( [M^dag M]/[V^dag V] )^{1/inv_pow} ) + = chi^dag ( [M^dag M]/[V^dag V] )^{-1/inv_pow} chi\ + = chi^dag ( [V^dag V]^{-1/2} [M^dag M] [V^dag V]^{-1/2} )^{-1/inv_pow} chi\ + = chi^dag [V^dag V]^{1/(2*inv_pow)} [M^dag M]^{-1/inv_pow} [V^dag V]^{1/(2*inv_pow)} chi\ + + S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi + + BIG WARNING: + Here V^dag V is referred to in this code as the "numerator" operator and M^dag M is the *denominator* operator. + this refers to their position in the pseudofermion action, which is the *inverse* of what appears in the determinant + Thus for DWF the numerator operator is the Pauli-Villars operator + + Here P/Q \sim R_{1/(2*inv_pow)} ~ (V^dagV)^{1/(2*inv_pow)} + Here N/D \sim R_{-1/inv_pow} ~ (M^dagM)^{-1/inv_pow} + */ + + template + class GeneralEvenOddRatioRationalPseudoFermionAction : public Action { + public: + + INHERIT_IMPL_TYPES(Impl); + + typedef RationalActionParams Params; + Params param; + + //For action evaluation + MultiShiftFunction ApproxPowerAction ; //rational approx for X^{1/inv_pow} + MultiShiftFunction ApproxNegPowerAction; //rational approx for X^{-1/inv_pow} + MultiShiftFunction ApproxHalfPowerAction; //rational approx for X^{1/(2*inv_pow)} + MultiShiftFunction ApproxNegHalfPowerAction; //rational approx for X^{-1/(2*inv_pow)} + + //For the MD integration + MultiShiftFunction ApproxPowerMD ; //rational approx for X^{1/inv_pow} + MultiShiftFunction ApproxNegPowerMD; //rational approx for X^{-1/inv_pow} + MultiShiftFunction ApproxHalfPowerMD; //rational approx for X^{1/(2*inv_pow)} + MultiShiftFunction ApproxNegHalfPowerMD; //rational approx for X^{-1/(2*inv_pow)} + + private: + + FermionOperator & NumOp;// the basic operator + FermionOperator & DenOp;// the basic operator + FermionField PhiEven; // the pseudo fermion field for this trajectory + FermionField PhiOdd; // the pseudo fermion field for this trajectory + + //Generate the approximation to x^{1/inv_pow} (->approx) and x^{-1/inv_pow} (-> approx_inv) by an approx_degree degree rational approximation + //CG_tolerance is used to issue a warning if the approximation error is larger than the tolerance of the CG and is otherwise just stored in the MultiShiftFunction for use by the multi-shift + static void generateApprox(MultiShiftFunction &approx, MultiShiftFunction &approx_inv, int inv_pow, int approx_degree, double CG_tolerance, AlgRemez &remez){ + std::cout< CG_tolerance) + std::cout< schurOp(numerator ? NumOp : DenOp); + ConjugateGradientMultiShift msCG(MaxIter, approx); + msCG(schurOp,in, out); + } + virtual void multiShiftInverse(bool numerator, const MultiShiftFunction &approx, const Integer MaxIter, const FermionField &in, std::vector &out_elems, FermionField &out){ + SchurDifferentiableOperator schurOp(numerator ? NumOp : DenOp); + ConjugateGradientMultiShift msCG(MaxIter, approx); + msCG(schurOp,in, out_elems, out); + } + //Allow derived classes to override the gauge import + virtual void ImportGauge(const GaugeField &U){ + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); + } + + public: + + GeneralEvenOddRatioRationalPseudoFermionAction(FermionOperator &_NumOp, + FermionOperator &_DenOp, + const Params & p + ) : + NumOp(_NumOp), + DenOp(_DenOp), + PhiOdd (_NumOp.FermionRedBlackGrid()), + PhiEven(_NumOp.FermionRedBlackGrid()), + param(p) + { + std::cout<Broadcast(0,r); + + if ( param.BoundsCheckFreq != 0 && (r % param.BoundsCheckFreq)==0 ) { + std::cout< MdagM(DenOp); + std::cout< MpvPhi_k (n_pv,NumOp.FermionRedBlackGrid()); + std::vector MpvMfMpvPhi_k(n_pv,NumOp.FermionRedBlackGrid()); + std::vector MfMpvPhi_k (n_f ,NumOp.FermionRedBlackGrid()); + + FermionField MpvPhi(NumOp.FermionRedBlackGrid()); + FermionField MfMpvPhi(NumOp.FermionRedBlackGrid()); + FermionField MpvMfMpvPhi(NumOp.FermionRedBlackGrid()); + FermionField Y(NumOp.FermionRedBlackGrid()); + + GaugeField tmp(NumOp.GaugeGrid()); + + ImportGauge(U); + + std::cout< MdagM(DenOp); + SchurDifferentiableOperator VdagV(NumOp); + + + RealD ak; + + dSdU = Zero(); + + // With these building blocks + // + // dS/dU = + // \sum_k -ak MfMpvPhi_k^dag [ dM^dag M + M^dag dM ] MfMpvPhi_k (1) + // + \sum_k -ak MpvMfMpvPhi_k^\dag [ dV^dag V + V^dag dV ] MpvPhi_k (2) + // -ak MpvPhi_k^dag [ dV^dag V + V^dag dV ] MpvMfMpvPhi_k (3) + + //(1) + std::cout< + Author: Peter Boyle + + 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 */ +#ifndef QCD_PSEUDOFERMION_GENERAL_EVEN_ODD_RATIONAL_RATIO_MIXED_PREC_H +#define QCD_PSEUDOFERMION_GENERAL_EVEN_ODD_RATIONAL_RATIO_MIXED_PREC_H + +NAMESPACE_BEGIN(Grid); + + ///////////////////////////////////////////////////////////////////////////////////////////////////////////// + // Generic rational approximation for ratios of operators utilizing the mixed precision multishift algorithm + // cf. GeneralEvenOddRational.h for details + ///////////////////////////////////////////////////////////////////////////////////////////////////////////// + + template + class GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction : public GeneralEvenOddRatioRationalPseudoFermionAction { + private: + typedef typename ImplD::FermionField FermionFieldD; + typedef typename ImplF::FermionField FermionFieldF; + + FermionOperator & NumOpD; + FermionOperator & DenOpD; + + FermionOperator & NumOpF; + FermionOperator & DenOpF; + + Integer ReliableUpdateFreq; + protected: + + //Allow derived classes to override the multishift CG + virtual void multiShiftInverse(bool numerator, const MultiShiftFunction &approx, const Integer MaxIter, const FermionFieldD &in, FermionFieldD &out){ + SchurDifferentiableOperator schurOpD(numerator ? NumOpD : DenOpD); + SchurDifferentiableOperator schurOpF(numerator ? NumOpF : DenOpF); + + ConjugateGradientMultiShiftMixedPrec msCG(MaxIter, approx, NumOpF.FermionRedBlackGrid(), schurOpF, ReliableUpdateFreq); + msCG(schurOpD, in, out); + } + virtual void multiShiftInverse(bool numerator, const MultiShiftFunction &approx, const Integer MaxIter, const FermionFieldD &in, std::vector &out_elems, FermionFieldD &out){ + SchurDifferentiableOperator schurOpD(numerator ? NumOpD : DenOpD); + SchurDifferentiableOperator schurOpF(numerator ? NumOpF : DenOpF); + + ConjugateGradientMultiShiftMixedPrec msCG(MaxIter, approx, NumOpF.FermionRedBlackGrid(), schurOpF, ReliableUpdateFreq); + msCG(schurOpD, in, out_elems, out); + } + //Allow derived classes to override the gauge import + virtual void ImportGauge(const typename ImplD::GaugeField &Ud){ + typename ImplF::GaugeField Uf(NumOpF.GaugeGrid()); + precisionChange(Uf, Ud); + + NumOpD.ImportGauge(Ud); + DenOpD.ImportGauge(Ud); + + NumOpF.ImportGauge(Uf); + DenOpF.ImportGauge(Uf); + } + + public: + GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction(FermionOperator &_NumOpD, FermionOperator &_DenOpD, + FermionOperator &_NumOpF, FermionOperator &_DenOpF, + const RationalActionParams & p, Integer _ReliableUpdateFreq + ) : GeneralEvenOddRatioRationalPseudoFermionAction(_NumOpD, _DenOpD, p), + ReliableUpdateFreq(_ReliableUpdateFreq), NumOpD(_NumOpD), DenOpD(_DenOpD), NumOpF(_NumOpF), DenOpF(_DenOpF){} + + virtual std::string action_name(){return "GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction";} + }; + +NAMESPACE_END(Grid); + +#endif diff --git a/Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h b/Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h index 6752ea19..1b36ae0f 100644 --- a/Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h +++ b/Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h @@ -40,257 +40,31 @@ NAMESPACE_BEGIN(Grid); // Here N/D \sim R_{-1/2} ~ (M^dagM)^{-1/2} template - class OneFlavourEvenOddRatioRationalPseudoFermionAction : public Action { + class OneFlavourEvenOddRatioRationalPseudoFermionAction : public GeneralEvenOddRatioRationalPseudoFermionAction { public: - - INHERIT_IMPL_TYPES(Impl); - typedef OneFlavourRationalParams Params; - Params param; - - MultiShiftFunction PowerHalf ; - MultiShiftFunction PowerNegHalf; - MultiShiftFunction PowerQuarter; - MultiShiftFunction PowerNegQuarter; - private: - - FermionOperator & NumOp;// the basic operator - FermionOperator & DenOp;// the basic operator - FermionField PhiEven; // the pseudo fermion field for this trajectory - FermionField PhiOdd; // the pseudo fermion field for this trajectory - FermionField Noise; // spare noise field for bounds check + 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: - OneFlavourEvenOddRatioRationalPseudoFermionAction(FermionOperator &_NumOp, - FermionOperator &_DenOp, - Params & p - ) : - NumOp(_NumOp), - DenOp(_DenOp), - PhiOdd (_NumOp.FermionRedBlackGrid()), - PhiEven(_NumOp.FermionRedBlackGrid()), - Noise(_NumOp.FermionRedBlackGrid()), - param(p) - { - AlgRemez remez(param.lo,param.hi,param.precision); + FermionOperator &_DenOp, + const Params & p + ) : + GeneralEvenOddRatioRationalPseudoFermionAction(_NumOp, _DenOp, transcribe(p)){} - // MdagM^(+- 1/2) - std::cout< sig^2 = 0.5. - // - // So eta should be of width sig = 1/sqrt(2). - - RealD scale = std::sqrt(0.5); - - FermionField eta(NumOp.FermionGrid()); - FermionField etaOdd (NumOp.FermionRedBlackGrid()); - FermionField etaEven(NumOp.FermionRedBlackGrid()); - FermionField tmp(NumOp.FermionRedBlackGrid()); - - gaussian(pRNG,eta); eta=eta*scale; - - pickCheckerboard(Even,etaEven,eta); - pickCheckerboard(Odd,etaOdd,eta); - - Noise = etaOdd; - NumOp.ImportGauge(U); - DenOp.ImportGauge(U); - - - // MdagM^1/4 eta - SchurDifferentiableOperator MdagM(DenOp); - ConjugateGradientMultiShift msCG_M(param.MaxIter,PowerQuarter); - msCG_M(MdagM,etaOdd,tmp); - - // VdagV^-1/4 MdagM^1/4 eta - SchurDifferentiableOperator VdagV(NumOp); - ConjugateGradientMultiShift msCG_V(param.MaxIter,PowerNegQuarter); - msCG_V(VdagV,tmp,PhiOdd); - - assert(NumOp.ConstEE() == 1); - assert(DenOp.ConstEE() == 1); - PhiEven = Zero(); - - }; - - ////////////////////////////////////////////////////// - // S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi - ////////////////////////////////////////////////////// - virtual RealD S(const GaugeField &U) { - - NumOp.ImportGauge(U); - DenOp.ImportGauge(U); - - FermionField X(NumOp.FermionRedBlackGrid()); - FermionField Y(NumOp.FermionRedBlackGrid()); - - // VdagV^1/4 Phi - SchurDifferentiableOperator VdagV(NumOp); - ConjugateGradientMultiShift msCG_V(param.MaxIter,PowerQuarter); - msCG_V(VdagV,PhiOdd,X); - - // MdagM^-1/4 VdagV^1/4 Phi - SchurDifferentiableOperator MdagM(DenOp); - ConjugateGradientMultiShift msCG_M(param.MaxIter,PowerNegQuarter); - msCG_M(MdagM,X,Y); - - // Randomly apply rational bounds checks. - auto grid = NumOp.FermionGrid(); - auto r=rand(); - grid->Broadcast(0,r); - if ( (r%param.BoundsCheckFreq)==0 ) { - FermionField gauss(NumOp.FermionRedBlackGrid()); - gauss = Noise; - HighBoundCheck(MdagM,gauss,param.hi); - InverseSqrtBoundsCheck(param.MaxIter,param.tolerance*100,MdagM,gauss,PowerNegHalf); - ChebyBoundsCheck(MdagM,Noise,param.lo,param.hi); - } - - // Phidag VdagV^1/4 MdagM^-1/4 MdagM^-1/4 VdagV^1/4 Phi - RealD action = norm2(Y); - - return action; - }; - - // S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi - // - // Here, M is some 5D operator and V is the Pauli-Villars field - // N and D makeup the rat. poly of the M term and P and & makeup the rat.poly of the denom term - // - // Need - // dS_f/dU = chi^dag d[P/Q] N/D P/Q chi - // + chi^dag P/Q d[N/D] P/Q chi - // + chi^dag P/Q N/D d[P/Q] chi - // - // P/Q is expressed as partial fraction expansion: - // - // a0 + \sum_k ak/(V^dagV + bk) - // - // d[P/Q] is then - // - // \sum_k -ak [V^dagV+bk]^{-1} [ dV^dag V + V^dag dV ] [V^dag V + bk]^{-1} - // - // and similar for N/D. - // - // Need - // MpvPhi_k = [Vdag V + bk]^{-1} chi - // MpvPhi = {a0 + \sum_k ak [Vdag V + bk]^{-1} }chi - // - // MfMpvPhi_k = [MdagM+bk]^{-1} MpvPhi - // MfMpvPhi = {a0 + \sum_k ak [Mdag M + bk]^{-1} } MpvPhi - // - // MpvMfMpvPhi_k = [Vdag V + bk]^{-1} MfMpvchi - // - - virtual void deriv(const GaugeField &U,GaugeField & dSdU) { - - const int n_f = PowerNegHalf.poles.size(); - const int n_pv = PowerQuarter.poles.size(); - - std::vector MpvPhi_k (n_pv,NumOp.FermionRedBlackGrid()); - std::vector MpvMfMpvPhi_k(n_pv,NumOp.FermionRedBlackGrid()); - std::vector MfMpvPhi_k (n_f ,NumOp.FermionRedBlackGrid()); - - FermionField MpvPhi(NumOp.FermionRedBlackGrid()); - FermionField MfMpvPhi(NumOp.FermionRedBlackGrid()); - FermionField MpvMfMpvPhi(NumOp.FermionRedBlackGrid()); - FermionField Y(NumOp.FermionRedBlackGrid()); - - GaugeField tmp(NumOp.GaugeGrid()); - - NumOp.ImportGauge(U); - DenOp.ImportGauge(U); - - SchurDifferentiableOperator VdagV(NumOp); - SchurDifferentiableOperator MdagM(DenOp); - - ConjugateGradientMultiShift msCG_V(param.MaxIter,PowerQuarter); - ConjugateGradientMultiShift msCG_M(param.MaxIter,PowerNegHalf); - - msCG_V(VdagV,PhiOdd,MpvPhi_k,MpvPhi); - msCG_M(MdagM,MpvPhi,MfMpvPhi_k,MfMpvPhi); - msCG_V(VdagV,MfMpvPhi,MpvMfMpvPhi_k,MpvMfMpvPhi); - - RealD ak; - - dSdU = Zero(); - - // With these building blocks - // - // dS/dU = - // \sum_k -ak MfMpvPhi_k^dag [ dM^dag M + M^dag dM ] MfMpvPhi_k (1) - // + \sum_k -ak MpvMfMpvPhi_k^\dag [ dV^dag V + V^dag dV ] MpvPhi_k (2) - // -ak MpvPhi_k^dag [ dV^dag V + V^dag dV ] MpvMfMpvPhi_k (3) - - //(1) - for(int k=0;k #include #include +#include +#include #include #include diff --git a/Grid/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h b/Grid/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h index 6d06e090..56dd6840 100644 --- a/Grid/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h +++ b/Grid/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h @@ -88,15 +88,9 @@ NAMESPACE_BEGIN(Grid); } - virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) { + const FermionField &getPhiOdd() const{ return PhiOdd; } - // P(phi) = e^{- phi^dag Vpc (MpcdagMpc)^-1 Vpcdag phi} - // - // NumOp == V - // DenOp == M - // - // Take phi_o = Vpcdag^{-1} Mpcdag eta_o ; eta_o = Mpcdag^{-1} Vpcdag Phi - // + virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) { // P(eta_o) = e^{- eta_o^dag eta_o} // // e^{x^2/2 sig^2} => sig^2 = 0.5. @@ -104,12 +98,22 @@ NAMESPACE_BEGIN(Grid); RealD scale = std::sqrt(0.5); FermionField eta (NumOp.FermionGrid()); + gaussian(pRNG,eta); eta = eta * scale; + + refresh(U,eta); + } + + void refresh(const GaugeField &U, const FermionField &eta) { + + // P(phi) = e^{- phi^dag Vpc (MpcdagMpc)^-1 Vpcdag phi} + // + // NumOp == V + // DenOp == M + // FermionField etaOdd (NumOp.FermionRedBlackGrid()); FermionField etaEven(NumOp.FermionRedBlackGrid()); FermionField tmp (NumOp.FermionRedBlackGrid()); - gaussian(pRNG,eta); - pickCheckerboard(Even,etaEven,eta); pickCheckerboard(Odd,etaOdd,eta); @@ -128,10 +132,6 @@ NAMESPACE_BEGIN(Grid); // Even det factors DenOp.MooeeDag(etaEven,tmp); NumOp.MooeeInvDag(tmp,PhiEven); - - PhiOdd =PhiOdd*scale; - PhiEven=PhiEven*scale; - }; ////////////////////////////////////////////////////// diff --git a/Grid/qcd/hmc/GenericHMCrunner.h b/Grid/qcd/hmc/GenericHMCrunner.h index 098f8f22..727b3e24 100644 --- a/Grid/qcd/hmc/GenericHMCrunner.h +++ b/Grid/qcd/hmc/GenericHMCrunner.h @@ -151,12 +151,22 @@ public: Resources.GetCheckPointer()->CheckpointRestore(Parameters.StartTrajectory, U, Resources.GetSerialRNG(), Resources.GetParallelRNG()); + } else if (Parameters.StartingType == "CheckpointStartReseed") { + // Same as CheckpointRestart but reseed the RNGs using the fixed integer seeding used for ColdStart and HotStart + // Useful for creating new evolution streams from an existing stream + + // WARNING: Unfortunately because the checkpointer doesn't presently allow us to separately restore the RNG and gauge fields we have to load + // an existing RNG checkpoint first; make sure one is available and named correctly + Resources.GetCheckPointer()->CheckpointRestore(Parameters.StartTrajectory, U, + Resources.GetSerialRNG(), + Resources.GetParallelRNG()); + Resources.SeedFixedIntegers(); } else { // others std::cout << GridLogError << "Unrecognized StartingType\n"; std::cout << GridLogError - << "Valid [HotStart, ColdStart, TepidStart, CheckpointStart]\n"; + << "Valid [HotStart, ColdStart, TepidStart, CheckpointStart, CheckpointStartReseed]\n"; exit(1); } } diff --git a/Grid/qcd/hmc/HMCModules.h b/Grid/qcd/hmc/HMCModules.h index 4c61a006..cf0edd26 100644 --- a/Grid/qcd/hmc/HMCModules.h +++ b/Grid/qcd/hmc/HMCModules.h @@ -80,7 +80,9 @@ public: std::cout << GridLogError << "Seeds not initialized" << std::endl; exit(1); } + std::cout << GridLogMessage << "Reseeding serial RNG with seed vector " << SerialSeeds << std::endl; sRNG_.SeedFixedIntegers(SerialSeeds); + std::cout << GridLogMessage << "Reseeding parallel RNG with seed vector " << ParallelSeeds << std::endl; pRNG_->SeedFixedIntegers(ParallelSeeds); } }; diff --git a/Grid/qcd/hmc/integrators/Integrator.h b/Grid/qcd/hmc/integrators/Integrator.h index 070cbea1..9563698c 100644 --- a/Grid/qcd/hmc/integrators/Integrator.h +++ b/Grid/qcd/hmc/integrators/Integrator.h @@ -334,15 +334,19 @@ public: void refresh(Field& U, GridSerialRNG & sRNG, GridParallelRNG& pRNG) { assert(P.Grid() == U.Grid()); - std::cout << GridLogIntegrator << "Integrator refresh\n"; + std::cout << GridLogIntegrator << "Integrator refresh" << std::endl; + std::cout << GridLogIntegrator << "Generating momentum" << std::endl; FieldImplementation::generate_momenta(P, sRNG, pRNG); // Update the smeared fields, can be implemented as observer // necessary to keep the fields updated even after a reject // of the Metropolis + std::cout << GridLogIntegrator << "Updating smeared fields" << std::endl; Smearer.set_Field(U); // Set the (eventual) representations gauge fields + + std::cout << GridLogIntegrator << "Updating representations" << std::endl; Representations.update(U); // The Smearer is attached to a pointer of the gauge field diff --git a/tests/forces/Test_gpdwf_force_1f_2f.cc b/tests/forces/Test_gpdwf_force_1f_2f.cc new file mode 100644 index 00000000..7e14eb08 --- /dev/null +++ b/tests/forces/Test_gpdwf_force_1f_2f.cc @@ -0,0 +1,446 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./forces/Test_gpdwf_force_1f_2f.cc + + Copyright (C) 2015 + +Author: Christopher Kelly +Author: paboyle + + 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 + +using namespace std; +using namespace Grid; + +//Here we test the G-parity action and force between the 1f (doubled-lattice) and 2f approaches + + +void copyConjGauge(LatticeGaugeFieldD &Umu_1f, const LatticeGaugeFieldD &Umu_2f, const int nu){ + GridBase* UGrid_2f = Umu_2f.Grid(); + GridBase* UGrid_1f = Umu_1f.Grid(); + + Replicate(Umu_2f,Umu_1f); + + int L_2f = UGrid_2f->FullDimensions()[nu]; + int L_1f = UGrid_1f->FullDimensions()[nu]; + assert(L_1f == 2 * L_2f); + + //Coordinate grid for reference + LatticeInteger xcoor_1f(UGrid_1f); + LatticeCoordinate(xcoor_1f,nu); + + //Copy-conjugate the gauge field + //First C-shift the lattice by Lx/2 + { + LatticeGaugeField Umu_shift = conjugate( Cshift(Umu_1f,nu,L_2f) ); + Umu_1f = where( xcoor_1f >= Integer(L_2f), Umu_shift, Umu_1f ); + + //We use the in built APBC + //Make the gauge field antiperiodic in nu-direction + //decltype(PeekIndex(Umu_1f,nu)) Unu(UGrid_1f); + //Unu = PeekIndex(Umu_1f,nu); + //Unu = where(xcoor_1f == Integer(2*L_2f-1), -Unu, Unu); + //PokeIndex(Umu_1f,Unu,nu); + } +} + +template +void convertFermion1f_from_2f(FermionField1f &out_1f, const FermionField2f &in_2f, const int nu, bool is_4d){ + GridBase* FGrid_1f = out_1f.Grid(); + GridBase* FGrid_2f = in_2f.Grid(); + + int nuoff = is_4d ? 0 : 1; //s in 0 direction + + int L_2f = FGrid_2f->FullDimensions()[nu+nuoff]; + int L_1f = FGrid_1f->FullDimensions()[nu+nuoff]; + assert(L_1f == 2 * L_2f); + + auto in_f0_2fgrid = PeekIndex(in_2f,0); //flavor 0 on 2f Grid + FermionField1f in_f0_1fgrid(FGrid_1f); + Replicate(in_f0_2fgrid, in_f0_1fgrid); //has flavor 0 on both halves + + auto in_f1_2fgrid = PeekIndex(in_2f,1); //flavor 1 on 2f Grid + FermionField1f in_f1_1fgrid(FGrid_1f); + Replicate(in_f1_2fgrid, in_f1_1fgrid); //has flavor 1 on both halves + + LatticeInteger xcoor_1f(FGrid_1f); + LatticeCoordinate(xcoor_1f,nu+nuoff); + + out_1f = where(xcoor_1f < L_2f, in_f0_1fgrid, in_f1_1fgrid); +} + +template +class RatioActionSetupBase{ +protected: + TwoFlavourEvenOddRatioPseudoFermionAction *pf_1f; + TwoFlavourEvenOddRatioPseudoFermionAction *pf_2f; + + GparityAction* action_2f; + GparityAction* action_PV_2f; + StandardAction* action_1f; + StandardAction* action_PV_1f; + + ConjugateGradient CG_1f; + ConjugateGradient CG_2f; + + RatioActionSetupBase(): CG_1f(1.0e-8,10000), CG_2f(1.0e-8,10000){} + + void setupPseudofermion(){ + pf_1f = new TwoFlavourEvenOddRatioPseudoFermionAction(*action_PV_1f, *action_1f, CG_1f, CG_1f); + pf_2f = new TwoFlavourEvenOddRatioPseudoFermionAction(*action_PV_2f, *action_2f, CG_2f, CG_2f); + } + +public: + GparityAction & action2f(){ return *action_2f; } + StandardAction & action1f(){ return *action_1f; } + + void refreshAction(LatticeGaugeField &Umu_2f, typename GparityAction::FermionField &eta_2f, + LatticeGaugeField &Umu_1f, typename StandardAction::FermionField &eta_1f){ + pf_1f->refresh(Umu_1f, eta_1f); + pf_2f->refresh(Umu_2f, eta_2f); + + //Compare PhiOdd + RealD norm_1f = norm2(pf_1f->getPhiOdd()); + RealD norm_2f = norm2(pf_2f->getPhiOdd()); + + std::cout << "Test PhiOdd 2f: " << norm_2f << " 1f: " << norm_1f << std::endl; + } + + void computeAction(RealD &S_2f, RealD &S_1f, LatticeGaugeField &Umu_2f, LatticeGaugeField &Umu_1f){ + S_1f = pf_1f->S(Umu_1f); + S_2f = pf_2f->S(Umu_2f); + } + + void computeDeriv(LatticeGaugeField &deriv_2f, LatticeGaugeField &deriv_1f, LatticeGaugeField &Umu_2f, LatticeGaugeField &Umu_1f){ + pf_1f->deriv(Umu_1f, deriv_1f); + pf_2f->deriv(Umu_2f, deriv_2f); + } + +}; + + + + +template +struct setupAction{}; + +template<> +struct setupAction: public RatioActionSetupBase{ + typedef GparityWilsonTMFermionD GparityAction; + typedef WilsonTMFermionD StandardAction; + + setupAction(GridCartesian* UGrid_2f, GridRedBlackCartesian* UrbGrid_2f, GridCartesian* FGrid_2f, GridRedBlackCartesian* FrbGrid_2f, + GridCartesian* UGrid_1f, GridRedBlackCartesian* UrbGrid_1f, GridCartesian* FGrid_1f, GridRedBlackCartesian* FrbGrid_1f, + LatticeGaugeField &Umu_2f, LatticeGaugeField &Umu_1f, int nu): RatioActionSetupBase(){ + RealD mass=-1.8; + //Use same DSDR twists as https://arxiv.org/pdf/1208.4412.pdf + RealD epsilon_f = 0.02; //numerator (in determinant) + RealD epsilon_b = 0.5; + + std::vector twists(Nd,0); + twists[nu] = 1; //GPBC in y + twists[3] = 1; //APBC + GparityAction::ImplParams params_2f; params_2f.twists = twists; + action_2f = new GparityWilsonTMFermionD(Umu_2f,*UGrid_2f,*UrbGrid_2f, mass, epsilon_f, params_2f); + action_PV_2f = new GparityWilsonTMFermionD(Umu_2f,*UGrid_2f,*UrbGrid_2f, mass, epsilon_b, params_2f); + + DomainWallFermionD::ImplParams params_1f; + params_1f.boundary_phases[nu] = -1; + params_1f.boundary_phases[3] = -1; + + action_1f = new WilsonTMFermionD(Umu_1f,*UGrid_1f,*UrbGrid_1f, mass, epsilon_f, params_1f); + action_PV_1f = new WilsonTMFermionD(Umu_1f,*UGrid_1f,*UrbGrid_1f, mass, epsilon_b, params_1f); + + setupPseudofermion(); + } + + static bool is4d(){ return true; } +}; + + +template<> +struct setupAction: public RatioActionSetupBase{ + typedef GparityDomainWallFermionD GparityAction; + typedef DomainWallFermionD StandardAction; + + setupAction(GridCartesian* UGrid_2f, GridRedBlackCartesian* UrbGrid_2f, GridCartesian* FGrid_2f, GridRedBlackCartesian* FrbGrid_2f, + GridCartesian* UGrid_1f, GridRedBlackCartesian* UrbGrid_1f, GridCartesian* FGrid_1f, GridRedBlackCartesian* FrbGrid_1f, + LatticeGaugeField &Umu_2f, LatticeGaugeField &Umu_1f, int nu): RatioActionSetupBase(){ + RealD mass=0.01; + RealD M5=1.8; + + std::vector twists(Nd,0); + twists[nu] = 1; //GPBC in y + twists[3] = 1; //APBC + GparityDomainWallFermionD::ImplParams params_2f; params_2f.twists = twists; + action_2f = new GparityDomainWallFermionD(Umu_2f,*FGrid_2f,*FrbGrid_2f,*UGrid_2f,*UrbGrid_2f,mass,M5,params_2f); + action_PV_2f = new GparityDomainWallFermionD(Umu_2f,*FGrid_2f,*FrbGrid_2f,*UGrid_2f,*UrbGrid_2f,1.0,M5,params_2f); + + DomainWallFermionD::ImplParams params_1f; + params_1f.boundary_phases[nu] = -1; + params_1f.boundary_phases[3] = -1; + + action_1f = new DomainWallFermionD(Umu_1f,*FGrid_1f,*FrbGrid_1f,*UGrid_1f,*UrbGrid_1f,mass,M5,params_1f); + action_PV_1f = new DomainWallFermionD(Umu_1f,*FGrid_1f,*FrbGrid_1f,*UGrid_1f,*UrbGrid_1f,1.0,M5,params_1f); + + setupPseudofermion(); + } + + static bool is4d(){ return false; } +}; + + + + + +//For EOFA we need a different pseudofermion type +template<> +struct setupAction{ + typedef GparityDomainWallEOFAFermionD GparityAction; + typedef DomainWallEOFAFermionD StandardAction; + + ExactOneFlavourRatioPseudoFermionAction *pf_1f; + ExactOneFlavourRatioPseudoFermionAction *pf_2f; + + GparityAction* action_2f; + GparityAction* action_PV_2f; + StandardAction* action_1f; + StandardAction* action_PV_1f; + + ConjugateGradient CG_1f; + ConjugateGradient CG_2f; + +public: + GparityAction & action2f(){ return *action_2f; } + StandardAction & action1f(){ return *action_1f; } + + void refreshAction(LatticeGaugeField &Umu_2f, typename GparityAction::FermionField &eta_2f, + LatticeGaugeField &Umu_1f, typename StandardAction::FermionField &eta_1f){ + pf_1f->refresh(Umu_1f, eta_1f); + pf_2f->refresh(Umu_2f, eta_2f); + + //Compare PhiOdd + RealD norm_1f = norm2(pf_1f->getPhi()); + RealD norm_2f = norm2(pf_2f->getPhi()); + + std::cout << "Test Phi 2f: " << norm_2f << " 1f: " << norm_1f << std::endl; + } + + void computeAction(RealD &S_2f, RealD &S_1f, LatticeGaugeField &Umu_2f, LatticeGaugeField &Umu_1f){ + S_1f = pf_1f->S(Umu_1f); + S_2f = pf_2f->S(Umu_2f); + } + + void computeDeriv(LatticeGaugeField &deriv_2f, LatticeGaugeField &deriv_1f, LatticeGaugeField &Umu_2f, LatticeGaugeField &Umu_1f){ + pf_1f->deriv(Umu_1f, deriv_1f); + pf_2f->deriv(Umu_2f, deriv_2f); + } + + + setupAction(GridCartesian* UGrid_2f, GridRedBlackCartesian* UrbGrid_2f, GridCartesian* FGrid_2f, GridRedBlackCartesian* FrbGrid_2f, + GridCartesian* UGrid_1f, GridRedBlackCartesian* UrbGrid_1f, GridCartesian* FGrid_1f, GridRedBlackCartesian* FrbGrid_1f, + LatticeGaugeField &Umu_2f, LatticeGaugeField &Umu_1f, int nu): CG_1f(1.0e-8,10000), CG_2f(1.0e-8,10000){ + RealD mass=0.01; + RealD M5=1.8; + + std::vector twists(Nd,0); + twists[nu] = 1; //GPBC in y + twists[3] = 1; //APBC + GparityAction::ImplParams params_2f; params_2f.twists = twists; + action_2f = new GparityAction(Umu_2f,*FGrid_2f,*FrbGrid_2f,*UGrid_2f,*UrbGrid_2f, mass, mass, 1.0, 0.0, -1, M5, params_2f); + action_PV_2f = new GparityAction(Umu_2f,*FGrid_2f,*FrbGrid_2f,*UGrid_2f,*UrbGrid_2f, 1.0, mass, 1.0, -1.0, 1, M5, params_2f); //cf Test_dwf_gpforce_eofa.cc + + StandardAction::ImplParams params_1f; + params_1f.boundary_phases[nu] = -1; + params_1f.boundary_phases[3] = -1; + + action_1f = new StandardAction(Umu_1f,*FGrid_1f,*FrbGrid_1f,*UGrid_1f,*UrbGrid_1f, mass, mass, 1.0, 0.0, -1, M5, params_1f); + action_PV_1f = new StandardAction(Umu_1f,*FGrid_1f,*FrbGrid_1f,*UGrid_1f,*UrbGrid_1f, 1.0, mass, 1.0, -1.0, 1, M5, params_1f); + + OneFlavourRationalParams RationalParams(0.95, 100.0, 5000, 1.0e-12, 12); + + pf_1f = new ExactOneFlavourRatioPseudoFermionAction(*action_1f, *action_PV_1f, CG_1f, CG_1f, CG_1f, CG_1f, CG_1f, RationalParams, true); + pf_2f = new ExactOneFlavourRatioPseudoFermionAction(*action_2f, *action_PV_2f, CG_2f, CG_2f, CG_2f, CG_2f, CG_2f, RationalParams, true); + } + + static bool is4d(){ return false; } +}; + + +template +void runTest(int argc, char** argv){ + Grid_init(&argc,&argv); + + const int nu = 1; + Coordinate latt_2f = GridDefaultLatt(); + Coordinate latt_1f = latt_2f; + latt_1f[nu] *= 2; + + Coordinate simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + Coordinate mpi_layout = GridDefaultMpi(); + + const int Ls=8; + + GridCartesian * UGrid_1f = SpaceTimeGrid::makeFourDimGrid(latt_1f, simd_layout, mpi_layout); + GridRedBlackCartesian * UrbGrid_1f = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid_1f); + GridCartesian * FGrid_1f = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid_1f); + GridRedBlackCartesian * FrbGrid_1f = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid_1f); + + + GridCartesian * UGrid_2f = SpaceTimeGrid::makeFourDimGrid(latt_2f, simd_layout, mpi_layout); + GridRedBlackCartesian * UrbGrid_2f = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid_2f); + GridCartesian * FGrid_2f = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid_2f); + GridRedBlackCartesian * FrbGrid_2f = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid_2f); + + std::vector seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + GridParallelRNG RNG5_2f(FGrid_2f); RNG5_2f.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4_2f(UGrid_2f); RNG4_2f.SeedFixedIntegers(seeds4); + + LatticeGaugeField Umu_2f(UGrid_2f); + SU::HotConfiguration(RNG4_2f,Umu_2f); + + LatticeGaugeField Umu_1f(UGrid_1f); + copyConjGauge(Umu_1f, Umu_2f, nu); + + typedef typename GparityAction::FermionField GparityFermionField; + typedef typename StandardAction::FermionField StandardFermionField; + + setupAction setup(UGrid_2f, UrbGrid_2f, FGrid_2f, FrbGrid_2f, + UGrid_1f, UrbGrid_1f, FGrid_1f, FrbGrid_1f, + Umu_2f, Umu_1f, nu); + GridBase* FGrid_2f_a = setup.action2f().FermionGrid(); + GridBase* FGrid_1f_a = setup.action1f().FermionGrid(); + GridBase* FrbGrid_2f_a = setup.action2f().FermionRedBlackGrid(); + GridBase* FrbGrid_1f_a = setup.action1f().FermionRedBlackGrid(); + bool is_4d = setup.is4d(); + + //Check components by doing an inversion + { + setup.action2f().ImportGauge(Umu_2f); + setup.action1f().ImportGauge(Umu_1f); + + GparityFermionField src_2f(FGrid_2f_a); + gaussian(is_4d ? RNG4_2f : RNG5_2f, src_2f); + + StandardFermionField src_1f(FGrid_1f_a); + convertFermion1f_from_2f(src_1f, src_2f, nu, is_4d); + + StandardFermionField src_o_1f(FrbGrid_1f_a); + StandardFermionField result_o_1f(FrbGrid_1f_a); + pickCheckerboard(Odd,src_o_1f,src_1f); + result_o_1f=Zero(); + + SchurDiagMooeeOperator HermOpEO_1f(setup.action1f()); + ConjugateGradient CG_1f(1.0e-8,10000); + CG_1f(HermOpEO_1f,src_o_1f,result_o_1f); + + + GparityFermionField src_o_2f(FrbGrid_2f_a); + GparityFermionField result_o_2f(FrbGrid_2f_a); + pickCheckerboard(Odd,src_o_2f,src_2f); + result_o_2f=Zero(); + + SchurDiagMooeeOperator HermOpEO_2f(setup.action2f()); + ConjugateGradient CG_2f(1.0e-8,10000); + CG_2f(HermOpEO_2f,src_o_2f,result_o_2f); + + RealD norm_1f = norm2(result_o_1f); + RealD norm_2f = norm2(result_o_2f); + + std::cout << "Test fermion inversion 2f: " << norm_2f << " 1f: " << norm_1f << std::endl; + } + + //Generate eta + RealD scale = std::sqrt(0.5); + + GparityFermionField eta_2f(FGrid_2f_a); + gaussian(is_4d ? RNG4_2f : RNG5_2f,eta_2f); eta_2f = eta_2f * scale; + + StandardFermionField eta_1f(FGrid_1f_a); + convertFermion1f_from_2f(eta_1f, eta_2f, nu, is_4d); + + setup.refreshAction(Umu_2f, eta_2f, Umu_1f, eta_1f); + + //Initial action is just |eta^2| + RealD S_1f, S_2f; + + setup.computeAction(S_2f, S_1f, Umu_2f, Umu_1f); + + std::cout << "Test Initial action 2f: " << S_2f << " 1f: " << S_1f << " diff: " << S_2f - S_1f << std::endl; + + //Do a random gauge field refresh + SU::HotConfiguration(RNG4_2f,Umu_2f); + copyConjGauge(Umu_1f, Umu_2f, nu); + + //Compute the action again + setup.computeAction(S_2f, S_1f, Umu_2f, Umu_1f); + + std::cout << "Test Action after gauge field randomize 2f: " << S_2f << " 1f: " << S_1f << " diff: " << S_2f - S_1f << std::endl; + + //Compute the derivative and test the conjugate relation + LatticeGaugeField deriv_2f(UGrid_2f); + LatticeGaugeField deriv_1f(UGrid_1f); + setup.computeDeriv(deriv_2f, deriv_1f, Umu_2f, Umu_1f); + + //Have to combine the two forces on the 1f by symmetrizing under the complex conjugate + { + RealD norm2_pre = norm2(deriv_1f); + LatticeGaugeField deriv_1f_shift = conjugate( Cshift(deriv_1f, nu, latt_2f[nu]) ); + deriv_1f = deriv_1f + deriv_1f_shift; + std::cout << "Test combine/symmetrize forces on 1f lattice, dS/dU : " << norm2_pre << " -> " << norm2(deriv_1f) << std::endl; + } + + LatticeGaugeField deriv_1f_from_2f(UGrid_1f); + copyConjGauge(deriv_1f_from_2f, deriv_2f, nu); + std::cout << "Test copy-conj 2f dS/dU to obtain equivalent 1f force : " << norm2(deriv_2f) << " -> " << norm2(deriv_1f_from_2f) << std::endl; + + LatticeGaugeField diff_deriv_1f = deriv_1f - deriv_1f_from_2f; + + std::cout << "Test dS/dU 1f constructed from 2f derivative: " << norm2(deriv_1f_from_2f) << " dS/dU 1f actual: " << norm2(deriv_1f) << " Norm of difference: " << norm2(diff_deriv_1f) << std::endl; + + std::cout<< GridLogMessage << "Done" <(argc, argv); + }else if(action == "EOFA"){ + runTest(argc, argv); + }else if(action == "DSDR"){ + runTest(argc,argv); + }else{ + assert(0); + } +} diff --git a/tests/forces/Test_mobius_force_eofa.cc b/tests/forces/Test_mobius_force_eofa.cc index eea3e3f4..1d25771a 100644 --- a/tests/forces/Test_mobius_force_eofa.cc +++ b/tests/forces/Test_mobius_force_eofa.cc @@ -89,7 +89,49 @@ int main (int argc, char** argv) ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, CG, CG, CG, CG, Params, false); GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds4); + + //Check the rational approximation + { + RealD scale = std::sqrt(0.5); + LatticeFermion eta (Lop.FermionGrid()); + gaussian(RNG5,eta); eta = eta * scale; + + Meofa.refresh(U, eta); + + //Phi = M^{-1/2} eta + //M is Hermitian + //(Phi, M Phi) = eta^\dagger M^{-1/2} M M^{-1/2} eta = eta^\dagger eta + LatticeFermion phi = Meofa.getPhi(); + LatticeFermion Mphi(FGrid); + + Meofa.Meofa(U, phi, Mphi); + std::cout << "Computing inner product" << std::endl; + ComplexD inner = innerProduct(phi, Mphi); + ComplexD test = inner - norm2(eta); + + std::cout << "(phi, Mphi) - (eta,eta): " << test << " expect 0" << std::endl; + + assert(test.real() < 1e-8); + assert(test.imag() < 1e-8); + + //Another test is to use heatbath twice to apply M^{-1/2} to Phi then apply M + // M Phi' + //= M M^{-1/2} Phi + //= M M^{-1/2} M^{-1/2} eta + //= eta + Meofa.refresh(U, phi); + LatticeFermion phi2 = Meofa.getPhi(); + LatticeFermion test2(FGrid); + Meofa.Meofa(U, phi2, test2); + test2 = test2 - eta; + RealD test2_norm = norm2(test2); + std::cout << "|M M^{-1/2} M^{-1/2} eta - eta|^2 = " << test2_norm << " expect 0" << std::endl; + assert( test2_norm < 1e-8 ); + } + + Meofa.refresh(U, sRNG, RNG5 ); + RealD S = Meofa.S(U); // pdag M p // get the deriv of phidag M phi with respect to "U" diff --git a/tests/forces/Test_mobius_gparity_eofa_mixed.cc b/tests/forces/Test_mobius_gparity_eofa_mixed.cc new file mode 100644 index 00000000..d490e838 --- /dev/null +++ b/tests/forces/Test_mobius_gparity_eofa_mixed.cc @@ -0,0 +1,233 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/forces/Test_mobius_gparity_eofa_mixed.cc + +Copyright (C) 2017 + +Author: Christopher Kelly +Author: Peter Boyle +Author: David Murphy + +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 + +using namespace std; +using namespace Grid; + ; + +typedef GparityWilsonImplD FermionImplPolicyD; +typedef GparityMobiusEOFAFermionD FermionActionD; +typedef typename FermionActionD::FermionField FermionFieldD; + +typedef GparityWilsonImplF FermionImplPolicyF; +typedef GparityMobiusEOFAFermionF FermionActionF; +typedef typename FermionActionF::FermionField FermionFieldF; + +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.) + { + }; + + void operator()(LinearOperatorBase &LinOpU, const FieldD &src, FieldD &psi) { + + std::cout << GridLogMessage << " Mixed precision CG wrapper operator() "<(&LinOpU); + assert(&(SchurOpU->_Mat)==&(LinOpD._Mat)); + + precisionChange(FermOpF.Umu, FermOpD.Umu); + + pickCheckerboard(Even,FermOpF.UmuEven,FermOpF.Umu); + pickCheckerboard(Odd ,FermOpF.UmuOdd ,FermOpF.Umu); + + //////////////////////////////////////////////////////////////////////////////////// + // Make a mixed precision conjugate gradient + //////////////////////////////////////////////////////////////////////////////////// + MixedPrecisionConjugateGradient MPCG(Tolerance,MaxInnerIterations,MaxOuterIterations,SinglePrecGrid5,LinOpF,LinOpD); + MPCG.InnerTolerance = InnerTolerance; + std::cout << GridLogMessage << "Calling mixed precision Conjugate Gradient" < seeds4({1,2,3,5}); + std::vector seeds5({5,6,7,8}); + GridParallelRNG RNG5(FGridD); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGridD); RNG4.SeedFixedIntegers(seeds4); + + int threads = GridThread::GetThreads(); + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + LatticeGaugeFieldD Ud(UGridD); + SU::HotConfiguration(RNG4,Ud); + + LatticeGaugeFieldF Uf(UGridF); + precisionChange(Uf, Ud); + + RealD b = 2.5; + RealD c = 1.5; + RealD mf = 0.01; + RealD mb = 1.0; + RealD M5 = 1.8; + FermionActionD::ImplParams params; + params.twists[0] = 1; //GPBC in X + params.twists[Nd-1] = 1; //APRD in T + + std::vector gtwists(4,0); + gtwists[0] = 1; + + ConjugateGimplD::setDirections(gtwists); + + FermionActionD LopD(Ud, *FGridD, *FrbGridD, *UGridD, *UrbGridD, mf, mf, mb, 0.0, -1, M5, b, c, params); + FermionActionD RopD(Ud, *FGridD, *FrbGridD, *UGridD, *UrbGridD, mb, mf, mb, -1.0, 1, M5, b, c, params); + + FermionActionF LopF(Uf, *FGridF, *FrbGridF, *UGridF, *UrbGridF, mf, mf, mb, 0.0, -1, M5, b, c, params); + FermionActionF RopF(Uf, *FGridF, *FrbGridF, *UGridF, *UrbGridF, mb, mf, mb, -1.0, 1, M5, b, c, params); + + + OneFlavourRationalParams OFRp(0.95, 100.0, 5000, 1.0e-12, 12); + ConjugateGradient CG(1.0e-10, 10000); + + + typedef SchurDiagMooeeOperator EOFAschuropD; + typedef SchurDiagMooeeOperator EOFAschuropF; + + EOFAschuropD linopL_D(LopD); + EOFAschuropD linopR_D(RopD); + + EOFAschuropF linopL_F(LopF); + EOFAschuropF linopR_F(RopF); + + typedef MixedPrecisionConjugateGradientOperatorFunction EOFA_mxCG; + + EOFA_mxCG MCG_L(1e-10, 10000, 1000, UGridF, FrbGridF, LopF, LopD, linopL_F, linopL_D); + MCG_L.InnerTolerance = 1e-5; + + EOFA_mxCG MCG_R(1e-10, 10000, 1000, UGridF, FrbGridF, RopF, RopD, linopR_F, linopR_D); + MCG_R.InnerTolerance = 1e-5; + + ExactOneFlavourRatioPseudoFermionAction MeofaD(LopD, RopD, CG, CG, CG, CG, CG, OFRp, true); + ExactOneFlavourRatioMixedPrecHeatbathPseudoFermionAction MeofaMx(LopF, RopF, LopD, RopD, MCG_L, MCG_R, MCG_L, MCG_R, MCG_L, MCG_R, OFRp, true); + + FermionFieldD eta(FGridD); + gaussian(RNG5, eta); + + MeofaD.refresh(Ud, eta); + MeofaMx.refresh(Ud, eta); + + FermionFieldD diff_phi(FGridD); + diff_phi = MeofaD.getPhi() - MeofaMx.getPhi(); + + RealD n = norm2(diff_phi); + + std::cout << GridLogMessage << "Phi(double)=" << norm2(MeofaD.getPhi()) << " Phi(mixed)=" << norm2(MeofaMx.getPhi()) << " diff=" << n << std::endl; + + assert(n < 1e-8); + + RealD Sd = MeofaD.S(Ud); + RealD Smx = MeofaMx.S(Ud); + + std::cout << GridLogMessage << "Initial action double=" << Sd << " mixed=" << Smx << " diff=" << Sd-Smx << std::endl; + + assert(fabs(Sd-Smx) < 1e-6); + + SU::HotConfiguration(RNG4,Ud); + precisionChange(Uf, Ud); + + Sd = MeofaD.S(Ud); + Smx = MeofaMx.S(Ud); + + std::cout << GridLogMessage << "After randomizing U, action double=" << Sd << " mixed=" << Smx << " diff=" << Sd-Smx << std::endl; + + assert(fabs(Sd-Smx) < 1e-6); + + std::cout << GridLogMessage << "Done" << std::endl; + Grid_finalize(); +} diff --git a/tests/hmc/Test_action_dwf_gparity2fvs1f.cc b/tests/hmc/Test_action_dwf_gparity2fvs1f.cc new file mode 100644 index 00000000..830bcead --- /dev/null +++ b/tests/hmc/Test_action_dwf_gparity2fvs1f.cc @@ -0,0 +1,257 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: tests/hmc/Test_action_dwf_gparity2fvs1f.cc + + Copyright (C) 2015 + + Author: Christopher Kelly + Author: paboyle + + 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 + +using namespace Grid; + + + +template +void copy2fTo1fFermionField(FermionField1f &out, const FermionField2f &in, int gpdir){ + auto f0_halfgrid = PeekIndex(in,0); //on 2f Grid + FermionField1f f0_fullgrid_dbl(out.Grid()); + Replicate(f0_halfgrid, f0_fullgrid_dbl); //double it up to live on the 1f Grid + + auto f1_halfgrid = PeekIndex(in,1); + FermionField1f f1_fullgrid_dbl(out.Grid()); + Replicate(f1_halfgrid, f1_fullgrid_dbl); + + const Coordinate &dim_2f = in.Grid()->GlobalDimensions(); + const Coordinate &dim_1f = out.Grid()->GlobalDimensions(); + + //We have to be careful for 5d fields; the s-direction is placed before the x,y,z,t and so we need to shift gpdir by 1 + std::cout << "gpdir " << gpdir << std::endl; + + gpdir+=1; + std::cout << "gpdir for 5D fields " << gpdir << std::endl; + + std::cout << "dim_2f " << dim_2f << std::endl; + std::cout << "dim_1f " << dim_1f << std::endl; + + assert(dim_1f[gpdir] == 2*dim_2f[gpdir]); + + LatticeInteger xcoor_1f(out.Grid()); //5d lattice integer + LatticeCoordinate(xcoor_1f,gpdir); + + int L = dim_2f[gpdir]; + + out = where(xcoor_1f < L, f0_fullgrid_dbl, f1_fullgrid_dbl); +} + +//Both have the same field type +void copy2fTo1fGaugeField(LatticeGaugeField &out, const LatticeGaugeField &in, int gpdir){ + LatticeGaugeField U_dbl(out.Grid()); + Replicate(in, U_dbl); + + LatticeGaugeField Uconj_dbl = conjugate( U_dbl ); + + const Coordinate &dim_2f = in.Grid()->GlobalDimensions(); + + LatticeInteger xcoor_1f(out.Grid()); + LatticeCoordinate(xcoor_1f,gpdir); + + int L = dim_2f[gpdir]; + + out = where(xcoor_1f < L, U_dbl, Uconj_dbl); +} + + +std::ostream & operator<<(std::ostream &os, const Coordinate &x){ + os << "("; + for(int i=0;i seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + GridParallelRNG RNG5_2f(FGrid_2f); RNG5_2f.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4_2f(UGrid_2f); RNG4_2f.SeedFixedIntegers(seeds4); + + std::cout << "Generating hot 2f gauge configuration" << std::endl; + LatticeGaugeField Umu_2f(UGrid_2f); + SU::HotConfiguration(RNG4_2f,Umu_2f); + + std::cout << "Copying 2f->1f gauge field" << std::endl; + LatticeGaugeField Umu_1f(UGrid_1f); + copy2fTo1fGaugeField(Umu_1f, Umu_2f, mu); + + typedef GparityWilsonImplR FermionImplPolicy2f; + typedef GparityDomainWallFermionR FermionAction2f; + typedef typename FermionAction2f::FermionField FermionField2f; + + typedef WilsonImplR FermionImplPolicy1f; + typedef DomainWallFermionR FermionAction1f; + typedef typename FermionAction1f::FermionField FermionField1f; + + std::cout << "Generating eta 2f" << std::endl; + FermionField2f eta_2f(FGrid_2f); + gaussian(RNG5_2f, eta_2f); + + RealD scale = std::sqrt(0.5); + eta_2f=eta_2f*scale; + + std::cout << "Copying 2f->1f eta" << std::endl; + FermionField1f eta_1f(FGrid_1f); + copy2fTo1fFermionField(eta_1f, eta_2f, mu); + + Real beta = 2.13; + Real light_mass = 0.01; + Real strange_mass = 0.032; + Real pv_mass = 1.0; + RealD M5 = 1.8; + + //Setup the Dirac operators + std::cout << "Initializing Dirac operators" << std::endl; + + FermionAction2f::ImplParams Params_2f; + Params_2f.twists[mu] = 1; + Params_2f.twists[Nd-1] = 1; //APBC in time direction + + //note 'Num' and 'Den' here refer to the determinant ratio, not the operator ratio in the pseudofermion action where the two are inverted + //to my mind the Pauli Villars and 'denominator' are synonymous but the Grid convention has this as the 'Numerator' operator in the RHMC implementation + FermionAction2f NumOp_2f(Umu_2f,*FGrid_2f,*FrbGrid_2f,*UGrid_2f, *UrbGrid_2f, light_mass,M5,Params_2f); + FermionAction2f DenOp_2f(Umu_2f,*FGrid_2f,*FrbGrid_2f,*UGrid_2f, *UrbGrid_2f, pv_mass, M5,Params_2f); + + FermionAction1f::ImplParams Params_1f; + Params_1f.boundary_phases[mu] = -1; //antiperiodic in doubled lattice in GP direction + Params_1f.boundary_phases[Nd-1] = -1; + + FermionAction1f NumOp_1f(Umu_1f,*FGrid_1f,*FrbGrid_1f,*UGrid_1f, *UrbGrid_1f, light_mass,M5,Params_1f); + FermionAction1f DenOp_1f(Umu_1f,*FGrid_1f,*FrbGrid_1f,*UGrid_1f, *UrbGrid_1f, pv_mass, M5,Params_1f); + + //Test the replication routines by running a CG on eta + double StoppingCondition = 1e-10; + double MaxCGIterations = 30000; + ConjugateGradient CG_2f(StoppingCondition,MaxCGIterations); + ConjugateGradient CG_1f(StoppingCondition,MaxCGIterations); + + NumOp_1f.ImportGauge(Umu_1f); + NumOp_2f.ImportGauge(Umu_2f); + + FermionField1f test_1f(FGrid_1f); + FermionField2f test_2f(FGrid_2f); + + MdagMLinearOperator Linop_1f(NumOp_1f); + MdagMLinearOperator Linop_2f(NumOp_2f); + + CG_1f(Linop_1f, eta_1f, test_1f); + CG_2f(Linop_2f, eta_2f, test_2f); + RealD test_1f_norm = norm2(test_1f); + RealD test_2f_norm = norm2(test_2f); + + std::cout << "Verification of replication routines: " << test_1f_norm << " " << test_2f_norm << " " << test_1f_norm - test_2f_norm << std::endl; + + +#if 1 + typedef GeneralEvenOddRatioRationalPseudoFermionAction Action2f; + typedef GeneralEvenOddRatioRationalPseudoFermionAction Action1f; + + RationalActionParams rational_params; + rational_params.inv_pow = 2; + rational_params.lo = 1e-5; + rational_params.hi = 32; + rational_params.md_degree = 16; + rational_params.action_degree = 16; + + Action2f action_2f(DenOp_2f, NumOp_2f, rational_params); + Action1f action_1f(DenOp_1f, NumOp_1f, rational_params); +#else + typedef TwoFlavourEvenOddRatioPseudoFermionAction Action2f; + typedef TwoFlavourEvenOddRatioPseudoFermionAction Action1f; + + Action2f action_2f(DenOp_2f, NumOp_2f, CG_2f, CG_2f); + Action1f action_1f(DenOp_1f, NumOp_1f, CG_1f, CG_1f); +#endif + + + std::cout << "Action refresh" << std::endl; + action_2f.refresh(Umu_2f, eta_2f); + action_1f.refresh(Umu_1f, eta_1f); + + std::cout << "Action compute post heatbath" << std::endl; + RealD S_2f = action_2f.S(Umu_2f); + RealD S_1f = action_1f.S(Umu_1f); + + std::cout << "Action comparison post heatbath" << std::endl; + std::cout << S_2f << " " << S_1f << " " << S_2f-S_1f << std::endl; + + //Change the gauge field between refresh and action eval else the matrix and inverse matrices all cancel and we just get |eta|^2 + SU::HotConfiguration(RNG4_2f,Umu_2f); + copy2fTo1fGaugeField(Umu_1f, Umu_2f, mu); + + //Now compute the action with the new gauge field + std::cout << "Action compute post gauge field update" << std::endl; + S_2f = action_2f.S(Umu_2f); + S_1f = action_1f.S(Umu_1f); + + std::cout << "Action comparison post gauge field update" << std::endl; + std::cout << S_2f << " " << S_1f << " " << S_2f-S_1f << std::endl; + + Grid_finalize(); +} // main + + diff --git a/tests/hmc/Test_hmc_GparityIwasakiGauge.cc b/tests/hmc/Test_hmc_GparityIwasakiGauge.cc index 7f74d5d8..d4bfa0a5 100644 --- a/tests/hmc/Test_hmc_GparityIwasakiGauge.cc +++ b/tests/hmc/Test_hmc_GparityIwasakiGauge.cc @@ -58,7 +58,7 @@ int main(int argc, char **argv) { CheckpointerParameters CPparams; CPparams.config_prefix = "ckpoint_EODWF_lat"; CPparams.rng_prefix = "ckpoint_EODWF_rng"; - CPparams.saveInterval = 5; + CPparams.saveInterval = 1; CPparams.format = "IEEE64BIG"; TheHMC.Resources.LoadNerscCheckpointer(CPparams); @@ -79,7 +79,7 @@ int main(int argc, char **argv) { // that have a complex construction // standard RealD beta = 2.6 ; - const int nu = 3; + const int nu = 1; std::vector twists(Nd,0); twists[nu] = 1; ConjugateGimplD::setDirections(twists); diff --git a/tests/solver/Test_eofa_inv.cc b/tests/solver/Test_eofa_inv.cc new file mode 100644 index 00000000..564405c2 --- /dev/null +++ b/tests/solver/Test_eofa_inv.cc @@ -0,0 +1,125 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/solver/Test_eofa_inv.cc + +Copyright (C) 2017 + +Author: Christopher Kelly +Author: Peter Boyle +Author: David Murphy + +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 + +using namespace std; +using namespace Grid; + ; + +int main (int argc, char** argv) +{ + Grid_init(&argc, &argv); + + Coordinate latt_size = GridDefaultLatt(); + Coordinate simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + Coordinate mpi_layout = GridDefaultMpi(); + + const int Ls = 8; + + GridCartesian *UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()), GridDefaultMpi()); + GridRedBlackCartesian *UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridCartesian *FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls, UGrid); + GridRedBlackCartesian *FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, UGrid); + + // Want a different conf at every run + // First create an instance of an engine. + std::random_device rnd_device; + // Specify the engine and distribution. + std::mt19937 mersenne_engine(rnd_device()); + std::uniform_int_distribution dist(1, 100); + + auto gen = std::bind(dist, mersenne_engine); + std::vector seeds4(4); + generate(begin(seeds4), end(seeds4), gen); + + //std::vector seeds4({1,2,3,5}); + std::vector seeds5({5,6,7,8}); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + + int threads = GridThread::GetThreads(); + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + LatticeFermion phi (FGrid); gaussian(RNG5, phi); + LatticeFermion Mphi (FGrid); + LatticeFermion MphiPrime (FGrid); + + LatticeGaugeField U(UGrid); + SU::HotConfiguration(RNG4,U); + + //////////////////////////////////// + // Unmodified matrix element + //////////////////////////////////// + RealD b = 2.5; + RealD c = 1.5; + RealD mf = 0.01; + RealD mb = 1.0; + RealD M5 = 1.8; + MobiusEOFAFermionR Lop(U, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, mf, mb, 0.0, -1, M5, b, c); + MobiusEOFAFermionR Rop(U, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mb, mf, mb, -1.0, 1, M5, b, c); + OneFlavourRationalParams Params(0.95, 100.0, 5000, 1.0e-10, 12); + ConjugateGradient CG(1.0e-10, 5000); + ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, CG, CG, CG, CG, Params, false); + + GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds4); + + + //Random field + LatticeFermion eta(FGrid); + gaussian(RNG5,eta); + + //Check left inverse + LatticeFermion Meta(FGrid); + Meofa.Meofa(U, eta, Meta); + + LatticeFermion MinvMeta(FGrid); + Meofa.MeofaInv(U, Meta, MinvMeta); + + LatticeFermion diff = MinvMeta - eta; + + std::cout << GridLogMessage << "eta: " << norm2(eta) << " M*eta: " << norm2(Meta) << " M^{-1}*M*eta: " << norm2(MinvMeta) << " M^{-1}*M*eta - eta: " << norm2(diff) << " (expect 0)" << std::endl; + assert(norm2(diff) < 1e-8); + + //Check right inverse + LatticeFermion MinvEta(FGrid); + Meofa.MeofaInv(U, eta, MinvEta); + + LatticeFermion MMinvEta(FGrid); + Meofa.Meofa(U, MinvEta, MMinvEta); + + diff = MMinvEta - eta; + + std::cout << GridLogMessage << "eta: " << norm2(eta) << " M^{-1}*eta: " << norm2(MinvEta) << " M*M^{-1}*eta: " << norm2(MMinvEta) << " M*M^{-1}*eta - eta: " << norm2(diff) << " (expect 0)" << std::endl; + assert(norm2(diff) < 1e-8); + + std::cout << GridLogMessage << "Done" << std::endl; + Grid_finalize(); +}