From 91405de3f71e4eafc0dd3c7c1b490ce9c4612c48 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 22 Jun 2018 12:14:37 +0200 Subject: [PATCH] Hadrons: new solver exposing fermion matrix and generic source/solve import/export --- extras/Hadrons/Global.hpp | 8 +-- extras/Hadrons/Makefile.am | 3 +- extras/Hadrons/Modules/MAction/DWF.hpp | 2 +- extras/Hadrons/Modules/MAction/Wilson.hpp | 2 +- .../Hadrons/Modules/MAction/WilsonClover.hpp | 2 +- extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp | 2 +- extras/Hadrons/Modules/MFermion/FreeProp.hpp | 2 +- extras/Hadrons/Modules/MFermion/GaugeProp.hpp | 21 +++++-- extras/Hadrons/Modules/MSolver/RBPrecCG.hpp | 6 +- extras/Hadrons/Modules/MSource/SeqGamma.hpp | 2 +- extras/Hadrons/Solver.hpp | 62 +++++++++++++++++++ 11 files changed, 92 insertions(+), 20 deletions(-) create mode 100644 extras/Hadrons/Solver.hpp diff --git a/extras/Hadrons/Global.hpp b/extras/Hadrons/Global.hpp index 72fbdf80..bae34a8d 100644 --- a/extras/Hadrons/Global.hpp +++ b/extras/Hadrons/Global.hpp @@ -93,17 +93,15 @@ typedef typename SImpl::Field ScalarField##suffix;\ typedef typename SImpl::Field PropagatorField##suffix; #define SOLVER_TYPE_ALIASES(FImpl, suffix)\ -typedef std::function SolverFn##suffix; +typedef Solver Solver##suffix; #define SINK_TYPE_ALIASES(suffix)\ typedef std::function SinkFn##suffix; -#define FGS_TYPE_ALIASES(FImpl, suffix)\ +#define FG_TYPE_ALIASES(FImpl, suffix)\ FERM_TYPE_ALIASES(FImpl, suffix)\ -GAUGE_TYPE_ALIASES(FImpl, suffix)\ -SOLVER_TYPE_ALIASES(FImpl, suffix) +GAUGE_TYPE_ALIASES(FImpl, suffix) // logger class HadronsLogger: public Logger diff --git a/extras/Hadrons/Makefile.am b/extras/Hadrons/Makefile.am index 3b945124..eb0bc3ad 100644 --- a/extras/Hadrons/Makefile.am +++ b/extras/Hadrons/Makefile.am @@ -15,16 +15,17 @@ libHadrons_adir = $(pkgincludedir)/Hadrons nobase_libHadrons_a_HEADERS = \ $(modules_hpp) \ Application.hpp \ + EigenPack.hpp \ Environment.hpp \ Exceptions.hpp \ Factory.hpp \ GeneticScheduler.hpp \ Global.hpp \ Graph.hpp \ - EigenPack.hpp \ Module.hpp \ Modules.hpp \ ModuleFactory.hpp \ + Solver.hpp \ VirtualMachine.hpp HadronsXmlRun_SOURCES = HadronsXmlRun.cc diff --git a/extras/Hadrons/Modules/MAction/DWF.hpp b/extras/Hadrons/Modules/MAction/DWF.hpp index f0dda4b0..94a05cc8 100644 --- a/extras/Hadrons/Modules/MAction/DWF.hpp +++ b/extras/Hadrons/Modules/MAction/DWF.hpp @@ -56,7 +56,7 @@ template class TDWF: public Module { public: - FGS_TYPE_ALIASES(FImpl,); + FG_TYPE_ALIASES(FImpl,); public: // constructor TDWF(const std::string name); diff --git a/extras/Hadrons/Modules/MAction/Wilson.hpp b/extras/Hadrons/Modules/MAction/Wilson.hpp index 1a119571..a0eabd82 100644 --- a/extras/Hadrons/Modules/MAction/Wilson.hpp +++ b/extras/Hadrons/Modules/MAction/Wilson.hpp @@ -54,7 +54,7 @@ template class TWilson: public Module { public: - FGS_TYPE_ALIASES(FImpl,); + FG_TYPE_ALIASES(FImpl,); public: // constructor TWilson(const std::string name); diff --git a/extras/Hadrons/Modules/MAction/WilsonClover.hpp b/extras/Hadrons/Modules/MAction/WilsonClover.hpp index ead3accc..bb2c5b59 100644 --- a/extras/Hadrons/Modules/MAction/WilsonClover.hpp +++ b/extras/Hadrons/Modules/MAction/WilsonClover.hpp @@ -59,7 +59,7 @@ template class TWilsonClover: public Module { public: - FGS_TYPE_ALIASES(FImpl,); + FG_TYPE_ALIASES(FImpl,); public: // constructor TWilsonClover(const std::string name); diff --git a/extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp b/extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp index 52e65322..88e4e85c 100644 --- a/extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp +++ b/extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp @@ -57,7 +57,7 @@ template class TZMobiusDWF: public Module { public: - FGS_TYPE_ALIASES(FImpl,); + FG_TYPE_ALIASES(FImpl,); public: // constructor TZMobiusDWF(const std::string name); diff --git a/extras/Hadrons/Modules/MFermion/FreeProp.hpp b/extras/Hadrons/Modules/MFermion/FreeProp.hpp index b1038ddd..0786eb2b 100644 --- a/extras/Hadrons/Modules/MFermion/FreeProp.hpp +++ b/extras/Hadrons/Modules/MFermion/FreeProp.hpp @@ -57,7 +57,7 @@ template class TFreeProp: public Module { public: - FGS_TYPE_ALIASES(FImpl,); + FG_TYPE_ALIASES(FImpl,); public: // constructor TFreeProp(const std::string name); diff --git a/extras/Hadrons/Modules/MFermion/GaugeProp.hpp b/extras/Hadrons/Modules/MFermion/GaugeProp.hpp index ee21cba9..fdd23766 100644 --- a/extras/Hadrons/Modules/MFermion/GaugeProp.hpp +++ b/extras/Hadrons/Modules/MFermion/GaugeProp.hpp @@ -35,6 +35,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include BEGIN_HADRONS_NAMESPACE @@ -76,7 +77,8 @@ template class TGaugeProp: public Module { public: - FGS_TYPE_ALIASES(FImpl,); + FG_TYPE_ALIASES(FImpl,); + SOLVER_TYPE_ALIASES(FImpl,); public: // constructor TGaugeProp(const std::string name); @@ -92,7 +94,7 @@ protected: virtual void execute(void); private: unsigned int Ls_; - SolverFn *solver_{nullptr}; + Solver *solver_{nullptr}; }; MODULE_REGISTER_TMP(GaugeProp, TGaugeProp, MFermion); @@ -147,7 +149,8 @@ void TGaugeProp::execute(void) std::string propName = (Ls_ == 1) ? getName() : (getName() + "_5d"); auto &prop = envGet(PropagatorField, propName); auto &fullSrc = envGet(PropagatorField, par().source); - auto &solver = envGet(SolverFn, par().solver); + auto &solver = envGet(Solver, par().solver); + auto &mat = solver.getFMat(); envGetTmp(FermionField, source); envGetTmp(FermionField, sol); @@ -155,11 +158,12 @@ void TGaugeProp::execute(void) LOG(Message) << "Inverting using solver '" << par().solver << "' on source '" << par().source << "'" << std::endl; for (unsigned int s = 0; s < Ns; ++s) - for (unsigned int c = 0; c < FImpl::Dimension; ++c) + for (unsigned int c = 0; c < FImpl::Dimension; ++c) { LOG(Message) << "Inversion for spin= " << s << ", color= " << c << std::endl; // source conversion for 4D sources + LOG(Message) << "Import source" << std::endl; if (!env().isObject5d(par().source)) { if (Ls_ == 1) @@ -169,7 +173,7 @@ void TGaugeProp::execute(void) else { PropToFerm(tmp, fullSrc, s, c); - make_5D(tmp, source, Ls_); + mat.ImportPhysicalFermionSource(tmp, source); } } // source conversion for 5D sources @@ -184,14 +188,19 @@ void TGaugeProp::execute(void) PropToFerm(source, fullSrc, s, c); } } + LOG(Message) << "Prepare source" << std::endl; + sol = source; + mat.Dminus(sol, source); sol = zero; + LOG(Message) << "Solve" << std::endl; solver(sol, source); + LOG(Message) << "Export solution" << std::endl; FermToProp(prop, sol, s, c); // create 4D propagators from 5D one if necessary if (Ls_ > 1) { PropagatorField &p4d = envGet(PropagatorField, getName()); - make_4D(sol, tmp, Ls_); + mat.ExportPhysicalFermionSolution(sol, tmp); FermToProp(p4d, tmp, s, c); } } diff --git a/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp b/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp index 206d44d1..31be621f 100644 --- a/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp +++ b/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp @@ -32,6 +32,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include #include BEGIN_HADRONS_NAMESPACE @@ -55,7 +56,8 @@ template class TRBPrecCG: public Module { public: - FGS_TYPE_ALIASES(FImpl,); + FG_TYPE_ALIASES(FImpl,); + SOLVER_TYPE_ALIASES(FImpl,); typedef FermionEigenPack EPack; typedef CoarseFermionEigenPack CoarseEPack; typedef std::shared_ptr> GuesserPt; @@ -175,7 +177,7 @@ void TRBPrecCG::setup(void) schurSolver(mat, source, sol, *guesser); }; - envCreate(SolverFn, getName(), Ls, solver); + envCreate(Solver, getName(), Ls, solver, mat); } diff --git a/extras/Hadrons/Modules/MSource/SeqGamma.hpp b/extras/Hadrons/Modules/MSource/SeqGamma.hpp index 4fe4dd76..2d8e1c31 100644 --- a/extras/Hadrons/Modules/MSource/SeqGamma.hpp +++ b/extras/Hadrons/Modules/MSource/SeqGamma.hpp @@ -71,7 +71,7 @@ template class TSeqGamma: public Module { public: - FGS_TYPE_ALIASES(FImpl,); + FG_TYPE_ALIASES(FImpl,); public: // constructor TSeqGamma(const std::string name); diff --git a/extras/Hadrons/Solver.hpp b/extras/Hadrons/Solver.hpp new file mode 100644 index 00000000..08867634 --- /dev/null +++ b/extras/Hadrons/Solver.hpp @@ -0,0 +1,62 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/EigenPack.hpp + +Copyright (C) 2015-2018 + +Author: Antonin Portelli + +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 Hadrons_Solver_hpp_ +#define Hadrons_Solver_hpp_ + +#include + +BEGIN_HADRONS_NAMESPACE + +template +class Solver +{ +public: + typedef typename FImpl::FermionField FermionField; + typedef FermionOperator FMat; + typedef std::function SolverFn; +public: + Solver(SolverFn fn, FMat &mat): mat_(mat), fn_(fn) {} + + void operator()(FermionField &sol, const FermionField &src) + { + fn_(sol, src); + } + + FMat & getFMat(void) + { + return mat_; + } +private: + FMat &mat_; + SolverFn fn_; +}; + +END_HADRONS_NAMESPACE + +#endif // Hadrons_Solver_hpp_