Hadrons: new solver exposing fermion matrix and generic source/solve import/export

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<void(FermionField##suffix &,\
+typedef Solver<FImpl> Solver##suffix;
-                      const FermionField##suffix &)> SolverFn##suffix;
 
 #define SINK_TYPE_ALIASES(suffix)\
 typedef std::function<SlicedPropagator##suffix\
                       (const PropagatorField##suffix &)> SinkFn##suffix;
 
-#define FGS_TYPE_ALIASES(FImpl, suffix)\
+#define FG_TYPE_ALIASES(FImpl, suffix)\
 FERM_TYPE_ALIASES(FImpl, suffix)\
-GAUGE_TYPE_ALIASES(FImpl, suffix)\
+GAUGE_TYPE_ALIASES(FImpl, suffix)
-SOLVER_TYPE_ALIASES(FImpl, suffix)
 
 // logger
 class HadronsLogger: public Logger

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

extras/Hadrons/Modules/MAction/DWF.hpp

@@ -56,7 +56,7 @@ template <typename FImpl>
 class TDWF: public Module<DWFPar>
 {
 public:
-    FGS_TYPE_ALIASES(FImpl,);
+    FG_TYPE_ALIASES(FImpl,);
 public:
     // constructor
     TDWF(const std::string name);

extras/Hadrons/Modules/MAction/Wilson.hpp

@@ -54,7 +54,7 @@ template <typename FImpl>
 class TWilson: public Module<WilsonPar>
 {
 public:
-    FGS_TYPE_ALIASES(FImpl,);
+    FG_TYPE_ALIASES(FImpl,);
 public:
     // constructor
     TWilson(const std::string name);

extras/Hadrons/Modules/MAction/WilsonClover.hpp

@@ -59,7 +59,7 @@ template <typename FImpl>
 class TWilsonClover: public Module<WilsonCloverPar>
 {
 public:
-    FGS_TYPE_ALIASES(FImpl,);
+    FG_TYPE_ALIASES(FImpl,);
 public:
     // constructor
     TWilsonClover(const std::string name);

extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp

@@ -57,7 +57,7 @@ template <typename FImpl>
 class TZMobiusDWF: public Module<ZMobiusDWFPar>
 {
 public:
-    FGS_TYPE_ALIASES(FImpl,);
+    FG_TYPE_ALIASES(FImpl,);
 public:
     // constructor
     TZMobiusDWF(const std::string name);

extras/Hadrons/Modules/MFermion/FreeProp.hpp

@@ -57,7 +57,7 @@ template <typename FImpl>
 class TFreeProp: public Module<FreePropPar>
 {
 public:
-    FGS_TYPE_ALIASES(FImpl,);
+    FG_TYPE_ALIASES(FImpl,);
 public:
     // constructor
     TFreeProp(const std::string name);

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 <Grid/Hadrons/Global.hpp>
 #include <Grid/Hadrons/Module.hpp>
 #include <Grid/Hadrons/ModuleFactory.hpp>
+#include <Grid/Hadrons/Solver.hpp>
 
 BEGIN_HADRONS_NAMESPACE
 
@@ -76,7 +77,8 @@ template <typename FImpl>
 class TGaugeProp: public Module<GaugePropPar>
 {
 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<FIMPL>, MFermion);
@@ -147,7 +149,8 @@ void TGaugeProp<FImpl>::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<FImpl>::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<FImpl>::execute(void)
             else
             {
                 PropToFerm<FImpl>(tmp, fullSrc, s, c);
-                make_5D(tmp, source, Ls_);
+                mat.ImportPhysicalFermionSource(tmp, source);
             }
         }
         // source conversion for 5D sources
@@ -184,14 +188,19 @@ void TGaugeProp<FImpl>::execute(void)
                 PropToFerm<FImpl>(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<FImpl>(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<FImpl>(p4d, tmp, s, c);
         }
     }

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 <Grid/Hadrons/Global.hpp>
 #include <Grid/Hadrons/Module.hpp>
 #include <Grid/Hadrons/ModuleFactory.hpp>
+#include <Grid/Hadrons/Solver.hpp>
 #include <Grid/Hadrons/EigenPack.hpp>
 
 BEGIN_HADRONS_NAMESPACE
@@ -55,7 +56,8 @@ template <typename FImpl, int nBasis>
 class TRBPrecCG: public Module<RBPrecCGPar>
 {
 public:
-    FGS_TYPE_ALIASES(FImpl,);
+    FG_TYPE_ALIASES(FImpl,);
+    SOLVER_TYPE_ALIASES(FImpl,);
     typedef FermionEigenPack<FImpl>                       EPack;
     typedef CoarseFermionEigenPack<FImpl, nBasis>         CoarseEPack;
     typedef std::shared_ptr<Guesser<FermionField>>        GuesserPt;
@@ -175,7 +177,7 @@ void TRBPrecCG<FImpl, nBasis>::setup(void)
         
         schurSolver(mat, source, sol, *guesser);
     };
-    envCreate(SolverFn, getName(), Ls, solver);
+    envCreate(Solver, getName(), Ls, solver, mat);
 }
 
 

extras/Hadrons/Modules/MSource/SeqGamma.hpp

@@ -71,7 +71,7 @@ template <typename FImpl>
 class TSeqGamma: public Module<SeqGammaPar>
 {
 public:
-    FGS_TYPE_ALIASES(FImpl,);
+    FG_TYPE_ALIASES(FImpl,);
 public:
     // constructor
     TSeqGamma(const std::string name);

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 <antonin.portelli@me.com>
+
+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 <Grid/Hadrons/Global.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+template <typename FImpl>
+class Solver
+{
+public:
+    typedef typename FImpl::FermionField                      FermionField;
+    typedef FermionOperator<FImpl>                            FMat; 
+    typedef std::function<void(FermionField &, 
+                               const FermionField &)>         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_