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(); +}