diff --git a/Grid/qcd/action/ActionParams.h b/Grid/qcd/action/ActionParams.h index 0e6a11c6..6374256e 100644 --- a/Grid/qcd/action/ActionParams.h +++ b/Grid/qcd/action/ActionParams.h @@ -36,7 +36,8 @@ NAMESPACE_BEGIN(Grid); // These can move into a params header and be given MacroMagic serialisation struct GparityWilsonImplParams { - Coordinate twists; + Coordinate twists; //Here the first Nd-1 directions are treated as "spatial", and a twist value of 1 indicates G-parity BCs in that direction. + //mu=Nd-1 is assumed to be the time direction and a twist value of 1 indicates antiperiodic BCs GparityWilsonImplParams() : twists(Nd, 0) {}; }; @@ -85,6 +86,44 @@ struct StaggeredImplParams { precision(_precision), BoundsCheckFreq(_BoundsCheckFreq){}; }; + + + /*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, + RealD, hi, + int, MaxIter, + RealD, tolerance, + int, degree, + int, precision, + int, BoundsCheckFreq); + // constructor + RationalActionParams(int _inv_pow = 2, + RealD _lo = 0.0, + RealD _hi = 1.0, + int _maxit = 1000, + RealD tol = 1.0e-8, + int _degree = 10, + int _precision = 64, + int _BoundsCheckFreq=20) + : inv_pow(_inv_pow), + lo(_lo), + hi(_hi), + MaxIter(_maxit), + tolerance(tol), + degree(_degree), + precision(_precision), + BoundsCheckFreq(_BoundsCheckFreq){}; + }; + + NAMESPACE_END(Grid); diff --git a/Grid/qcd/action/pseudofermion/Bounds.h b/Grid/qcd/action/pseudofermion/Bounds.h index 535e1a49..c428d5ee 100644 --- a/Grid/qcd/action/pseudofermion/Bounds.h +++ b/Grid/qcd/action/pseudofermion/Bounds.h @@ -48,5 +48,56 @@ NAMESPACE_BEGIN(Grid); assert( (std::sqrt(Nd/Nx) 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 + 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; + + MultiShiftFunction ApproxPower ; //rational approx for X^{1/inv_pow} + MultiShiftFunction ApproxNegPower; //rational approx for X^{-1/inv_pow} + MultiShiftFunction ApproxHalfPower; //rational approx for X^{1/(2*inv_pow)} + MultiShiftFunction ApproxNegHalfPower; //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 + + public: + + GeneralEvenOddRatioRationalPseudoFermionAction(FermionOperator &_NumOp, + FermionOperator &_DenOp, + Params & p + ) : + NumOp(_NumOp), + DenOp(_DenOp), + PhiOdd (_NumOp.FermionRedBlackGrid()), + PhiEven(_NumOp.FermionRedBlackGrid()), + param(p) + { + AlgRemez remez(param.lo,param.hi,param.precision); + + int inv_pow = param.inv_pow; + int _2_inv_pow = 2*inv_pow; + + // MdagM^(+- 1/inv_pow) + 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); + + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); + + + // MdagM^1/(2*inv_pow) eta + SchurDifferentiableOperator MdagM(DenOp); + ConjugateGradientMultiShift msCG_M(param.MaxIter,ApproxHalfPower); + msCG_M(MdagM,etaOdd,tmp); + + // VdagV^-1/(2*inv_pow) MdagM^1/(2*inv_pow) eta + SchurDifferentiableOperator VdagV(NumOp); + ConjugateGradientMultiShift msCG_V(param.MaxIter,ApproxNegHalfPower); + 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/(2*inv_pow) Phi + SchurDifferentiableOperator VdagV(NumOp); + ConjugateGradientMultiShift msCG_V(param.MaxIter,ApproxHalfPower); + msCG_V(VdagV,PhiOdd,X); + + // MdagM^-1/(2*inv_pow) VdagV^1/(2*inv_pow) Phi + SchurDifferentiableOperator MdagM(DenOp); + ConjugateGradientMultiShift msCG_M(param.MaxIter,ApproxNegHalfPower); + msCG_M(MdagM,X,Y); + + // Randomly apply rational bounds checks. + if ( (rand()%param.BoundsCheckFreq)==0 ) { + FermionField gauss(NumOp.FermionRedBlackGrid()); + gauss = PhiOdd; + HighBoundCheck(MdagM,gauss,param.hi); + + InversePowerBoundsCheck(param.inv_pow,param.MaxIter,param.tolerance*100,MdagM,gauss,ApproxNegPower); + } + + // Phidag VdagV^1/(2*inv_pow) MdagM^-1/(2*inv_pow) MdagM^-1/(2*inv_pow) VdagV^1/(2*inv_pow) 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 = ApproxNegPower.poles.size(); + const int n_pv = ApproxHalfPower.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,ApproxHalfPower); + ConjugateGradientMultiShift msCG_M(param.MaxIter,ApproxNegPower); + + 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 #endif diff --git a/Grid/qcd/hmc/HMC.h b/Grid/qcd/hmc/HMC.h index f168b69a..9bf5261b 100644 --- a/Grid/qcd/hmc/HMC.h +++ b/Grid/qcd/hmc/HMC.h @@ -207,7 +207,7 @@ public: DeltaH = evolve_hmc_step(Ucopy); // Metropolis-Hastings test bool accept = true; - if (traj >= Params.StartTrajectory + Params.NoMetropolisUntil) { + if (Params.MetropolisTest && traj >= Params.StartTrajectory + Params.NoMetropolisUntil) { accept = metropolis_test(DeltaH); } else { std::cout << GridLogMessage << "Skipping Metropolis test" << std::endl; diff --git a/tests/hmc/Test_rhmc_EOWilsonRatioPowQuarter.cc b/tests/hmc/Test_rhmc_EOWilsonRatioPowQuarter.cc new file mode 100644 index 00000000..3618e422 --- /dev/null +++ b/tests/hmc/Test_rhmc_EOWilsonRatioPowQuarter.cc @@ -0,0 +1,139 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_rhmc_EOWilsonRatio.cc + + Copyright (C) 2015 + +Author: Peter Boyle +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 + +//This test is for the Wilson action with the determinant det( M^dag M)^1/4 +//testing the generic RHMC + +int main(int argc, char **argv) { + using namespace Grid; + ; + + Grid_init(&argc, &argv); + int threads = GridThread::GetThreads(); + // here make a routine to print all the relevant information on the run + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + // Typedefs to simplify notation + typedef GenericHMCRunner HMCWrapper; // Uses the default minimum norm + typedef WilsonImplR FermionImplPolicy; + typedef WilsonFermionR FermionAction; + typedef typename FermionAction::FermionField FermionField; + + + //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + HMCWrapper TheHMC; + + // Grid from the command line + TheHMC.Resources.AddFourDimGrid("gauge"); + + // Checkpointer definition + CheckpointerParameters CPparams; + CPparams.config_prefix = "ckpoint_lat"; + CPparams.rng_prefix = "ckpoint_rng"; + CPparams.saveInterval = 5; + CPparams.format = "IEEE64BIG"; + + TheHMC.Resources.LoadNerscCheckpointer(CPparams); + + RNGModuleParameters RNGpar; + RNGpar.serial_seeds = "1 2 3 4 5"; + RNGpar.parallel_seeds = "6 7 8 9 10"; + TheHMC.Resources.SetRNGSeeds(RNGpar); + + // Construct observables + typedef PlaquetteMod PlaqObs; + TheHMC.Resources.AddObservable(); + ////////////////////////////////////////////// + + ///////////////////////////////////////////////////////////// + // Collect actions, here use more encapsulation + // need wrappers of the fermionic classes + // that have a complex construction + // standard + RealD beta = 5.6 ; + WilsonGaugeActionR Waction(beta); + + auto GridPtr = TheHMC.Resources.GetCartesian(); + auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); + + // temporarily need a gauge field + LatticeGaugeField U(GridPtr); + + Real mass = -0.77; + Real pv = 0.0; + + // Can we define an overloaded operator that does not need U and initialises + // it with zeroes? + FermionAction DenOp(U, *GridPtr, *GridRBPtr, mass); + FermionAction NumOp(U, *GridPtr, *GridRBPtr, pv); + + + // 1/2+1/2 flavour + // RationalActionParams(int _inv_pow = 2, + // RealD _lo = 0.0, + // RealD _hi = 1.0, + // int _maxit = 1000, + // RealD tol = 1.0e-8, + // int _degree = 10, + // int _precision = 64, + // int _BoundsCheckFreq=20) + + + int inv_pow = 4; + RationalActionParams Params(inv_pow,1.0e-2,64.0,1000,1.0e-6,14,64,1); + + GeneralEvenOddRatioRationalPseudoFermionAction RHMC(NumOp,DenOp,Params); + + // Collect actions + ActionLevel Level1(1); + Level1.push_back(&RHMC); + + ActionLevel Level2(4); + Level2.push_back(&Waction); + + TheHMC.TheAction.push_back(Level1); + TheHMC.TheAction.push_back(Level2); + ///////////////////////////////////////////////////////////// + + // HMC parameters are serialisable + TheHMC.Parameters.MD.MDsteps = 20; + TheHMC.Parameters.MD.trajL = 1.0; + + TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file + TheHMC.Run(); + + Grid_finalize(); + +} // main + + + + +