Hadrons: first prototype with working inversions

programs/Hadrons/AWilson.cc

@@ -0,0 +1,56 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/AWilson.cc
+
+Copyright (C) 2016
+
+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.
+*******************************************************************************/
+
+#include <Hadrons/AWilson.hpp>
+#include <Hadrons/Environment.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+
+/******************************************************************************
+*                          AWilson implementation                             *
+******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+AWilson::AWilson(const std::string name)
+: FermionAction(name)
+{}
+
+// parse parameters ////////////////////////////////////////////////////////////
+void AWilson::parseParameters(XmlReader &reader, const std::string name)
+{
+   read(reader, name, par_);
+}
+
+// create operator /////////////////////////////////////////////////////////////
+void AWilson::create(Environment &env)
+{
+    auto &U      = *env.getGauge();
+    auto &grid   = *env.getGrid();
+    auto &gridRb = *env.getRbGrid();
+    
+    setFMat(new WilsonFermionR(U, grid, gridRb, par_.mass));
+}

programs/Hadrons/AWilson.hpp

@@ -0,0 +1,65 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/AWilson.hpp
+
+Copyright (C) 2016
+
+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.
+*******************************************************************************/
+
+#ifndef Hadrons_AWilson_hpp_
+#define Hadrons_AWilson_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/FermionAction.hpp>
+#include <Hadrons/FermionActionFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                           Wilson fermions                                  *
+ ******************************************************************************/
+class AWilson: public FermionAction
+{
+public:
+    class Par: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Par, double, mass);
+    };
+public:
+    // constructor
+    AWilson(const std::string name);
+    // destructor
+    virtual ~AWilson(void) = default;
+    // parse parameters
+    virtual void parseParameters(XmlReader &reader, const std::string name);
+    // create operator
+    virtual void create(Environment &env);
+private:
+    Par par_;
+};
+
+ACTION_REGISTER(AWilson);
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_AWilson_hpp_

programs/Hadrons/Application.cc

@@ -38,10 +38,17 @@ using namespace Hadrons;
 Application::Application(const std::string parameterFileName)
 : parameterFileName_(parameterFileName)
 , env_(Environment::getInstance())
+, actionFactory_(FermionActionFactory::getInstance())
 , modFactory_(ModuleFactory::getInstance())
 {
+    LOG(Message) << "Fermion actions available:" << std::endl;
+    auto list = actionFactory_.getBuilderList();
+    for (auto &m: list)
+    {
+        LOG(Message) << "  " << m << std::endl;
+    }
     LOG(Message) << "Modules available:" << std::endl;
-    auto list = modFactory_.getModuleList();
+    list = modFactory_.getBuilderList();
     for (auto &m: list)
     {
         LOG(Message) << "  " << m << std::endl;
@@ -61,10 +68,10 @@ void Application::run(void)
 }
 
 // parse parameter file ////////////////////////////////////////////////////////
-class ModuleId: Serializable
+class ObjectId: Serializable
 {
 public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(ModuleId,
+    GRID_SERIALIZABLE_CLASS_MEMBERS(ObjectId,
                                     std::string, name,
                                     std::string, type);
 };
@@ -72,10 +79,22 @@ public:
 void Application::parseParameterFile(void)
 {
     XmlReader reader(parameterFileName_);
-    ModuleId  id;
+    ObjectId  id;
     
     LOG(Message) << "Reading '" << parameterFileName_ << "'..." << std::endl;
     read(reader, "parameters", par_);
+    push(reader, "actions");
+    push(reader, "action");
+    do
+    {
+        read(reader, "id", id);
+        env_.addFermionAction(actionFactory_.create(id.type, id.name));
+        auto &action = *env_.getFermionAction(id.name);
+        action.parseParameters(reader, "options");
+        action.create(env_);
+    } while (reader.nextElement("action"));
+    pop(reader);
+    pop(reader);
     push(reader, "modules");
     push(reader, "module");
     do
@@ -92,6 +111,7 @@ void Application::parseParameterFile(void)
     } while (reader.nextElement("module"));
     pop(reader);
     pop(reader);
+    env_.setSeed(strToVec<int>(par_.seed));
 }
 
 // schedule computation ////////////////////////////////////////////////////////
@@ -122,33 +142,38 @@ void Application::schedule(void)
     unsigned int k = 0;
     
     std::vector<Graph<std::string>> con = moduleGraph.getConnectedComponents();
-    LOG(Message) << "Program:" << std::endl;
+    
     for (unsigned int i = 0; i < con.size(); ++i)
     {
-        std::vector<std::vector<std::string>> t = con[i].allTopoSort();
+//        std::vector<std::vector<std::string>> t = con[i].allTopoSort();
-        int                                   memPeak, minMemPeak = -1;
+//        int                                   memPeak, minMemPeak = -1;
-        unsigned int                          bestInd;
+//        unsigned int                          bestInd;
-        bool                                  msg;
+//        bool                                  msg;
-
+//
-        env_.dryRun(true);
+//        LOG(Message) << "analyzing " << t.size() << " possible programs..."
-        for (unsigned int p = 0; p < t.size(); ++p)
+//                     << std::endl;
+//        env_.dryRun(true);
+//        for (unsigned int p = 0; p < t.size(); ++p)
+//        {
+//            msg = HadronsLogMessage.isActive();
+//            HadronsLogMessage.Active(false);
+//        
+//            memPeak = execute(t[p]);
+//            if ((memPeak < minMemPeak) or (minMemPeak < 0))
+//            {
+//                minMemPeak = memPeak;
+//                bestInd = p;
+//            }
+//            HadronsLogMessage.Active(msg);
+//            env_.freeAll();
+//        }
+//        env_.dryRun(false);
+        std::vector<std::string> t = con[i].topoSort();
+        LOG(Message) << "Program " << i + 1 << ":" << std::endl;
+        for (unsigned int j = 0; j < t.size(); ++j)
         {
-            msg = HadronsLogMessage.isActive();
+            program_.push_back(t[j]);
-            HadronsLogMessage.Active(false);
+            LOG(Message) << std::setw(4) << std::right << k + 1 << ": "
-            memPeak = execute(t[p]);
-            if ((memPeak < minMemPeak) or (minMemPeak < 0))
-            {
-                minMemPeak = memPeak;
-                bestInd = p;
-            }
-            HadronsLogMessage.Active(msg);
-            env_.freeAll();
-        }
-        env_.dryRun(false);
-        for (unsigned int j = 0; j < t[bestInd].size(); ++j)
-        {
-            program_.push_back(t[bestInd][j]);
-            LOG(Message) << std::setw(4) << std::right << k << ": "
                          << program_[k] << std::endl;
             k++;
         }
@@ -164,7 +189,7 @@ void Application::configLoop(void)
     {
         LOG(Message) << "Starting measurement for trajectory " << t
                      << std::endl;
-        env_.loadUnitGauge();
+        env_.loadRandomGauge();
         execute(program_);
         env_.freeAll();
     }
@@ -172,8 +197,25 @@ void Application::configLoop(void)
 
 unsigned int Application::execute(const std::vector<std::string> &program)
 {
-    unsigned int memPeak = 0;
+    unsigned int                          memPeak = 0;
+    std::vector<std::vector<std::string>> freeProg;
     
+    freeProg.resize(program.size());
+    for (auto &n: associatedModule_)
+    {
+        auto pred = [&n, this](const std::string &s)
+        {
+            auto &in = input_[s];
+            auto it  = std::find(in.begin(), in.end(), n.first);
+            
+            return (it != in.end()) or (s == n.second);
+        };
+        auto it = std::find_if(program.rbegin(), program.rend(), pred);
+        if (it != program.rend())
+        {
+            freeProg[program.rend() - it - 1].push_back(n.first);
+        }
+    }
     for (unsigned int i = 0; i < program.size(); ++i)
     {
         LOG(Message) << "Measurement step (" << i+1 << "/" << program.size()
@@ -184,21 +226,9 @@ unsigned int Application::execute(const std::vector<std::string> &program)
         {
             memPeak = env_.nProp();
         }
-        for (auto &n: associatedModule_)
+        for (auto &n: freeProg[i])
         {
-            bool canFree = true;
+            env_.free(n);
-            
-            for (unsigned int j = i + 1; j < program.size(); ++j)
-            {
-                auto &in = input_[program[j]];
-                auto it  = std::find(in.begin(), in.end(), n.first);
-                canFree  = canFree and (it == in.end());
-            }
-            if (canFree and env_.propExists(n.first))
-            {
-                LOG(Message) << "freeing '" << n.first << "'" << std::endl;
-                env_.free(n.first);
-            }
         }
     }
     

programs/Hadrons/Application.hpp

@@ -29,6 +29,7 @@ directory.
 #define Hadrons_Application_hpp_
 
 #include <Hadrons/Global.hpp>
+#include <Hadrons/FermionActionFactory.hpp>
 #include <Hadrons/Environment.hpp>
 #include <Hadrons/ModuleFactory.hpp>
 
@@ -59,7 +60,8 @@ class GlobalPar: Serializable
 {
 public:
     GRID_SERIALIZABLE_CLASS_MEMBERS(GlobalPar,
-                                    ConfigPar, configs);
+                                    ConfigPar,   configs,
+                                    std::string, seed);
 };
 
 /******************************************************************************
@@ -88,11 +90,13 @@ private:
     std::string                                     parameterFileName_;
     GlobalPar                                       par_;
     Environment                                     &env_;
+    FermionActionFactory                            &actionFactory_;
     ModuleFactory                                   &modFactory_;
     std::map<std::string, std::unique_ptr<Module>>  module_;
     std::map<std::string, std::string>              associatedModule_;
     std::map<std::string, std::vector<std::string>> input_;
     std::vector<std::string>                        program_;
+    std::vector<std::vector<std::string>>           freeProg_;
 };
 
 END_HADRONS_NAMESPACE

programs/Hadrons/CMeson.cc

@@ -59,10 +59,6 @@ std::vector<std::string> CMeson::getOutput(void)
     return output;
 }
 
-// memory footprint ////////////////////////////////////////////////////////////
-void CMeson::allocate(Environment &env)
-{}
-
 // execution ///////////////////////////////////////////////////////////////////
 void CMeson::execute(Environment &env)
 {

programs/Hadrons/CMeson.hpp

@@ -58,8 +58,6 @@ public:
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
-    // allocation
-    virtual void allocate(Environment &env);
     // execution
     virtual void execute(Environment &env);
 private:

programs/Hadrons/Environment.cc

@@ -37,15 +37,13 @@ using namespace Hadrons;
 // constructor /////////////////////////////////////////////////////////////////
 Environment::Environment(void)
 {
-    std::vector<int> seed4d({1,2,3,4});
-    
     grid4d_.reset(SpaceTimeGrid::makeFourDimGrid(
         GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()),
         GridDefaultMpi()));
     gridRb4d_.reset(SpaceTimeGrid::makeFourDimRedBlackGrid(grid4d_.get()));
     rng4d_.reset(new GridParallelRNG(grid4d_.get()));
-    rng4d_->SeedFixedIntegers(seed4d);
     gauge_.reset(new LatticeGaugeField(grid4d_.get()));
+    loadUnitGauge();
 }
 
 // dry run /////////////////////////////////////////////////////////////////////
@@ -54,27 +52,24 @@ void Environment::dryRun(const bool isDry)
     dryRun_ = isDry;
 }
 
-bool Environment::isDryRun(void)
+bool Environment::isDryRun(void) const
 {
     return dryRun_;
 }
 
 // grids ///////////////////////////////////////////////////////////////////////
-GridCartesian * Environment::get4dGrid(void)
+GridCartesian * Environment::getGrid(const unsigned int Ls) const
-{
-    return grid4d_.get();
-}
-
-GridRedBlackCartesian * Environment::getRb4dGrid(void)
-{
-    return gridRb4d_.get();
-}
-
-GridCartesian * Environment::get5dGrid(const unsigned int Ls)
 {
     try
     {
-        return grid5d_.at(Ls).get();
+        if (Ls == 1)
+        {
+            return grid4d_.get();
+        }
+        else
+        {
+            return grid5d_.at(Ls).get();
+        }
     }
     catch(std::out_of_range &)
     {
@@ -82,11 +77,18 @@ GridCartesian * Environment::get5dGrid(const unsigned int Ls)
     }
 }
 
-GridRedBlackCartesian * Environment::getRb5dGrid(const unsigned int Ls)
+GridRedBlackCartesian * Environment::getRbGrid(const unsigned int Ls) const
 {
     try
     {
-        return gridRb5d_.at(Ls).get();
+        if (Ls == 1)
+        {
+            return gridRb4d_.get();
+        }
+        else
+        {
+            return gridRb5d_.at(Ls).get();
+        }
     }
     catch(std::out_of_range &)
     {
@@ -94,6 +96,52 @@ GridRedBlackCartesian * Environment::getRb5dGrid(const unsigned int Ls)
     }
 }
 
+// fermion actions /////////////////////////////////////////////////////////////
+void Environment::addFermionAction(FActionPt action)
+{
+    fAction_[action->getName()] = std::move(action);
+}
+
+FermionAction * Environment::getFermionAction(const std::string name) const
+{
+    try
+    {
+        return fAction_.at(name).get();
+    }
+    catch(std::out_of_range &)
+    {
+        try
+        {
+            return fAction_.at(solverAction_.at(name)).get();
+        }
+        catch (std::out_of_range &)
+        {
+            HADRON_ERROR("no action with name '" << name << "'");
+        }
+    }
+}
+
+// solvers /////////////////////////////////////////////////////////////////////
+void Environment::addSolver(const std::string name, Solver s,
+                            const std::string actionName)
+{
+    solver_[name]       = s;
+    solverAction_[name] = actionName;
+}
+
+void Environment::callSolver(const std::string name, LatticeFermion &sol,
+                             const LatticeFermion &source) const
+{
+    try
+    {
+        solver_.at(name)(sol, source);
+    }
+    catch(std::out_of_range &)
+    {
+        HADRON_ERROR("no solver with name '" << name << "'");
+    }
+}
+
 // quark propagators ///////////////////////////////////////////////////////////
 void Environment::addProp(const std::string name, const unsigned int Ls)
 {
@@ -154,11 +202,46 @@ void Environment::freeProp(const std::string name)
     }
 }
 
-LatticePropagator * Environment::getProp(const std::string name)
+bool Environment::isProp5d(const std::string name) const
 {
     if (propExists(name))
     {
-        return prop_[name].get();
+        return (getProp(name)->_grid->GlobalDimensions().size() == Nd + 1);
+    }
+    else
+    {
+        HADRON_ERROR("propagator '" + name + "' unknown");
+        
+        return false;
+    }
+}
+
+unsigned int Environment::getPropLs(const std::string name) const
+{
+    if (propExists(name))
+    {
+        if (isProp5d(name))
+        {
+            return getProp(name)->_grid->GlobalDimensions()[0];
+        }
+        else
+        {
+            return 1;
+        }
+    }
+    else
+    {
+        HADRON_ERROR("propagator '" + name + "' unknown");
+        
+        return 0;
+    }
+}
+
+LatticePropagator * Environment::getProp(const std::string name) const
+{
+    if (propExists(name))
+    {
+        return prop_.at(name).get();
     }
     else
     {
@@ -168,21 +251,12 @@ LatticePropagator * Environment::getProp(const std::string name)
     }
 }
 
-bool Environment::propExists(const std::string name)
+bool Environment::propExists(const std::string name) const
 {
-    auto it = prop_.find(name);
+    return (prop_.find(name) != prop_.end());
-    
-    if (it == prop_.end())
-    {
-        return false;
-    }
-    else
-    {
-        return true;
-    }
 }
 
-unsigned int Environment::nProp(void)
+unsigned int Environment::nProp(void) const
 {
     unsigned int size = 0;
     
@@ -195,7 +269,7 @@ unsigned int Environment::nProp(void)
 }
 
 // gauge configuration /////////////////////////////////////////////////////////
-LatticeGaugeField * Environment::getGauge(void)
+LatticeGaugeField * Environment::getGauge(void) const
 {
     return gauge_.get();
 }
@@ -210,11 +284,23 @@ void Environment::loadRandomGauge(void)
     SU3::HotConfiguration(*rng4d_, *gauge_);
 }
 
+// random number generator /////////////////////////////////////////////////////
+void Environment::setSeed(const std::vector<int> &seed)
+{
+    rng4d_->SeedFixedIntegers(seed);
+}
+
+GridParallelRNG * Environment::get4dRng(void) const
+{
+    return rng4d_.get();
+}
+
 // general free ////////////////////////////////////////////////////////////////
 void Environment::free(const std::string name)
 {
     if (propExists(name))
     {
+        LOG(Message) << "freeing '" << name << "'" << std::endl;
         freeProp(name);
     }
 }

programs/Hadrons/Environment.hpp

@@ -29,6 +29,7 @@ directory.
 #define Hadrons_Environment_hpp_
 
 #include <Hadrons/Global.hpp>
+#include <Hadrons/FermionAction.hpp>
 
 BEGIN_HADRONS_NAMESPACE
 
@@ -39,31 +40,46 @@ class Environment
 {
     SINGLETON(Environment);
 public:
-    typedef std::unique_ptr<GridCartesian>         GridPt;
+    typedef std::function<void(LatticeFermion &,
-    typedef std::unique_ptr<GridRedBlackCartesian> GridRbPt;
+                               const LatticeFermion &)> Solver;
-    typedef std::unique_ptr<GridParallelRNG>       RngPt;
+    typedef std::unique_ptr<GridCartesian>              GridPt;
-    typedef std::unique_ptr<LatticePropagator>     PropPt;
+    typedef std::unique_ptr<GridRedBlackCartesian>      GridRbPt;
-    typedef std::unique_ptr<LatticeGaugeField>     GaugePt;
+    typedef std::unique_ptr<GridParallelRNG>            RngPt;
+    typedef std::unique_ptr<FermionAction>              FActionPt;
+    typedef std::unique_ptr<LatticePropagator>          PropPt;
+    typedef std::unique_ptr<LatticeGaugeField>          GaugePt;
 public:
     // dry run
     void                    dryRun(const bool isDry);
-    bool                    isDryRun(void);
+    bool                    isDryRun(void) const;
     // grids
-    GridCartesian *         get4dGrid(void);
+    GridCartesian *         getGrid(const unsigned int Ls = 1) const;
-    GridRedBlackCartesian * getRb4dGrid(void);
+    GridRedBlackCartesian * getRbGrid(const unsigned int Ls = 1) const;
-    GridCartesian *         get5dGrid(const unsigned int Ls);
+    // fermion actions
-    GridRedBlackCartesian * getRb5dGrid(const unsigned int Ls);
+    void                    addFermionAction(FActionPt action);
+    FermionAction *         getFermionAction(const std::string name) const;
+    // solvers
+    void                    addSolver(const std::string name, Solver s,
+                                      const std::string actionName);
+    void                    callSolver(const std::string name,
+                                       LatticeFermion &sol,
+                                       const LatticeFermion &src) const;
     // quark propagators
     void                    addProp(const std::string name,
                                     const unsigned int Ls = 1);
     void                    freeProp(const std::string name);
-    LatticePropagator *     getProp(const std::string name);
+    bool                    isProp5d(const std::string name) const;
-    bool                    propExists(const std::string name);
+    unsigned int            getPropLs(const std::string name) const;
-    unsigned int            nProp(void);
+    LatticePropagator *     getProp(const std::string name) const;
+    bool                    propExists(const std::string name) const;
+    unsigned int            nProp(void) const;
     // gauge configuration
-    LatticeGaugeField *     getGauge(void);
+    LatticeGaugeField *     getGauge(void) const;
     void                    loadUnitGauge(void);
     void                    loadRandomGauge(void);
+    // random number generator
+    void                    setSeed(const std::vector<int> &seed);
+    GridParallelRNG *       get4dRng(void) const;
     // general free
     void                    free(const std::string name);
     void                    freeAll(void);
@@ -74,6 +90,9 @@ private:
     GridRbPt                            gridRb4d_;
     std::map<unsigned int, GridRbPt>    gridRb5d_;
     RngPt                               rng4d_;
+    std::map<std::string, FActionPt>    fAction_;
+    std::map<std::string, Solver>       solver_;
+    std::map<std::string, std::string>  solverAction_;
     std::map<std::string, PropPt>       prop_;
     std::map<std::string, unsigned int> propSize_;
     GaugePt                             gauge_;

programs/Hadrons/Factory.hpp

@@ -0,0 +1,104 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/Factory.hpp
+
+Copyright (C) 2016
+
+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.
+*******************************************************************************/
+
+#ifndef Hadrons_Factory_hpp_
+#define Hadrons_Factory_hpp_
+
+#include <Hadrons/Global.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                        abstract factory class                              *
+ ******************************************************************************/
+template <typename T>
+class Factory
+{
+public:
+    typedef std::function<std::unique_ptr<T>(const std::string)> Func;
+public:
+    // constructor
+    Factory(void) = default;
+    // destructor
+    virtual ~Factory(void) = default;
+    // registration
+    void registerBuilder(const std::string type, const Func &f);
+    // get module list
+    std::vector<std::string> getBuilderList(void) const;
+    // factory
+    std::unique_ptr<T> create(const std::string type,
+                              const std::string name) const;
+private:
+    std::map<std::string, Func> builder_;
+};
+
+/******************************************************************************
+ *                         template implementation                            *
+ ******************************************************************************/
+// registration ////////////////////////////////////////////////////////////////
+template <typename T>
+void Factory<T>::registerBuilder(const std::string type, const Func &f)
+{
+    builder_[type] = f;
+}
+
+// get module list /////////////////////////////////////////////////////////////
+template <typename T>
+std::vector<std::string> Factory<T>::getBuilderList(void) const
+{
+    std::vector<std::string> list;
+    
+    for (auto &b: builder_)
+    {
+        list.push_back(b.first);
+    }
+    
+    return list;
+}
+
+// factory /////////////////////////////////////////////////////////////////////
+template <typename T>
+std::unique_ptr<T> Factory<T>::create(const std::string type,
+                                      const std::string name) const
+{
+    Func func;
+    
+    try
+    {
+        func = builder_.at(type);
+    }
+    catch (std::out_of_range)
+    {
+        HADRON_ERROR("object of type '" + type + "' unknown");
+    }
+    
+    return func(name);
+}
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Factory_hpp_

programs/Hadrons/FermionAction.cc

@@ -1,9 +1,9 @@
 /*******************************************************************************
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: programs/Hadrons/ModuleFactory.cc
+Source file: programs/Hadrons/FermionAction.cc
 
-Copyright (C) 2015
+Copyright (C) 2016
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 
@@ -25,53 +25,36 @@ See the full license in the file "LICENSE" in the top level distribution
 directory.
 *******************************************************************************/
 
-#include <Hadrons/ModuleFactory.hpp>
+#include <Hadrons/FermionAction.hpp>
 
 using namespace Grid;
 using namespace Hadrons;
 
 /******************************************************************************
- *                      ModuleFactory implementation                          *
+ *                      FermionAction implementation                          *
  ******************************************************************************/
 // constructor /////////////////////////////////////////////////////////////////
-ModuleFactory::ModuleFactory(void)
+FermionAction::FermionAction(const std::string name)
+: name_(name)
+{}
+
+// access //////////////////////////////////////////////////////////////////////
+std::string FermionAction::getName(void) const
 {
+    return name_;
 }
 
-// registration ////////////////////////////////////////////////////////////////
+unsigned int FermionAction::getLs(void) const
-void ModuleFactory::registerModule(const std::string type,
-                                   const FactoryFunc &f)
 {
-    factory_[type] = f;
+    return 1;
 }
 
-// get module list /////////////////////////////////////////////////////////////
+void FermionAction::setFMat(FMat *fMat)
-std::vector<std::string> ModuleFactory::getModuleList(void) const
 {
-    std::vector<std::string> list;
+    fMat_.reset(fMat);
-    
-    for (auto &f: factory_)
-    {
-        list.push_back(f.first);
-    }
-    
-    return list;
 }
 
-// factory /////////////////////////////////////////////////////////////////////
+FermionAction::FMat * FermionAction::getFMat(void)
-std::unique_ptr<Module> ModuleFactory::create(const std::string type,
-                                              const std::string name) const
 {
-    FactoryFunc func;
+    return fMat_.get();
-    
-    try
-    {
-        func = factory_.at(type);
-    }
-    catch (std::out_of_range)
-    {
-        HADRON_ERROR("module type '" + type + "' unknown");
-    }
-    
-    return func(name);
 }

programs/Hadrons/FermionAction.hpp

@@ -0,0 +1,83 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/FermionAction.hpp
+
+Copyright (C) 2016
+
+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.
+*******************************************************************************/
+
+#ifndef Hadrons_FermionAction_hpp_
+#define Hadrons_FermionAction_hpp_
+
+#include <Hadrons/Global.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+// action registration macro
+#define ACTION_REGISTER(action)\
+class action##ActionRegistrar\
+{\
+public:\
+    action##ActionRegistrar(void)\
+    {\
+        FermionActionFactory &actionFac = FermionActionFactory::getInstance();\
+        actionFac.registerBuilder(#action, [&](const std::string name)\
+                                  {\
+                                      return std::unique_ptr<action>(\
+                                          new action(name));\
+                                  });\
+    }\
+};\
+static action##ActionRegistrar action##ActionRegistrarInstance;
+
+/******************************************************************************
+ *                             FermionAction                                  *
+ ******************************************************************************/
+class Environment;
+
+class FermionAction
+{
+public:
+    typedef FermionOperator<WilsonImplR> FMat;
+    typedef std::unique_ptr<FMat>        FMatPt;
+public:
+    // constructor
+    FermionAction(const std::string name);
+    // destructor
+    virtual ~FermionAction(void) = default;
+    // access
+            std::string  getName(void) const;
+    virtual unsigned int getLs(void) const;
+            void         setFMat(FMat *fMat);
+            FMat *       getFMat(void);
+    // parse parameters
+    virtual void parseParameters(XmlReader &reader, const std::string name) = 0;
+    // create operator
+    virtual void create(Environment &env) = 0;
+private:
+    std::string name_;
+    FMatPt      fMat_;
+};
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_FermionAction_hpp_

programs/Hadrons/FermionActionFactory.hpp

@@ -0,0 +1,47 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/FermionActionFactory.hpp
+
+Copyright (C) 2016
+
+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.
+*******************************************************************************/
+
+#ifndef Hadrons_FermionActionFactory_hpp_
+#define Hadrons_FermionActionFactory_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Factory.hpp>
+#include <Hadrons/FermionAction.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         FermionActionFactory                               *
+ ******************************************************************************/
+class FermionActionFactory: public Factory<FermionAction>
+{
+    SINGLETON_DEFCTOR(FermionActionFactory)
+};
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_FermionActionFactory_hpp_

programs/Hadrons/Global.hpp

@@ -49,7 +49,8 @@ public:
 
 #define LOG(channel) std::cout << HadronsLog##channel
 #define HADRON_ERROR(msg)\
-LOG(Error) << msg << std::endl;\
+LOG(Error) << msg << " (" << __FUNCTION__ << " at " << __FILE__ << ":"\
+           << __LINE__ << ")" << std::endl;\
 abort();
 
 #define DEBUG_VAR(var) LOG(Debug) << #var << "= " << (var) << std::endl;
@@ -72,6 +73,18 @@ public:\
 private:\
     name(void);
 
+#define SINGLETON_DEFCTOR(name)\
+public:\
+    name(const name &e) = delete;\
+    void operator=(const name &e) = delete;\
+    static name & getInstance(void)\
+    {\
+        static name e;\
+        return e;\
+    }\
+private:\
+    name(void) = default;
+
 END_HADRONS_NAMESPACE
 
 #endif // Hadrons_Global_hpp_

programs/Hadrons/MQuark.cc

@@ -47,7 +47,9 @@ void MQuark::parseParameters(XmlReader &reader, const std::string name)
 // dependency relation
 std::vector<std::string> MQuark::getInput(void)
 {
-    return std::vector<std::string>();
+    std::vector<std::string> in = {par_.source, par_.solver};
+    
+    return in;
 }
 
 std::vector<std::string> MQuark::getOutput(void)
@@ -57,14 +59,20 @@ std::vector<std::string> MQuark::getOutput(void)
     return out;
 }
 
+// setup ///////////////////////////////////////////////////////////////////////
+void MQuark::setup(Environment &env)
+{
+    Ls_ = env.getFermionAction(par_.solver)->getLs();
+}
+
 // allocation //////////////////////////////////////////////////////////////////
 void MQuark::allocate(Environment &env)
 {
     env.addProp(getName());
     quark_ = env.getProp(getName());
-    if (par_.Ls > 1)
+    if (Ls_ > 1)
     {
-        env.addProp(getName() + "_5d", par_.Ls);
+        env.addProp(getName() + "_5d", Ls_);
         quark5d_ = env.getProp(getName() + "_5d");
     }
 }
@@ -72,22 +80,53 @@ void MQuark::allocate(Environment &env)
 // execution
 void MQuark::execute(Environment &env)
 {
+    LatticePropagator *fullSource;
+    LatticeFermion    source(env.getGrid(Ls_)), sol(env.getGrid(Ls_));
+    
     LOG(Message) << "computing quark propagator '" << getName() << "'"
                  << std::endl;
+    if (!env.isProp5d(par_.source))
+    {
+        if (Ls_ == 1)
+        {
+            fullSource = env.getProp(par_.source);
+        }
+        else
+        {
+            HADRON_ERROR("MQuark not implemented with 5D actions");
+        }
+    }
+    else
+    {
+        if (Ls_ == 1)
+        {
+            HADRON_ERROR("MQuark not implemented with 5D actions");
+        }
+        else if (Ls_ != env.getPropLs(par_.source))
+        {
+            HADRON_ERROR("MQuark not implemented with 5D actions");
+        }
+        else
+        {
+            fullSource = env.getProp(par_.source);
+        }
+    }
     
-    GridCartesian         *g4d   = env.get4dGrid(),
+    LOG(Message) << "inverting using solver '" << par_.solver
-                          *g5d   = env.get5dGrid(par_.Ls);
+                 << "' on source '" << par_.source << "'" << std::endl;
-    GridRedBlackCartesian *gRb4d = env.getRb4dGrid(),
+    for (unsigned int s = 0; s < Ns; ++s)
-                          *gRb5d = env.getRb5dGrid(par_.Ls);
+    for (unsigned int c = 0; c < Nc; ++c)
-    LatticeGaugeField     &Umu   = *env.getGauge();
+    {
-    LatticeFermion        src(g5d); src=zero;
+        PropToFerm(source, *fullSource, s, c);
-    LatticeFermion        result(g5d); result=zero;
+        sol = zero;
-    
+        env.callSolver(par_.solver, sol, source);
-    RealD mass=0.1;
+        if (Ls_ == 1)
-    RealD M5=1.8;
+        {
-    DomainWallFermionR Ddwf(Umu, *g5d, *gRb5d, *g4d, *gRb4d, mass, M5);
+            FermToProp(*quark_, sol, s, c);
-    
+        }
-    ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
+        else
-    SchurRedBlackDiagMooeeSolve<LatticeFermion> SchurSolver(CG);
+        {
-    SchurSolver(Ddwf,src,result);
+            HADRON_ERROR("MQuark not implemented with 5D actions");
+        }
+    }
 }

programs/Hadrons/MQuark.hpp

@@ -43,7 +43,8 @@ public:
     class Par: Serializable
     {
     public:
-        GRID_SERIALIZABLE_CLASS_MEMBERS(Par, unsigned int, Ls);
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Par, std::string , source,
+                                             std::string , solver);
     };
 public:
     // constructor
@@ -55,13 +56,17 @@ public:
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(Environment &env);
     // allocation
     virtual void allocate(Environment &env);
     // execution
     virtual void execute(Environment &env);
 private:
-    Par par_;
+    Par                 par_;
-    LatticePropagator *quark_{nullptr}, *quark5d_{nullptr};
+    unsigned int        Ls_;
+    LatticePropagator   *source_{nullptr}, *quark_{nullptr}, *quark5d_{nullptr};
+    Environment::Solver *solver_{nullptr};
 };
 
 MODULE_REGISTER(MQuark);

programs/Hadrons/MSource.cc

@@ -1,9 +1,141 @@
-//
+/*******************************************************************************
-//  MSource.cpp
+Grid physics library, www.github.com/paboyle/Grid 
-//  Grid
-//
-//  Created by Antonin Portelli on 16/04/2016.
-//
-//
 
-#include "MSource.hpp"
+Source file: programs/Hadrons/MSource.cc
+
+Copyright (C) 2016
+
+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.
+*******************************************************************************/
+
+#include <Hadrons/MSource.hpp>
+
+#define ERROR_SUF " (source '" << getName() << "')"
+using namespace Grid;
+using namespace Hadrons;
+
+/******************************************************************************
+*                         MSource implementation                              *
+******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+MSource::MSource(const std::string name)
+: Module(name)
+{}
+
+// parse parameters
+void MSource::parseParameters(XmlReader &reader, const std::string name)
+{
+    read(reader, name, par_);
+}
+
+// dependency relation
+std::vector<std::string> MSource::getInput(void)
+{
+    return std::vector<std::string>();
+}
+
+std::vector<std::string> MSource::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// allocation //////////////////////////////////////////////////////////////////
+void MSource::allocate(Environment &env)
+{
+    switch (par_.sourceType)
+    {
+        // 4D sources
+        case Grid::SourceType::point:
+        case Grid::SourceType::z2Band:
+            env.addProp(getName());
+            src_ = env.getProp(getName());
+            break;
+        // error
+        default:
+            HADRON_ERROR("no allocation implemented for source type '"
+                         << par_.sourceType << "'" << ERROR_SUF);
+            break;
+    }
+}
+
+// execution
+#define ARG_CHECK(n)\
+if (par_.arguments.size() != (n))\
+{\
+    HADRON_ERROR("source type '" << par_.sourceType << "' expect "\
+                 << (n) << " arguments (got "\
+                 << par_.arguments.size() << ")" << ERROR_SUF);\
+}
+
+void MSource::execute(Environment &env)
+{
+    LOG(Message) << "generating source '" << getName() << "' of type '"
+                 << par_.sourceType << "'" << std::endl;
+    switch (par_.sourceType)
+    {
+        // point source
+        case Grid::SourceType::point:
+        {
+            ARG_CHECK(1);
+            
+            std::vector<int> origin = strToVec<int>(par_.arguments[0]);
+            SpinColourMatrix id(1.);
+            
+            if (origin.size() != Nd)
+            {
+                HADRON_ERROR("point source origin dimension different from "
+                             << Nd << ERROR_SUF);
+            }
+            *src_ = zero;
+            pokeSite(id, *src_, origin);
+            
+            break;
+        }
+        // z2Band source
+        case Grid::SourceType::z2Band:
+        {
+            ARG_CHECK(2);
+            
+            int                        ta = std::stoi(par_.arguments[0]);
+            int                        tb = std::stoi(par_.arguments[1]);
+            Lattice<iScalar<vInteger>> t(env.getGrid());
+            LatticeComplex             eta(env.getGrid());
+            LatticeFermion             phi(env.getGrid());
+            ComplexD                   shift(1., 1.);
+            
+            LatticeCoordinate(t, Tp);
+            bernoulli(*env.get4dRng(), eta);
+            eta   = (2.*eta - shift)*(1./::sqrt(2.));
+            *src_ = 1.;
+            *src_ = where((t >= ta) and (t <= tb), (*src_)*eta, (*src_)*0.);
+            
+            break;
+        }
+        // error
+        default:
+        {
+            HADRON_ERROR("no definition implemented for source type '"
+                         << par_.sourceType << "'" << ERROR_SUF);
+            break;
+        }
+    }
+}

programs/Hadrons/MSource.hpp

@@ -1,14 +1,103 @@
-//
+/*******************************************************************************
-//  MSource.hpp
+Grid physics library, www.github.com/paboyle/Grid 
-//  Grid
-//
-//  Created by Antonin Portelli on 16/04/2016.
-//
-//
 
-#ifndef MSource_hpp
+Source file: programs/Hadrons/MSource.hpp
-#define MSource_hpp
 
-#include <stdio.h>
+Copyright (C) 2016
 
-#endif /* MSource_hpp */
+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.
+*******************************************************************************/
+
+/************
+ * Sources  *
+ ************
+ 
+ Description of all source types.
+ Convention: the discrete Heavyside function verifies theta(0) = 1.
+ 
+ point: Point source
+ -------------------
+ * src(x) = delta_x,o
+ 
+ * arguments: o
+    - o: origin, space-separated integer sequence (e.g. "0 1 1 0")
+ 
+ z2Band: Z_2 stochastic source
+ -----------------------------
+ * src(x) = eta_x * theta(x_0 - ta) * theta(tb - x_0)
+ 
+ * arguments: ta tb
+    - ta: begin timeslice (integer)
+    - tb: end timesilce (integer)
+ 
+ */
+
+
+#ifndef Hadrons_MSource_hpp_
+#define Hadrons_MSource_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+
+namespace Grid{
+    GRID_SERIALIZABLE_ENUM(SourceType, undef,
+                                       point,  1,
+                                       z2Band, 2);
+}
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                            Source module                                   *
+ ******************************************************************************/
+class MSource: public Module
+{
+public:
+    class Par: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Par, SourceType, sourceType,
+                                        std::vector<std::string>, arguments);
+    };
+public:
+    // constructor
+    MSource(const std::string name);
+    // destructor
+    virtual ~MSource(void) = default;
+    // parse parameters
+    virtual void parseParameters(XmlReader &reader, const std::string name);
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // allocation
+    virtual void allocate(Environment &env);
+    // execution
+    virtual void execute(Environment &env);
+private:
+    Par               par_;
+    LatticePropagator *src_{nullptr};
+};
+
+MODULE_REGISTER(MSource);
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MSource_hpp_

programs/Hadrons/Makefile.am

@@ -1,19 +1,31 @@
 AM_CXXFLAGS = -I$(top_srcdir)/programs -I$(top_srcdir)/lib
-AM_LDFLAGS = -L$(top_builddir)/lib
+AM_LDFLAGS  = -L$(top_builddir)/lib
 
 bin_PROGRAMS = Hadrons
 
+## general sources
 Hadrons_SOURCES =   \
     Application.cc  \
     Environment.cc  \
+    FermionAction.cc\
     Global.cc       \
     Hadrons.cc      \
-    Module.cc       \
+    Module.cc
-    ModuleFactory.cc
 
+## fermion actions
 Hadrons_SOURCES +=  \
-    MQuark.cc
+    AWilson.cc
 
+## general modules
+Hadrons_SOURCES +=  \
+    MQuark.cc       \
+    MSource.cc
+
+## solver modules
+Hadrons_SOURCES +=  \
+    SRBPrecCG.cc
+
+## contraction modules
 Hadrons_SOURCES +=  \
     CMeson.cc
 

programs/Hadrons/Module.cc

@@ -46,6 +46,7 @@ std::string Module::getName(void) const
 
 void Module::operator()(Environment &env)
 {
+    setup(env);
     allocate(env);
     if (!env.isDryRun())
     {

programs/Hadrons/Module.hpp

@@ -35,19 +35,19 @@ BEGIN_HADRONS_NAMESPACE
 
 // module registration macro
 #define MODULE_REGISTER(mod)\
-class mod##Registrar\
+class mod##ModuleRegistrar\
 {\
 public:\
-    mod##Registrar(void)\
+    mod##ModuleRegistrar(void)\
     {\
         ModuleFactory &modFac = ModuleFactory::getInstance();\
-        modFac.registerModule(#mod, [&](const std::string name)\
+        modFac.registerBuilder(#mod, [&](const std::string name)\
                               {\
                                   return std::unique_ptr<mod>(new mod(name));\
                               });\
     }\
 };\
-static mod##Registrar mod##RegistrarInstance;
+static mod##ModuleRegistrar mod##ModuleRegistrarInstance;
 
 /******************************************************************************
  *                                 Module                                     *
@@ -66,8 +66,10 @@ public:
     // dependency relation
     virtual std::vector<std::string> getInput(void) = 0;
     virtual std::vector<std::string> getOutput(void) = 0;
+    // setup
+    virtual void setup(Environment &env) {};
     // allocation
-    virtual void allocate(Environment &env) = 0;
+    virtual void allocate(Environment &env) {};
     // execution
     void operator()(Environment &env);
     virtual void execute(Environment &env) = 0;

programs/Hadrons/ModuleFactory.hpp

@@ -29,6 +29,7 @@ directory.
 #define Hadrons_ModuleFactory_hpp_
 
 #include <Hadrons/Global.hpp>
+#include <Hadrons/Factory.hpp>
 #include <Hadrons/Module.hpp>
 
 BEGIN_HADRONS_NAMESPACE
@@ -36,22 +37,9 @@ BEGIN_HADRONS_NAMESPACE
 /******************************************************************************
  *                            ModuleFactory                                   *
  ******************************************************************************/
-class ModuleFactory
+class ModuleFactory: public Factory<Module>
 {
-    SINGLETON(ModuleFactory)
+    SINGLETON_DEFCTOR(ModuleFactory)
-public:
-    typedef std::function<std::unique_ptr<Module>(const std::string )>
-        FactoryFunc;
-public:
-    // registration
-    void registerModule(const std::string type, const FactoryFunc &f);
-    // get module list
-    std::vector<std::string> getModuleList(void) const;
-    // factory
-    std::unique_ptr<Module> create(const std::string type,
-                                   const std::string name) const;
-private:
-    std::map<std::string, FactoryFunc> factory_;
 };
 
 END_HADRONS_NAMESPACE

programs/Hadrons/SRBPrecCG.cc

@@ -0,0 +1,74 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/SRBPrecCG.cc
+
+Copyright (C) 2016
+
+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.
+*******************************************************************************/
+
+#include <Hadrons/SRBPrecCG.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+
+/******************************************************************************
+*                       SRBPrecCG implementation                              *
+******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+SRBPrecCG::SRBPrecCG(const std::string name)
+: Module(name)
+{}
+
+// parse parameters
+void SRBPrecCG::parseParameters(XmlReader &reader, const std::string name)
+{
+   read(reader, name, par_);
+}
+
+// dependency relation
+std::vector<std::string> SRBPrecCG::getInput(void)
+{
+    return std::vector<std::string>();
+}
+
+std::vector<std::string> SRBPrecCG::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+void SRBPrecCG::execute(Environment &env)
+{
+    auto &mat   = *(env.getFermionAction(par_.action)->getFMat());
+    auto solver = [&mat, this](LatticeFermion &sol,
+                               const LatticeFermion &source)
+    {
+        ConjugateGradient<LatticeFermion>           cg(par_.residual, 10000);
+        SchurRedBlackDiagMooeeSolve<LatticeFermion> schurSolver(cg);
+        
+        schurSolver(mat, source, sol);
+    };
+    
+    env.addSolver(getName(), solver, par_.action);
+}

programs/Hadrons/SRBPrecCG.hpp

@@ -0,0 +1,69 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/SRBPrecCG.hpp
+
+Copyright (C) 2016
+
+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.
+*******************************************************************************/
+
+#ifndef Hadrons_SRBPrecCG_hpp_
+#define Hadrons_SRBPrecCG_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         SRBPrecCG                                 *
+ ******************************************************************************/
+class SRBPrecCG: public Module
+{
+public:
+    class Par: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Par, std::string, action,
+                                             double     , residual);
+    };
+public:
+    // constructor
+    SRBPrecCG(const std::string name);
+    // destructor
+    virtual ~SRBPrecCG(void) = default;
+    // parse parameters
+    virtual void parseParameters(XmlReader &reader, const std::string name);
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // execution
+    virtual void execute(Environment &env);
+private:
+    Par par_;
+};
+
+MODULE_REGISTER(SRBPrecCG);
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_SRBPrecCG_hpp_