diff --git a/lib/qcd/action/pseudofermion/ExactOneFlavourRatio.h b/lib/qcd/action/pseudofermion/ExactOneFlavourRatio.h new file mode 100644 index 00000000..38c7380c --- /dev/null +++ b/lib/qcd/action/pseudofermion/ExactOneFlavourRatio.h @@ -0,0 +1,230 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/action/pseudofermion/ExactOneFlavourRatio.h + +Copyright (C) 2017 + +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 */ + +///////////////////////////////////////////////////////////////// +// Implementation of exact one flavour algorithm (EOFA) // +// using fermion classes defined in: // +// Grid/qcd/action/fermion/DomainWallEOFAFermion.h (Shamir) // +// Grid/qcd/action/fermion/MobiusEOFAFermion.h (Mobius) // +// arXiv: 1403.1683, 1706.05843 // +///////////////////////////////////////////////////////////////// + +#ifndef QCD_PSEUDOFERMION_EXACT_ONE_FLAVOUR_RATIO_H +#define QCD_PSEUDOFERMION_EXACT_ONE_FLAVOUR_RATIO_H + +namespace Grid{ +namespace QCD{ + + /////////////////////////////////////////////////////////////// + // Exact one flavour implementation of DWF determinant ratio // + /////////////////////////////////////////////////////////////// + + template + class ExactOneFlavourRatioPseudoFermionAction : public Action + { + public: + INHERIT_IMPL_TYPES(Impl); + typedef OneFlavourRationalParams Params; + Params param; + MultiShiftFunction PowerNegHalf; + + private: + bool use_heatbath_forecasting; + AbstractEOFAFermion& Lop; // the basic LH operator + AbstractEOFAFermion& Rop; // the basic RH operator + SchurRedBlackDiagMooeeSolve Solver; + FermionField Phi; // the pseudofermion field for this trajectory + + public: + ExactOneFlavourRatioPseudoFermionAction(AbstractEOFAFermion& _Lop, AbstractEOFAFermion& _Rop, + OperatorFunction& S, Params& p, bool use_fc=false) : Lop(_Lop), Rop(_Rop), Solver(S), + Phi(_Lop.FermionGrid()), param(p), use_heatbath_forecasting(use_fc) + { + AlgRemez remez(param.lo, param.hi, param.precision); + + // MdagM^(+- 1/2) + std::cout << GridLogMessage << "Generating degree " << param.degree << " for x^(-1/2)" << std::endl; + remez.generateApprox(param.degree, 1, 2); + PowerNegHalf.Init(remez, param.tolerance, true); + }; + + virtual std::string action_name() { return "ExactOneFlavourRatioPseudoFermionAction"; } + + virtual std::string LogParameters() { + std::stringstream sstream; + sstream << GridLogMessage << "[" << action_name() << "] Low :" << param.lo << std::endl; + sstream << GridLogMessage << "[" << action_name() << "] High :" << param.hi << std::endl; + sstream << GridLogMessage << "[" << action_name() << "] Max iterations :" << param.MaxIter << std::endl; + sstream << GridLogMessage << "[" << action_name() << "] Tolerance :" << param.tolerance << std::endl; + sstream << GridLogMessage << "[" << action_name() << "] Degree :" << param.degree << std::endl; + sstream << GridLogMessage << "[" << action_name() << "] Precision :" << param.precision << std::endl; + return sstream.str(); + } + + // Spin projection + void spProj(const FermionField& in, FermionField& out, int sign, int Ls) + { + if(sign == 1){ for(int s=0; s tmp(2, Lop.FermionGrid()); + + // Use chronological inverter to forecast solutions across poles + std::vector prev_solns; + 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; + printf("Heatbath source vector: <\eta|\eta> = %1.15e\n", norm2(eta)); + + // \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()); + + // S = <\Phi|\Phi> + RealD action(norm2(Phi)); + + // 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); + G5R5(tmp[1], tmp[0]); + tmp[0] = zero; + Solver(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); + 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> + spProj(Phi, spProj_Phi, 1, Rop.Ls); + Rop.Omega(spProj_Phi, tmp[0], 1, 0); + G5R5(tmp[1], tmp[0]); + tmp[0] = zero; + Solver(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(); + + return action; + }; + + // EOFA pseudofermion force: see Eqns. (34)-(36) of arXiv:1706.05843 + virtual void deriv(const GaugeField& U, GaugeField& dSdU) + { + }; + }; +}} + +#endif diff --git a/lib/qcd/action/pseudofermion/PseudoFermion.h b/lib/qcd/action/pseudofermion/PseudoFermion.h index bccca3d4..133ebb7d 100644 --- a/lib/qcd/action/pseudofermion/PseudoFermion.h +++ b/lib/qcd/action/pseudofermion/PseudoFermion.h @@ -38,5 +38,6 @@ directory #include #include #include +#include #endif diff --git a/tests/debug/Test_heatbath_dwf_eofa.cc b/tests/debug/Test_heatbath_dwf_eofa.cc new file mode 100644 index 00000000..b77bc982 --- /dev/null +++ b/tests/debug/Test_heatbath_dwf_eofa.cc @@ -0,0 +1,102 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/debug/Test_heatbath_dwf_eofa.cc + + Copyright (C) 2017 + +Author: Peter Boyle +Author: paboyle +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 */ + +////////////////////////////////////////////////////////////////////////////////////////// +// This program sets up the initial pseudofermion field |Phi> = Meofa^{-1/2}*|eta>, and +// then uses this Phi to compute the action . +// If all is working, one should find that = . +////////////////////////////////////////////////////////////////////////////////////////// + +#include + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +// Parameters for test +const std::vector grid_dim = { 8, 8, 8, 8 }; +const int Ls = 8; +const int Npoles = 12; +const RealD mf = 0.01; +const RealD mpv = 1.0; +const RealD M5 = 1.8; + +int main(int argc, char** argv) +{ + Grid_init(&argc, &argv); + + int threads = GridThread::GetThreads(); + std::cout << GridLogMessage << "Grid is set up to use " << threads << " threads" << std::endl; + + // Initialize spacetime grid + std::cout << GridLogMessage << "Lattice dimensions: " << grid_dim << " Ls: " << Ls << std::endl; + GridCartesian* UGrid = SpaceTimeGrid::makeFourDimGrid(grid_dim, + GridDefaultSimd(Nd,vComplex::Nsimd()), GridDefaultMpi()); + GridRedBlackCartesian* UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridCartesian* FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls, UGrid); + GridRedBlackCartesian* FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, UGrid); + + // Set up RNGs + std::vector seeds4({1, 2, 3, 4}); + std::vector seeds5({5, 6, 7, 8}); + GridParallelRNG RNG5(FGrid); + RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); + RNG4.SeedFixedIntegers(seeds4); + + // Random gauge field + LatticeGaugeField Umu(UGrid); + SU3::HotConfiguration(RNG4, Umu); + + DomainWallEOFAFermionR Lop(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, mf, mpv, 0.0, -1, M5); + DomainWallEOFAFermionR Rop(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mpv, mf, mpv, -1.0, 1, M5); + + // Construct the action and test the heatbath (zero initial guess) + { + OneFlavourRationalParams Params(0.95, 100.0, 5000, 1.0e-12, Npoles); + ConjugateGradient CG(1.0e-12, 5000); + ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, false); + + Meofa.refresh(Umu, RNG5); + printf(" = %1.15e\n", Meofa.S(Umu)); + } + + // Construct the action and test the heatbath (forecasted initial guesses) + { + OneFlavourRationalParams Params(0.95, 100.0, 5000, 1.0e-12, Npoles); + ConjugateGradient CG(1.0e-12, 5000); + ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, true); + + Meofa.refresh(Umu, RNG5); + printf(" = %1.15e\n", Meofa.S(Umu)); + } + + return 0; +} diff --git a/tests/debug/Test_heatbath_dwf_eofa_gparity.cc b/tests/debug/Test_heatbath_dwf_eofa_gparity.cc new file mode 100644 index 00000000..5c9d4923 --- /dev/null +++ b/tests/debug/Test_heatbath_dwf_eofa_gparity.cc @@ -0,0 +1,108 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/debug/Test_heatbath_dwf_eofa.cc + + Copyright (C) 2017 + +Author: Peter Boyle +Author: paboyle +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 */ + +////////////////////////////////////////////////////////////////////////////////////////// +// This program sets up the initial pseudofermion field |Phi> = Meofa^{-1/2}*|eta>, and +// then uses this Phi to compute the action . +// If all is working, one should find that = . +////////////////////////////////////////////////////////////////////////////////////////// + +#include + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +typedef GparityWilsonImplR FermionImplPolicy; +typedef GparityDomainWallEOFAFermionR FermionAction; +typedef typename FermionAction::FermionField FermionField; + +// Parameters for test +const std::vector grid_dim = { 8, 8, 8, 8 }; +const int Ls = 8; +const int Npoles = 12; +const RealD mf = 0.01; +const RealD mpv = 1.0; +const RealD M5 = 1.8; + +int main(int argc, char** argv) +{ + Grid_init(&argc, &argv); + + int threads = GridThread::GetThreads(); + std::cout << GridLogMessage << "Grid is set up to use " << threads << " threads" << std::endl; + + // Initialize spacetime grid + std::cout << GridLogMessage << "Lattice dimensions: " << grid_dim << " Ls: " << Ls << std::endl; + GridCartesian* UGrid = SpaceTimeGrid::makeFourDimGrid(grid_dim, + GridDefaultSimd(Nd,vComplex::Nsimd()), GridDefaultMpi()); + GridRedBlackCartesian* UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridCartesian* FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls, UGrid); + GridRedBlackCartesian* FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, UGrid); + + // Set up RNGs + std::vector seeds4({1, 2, 3, 4}); + std::vector seeds5({5, 6, 7, 8}); + GridParallelRNG RNG5(FGrid); + RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); + RNG4.SeedFixedIntegers(seeds4); + + // Random gauge field + LatticeGaugeField Umu(UGrid); + SU3::HotConfiguration(RNG4, Umu); + + // GparityDomainWallFermionR::ImplParams params; + FermionAction::ImplParams params; + FermionAction Lop(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, mf, mpv, 0.0, -1, M5, params); + FermionAction Rop(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mpv, mf, mpv, -1.0, 1, M5, params); + + // Construct the action and test the heatbath (zero initial guess) + { + OneFlavourRationalParams Params(0.95, 100.0, 5000, 1.0e-12, Npoles); + ConjugateGradient CG(1.0e-12, 5000); + ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, false); + + Meofa.refresh(Umu, RNG5); + printf(" = %1.15e\n", Meofa.S(Umu)); + } + + // Construct the action and test the heatbath (forecasted initial guesses) + { + OneFlavourRationalParams Params(0.95, 100.0, 5000, 1.0e-12, Npoles); + ConjugateGradient CG(1.0e-12, 5000); + ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, true); + + Meofa.refresh(Umu, RNG5); + printf(" = %1.15e\n", Meofa.S(Umu)); + } + + return 0; +}