Merge branch 'develop' of https://github.com/paboyle/Grid into develop

Makefile.am

@@ -1,5 +1,5 @@
 # additional include paths necessary to compile the C++ library
-SUBDIRS = lib benchmarks tests
+SUBDIRS = lib benchmarks tests extras
 
 include $(top_srcdir)/doxygen.inc
 

configure.ac

@@ -381,10 +381,13 @@ AC_CONFIG_FILES(tests/IO/Makefile)
 AC_CONFIG_FILES(tests/core/Makefile)
 AC_CONFIG_FILES(tests/debug/Makefile)
 AC_CONFIG_FILES(tests/forces/Makefile)
+AC_CONFIG_FILES(tests/hadrons/Makefile)
 AC_CONFIG_FILES(tests/hmc/Makefile)
 AC_CONFIG_FILES(tests/solver/Makefile)
 AC_CONFIG_FILES(tests/qdpxx/Makefile)
 AC_CONFIG_FILES(benchmarks/Makefile)
+AC_CONFIG_FILES(extras/Makefile)
+AC_CONFIG_FILES(extras/Hadrons/Makefile)
 AC_OUTPUT
 
 echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

extras/Hadrons/Application.cc

@@ -0,0 +1,312 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Application.cc
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#include <Grid/Hadrons/Application.hpp>
+#include <Grid/Hadrons/GeneticScheduler.hpp>
+
+using namespace Grid;
+using namespace QCD;
+using namespace Hadrons;
+
+#define BIG_SEP "==============="
+#define SEP     "---------------"
+
+/******************************************************************************
+ *                       Application implementation                           *
+ ******************************************************************************/
+// constructors ////////////////////////////////////////////////////////////////
+Application::Application(void)
+: env_(Environment::getInstance())
+{
+    LOG(Message) << "Modules available:" << std::endl;
+    auto list = ModuleFactory::getInstance().getBuilderList();
+    for (auto &m: list)
+    {
+        LOG(Message) << "  " << m << std::endl;
+    }
+    auto dim = GridDefaultLatt(), mpi = GridDefaultMpi(), loc(dim);
+    locVol_ = 1;
+    for (unsigned int d = 0; d < dim.size(); ++d)
+    {
+        loc[d]  /= mpi[d];
+        locVol_ *= loc[d];
+    }
+    LOG(Message) << "Global lattice: " << dim << std::endl;
+    LOG(Message) << "MPI partition : " << mpi << std::endl;
+    LOG(Message) << "Local lattice : " << loc << std::endl;
+}
+
+Application::Application(const Application::GlobalPar &par)
+: Application()
+{
+    setPar(par);
+}
+
+Application::Application(const std::string parameterFileName)
+: Application()
+{
+    parameterFileName_ = parameterFileName;
+}
+
+// access //////////////////////////////////////////////////////////////////////
+void Application::setPar(const Application::GlobalPar &par)
+{
+    par_ = par;
+    env_.setSeed(strToVec<int>(par_.seed));
+}
+
+const Application::GlobalPar & Application::getPar(void)
+{
+    return par_;
+}
+
+// execute /////////////////////////////////////////////////////////////////////
+void Application::run(void)
+{
+    if (!parameterFileName_.empty() and (env_.getNModule() == 0))
+    {
+        parseParameterFile(parameterFileName_);
+    }
+    if (!scheduled_)
+    {
+        schedule();
+    }
+    printSchedule();
+    configLoop();
+}
+
+// parse parameter file ////////////////////////////////////////////////////////
+class ObjectId: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(ObjectId,
+                                    std::string, name,
+                                    std::string, type);
+};
+
+void Application::parseParameterFile(const std::string parameterFileName)
+{
+    XmlReader reader(parameterFileName);
+    GlobalPar par;
+    ObjectId  id;
+    
+    LOG(Message) << "Building application from '" << parameterFileName << "'..." << std::endl;
+    read(reader, "parameters", par);
+    setPar(par);
+    push(reader, "modules");
+    push(reader, "module");
+    do
+    {
+        read(reader, "id", id);
+        env_.createModule(id.name, id.type, reader);
+    } while (reader.nextElement("module"));
+    pop(reader);
+    pop(reader);
+}
+
+void Application::saveParameterFile(const std::string parameterFileName)
+{
+    XmlWriter          writer(parameterFileName);
+    ObjectId           id;
+    const unsigned int nMod = env_.getNModule();
+    
+    LOG(Message) << "Saving application to '" << parameterFileName << "'..." << std::endl;
+    write(writer, "parameters", getPar());
+    push(writer, "modules");
+    for (unsigned int i = 0; i < nMod; ++i)
+    {
+        push(writer, "module");
+        id.name = env_.getModuleName(i);
+        id.type = env_.getModule(i)->getRegisteredName();
+        write(writer, "id", id);
+        env_.getModule(i)->saveParameters(writer, "options");
+        pop(writer);
+    }
+    pop(writer);
+    pop(writer);
+}
+
+// schedule computation ////////////////////////////////////////////////////////
+#define MEM_MSG(size)\
+sizeString((size)*locVol_) << " (" << sizeString(size)  << "/site)"
+
+#define DEFINE_MEMPEAK \
+auto memPeak = [this](const std::vector<unsigned int> &program)\
+{\
+    unsigned int memPeak;\
+    bool         msg;\
+    \
+    msg = HadronsLogMessage.isActive();\
+    HadronsLogMessage.Active(false);\
+    env_.dryRun(true);\
+    memPeak = env_.executeProgram(program);\
+    env_.dryRun(false);\
+    env_.freeAll();\
+    HadronsLogMessage.Active(true);\
+    \
+    return memPeak;\
+}
+
+void Application::schedule(void)
+{
+    DEFINE_MEMPEAK;
+    
+    // build module dependency graph
+    LOG(Message) << "Building module graph..." << std::endl;
+    auto graph = env_.makeModuleGraph();
+    auto con = graph.getConnectedComponents();
+    
+    // constrained topological sort using a genetic algorithm
+    LOG(Message) << "Scheduling computation..." << std::endl;
+    LOG(Message) << "               #module= " << graph.size() << std::endl;
+    LOG(Message) << "       population size= " << par_.genetic.popSize << std::endl;
+    LOG(Message) << "       max. generation= " << par_.genetic.maxGen << std::endl;
+    LOG(Message) << "  max. cst. generation= " << par_.genetic.maxCstGen << std::endl;
+    LOG(Message) << "         mutation rate= " << par_.genetic.mutationRate << std::endl;
+    
+    unsigned int                               k = 0, gen, prevPeak, nCstPeak = 0;
+    std::random_device                         rd;
+    GeneticScheduler<unsigned int>::Parameters par;
+    
+    par.popSize      = par_.genetic.popSize;
+    par.mutationRate = par_.genetic.mutationRate;
+    par.seed         = rd();
+    memPeak_         = 0;
+    CartesianCommunicator::BroadcastWorld(0, &(par.seed), sizeof(par.seed));
+    for (unsigned int i = 0; i < con.size(); ++i)
+    {
+        GeneticScheduler<unsigned int> scheduler(con[i], memPeak, par);
+        
+        gen = 0;
+        do
+        {
+            LOG(Debug) << "Generation " << gen << ":" << std::endl;
+            scheduler.nextGeneration();
+            if (gen != 0)
+            {
+                if (prevPeak == scheduler.getMinValue())
+                {
+                    nCstPeak++;
+                }
+                else
+                {
+                    nCstPeak = 0;
+                }
+            }
+            
+            prevPeak = scheduler.getMinValue();
+            if (gen % 10 == 0)
+            {
+                LOG(Iterative) << "Generation " << gen << ": "
+                               << MEM_MSG(scheduler.getMinValue()) << std::endl;
+            }
+            
+            gen++;
+        } while ((gen < par_.genetic.maxGen)
+                 and (nCstPeak < par_.genetic.maxCstGen));
+        auto &t = scheduler.getMinSchedule();
+        if (scheduler.getMinValue() > memPeak_)
+        {
+            memPeak_ = scheduler.getMinValue();
+        }
+        for (unsigned int j = 0; j < t.size(); ++j)
+        {
+            program_.push_back(t[j]);
+        }
+    }
+    scheduled_ = true;
+}
+
+void Application::saveSchedule(const std::string filename)
+{
+    TextWriter               writer(filename);
+    std::vector<std::string> program;
+    
+    if (!scheduled_)
+    {
+        HADRON_ERROR("Computation not scheduled");
+    }
+    LOG(Message) << "Saving current schedule to '" << filename << "'..."
+                 << std::endl;
+    for (auto address: program_)
+    {
+        program.push_back(env_.getModuleName(address));
+    }
+    write(writer, "schedule", program);
+}
+
+void Application::loadSchedule(const std::string filename)
+{
+    DEFINE_MEMPEAK;
+    
+    TextReader               reader(filename);
+    std::vector<std::string> program;
+    
+    LOG(Message) << "Loading schedule from '" << filename << "'..."
+                 << std::endl;
+    read(reader, "schedule", program);
+    program_.clear();
+    for (auto &name: program)
+    {
+        program_.push_back(env_.getModuleAddress(name));
+    }
+    scheduled_ = true;
+    memPeak_   = memPeak(program_);
+}
+
+void Application::printSchedule(void)
+{
+    if (!scheduled_)
+    {
+        HADRON_ERROR("Computation not scheduled");
+    }
+    LOG(Message) << "Schedule (memory peak: " << MEM_MSG(memPeak_) << "):"
+                 << std::endl;
+    for (unsigned int i = 0; i < program_.size(); ++i)
+    {
+        LOG(Message) << std::setw(4) << i + 1 << ": "
+                     << env_.getModuleName(program_[i]) << std::endl;
+    }
+}
+
+// loop on configurations //////////////////////////////////////////////////////
+void Application::configLoop(void)
+{
+    auto range = par_.trajCounter;
+    
+    for (unsigned int t = range.start; t < range.end; t += range.step)
+    {
+        LOG(Message) << BIG_SEP << " Starting measurement for trajectory " << t
+                     << " " << BIG_SEP << std::endl;
+        env_.setTrajectory(t);
+        env_.executeProgram(program_);
+        env_.freeAll();
+    }
+    LOG(Message) << BIG_SEP << " End of measurement " << BIG_SEP << std::endl;
+}

extras/Hadrons/Application.hpp

@@ -0,0 +1,130 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Application.hpp
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_Application_hpp_
+#define Hadrons_Application_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Environment.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+#include <Grid/Hadrons/Modules.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         Main program manager                               *
+ ******************************************************************************/
+class Application
+{
+public:
+    class TrajRange: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(TrajRange,
+                                        unsigned int, start,
+                                        unsigned int, end,
+                                        unsigned int, step);
+    };
+    class GeneticPar: Serializable
+    {
+    public:
+        GeneticPar(void):
+            popSize{20}, maxGen{1000}, maxCstGen{100}, mutationRate{.1} {};
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(GeneticPar,
+                                        unsigned int, popSize,
+                                        unsigned int, maxGen,
+                                        unsigned int, maxCstGen,
+                                        double      , mutationRate);
+    };
+    class GlobalPar: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(GlobalPar,
+                                        TrajRange,   trajCounter,
+                                        GeneticPar,  genetic,
+                                        std::string, seed);
+    };
+public:
+    // constructors
+    Application(void);
+    Application(const GlobalPar &par);
+    Application(const std::string parameterFileName);
+    // destructor
+    virtual ~Application(void) = default;
+    // access
+    void              setPar(const GlobalPar &par);
+    const GlobalPar & getPar(void);
+    // module creation
+    template <typename M>
+    void createModule(const std::string name);
+    template <typename M>
+    void createModule(const std::string name, const typename M::Par &par);
+    // execute
+    void run(void);
+    // XML parameter file I/O
+    void parseParameterFile(const std::string parameterFileName);
+    void saveParameterFile(const std::string parameterFileName);
+    // schedule computation
+    void schedule(void);
+    void saveSchedule(const std::string filename);
+    void loadSchedule(const std::string filename);
+    void printSchedule(void);
+    // loop on configurations
+    void configLoop(void);
+private:
+    long unsigned int         locVol_;
+    std::string               parameterFileName_{""};
+    GlobalPar                 par_;
+    Environment               &env_;
+    std::vector<unsigned int> program_;
+    Environment::Size         memPeak_;
+    bool                      scheduled_{false};
+};
+
+/******************************************************************************
+ *                     Application template implementation                    *
+ ******************************************************************************/
+// module creation /////////////////////////////////////////////////////////////
+template <typename M>
+void Application::createModule(const std::string name)
+{
+    env_.createModule<M>(name);
+}
+
+template <typename M>
+void Application::createModule(const std::string name,
+                               const typename M::Par &par)
+{
+    env_.createModule<M>(name, par);
+}
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Application_hpp_

extras/Hadrons/Environment.cc

@@ -0,0 +1,702 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Environment.cc
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#include <Grid/Hadrons/Environment.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+using namespace Grid;
+using namespace QCD;
+using namespace Hadrons;
+
+/******************************************************************************
+ *                       Environment implementation                           *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+Environment::Environment(void)
+{
+    grid4d_.reset(SpaceTimeGrid::makeFourDimGrid(
+        GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()),
+        GridDefaultMpi()));
+    gridRb4d_.reset(SpaceTimeGrid::makeFourDimRedBlackGrid(grid4d_.get()));
+    auto loc = getGrid()->LocalDimensions();
+    locVol_ = 1;
+    for (unsigned int d = 0; d < loc.size(); ++d)
+    {
+        locVol_ *= loc[d];
+    }
+    rng4d_.reset(new GridParallelRNG(grid4d_.get()));
+}
+
+// dry run /////////////////////////////////////////////////////////////////////
+void Environment::dryRun(const bool isDry)
+{
+    dryRun_ = isDry;
+}
+
+bool Environment::isDryRun(void) const
+{
+    return dryRun_;
+}
+
+// trajectory number ///////////////////////////////////////////////////////////
+void Environment::setTrajectory(const unsigned int traj)
+{
+    traj_ = traj;
+}
+
+unsigned int Environment::getTrajectory(void) const
+{
+    return traj_;
+}
+
+// grids ///////////////////////////////////////////////////////////////////////
+void Environment::createGrid(const unsigned int Ls)
+{
+    if (grid5d_.find(Ls) == grid5d_.end())
+    {
+        auto g = getGrid();
+        
+        grid5d_[Ls].reset(SpaceTimeGrid::makeFiveDimGrid(Ls, g));
+        gridRb5d_[Ls].reset(SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, g));
+    }
+}
+
+GridCartesian * Environment::getGrid(const unsigned int Ls) const
+{
+    try
+    {
+        if (Ls == 1)
+        {
+            return grid4d_.get();
+        }
+        else
+        {
+            return grid5d_.at(Ls).get();
+        }
+    }
+    catch(std::out_of_range &)
+    {
+        HADRON_ERROR("no grid with Ls= " << Ls);
+    }
+}
+
+GridRedBlackCartesian * Environment::getRbGrid(const unsigned int Ls) const
+{
+    try
+    {
+        if (Ls == 1)
+        {
+            return gridRb4d_.get();
+        }
+        else
+        {
+            return gridRb5d_.at(Ls).get();
+        }
+    }
+    catch(std::out_of_range &)
+    {
+        HADRON_ERROR("no red-black 5D grid with Ls= " << Ls);
+    }
+}
+
+// random number generator /////////////////////////////////////////////////////
+void Environment::setSeed(const std::vector<int> &seed)
+{
+    rng4d_->SeedFixedIntegers(seed);
+}
+
+GridParallelRNG * Environment::get4dRng(void) const
+{
+    return rng4d_.get();
+}
+
+// module management ///////////////////////////////////////////////////////////
+void Environment::pushModule(Environment::ModPt &pt)
+{
+    std::string name = pt->getName();
+    
+    if (!hasModule(name))
+    {
+        std::vector<unsigned int> inputAddress;
+        unsigned int              address;
+        ModuleInfo                m;
+        
+        m.data = std::move(pt);
+        m.type = typeIdPt(*m.data.get());
+        m.name = name;
+        auto input  = m.data->getInput();
+        for (auto &in: input)
+        {
+            if (!hasObject(in))
+            {
+                addObject(in , -1);
+            }
+            m.input.push_back(objectAddress_[in]);
+        }
+        auto output = m.data->getOutput();
+        module_.push_back(std::move(m));
+        address              = static_cast<unsigned int>(module_.size() - 1);
+        moduleAddress_[name] = address;
+        for (auto &out: output)
+        {
+            if (!hasObject(out))
+            {
+                addObject(out, address);
+            }
+            else
+            {
+                if (object_[objectAddress_[out]].module < 0)
+                {
+                    object_[objectAddress_[out]].module = address;
+                }
+                else
+                {
+                    HADRON_ERROR("object '" + out
+                                 + "' is already produced by module '"
+                                 + module_[object_[getObjectAddress(out)].module].name
+                                 + "' (while pushing module '" + name + "')");
+                }
+            }
+        }
+    }
+    else
+    {
+        HADRON_ERROR("module '" + name + "' already exists");
+    }
+}
+
+unsigned int Environment::getNModule(void) const
+{
+    return module_.size();
+}
+
+void Environment::createModule(const std::string name, const std::string type,
+                               XmlReader &reader)
+{
+    auto &factory = ModuleFactory::getInstance();
+    auto pt       = factory.create(type, name);
+    
+    pt->parseParameters(reader, "options");
+    pushModule(pt);
+}
+
+ModuleBase * Environment::getModule(const unsigned int address) const
+{
+    if (hasModule(address))
+    {
+        return module_[address].data.get();
+    }
+    else
+    {
+        HADRON_ERROR("no module with address " + std::to_string(address));
+    }
+}
+
+ModuleBase * Environment::getModule(const std::string name) const
+{
+    return getModule(getModuleAddress(name));
+}
+
+unsigned int Environment::getModuleAddress(const std::string name) const
+{
+    if (hasModule(name))
+    {
+        return moduleAddress_.at(name);
+    }
+    else
+    {
+        HADRON_ERROR("no module with name '" + name + "'");
+    }
+}
+
+std::string Environment::getModuleName(const unsigned int address) const
+{
+    if (hasModule(address))
+    {
+        return module_[address].name;
+    }
+    else
+    {
+        HADRON_ERROR("no module with address " + std::to_string(address));
+    }
+}
+
+std::string Environment::getModuleType(const unsigned int address) const
+{
+    if (hasModule(address))
+    {
+        return typeName(module_[address].type);
+    }
+    else
+    {
+        HADRON_ERROR("no module with address " + std::to_string(address));
+    }
+}
+
+std::string Environment::getModuleType(const std::string name) const
+{
+    return getModuleType(getModuleAddress(name));
+}
+
+bool Environment::hasModule(const unsigned int address) const
+{
+    return (address < module_.size());
+}
+
+bool Environment::hasModule(const std::string name) const
+{
+    return (moduleAddress_.find(name) != moduleAddress_.end());
+}
+
+Graph<unsigned int> Environment::makeModuleGraph(void) const
+{
+    Graph<unsigned int> moduleGraph;
+    
+    for (unsigned int i = 0; i < module_.size(); ++i)
+    {
+        moduleGraph.addVertex(i);
+        for (auto &j: module_[i].input)
+        {
+            moduleGraph.addEdge(object_[j].module, i);
+        }
+    }
+    
+    return moduleGraph;
+}
+
+#define BIG_SEP "==============="
+#define SEP     "---------------"
+#define MEM_MSG(size)\
+sizeString((size)*locVol_) << " (" << sizeString(size)  << "/site)"
+
+Environment::Size
+Environment::executeProgram(const std::vector<unsigned int> &p)
+{
+    Size                                memPeak = 0, sizeBefore, sizeAfter;
+    std::vector<std::set<unsigned int>> freeProg;
+    bool                                continueCollect, nothingFreed;
+    
+    // build garbage collection schedule
+    freeProg.resize(p.size());
+    for (unsigned int i = 0; i < object_.size(); ++i)
+    {
+        auto pred = [i, this](const unsigned int j)
+        {
+            auto &in = module_[j].input;
+            auto it  = std::find(in.begin(), in.end(), i);
+            
+            return (it != in.end()) or (j == object_[i].module);
+        };
+        auto it = std::find_if(p.rbegin(), p.rend(), pred);
+        if (it != p.rend())
+        {
+            freeProg[p.rend() - it - 1].insert(i);
+        }
+    }
+    
+    // program execution
+    for (unsigned int i = 0; i < p.size(); ++i)
+    {
+        // execute module
+        if (!isDryRun())
+        {
+            LOG(Message) << SEP << " Measurement step " << i+1 << "/"
+                         << p.size() << " (module '" << module_[p[i]].name
+                         << "') " << SEP << std::endl;
+        }
+        (*module_[p[i]].data)();
+        sizeBefore = getTotalSize();
+        // print used memory after execution
+        if (!isDryRun())
+        {
+            LOG(Message) << "Allocated objects: " << MEM_MSG(sizeBefore)
+                         << std::endl;
+        }
+        if (sizeBefore > memPeak)
+        {
+            memPeak = sizeBefore;
+        }
+        // garbage collection for step i
+        if (!isDryRun())
+        {
+            LOG(Message) << "Garbage collection..." << std::endl;
+        }
+        nothingFreed = true;
+        do
+        {
+            continueCollect = false;
+            auto toFree = freeProg[i];
+            for (auto &j: toFree)
+            {
+                // continue garbage collection while there are still
+                // objects without owners
+                continueCollect = continueCollect or !hasOwners(j);
+                if(freeObject(j))
+                {
+                    // if an object has been freed, remove it from
+                    // the garbage collection schedule
+                    freeProg[i].erase(j);
+                    nothingFreed = false;
+                }
+            }
+        } while (continueCollect);
+        // any remaining objects in step i garbage collection schedule
+        // is scheduled for step i + 1
+        if (i + 1 < p.size())
+        {
+            for (auto &j: freeProg[i])
+            {
+                freeProg[i + 1].insert(j);
+            }
+        }
+        // print used memory after garbage collection if necessary
+        if (!isDryRun())
+        {
+            sizeAfter = getTotalSize();
+            if (sizeBefore != sizeAfter)
+            {
+                LOG(Message) << "Allocated objects: " << MEM_MSG(sizeAfter)
+                             << std::endl;
+            }
+            else
+            {
+                LOG(Message) << "Nothing to free" << std::endl;
+            }
+        }
+    }
+    
+    return memPeak;
+}
+
+Environment::Size Environment::executeProgram(const std::vector<std::string> &p)
+{
+    std::vector<unsigned int> pAddress;
+    
+    for (auto &n: p)
+    {
+        pAddress.push_back(getModuleAddress(n));
+    }
+    
+    return executeProgram(pAddress);
+}
+
+// general memory management ///////////////////////////////////////////////////
+void Environment::addObject(const std::string name, const int moduleAddress)
+{
+    ObjInfo info;
+    
+    info.name   = name;
+    info.module = moduleAddress;
+    object_.push_back(std::move(info));
+    objectAddress_[name] = static_cast<unsigned int>(object_.size() - 1);
+}
+
+void Environment::registerObject(const unsigned int address,
+                                 const unsigned int size, const unsigned int Ls)
+{
+    if (!hasRegisteredObject(address))
+    {
+        if (hasObject(address))
+        {
+            object_[address].size         = size;
+            object_[address].Ls           = Ls;
+            object_[address].isRegistered = true;
+        }
+        else
+        {
+            HADRON_ERROR("no object with address " + std::to_string(address));
+        }
+    }
+    else
+    {
+        HADRON_ERROR("object with address " + std::to_string(address)
+                     + " already registered");
+    }
+}
+
+void Environment::registerObject(const std::string name,
+                                 const unsigned int size, const unsigned int Ls)
+{
+    registerObject(getObjectAddress(name), size, Ls);
+}
+
+unsigned int Environment::getObjectAddress(const std::string name) const
+{
+    if (hasObject(name))
+    {
+        return objectAddress_.at(name);
+    }
+    else
+    {
+        HADRON_ERROR("no object with name '" + name + "'");
+    }
+}
+
+std::string Environment::getObjectName(const unsigned int address) const
+{
+    if (hasObject(address))
+    {
+        return object_[address].name;
+    }
+    else
+    {
+        HADRON_ERROR("no object with address " + std::to_string(address));
+    }
+}
+
+std::string Environment::getObjectType(const unsigned int address) const
+{
+    if (hasRegisteredObject(address))
+    {
+        return typeName(object_[address].type);
+    }
+    else if (hasObject(address))
+    {
+        HADRON_ERROR("object with address " + std::to_string(address)
+                     + " exists but is not registered");
+    }
+    else
+    {
+        HADRON_ERROR("no object with address " + std::to_string(address));
+    }
+}
+
+std::string Environment::getObjectType(const std::string name) const
+{
+    return getObjectType(getObjectAddress(name));
+}
+
+Environment::Size Environment::getObjectSize(const unsigned int address) const
+{
+    if (hasRegisteredObject(address))
+    {
+        return object_[address].size;
+    }
+    else if (hasObject(address))
+    {
+        HADRON_ERROR("object with address " + std::to_string(address)
+                     + " exists but is not registered");
+    }
+    else
+    {
+        HADRON_ERROR("no object with address " + std::to_string(address));
+    }
+}
+
+Environment::Size Environment::getObjectSize(const std::string name) const
+{
+    return getObjectSize(getObjectAddress(name));
+}
+
+unsigned int Environment::getObjectLs(const unsigned int address) const
+{
+    if (hasRegisteredObject(address))
+    {
+        return object_[address].Ls;
+    }
+    else if (hasObject(address))
+    {
+        HADRON_ERROR("object with address " + std::to_string(address)
+                     + " exists but is not registered");
+    }
+    else
+    {
+        HADRON_ERROR("no object with address " + std::to_string(address));
+    }
+}
+
+unsigned int Environment::getObjectLs(const std::string name) const
+{
+    return getObjectLs(getObjectAddress(name));
+}
+
+bool Environment::hasObject(const unsigned int address) const
+{
+    return (address < object_.size());
+}
+
+bool Environment::hasObject(const std::string name) const
+{
+    auto it = objectAddress_.find(name);
+    
+    return ((it != objectAddress_.end()) and hasObject(it->second));
+}
+
+bool Environment::hasRegisteredObject(const unsigned int address) const
+{
+    if (hasObject(address))
+    {
+        return object_[address].isRegistered;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+bool Environment::hasRegisteredObject(const std::string name) const
+{
+    if (hasObject(name))
+    {
+        return hasRegisteredObject(getObjectAddress(name));
+    }
+    else
+    {
+        return false;
+    }
+}
+
+bool Environment::isObject5d(const unsigned int address) const
+{
+    return (getObjectLs(address) > 1);
+}
+
+bool Environment::isObject5d(const std::string name) const
+{
+    return (getObjectLs(name) > 1);
+}
+
+Environment::Size Environment::getTotalSize(void) const
+{
+    Environment::Size size = 0;
+    
+    for (auto &o: object_)
+    {
+        if (o.isRegistered)
+        {
+            size += o.size;
+        }
+    }
+    
+    return size;
+}
+
+void Environment::addOwnership(const unsigned int owner,
+                               const unsigned int property)
+{
+    if (hasObject(property))
+    {
+        object_[property].owners.insert(owner);
+    }
+    else
+    {
+        HADRON_ERROR("no object with address " + std::to_string(property));
+    }
+    if (hasObject(owner))
+    {
+        object_[owner].properties.insert(property);
+    }
+    else
+    {
+        HADRON_ERROR("no object with address " + std::to_string(owner));
+    }
+}
+
+void Environment::addOwnership(const std::string owner,
+                               const std::string property)
+{
+    addOwnership(getObjectAddress(owner), getObjectAddress(property));
+}
+
+bool Environment::hasOwners(const unsigned int address) const
+{
+    
+    if (hasObject(address))
+    {
+        return (!object_[address].owners.empty());
+    }
+    else
+    {
+        HADRON_ERROR("no object with address " + std::to_string(address));
+    }
+}
+
+bool Environment::hasOwners(const std::string name) const
+{
+    return hasOwners(getObjectAddress(name));
+}
+
+bool Environment::freeObject(const unsigned int address)
+{
+    if (!hasOwners(address))
+    {
+        if (!isDryRun() and object_[address].isRegistered)
+        {
+            LOG(Message) << "Destroying object '" << object_[address].name
+                         << "'" << std::endl;
+        }
+        for (auto &p: object_[address].properties)
+        {
+            object_[p].owners.erase(address);
+        }
+        object_[address].size         = 0;
+        object_[address].Ls           = 0;
+        object_[address].isRegistered = false;
+        object_[address].type         = nullptr;
+        object_[address].owners.clear();
+        object_[address].properties.clear();
+        object_[address].data.reset(nullptr);
+        
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+bool Environment::freeObject(const std::string name)
+{
+    return freeObject(getObjectAddress(name));
+}
+
+void Environment::freeAll(void)
+{
+    for (unsigned int i = 0; i < object_.size(); ++i)
+    {
+        freeObject(i);
+    }
+}
+
+void Environment::printContent(void)
+{
+    LOG(Message) << "Modules: " << std::endl;
+    for (unsigned int i = 0; i < module_.size(); ++i)
+    {
+        LOG(Message) << std::setw(4) << i << ": "
+                     << getModuleName(i) << std::endl;
+    }
+    LOG(Message) << "Objects: " << std::endl;
+    for (unsigned int i = 0; i < object_.size(); ++i)
+    {
+        LOG(Message) << std::setw(4) << i << ": "
+                     << getObjectName(i) << std::endl;
+    }
+}

extras/Hadrons/Environment.hpp

@@ -0,0 +1,381 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Environment.hpp
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_Environment_hpp_
+#define Hadrons_Environment_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Graph.hpp>
+
+#ifndef SITE_SIZE_TYPE
+#define SITE_SIZE_TYPE unsigned int
+#endif
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         Global environment                                 *
+ ******************************************************************************/
+// forward declaration of Module
+class ModuleBase;
+
+class Object
+{
+public:
+    Object(void) = default;
+    virtual ~Object(void) = default;
+};
+
+template <typename T>
+class Holder: public Object
+{
+public:
+    Holder(void) = default;
+    Holder(T *pt);
+    virtual ~Holder(void) = default;
+    T &       get(void) const;
+    T *       getPt(void) const;
+    void      reset(T *pt);
+private:
+    std::unique_ptr<T> objPt_{nullptr};
+};
+
+class Environment
+{
+    SINGLETON(Environment);
+public:
+    typedef SITE_SIZE_TYPE                         Size;
+    typedef std::unique_ptr<ModuleBase>            ModPt;
+    typedef std::unique_ptr<GridCartesian>         GridPt;
+    typedef std::unique_ptr<GridRedBlackCartesian> GridRbPt;
+    typedef std::unique_ptr<GridParallelRNG>       RngPt;
+    typedef std::unique_ptr<LatticeBase>           LatticePt;
+private:
+    struct ModuleInfo
+    {
+        const std::type_info      *type{nullptr};
+        std::string               name;
+        ModPt                     data{nullptr};
+        std::vector<unsigned int> input;
+    };
+    struct ObjInfo
+    {
+        Size                    size{0};
+        unsigned int            Ls{0};
+        bool                    isRegistered{false};
+        const std::type_info    *type{nullptr};
+        std::string             name;
+        int                     module{-1};
+        std::set<unsigned int>  owners, properties;
+        std::unique_ptr<Object> data{nullptr};
+    };
+public:
+    // dry run
+    void                    dryRun(const bool isDry);
+    bool                    isDryRun(void) const;
+    // trajectory number
+    void                    setTrajectory(const unsigned int traj);
+    unsigned int            getTrajectory(void) const;
+    // grids
+    void                    createGrid(const unsigned int Ls);
+    GridCartesian *         getGrid(const unsigned int Ls = 1) const;
+    GridRedBlackCartesian * getRbGrid(const unsigned int Ls = 1) const;
+    // random number generator
+    void                    setSeed(const std::vector<int> &seed);
+    GridParallelRNG *       get4dRng(void) const;
+    // module management
+    void                    pushModule(ModPt &pt);
+    template <typename M>
+    void                    createModule(const std::string name);
+    template <typename M>
+    void                    createModule(const std::string name,
+                                         const typename M::Par &par);
+    void                    createModule(const std::string name,
+                                         const std::string type,
+                                         XmlReader &reader);
+    unsigned int            getNModule(void) const;
+    ModuleBase *            getModule(const unsigned int address) const;
+    ModuleBase *            getModule(const std::string name) const;
+    template <typename M>
+    M *                     getModule(const unsigned int address) const;
+    template <typename M>
+    M *                     getModule(const std::string name) const;
+    unsigned int            getModuleAddress(const std::string name) const;
+    std::string             getModuleName(const unsigned int address) const;
+    std::string             getModuleType(const unsigned int address) const;
+    std::string             getModuleType(const std::string name) const;
+    bool                    hasModule(const unsigned int address) const;
+    bool                    hasModule(const std::string name) const;
+    Graph<unsigned int>     makeModuleGraph(void) const;
+    Size                    executeProgram(const std::vector<unsigned int> &p);
+    Size                    executeProgram(const std::vector<std::string> &p);
+    // general memory management
+    void                    addObject(const std::string name,
+                                      const int moduleAddress);
+    void                    registerObject(const unsigned int address,
+                                           const unsigned int size,
+                                           const unsigned int Ls = 1);
+    void                    registerObject(const std::string name,
+                                           const unsigned int size,
+                                           const unsigned int Ls = 1);
+    template <typename T>
+    unsigned int            lattice4dSize(void) const;
+    template <typename T>
+    void                    registerLattice(const unsigned int address,
+                                            const unsigned int Ls = 1);
+    template <typename T>
+    void                    registerLattice(const std::string name,
+                                            const unsigned int Ls = 1);
+    template <typename T>
+    void                    setObject(const unsigned int address, T *object);
+    template <typename T>
+    void                    setObject(const std::string name, T *object);
+    template <typename T>
+    T *                     getObject(const unsigned int address) const;
+    template <typename T>
+    T *                     getObject(const std::string name) const;
+    template <typename T>
+    T *                     createLattice(const unsigned int address);
+    template <typename T>
+    T *                     createLattice(const std::string name);
+    unsigned int            getObjectAddress(const std::string name) const;
+    std::string             getObjectName(const unsigned int address) const;
+    std::string             getObjectType(const unsigned int address) const;
+    std::string             getObjectType(const std::string name) const;
+    Size                    getObjectSize(const unsigned int address) const;
+    Size                    getObjectSize(const std::string name) const;
+    unsigned int            getObjectLs(const unsigned int address) const;
+    unsigned int            getObjectLs(const std::string name) const;
+    bool                    hasObject(const unsigned int address) const;
+    bool                    hasObject(const std::string name) const;
+    bool                    hasRegisteredObject(const unsigned int address) const;
+    bool                    hasRegisteredObject(const std::string name) const;
+    bool                    isObject5d(const unsigned int address) const;
+    bool                    isObject5d(const std::string name) const;
+    Environment::Size       getTotalSize(void) const;
+    void                    addOwnership(const unsigned int owner,
+                                         const unsigned int property);
+    void                    addOwnership(const std::string owner,
+                                         const std::string property);
+    bool                    hasOwners(const unsigned int address) const;
+    bool                    hasOwners(const std::string name) const;
+    bool                    freeObject(const unsigned int address);
+    bool                    freeObject(const std::string name);
+    void                    freeAll(void);
+    void                    printContent(void);
+private:
+    // general
+    bool                                   dryRun_{false};
+    unsigned int                           traj_, locVol_;
+    // grids
+    GridPt                                 grid4d_;
+    std::map<unsigned int, GridPt>         grid5d_;
+    GridRbPt                               gridRb4d_;
+    std::map<unsigned int, GridRbPt>       gridRb5d_;
+    // random number generator
+    RngPt                                  rng4d_;
+    // module and related maps
+    std::vector<ModuleInfo>                module_;
+    std::map<std::string, unsigned int>    moduleAddress_;
+    // lattice store
+    std::map<unsigned int, LatticePt>      lattice_;
+    // object store
+    std::vector<ObjInfo>                   object_;
+    std::map<std::string, unsigned int>    objectAddress_;
+};
+
+/******************************************************************************
+ *                       Holder template implementation                       *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename T>
+Holder<T>::Holder(T *pt)
+: objPt_(pt)
+{}
+
+// access //////////////////////////////////////////////////////////////////////
+template <typename T>
+T & Holder<T>::get(void) const
+{
+    return &objPt_.get();
+}
+
+template <typename T>
+T * Holder<T>::getPt(void) const
+{
+    return objPt_.get();
+}
+
+template <typename T>
+void Holder<T>::reset(T *pt)
+{
+    objPt_.reset(pt);
+}
+
+/******************************************************************************
+ *                     Environment template implementation                    *
+ ******************************************************************************/
+// module management ///////////////////////////////////////////////////////////
+template <typename M>
+void Environment::createModule(const std::string name)
+{
+    ModPt pt(new M(name));
+    
+    pushModule(pt);
+}
+
+template <typename M>
+void Environment::createModule(const std::string name,
+                               const typename M::Par &par)
+{
+    ModPt pt(new M(name));
+    
+    static_cast<M *>(pt.get())->setPar(par);
+    pushModule(pt);
+}
+
+template <typename M>
+M * Environment::getModule(const unsigned int address) const
+{
+    if (auto *pt = dynamic_cast<M *>(getModule(address)))
+    {
+        return pt;
+    }
+    else
+    {
+        HADRON_ERROR("module '" + module_[address].name
+                     + "' does not have type " + typeid(M).name()
+                     + "(object type: " + getModuleType(address) + ")");
+    }
+}
+
+template <typename M>
+M * Environment::getModule(const std::string name) const
+{
+    return getModule<M>(getModuleAddress(name));
+}
+
+template <typename T>
+unsigned int Environment::lattice4dSize(void) const
+{
+    return sizeof(typename T::vector_object)/getGrid()->Nsimd();
+}
+
+template <typename T>
+void Environment::registerLattice(const unsigned int address,
+                                  const unsigned int Ls)
+{
+    createGrid(Ls);
+    registerObject(address, Ls*lattice4dSize<T>(), Ls);
+}
+
+template <typename T>
+void Environment::registerLattice(const std::string name, const unsigned int Ls)
+{
+    createGrid(Ls);
+    registerObject(name, Ls*lattice4dSize<T>(), Ls);
+}
+
+template <typename T>
+void Environment::setObject(const unsigned int address, T *object)
+{
+    if (hasRegisteredObject(address))
+    {
+        object_[address].data.reset(new Holder<T>(object));
+        object_[address].type = &typeid(T);
+    }
+    else if (hasObject(address))
+    {
+        HADRON_ERROR("object with address " + std::to_string(address) +
+                     " exists but is not registered");
+    }
+    else
+    {
+        HADRON_ERROR("no object with address " + std::to_string(address));
+    }
+}
+
+template <typename T>
+void Environment::setObject(const std::string name, T *object)
+{
+    setObject(getObjectAddress(name), object);
+}
+
+template <typename T>
+T * Environment::getObject(const unsigned int address) const
+{
+    if (hasRegisteredObject(address))
+    {
+        if (auto h = dynamic_cast<Holder<T> *>(object_[address].data.get()))
+        {
+            return h->getPt();
+        }
+        else
+        {
+            HADRON_ERROR("object with address " + std::to_string(address) +
+                         " does not have type '" + typeid(T).name() +
+                         "' (has type '" + getObjectType(address) + "')");
+        }
+    }
+    else if (hasObject(address))
+    {
+        HADRON_ERROR("object with address " + std::to_string(address) +
+                     " exists but is not registered");
+    }
+    else
+    {
+        HADRON_ERROR("no object with address " + std::to_string(address));
+    }
+}
+
+template <typename T>
+T * Environment::getObject(const std::string name) const
+{
+    return getObject<T>(getObjectAddress(name));
+}
+
+template <typename T>
+T * Environment::createLattice(const unsigned int address)
+{
+    GridCartesian *g = getGrid(getObjectLs(address));
+    
+    setObject(address, new T(g));
+    
+    return getObject<T>(address);
+}
+
+template <typename T>
+T * Environment::createLattice(const std::string name)
+{
+    return createLattice<T>(getObjectAddress(name));
+}
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Environment_hpp_

extras/Hadrons/Factory.hpp

@@ -0,0 +1,106 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Factory.hpp
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_Factory_hpp_
+#define Hadrons_Factory_hpp_
+
+#include <Grid/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 builder 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_

extras/Hadrons/GeneticScheduler.hpp

@@ -0,0 +1,329 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/GeneticScheduler.hpp
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_GeneticScheduler_hpp_
+#define Hadrons_GeneticScheduler_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Graph.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                   Scheduler based on a genetic algorithm                   *
+ ******************************************************************************/
+template <typename T>
+class GeneticScheduler
+{
+public:
+    typedef std::vector<T>                   Gene;
+    typedef std::pair<Gene *, Gene *>        GenePair;
+    typedef std::function<int(const Gene &)> ObjFunc;
+    struct Parameters
+    {
+        double       mutationRate;
+        unsigned int popSize, seed;
+    };
+public:
+    // constructor
+    GeneticScheduler(Graph<T> &graph, const ObjFunc &func,
+                     const Parameters &par);
+    // destructor
+    virtual ~GeneticScheduler(void) = default;
+    // access
+    const Gene & getMinSchedule(void);
+    int          getMinValue(void);
+    // breed a new generation
+    void nextGeneration(void);
+    // heuristic benchmarks
+    void benchmarkCrossover(const unsigned int nIt);
+    // print population
+    friend std::ostream & operator<<(std::ostream &out,
+                                     const GeneticScheduler<T> &s)
+    {
+        out << "[";
+        for (auto &p: s.population_)
+        {
+            out << p.first << ", ";
+        }
+        out << "\b\b]";
+        
+        return out;
+    }
+private:
+    // evolution steps
+    void initPopulation(void);
+    void doCrossover(void);
+    void doMutation(void);
+    // genetic operators
+    GenePair selectPair(void);
+    void     crossover(Gene &c1, Gene &c2, const Gene &p1, const Gene &p2);
+    void     mutation(Gene &m, const Gene &c);
+    
+private:
+    Graph<T>                 &graph_;
+    const ObjFunc            &func_;
+    const Parameters         par_;
+    std::multimap<int, Gene> population_;
+    std::mt19937             gen_;
+};
+
+/******************************************************************************
+ *                       template implementation                              *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename T>
+GeneticScheduler<T>::GeneticScheduler(Graph<T> &graph, const ObjFunc &func,
+                                      const Parameters &par)
+: graph_(graph)
+, func_(func)
+, par_(par)
+{
+    gen_.seed(par_.seed);
+}
+
+// access //////////////////////////////////////////////////////////////////////
+template <typename T>
+const typename GeneticScheduler<T>::Gene &
+GeneticScheduler<T>::getMinSchedule(void)
+{
+    return population_.begin()->second;
+}
+
+template <typename T>
+int GeneticScheduler<T>::getMinValue(void)
+{
+    return population_.begin()->first;
+}
+
+// breed a new generation //////////////////////////////////////////////////////
+template <typename T>
+void GeneticScheduler<T>::nextGeneration(void)
+{
+    // random initialization of the population if necessary
+    if (population_.size() != par_.popSize)
+    {
+        initPopulation();
+    }
+    LOG(Debug) << "Starting population:\n" << *this << std::endl;
+    
+    // random mutations
+    PARALLEL_FOR_LOOP
+    for (unsigned int i = 0; i < par_.popSize; ++i)
+    {
+        doMutation();
+    }
+    LOG(Debug) << "After mutations:\n" << *this << std::endl;
+    
+    // mating
+    PARALLEL_FOR_LOOP
+    for (unsigned int i = 0; i < par_.popSize/2; ++i)
+    {
+        doCrossover();
+    }
+    LOG(Debug) << "After mating:\n" << *this << std::endl;
+    
+    // grim reaper
+    auto it = population_.begin();
+    
+    std::advance(it, par_.popSize);
+    population_.erase(it, population_.end());
+    LOG(Debug) << "After grim reaper:\n" << *this << std::endl;
+}
+
+// evolution steps /////////////////////////////////////////////////////////////
+template <typename T>
+void GeneticScheduler<T>::initPopulation(void)
+{
+    population_.clear();
+    for (unsigned int i = 0; i < par_.popSize; ++i)
+    {
+        auto p = graph_.topoSort(gen_);
+        
+        population_.emplace(func_(p), p);
+    }
+}
+
+template <typename T>
+void GeneticScheduler<T>::doCrossover(void)
+{
+    auto p = selectPair();
+    Gene &p1 = *(p.first), &p2 = *(p.second);
+    Gene c1, c2;
+    
+    crossover(c1, c2, p1, p2);
+    PARALLEL_CRITICAL
+    {
+        population_.emplace(func_(c1), c1);
+        population_.emplace(func_(c2), c2);
+    }
+}
+
+template <typename T>
+void GeneticScheduler<T>::doMutation(void)
+{
+    std::uniform_real_distribution<double>      mdis(0., 1.);
+    std::uniform_int_distribution<unsigned int> pdis(0, population_.size() - 1);
+    
+    if (mdis(gen_) < par_.mutationRate)
+    {
+        Gene m;
+        auto it = population_.begin();
+        
+        std::advance(it, pdis(gen_));
+        mutation(m, it->second);
+        PARALLEL_CRITICAL
+        {
+            population_.emplace(func_(m), m);
+        }
+    }
+}
+
+// genetic operators ///////////////////////////////////////////////////////////
+template <typename T>
+typename GeneticScheduler<T>::GenePair GeneticScheduler<T>::selectPair(void)
+{
+    std::vector<double> prob;
+    unsigned int        ind;
+    Gene                *p1, *p2;
+    
+    for (auto &c: population_)
+    {
+        prob.push_back(1./c.first);
+    }
+    do
+    {
+        double probCpy;
+        
+        std::discrete_distribution<unsigned int> dis1(prob.begin(), prob.end());
+        auto rIt = population_.begin();
+        ind = dis1(gen_);
+        std::advance(rIt, ind);
+        p1 = &(rIt->second);
+        probCpy   = prob[ind];
+        prob[ind] = 0.;
+        std::discrete_distribution<unsigned int> dis2(prob.begin(), prob.end());
+        rIt = population_.begin();
+        std::advance(rIt, dis2(gen_));
+        p2 = &(rIt->second);
+        prob[ind] = probCpy;
+    } while (p1 == p2);
+    
+    return std::make_pair(p1, p2);
+}
+
+template <typename T>
+void GeneticScheduler<T>::crossover(Gene &c1, Gene &c2, const Gene &p1,
+                                    const Gene &p2)
+{
+    Gene                                        buf;
+    std::uniform_int_distribution<unsigned int> dis(0, p1.size() - 1);
+    unsigned int                                cut = dis(gen_);
+    
+    c1.clear();
+    buf = p2;
+    for (unsigned int i = 0; i < cut; ++i)
+    {
+        c1.push_back(p1[i]);
+        buf.erase(std::find(buf.begin(), buf.end(), p1[i]));
+    }
+    for (unsigned int i = 0; i < buf.size(); ++i)
+    {
+        c1.push_back(buf[i]);
+    }
+    c2.clear();
+    buf = p2;
+    for (unsigned int i = cut; i < p1.size(); ++i)
+    {
+        buf.erase(std::find(buf.begin(), buf.end(), p1[i]));
+    }
+    for (unsigned int i = 0; i < buf.size(); ++i)
+    {
+        c2.push_back(buf[i]);
+    }
+    for (unsigned int i = cut; i < p1.size(); ++i)
+    {
+        c2.push_back(p1[i]);
+    }
+}
+
+template <typename T>
+void GeneticScheduler<T>::mutation(Gene &m, const Gene &c)
+{
+    Gene                                        buf;
+    std::uniform_int_distribution<unsigned int> dis(0, c.size() - 1);
+    unsigned int                                cut = dis(gen_);
+    Graph<T>                                    g1 = graph_, g2 = graph_;
+    
+    for (unsigned int i = 0; i < cut; ++i)
+    {
+        g1.removeVertex(c[i]);
+    }
+    for (unsigned int i = cut; i < c.size(); ++i)
+    {
+        g2.removeVertex(c[i]);
+    }
+    if (g1.size() > 0)
+    {
+        buf = g1.topoSort(gen_);
+    }
+    if (g2.size() > 0)
+    {
+        m = g2.topoSort(gen_);
+    }
+    for (unsigned int i = cut; i < c.size(); ++i)
+    {
+        m.push_back(buf[i - cut]);
+    }
+}
+
+template <typename T>
+void GeneticScheduler<T>::benchmarkCrossover(const unsigned int nIt)
+{
+    Gene   p1, p2, c1, c2;
+    double neg = 0., eq = 0., pos = 0., total;
+    int    improvement;
+    
+    LOG(Message) << "Benchmarking crossover..." << std::endl;
+    for (unsigned int i = 0; i < nIt; ++i)
+    {
+        p1 = graph_.topoSort(gen_);
+        p2 = graph_.topoSort(gen_);
+        crossover(c1, c2, p1, p2);
+        improvement = (func_(c1) + func_(c2) - func_(p1) - func_(p2))/2;
+        if (improvement < 0) neg++; else if (improvement == 0) eq++; else pos++;
+    }
+    total = neg + eq + pos;
+    LOG(Message) << "  -: " << neg/total << "  =: " << eq/total
+                 << "  +: " << pos/total << std::endl;
+}
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_GeneticScheduler_hpp_

extras/Hadrons/Global.cc

@@ -0,0 +1,82 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Global.cc
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#include <Grid/Hadrons/Global.hpp>
+
+using namespace Grid;
+using namespace QCD;
+using namespace Hadrons;
+
+HadronsLogger Hadrons::HadronsLogError(1,"Error");
+HadronsLogger Hadrons::HadronsLogWarning(1,"Warning");
+HadronsLogger Hadrons::HadronsLogMessage(1,"Message");
+HadronsLogger Hadrons::HadronsLogIterative(1,"Iterative");
+HadronsLogger Hadrons::HadronsLogDebug(1,"Debug");
+
+// pretty size formatting //////////////////////////////////////////////////////
+std::string Hadrons::sizeString(long unsigned int bytes)
+
+{
+    constexpr unsigned int bufSize = 256;
+    const char             *suffixes[7] = {"", "K", "M", "G", "T", "P", "E"};
+    char                   buf[256];
+    long unsigned int      s     = 0;
+    double                 count = bytes;
+    
+    while (count >= 1024 && s < 7)
+    {
+        s++;
+        count /= 1024;
+    }
+    if (count - floor(count) == 0.0)
+    {
+        snprintf(buf, bufSize, "%d %sB", (int)count, suffixes[s]);
+    }
+    else
+    {
+        snprintf(buf, bufSize, "%.1f %sB", count, suffixes[s]);
+    }
+    
+    return std::string(buf);
+}
+
+// type utilities //////////////////////////////////////////////////////////////
+constexpr unsigned int maxNameSize = 1024u;
+
+std::string Hadrons::typeName(const std::type_info *info)
+{
+    char        *buf;
+    std::string name;
+    
+    buf  = abi::__cxa_demangle(info->name(), nullptr, nullptr, nullptr);
+    name = buf;
+    free(buf);
+    
+    return name;
+}

extras/Hadrons/Global.hpp

@@ -0,0 +1,150 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Global.hpp
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_Global_hpp_
+#define Hadrons_Global_hpp_
+
+#include <set>
+#include <stack>
+#include <Grid/Grid.h>
+#include <cxxabi.h>
+
+#define BEGIN_HADRONS_NAMESPACE \
+namespace Grid {\
+using namespace QCD;\
+namespace Hadrons {\
+using Grid::operator<<;
+#define END_HADRONS_NAMESPACE }}
+
+#define BEGIN_MODULE_NAMESPACE(name)\
+namespace name {\
+using Grid::operator<<;
+#define END_MODULE_NAMESPACE }
+
+/* the 'using Grid::operator<<;' statement prevents a very nasty compilation
+ * error with GCC 5 (clang & GCC 6 compile fine without it).
+ */
+
+// FIXME: find a way to do that in a more general fashion
+#ifndef FIMPL
+#define FIMPL WilsonImplR
+#endif
+
+BEGIN_HADRONS_NAMESPACE
+
+// type aliases
+#define TYPE_ALIASES(FImpl, suffix)\
+typedef FermionOperator<FImpl>                       FMat##suffix;             \
+typedef typename FImpl::FermionField                 FermionField##suffix;     \
+typedef typename FImpl::PropagatorField              PropagatorField##suffix;  \
+typedef typename FImpl::SitePropagator               SitePropagator##suffix;   \
+typedef typename FImpl::DoubledGaugeField            DoubledGaugeField##suffix;\
+typedef std::function<void(FermionField##suffix &,                             \
+                      const FermionField##suffix &)> SolverFn##suffix;
+
+// logger
+class HadronsLogger: public Logger
+{
+public:
+    HadronsLogger(int on, std::string nm): Logger("Hadrons", on, nm,
+                                                  GridLogColours, "BLACK"){};
+};
+
+#define LOG(channel) std::cout << HadronsLog##channel
+#define HADRON_ERROR(msg)\
+LOG(Error) << msg << " (" << __FUNCTION__ << " at " << __FILE__ << ":"\
+           << __LINE__ << ")" << std::endl;\
+abort();
+
+#define DEBUG_VAR(var) LOG(Debug) << #var << "= " << (var) << std::endl;
+
+extern HadronsLogger HadronsLogError;
+extern HadronsLogger HadronsLogWarning;
+extern HadronsLogger HadronsLogMessage;
+extern HadronsLogger HadronsLogIterative;
+extern HadronsLogger HadronsLogDebug;
+
+// singleton pattern
+#define SINGLETON(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);
+
+#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;
+
+// pretty size formating
+std::string sizeString(long unsigned int bytes);
+
+// type utilities
+template <typename T>
+const std::type_info * typeIdPt(const T &x)
+{
+    return &typeid(x);
+}
+
+std::string typeName(const std::type_info *info);
+
+template <typename T>
+const std::type_info * typeIdPt(void)
+{
+    return &typeid(T);
+}
+
+template <typename T>
+std::string typeName(const T &x)
+{
+    return typeName(typeIdPt(x));
+}
+
+template <typename T>
+std::string typeName(void)
+{
+    return typeName(typeIdPt<T>());
+}
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Global_hpp_

extras/Hadrons/Graph.hpp

@@ -0,0 +1,760 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Graph.hpp
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_Graph_hpp_
+#define Hadrons_Graph_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                          Oriented graph class                              *
+ ******************************************************************************/
+// I/O for edges
+template <typename T>
+std::ostream & operator<<(std::ostream &out, const std::pair<T, T> &e)
+{
+    out << "\""  << e.first << "\" -> \"" << e.second << "\"";
+    
+    return out;
+}
+
+// main class
+template <typename T>
+class Graph
+{
+public:
+    typedef std::pair<T, T> Edge;
+public:
+    // constructor
+    Graph(void);
+    // destructor
+    virtual ~Graph(void) = default;
+    // access
+    void           addVertex(const T &value);
+    void           addEdge(const Edge &e);
+    void           addEdge(const T &start, const T &end);
+    std::vector<T> getVertices(void) const;
+    void           removeVertex(const T &value);
+    void           removeEdge(const Edge &e);
+    void           removeEdge(const T &start, const T &end);
+    unsigned int   size(void) const;
+    // tests
+    bool gotValue(const T &value) const;
+    // graph topological manipulations
+    std::vector<T>              getAdjacentVertices(const T &value) const;
+    std::vector<T>              getChildren(const T &value) const;
+    std::vector<T>              getParents(const T &value) const;
+    std::vector<T>              getRoots(void) const;
+    std::vector<Graph<T>>       getConnectedComponents(void) const;
+    std::vector<T>              topoSort(void);
+    template <typename Gen>
+    std::vector<T>              topoSort(Gen &gen);
+    std::vector<std::vector<T>> allTopoSort(void);
+    // I/O
+    friend std::ostream & operator<<(std::ostream &out, const Graph<T> &g)
+    {
+        out << "{";
+        for (auto &e: g.edgeSet_)
+        {
+            out << e << ", ";
+        }
+        if (g.edgeSet_.size() != 0)
+        {
+            out << "\b\b";
+        }
+        out << "}";
+        
+        return out;
+    }
+private:
+    // vertex marking
+    void      mark(const T &value, const bool doMark = true);
+    void      markAll(const bool doMark = true);
+    void      unmark(const T &value);
+    void      unmarkAll(void);
+    bool      isMarked(const T &value) const;
+    const T * getFirstMarked(const bool isMarked = true) const;
+    template <typename Gen>
+    const T * getRandomMarked(const bool isMarked, Gen &gen);
+    const T * getFirstUnmarked(void) const;
+    template <typename Gen>
+    const T * getRandomUnmarked(Gen &gen);
+    // prune marked/unmarked vertices
+    void removeMarked(const bool isMarked = true);
+    void removeUnmarked(void);
+    // depth-first search marking
+    void depthFirstSearch(void);
+    void depthFirstSearch(const T &root);
+private:
+    std::map<T, bool>  isMarked_;
+    std::set<Edge>     edgeSet_;
+};
+
+// build depedency matrix from topological sorts
+template <typename T>
+std::map<T, std::map<T, bool>>
+makeDependencyMatrix(const std::vector<std::vector<T>> &topSort);
+
+/******************************************************************************
+ *                       template implementation                              *
+ ******************************************************************************
+ * in all the following V is the number of vertex and E is the number of edge
+ * in the worst case E = V^2
+ */
+
+// constructor /////////////////////////////////////////////////////////////////
+template <typename T>
+Graph<T>::Graph(void)
+{}
+
+// access //////////////////////////////////////////////////////////////////////
+// complexity: log(V)
+template <typename T>
+void Graph<T>::addVertex(const T &value)
+{
+    isMarked_[value] = false;
+}
+
+// complexity: O(log(V))
+template <typename T>
+void Graph<T>::addEdge(const Edge &e)
+{
+    addVertex(e.first);
+    addVertex(e.second);
+    edgeSet_.insert(e);
+}
+
+// complexity: O(log(V))
+template <typename T>
+void Graph<T>::addEdge(const T &start, const T &end)
+{
+    addEdge(Edge(start, end));
+}
+
+template <typename T>
+std::vector<T> Graph<T>::getVertices(void) const
+{
+    std::vector<T> vertex;
+    
+    for (auto &v: isMarked_)
+    {
+        vertex.push_back(v.first);
+    }
+    
+    return vertex;
+}
+
+// complexity: O(V*log(V))
+template <typename T>
+void Graph<T>::removeVertex(const T &value)
+{
+    // remove vertex from the mark table
+    auto vIt = isMarked_.find(value);
+    
+    if (vIt != isMarked_.end())
+    {
+        isMarked_.erase(vIt);
+    }
+    else
+    {
+        HADRON_ERROR("vertex " << value << " does not exists");
+    }
+
+    // remove all edges containing the vertex
+    auto pred = [&value](const Edge &e)
+    {
+        return ((e.first == value) or (e.second == value));
+    };
+    auto eIt = find_if(edgeSet_.begin(), edgeSet_.end(), pred);
+    
+    while (eIt != edgeSet_.end())
+    {
+        edgeSet_.erase(eIt);
+        eIt = find_if(edgeSet_.begin(), edgeSet_.end(), pred);
+    }
+}
+
+// complexity: O(log(V))
+template <typename T>
+void Graph<T>::removeEdge(const Edge &e)
+{
+    auto eIt = edgeSet_.find(e);
+    
+    if (eIt != edgeSet_.end())
+    {
+        edgeSet_.erase(eIt);
+    }
+    else
+    {
+        HADRON_ERROR("edge "  << e << " does not exists");
+    }
+}
+
+// complexity: O(log(V))
+template <typename T>
+void Graph<T>::removeEdge(const T &start, const T &end)
+{
+    removeEdge(Edge(start, end));
+}
+
+// complexity: O(1)
+template <typename T>
+unsigned int Graph<T>::size(void) const
+{
+    return isMarked_.size();
+}
+
+// tests ///////////////////////////////////////////////////////////////////////
+// complexity: O(log(V))
+template <typename T>
+bool Graph<T>::gotValue(const T &value) const
+{
+    auto it = isMarked_.find(value);
+    
+    if (it == isMarked_.end())
+    {
+        return false;
+    }
+    else
+    {
+        return true;
+    }
+}
+
+// vertex marking //////////////////////////////////////////////////////////////
+// complexity: O(log(V))
+template <typename T>
+void Graph<T>::mark(const T &value, const bool doMark)
+{
+    if (gotValue(value))
+    {
+        isMarked_[value] = doMark;
+    }
+    else
+    {
+        HADRON_ERROR("vertex " << value << " does not exists");
+    }
+}
+
+// complexity: O(V*log(V))
+template <typename T>
+void Graph<T>::markAll(const bool doMark)
+{
+    for (auto &v: isMarked_)
+    {
+        mark(v.first, doMark);
+    }
+}
+
+// complexity: O(log(V))
+template <typename T>
+void Graph<T>::unmark(const T &value)
+{
+    mark(value, false);
+}
+
+// complexity: O(V*log(V))
+template <typename T>
+void Graph<T>::unmarkAll(void)
+{
+    markAll(false);
+}
+
+// complexity: O(log(V))
+template <typename T>
+bool Graph<T>::isMarked(const T &value) const
+{
+    if (gotValue(value))
+    {
+        return isMarked_.at(value);
+    }
+    else
+    {
+        HADRON_ERROR("vertex " << value << " does not exists");
+        
+        return false;
+    }
+}
+
+// complexity: O(log(V))
+template <typename T>
+const T * Graph<T>::getFirstMarked(const bool isMarked) const
+{
+    auto pred = [&isMarked](const std::pair<T, bool> &v)
+    {
+        return (v.second == isMarked);
+    };
+    auto vIt = std::find_if(isMarked_.begin(), isMarked_.end(), pred);
+    
+    if (vIt != isMarked_.end())
+    {
+        return &(vIt->first);
+    }
+    else
+    {
+        return nullptr;
+    }
+}
+
+// complexity: O(log(V))
+template <typename T>
+template <typename Gen>
+const T * Graph<T>::getRandomMarked(const bool isMarked, Gen &gen)
+{
+    auto pred = [&isMarked](const std::pair<T, bool> &v)
+    {
+        return (v.second == isMarked);
+    };
+    std::uniform_int_distribution<unsigned int> dis(0, size() - 1);
+    auto                                        rIt = isMarked_.begin();
+    
+    std::advance(rIt, dis(gen));
+    auto vIt = std::find_if(rIt, isMarked_.end(), pred);
+    if (vIt != isMarked_.end())
+    {
+        return &(vIt->first);
+    }
+    else
+    {
+        vIt = std::find_if(isMarked_.begin(), rIt, pred);
+        if (vIt != rIt)
+        {
+            return &(vIt->first);
+        }
+        else
+        {
+            return nullptr;
+        }
+    }
+}
+
+// complexity: O(log(V))
+template <typename T>
+const T * Graph<T>::getFirstUnmarked(void) const
+{
+    return getFirstMarked(false);
+}
+
+// complexity: O(log(V))
+template <typename T>
+template <typename Gen>
+const T * Graph<T>::getRandomUnmarked(Gen &gen)
+{
+    return getRandomMarked(false, gen);
+}
+
+// prune marked/unmarked vertices //////////////////////////////////////////////
+// complexity: O(V^2*log(V))
+template <typename T>
+void Graph<T>::removeMarked(const bool isMarked)
+{
+    auto isMarkedCopy = isMarked_;
+    
+    for (auto &v: isMarkedCopy)
+    {
+        if (v.second == isMarked)
+        {
+            removeVertex(v.first);
+        }
+    }
+}
+
+// complexity: O(V^2*log(V))
+template <typename T>
+void Graph<T>::removeUnmarked(void)
+{
+    removeMarked(false);
+}
+
+// depth-first search marking //////////////////////////////////////////////////
+// complexity: O(V*log(V))
+template <typename T>
+void Graph<T>::depthFirstSearch(void)
+{
+    depthFirstSearch(isMarked_.begin()->first);
+}
+
+// complexity: O(V*log(V))
+template <typename T>
+void Graph<T>::depthFirstSearch(const T &root)
+{
+    std::vector<T> adjacentVertex;
+    
+    mark(root);
+    adjacentVertex = getAdjacentVertices(root);
+    for (auto &v: adjacentVertex)
+    {
+        if (!isMarked(v))
+        {
+            depthFirstSearch(v);
+        }
+    }
+}
+
+// graph topological manipulations /////////////////////////////////////////////
+// complexity: O(V*log(V))
+template <typename T>
+std::vector<T> Graph<T>::getAdjacentVertices(const T &value) const
+{
+    std::vector<T> adjacentVertex;
+    
+    auto pred = [&value](const Edge &e)
+    {
+        return ((e.first == value) or (e.second == value));
+    };
+    auto eIt = find_if(edgeSet_.begin(), edgeSet_.end(), pred);
+    
+    while (eIt != edgeSet_.end())
+    {
+        if (eIt->first == value)
+        {
+            adjacentVertex.push_back((*eIt).second);
+        }
+        else if (eIt->second == value)
+        {
+            adjacentVertex.push_back((*eIt).first);
+        }
+        eIt = find_if(++eIt, edgeSet_.end(), pred);
+    }
+    
+    return adjacentVertex;
+}
+
+// complexity: O(V*log(V))
+template <typename T>
+std::vector<T> Graph<T>::getChildren(const T &value) const
+{
+    std::vector<T> child;
+    
+    auto pred = [&value](const Edge &e)
+    {
+        return (e.first == value);
+    };
+    auto eIt = find_if(edgeSet_.begin(), edgeSet_.end(), pred);
+    
+    while (eIt != edgeSet_.end())
+    {
+        child.push_back((*eIt).second);
+        eIt = find_if(++eIt, edgeSet_.end(), pred);
+    }
+    
+    return child;
+}
+
+// complexity: O(V*log(V))
+template <typename T>
+std::vector<T> Graph<T>::getParents(const T &value) const
+{
+    std::vector<T> parent;
+    
+    auto pred = [&value](const Edge &e)
+    {
+        return (e.second == value);
+    };
+    auto eIt = find_if(edgeSet_.begin(), edgeSet_.end(), pred);
+    
+    while (eIt != edgeSet_.end())
+    {
+        parent.push_back((*eIt).first);
+        eIt = find_if(++eIt, edgeSet_.end(), pred);
+    }
+    
+    return parent;
+}
+
+// complexity: O(V^2*log(V))
+template <typename T>
+std::vector<T> Graph<T>::getRoots(void) const
+{
+    std::vector<T> root;
+    
+    for (auto &v: isMarked_)
+    {
+        auto parent = getParents(v.first);
+        
+        if (parent.size() == 0)
+        {
+            root.push_back(v.first);
+        }
+    }
+    
+    return root;
+}
+
+// complexity: O(V^2*log(V))
+template <typename T>
+std::vector<Graph<T>> Graph<T>::getConnectedComponents(void) const
+{
+    std::vector<Graph<T>> res;
+    Graph<T>              copy(*this);
+    
+    while (copy.size() > 0)
+    {
+        copy.depthFirstSearch();
+        res.push_back(copy);
+        res.back().removeUnmarked();
+        res.back().unmarkAll();
+        copy.removeMarked();
+        copy.unmarkAll();
+    }
+    
+    return res;
+}
+
+// topological sort using a directed DFS algorithm
+// complexity: O(V*log(V))
+template <typename T>
+std::vector<T> Graph<T>::topoSort(void)
+{
+    std::stack<T>     buf;
+    std::vector<T>    res;
+    const T           *vPt;
+    std::map<T, bool> tmpMarked(isMarked_);
+
+    // visit function
+    std::function<void(const T &)> visit = [&](const T &v)
+    {
+        if (tmpMarked.at(v))
+        {
+            HADRON_ERROR("cannot topologically sort a cyclic graph");
+        }
+        if (!isMarked(v))
+        {
+            std::vector<T> child = getChildren(v);
+
+            tmpMarked[v] = true;
+            for (auto &c: child)
+            {
+                visit(c);
+            }
+            mark(v);
+            tmpMarked[v] = false;
+            buf.push(v);
+        }
+    };
+    
+    // reset temporary marks
+    for (auto &v: tmpMarked)
+    {
+        tmpMarked.at(v.first) = false;
+    }
+    
+    // loop on unmarked vertices
+    unmarkAll();
+    vPt = getFirstUnmarked();
+    while (vPt)
+    {
+        visit(*vPt);
+        vPt = getFirstUnmarked();
+    }
+    unmarkAll();
+    
+    // create result vector
+    while (!buf.empty())
+    {
+        res.push_back(buf.top());
+        buf.pop();
+    }
+    
+    return res;
+}
+
+// random version of the topological sort
+// complexity: O(V*log(V))
+template <typename T>
+template <typename Gen>
+std::vector<T> Graph<T>::topoSort(Gen &gen)
+{
+    std::stack<T>     buf;
+    std::vector<T>    res;
+    const T           *vPt;
+    std::map<T, bool> tmpMarked(isMarked_);
+    
+    // visit function
+    std::function<void(const T &)> visit = [&](const T &v)
+    {
+        if (tmpMarked.at(v))
+        {
+            HADRON_ERROR("cannot topologically sort a cyclic graph");
+        }
+        if (!isMarked(v))
+        {
+            std::vector<T> child = getChildren(v);
+            
+            tmpMarked[v] = true;
+            std::shuffle(child.begin(), child.end(), gen);
+            for (auto &c: child)
+            {
+                visit(c);
+            }
+            mark(v);
+            tmpMarked[v] = false;
+            buf.push(v);
+        }
+    };
+    
+    // reset temporary marks
+    for (auto &v: tmpMarked)
+    {
+        tmpMarked.at(v.first) = false;
+    }
+    
+    // loop on unmarked vertices
+    unmarkAll();
+    vPt = getRandomUnmarked(gen);
+    while (vPt)
+    {
+        visit(*vPt);
+        vPt = getRandomUnmarked(gen);
+    }
+    unmarkAll();
+    
+    // create result vector
+    while (!buf.empty())
+    {
+        res.push_back(buf.top());
+        buf.pop();
+    }
+    
+    return res;
+}
+
+// generate all possible topological sorts
+// Y. L. Varol & D. Rotem, Comput. J. 24(1), pp. 83–84, 1981
+// http://comjnl.oupjournals.org/cgi/doi/10.1093/comjnl/24.1.83
+// complexity: O(V*log(V)) (from the paper, but really ?)
+template <typename T>
+std::vector<std::vector<T>> Graph<T>::allTopoSort(void)
+{
+    std::vector<std::vector<T>>    res;
+    std::map<T, std::map<T, bool>> iMat;
+    
+    // create incidence matrix
+    for (auto &v1: isMarked_)
+    for (auto &v2: isMarked_)
+    {
+        iMat[v1.first][v2.first] = false;
+    }
+    for (auto &v: isMarked_)
+    {
+        auto cVec = getChildren(v.first);
+        
+        for (auto &c: cVec)
+        {
+            iMat[v.first][c] = true;
+        }
+    }
+    
+    // generate initial topological sort
+    res.push_back(topoSort());
+    
+    // generate all other topological sorts by permutation
+    std::vector<T>            p = res[0];
+    const unsigned int        n = size();
+    std::vector<unsigned int> loc(n);
+    unsigned int              i, k, k1;
+    T                         obj_k, obj_k1;
+    bool                      isFinal;
+    
+    for (unsigned int j = 0; j < n; ++j)
+    {
+        loc[j] = j;
+    }
+    i = 0;
+    while (i < n-1)
+    {
+        k      = loc[i];
+        k1     = k + 1;
+        obj_k  = p[k];
+        if (k1 >= n)
+        {
+            isFinal = true;
+            obj_k1  = obj_k;
+        }
+        else
+        {
+            isFinal = false;
+            obj_k1  = p[k1];
+        }
+        if (iMat[res[0][i]][obj_k1] or isFinal)
+        {
+            for (unsigned int l = k; l >= i + 1; --l)
+            {
+                p[l]   = p[l-1];
+            }
+            p[i]   = obj_k;
+            loc[i] = i;
+            i++;
+        }
+        else
+        {
+            p[k]   = obj_k1;
+            p[k1]  = obj_k;
+            loc[i] = k1;
+            i      = 0;
+            res.push_back(p);
+        }
+    }
+    
+    return res;
+}
+
+// build depedency matrix from topological sorts ///////////////////////////////
+// complexity: something like O(V^2*log(V!))
+template <typename T>
+std::map<T, std::map<T, bool>>
+makeDependencyMatrix(const std::vector<std::vector<T>> &topSort)
+{
+    std::map<T, std::map<T, bool>> m;
+    const std::vector<T>           &vList = topSort[0];
+    
+    for (auto &v1: vList)
+    for (auto &v2: vList)
+    {
+        bool dep = true;
+        
+        for (auto &t: topSort)
+        {
+            auto i1 = std::find(t.begin(), t.end(), v1);
+            auto i2 = std::find(t.begin(), t.end(), v2);
+            
+            dep = dep and (i1 - i2 > 0);
+            if (!dep) break;
+        }
+        m[v1][v2] = dep;
+    }
+    
+    return m;
+}
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Graph_hpp_

extras/Hadrons/HadronsXmlRun.cc

@@ -0,0 +1,80 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/HadronsXmlRun.cc
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#include <Grid/Hadrons/Application.hpp>
+
+using namespace Grid;
+using namespace QCD;
+using namespace Hadrons;
+
+int main(int argc, char *argv[])
+{
+    // parse command line
+    std::string parameterFileName, scheduleFileName = "";
+    
+    if (argc < 2)
+    {
+        std::cerr << "usage: " << argv[0] << " <parameter file> [<precomputed schedule>] [Grid options]";
+        std::cerr << std::endl;
+        std::exit(EXIT_FAILURE);
+    }
+    parameterFileName = argv[1];
+    if (argc > 2)
+    {
+        if (argv[2][0] != '-')
+        {
+            scheduleFileName = argv[2];
+        }
+    }
+    
+    // initialization
+    Grid_init(&argc, &argv);
+    HadronsLogError.Active(GridLogError.isActive());
+    HadronsLogWarning.Active(GridLogWarning.isActive());
+    HadronsLogMessage.Active(GridLogMessage.isActive());
+    HadronsLogIterative.Active(GridLogIterative.isActive());
+    HadronsLogDebug.Active(GridLogDebug.isActive());
+    LOG(Message) << "Grid initialized" << std::endl;
+    
+    // execution
+    Application application(parameterFileName);
+    
+    application.parseParameterFile(parameterFileName);
+    if (!scheduleFileName.empty())
+    {
+        application.loadSchedule(scheduleFileName);
+    }
+    application.run();
+    
+    // epilogue
+    LOG(Message) << "Grid is finalizing now" << std::endl;
+    Grid_finalize();
+    
+    return EXIT_SUCCESS;
+}

extras/Hadrons/HadronsXmlSchedule.cc

@@ -0,0 +1,72 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/HadronsXmlSchedule.cc
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#include <Grid/Hadrons/Application.hpp>
+
+using namespace Grid;
+using namespace QCD;
+using namespace Hadrons;
+
+int main(int argc, char *argv[])
+{
+    // parse command line
+    std::string parameterFileName, scheduleFileName;
+    
+    if (argc < 3)
+    {
+        std::cerr << "usage: " << argv[0] << " <parameter file> <schedule output> [Grid options]";
+        std::cerr << std::endl;
+        std::exit(EXIT_FAILURE);
+    }
+    parameterFileName = argv[1];
+    scheduleFileName  = argv[2];
+    
+    // initialization
+    Grid_init(&argc, &argv);
+    HadronsLogError.Active(GridLogError.isActive());
+    HadronsLogWarning.Active(GridLogWarning.isActive());
+    HadronsLogMessage.Active(GridLogMessage.isActive());
+    HadronsLogIterative.Active(GridLogIterative.isActive());
+    HadronsLogDebug.Active(GridLogDebug.isActive());
+    LOG(Message) << "Grid initialized" << std::endl;
+    
+    // execution
+    Application application;
+    
+    application.parseParameterFile(parameterFileName);
+    application.schedule();
+    application.printSchedule();
+    application.saveSchedule(scheduleFileName);
+    
+    // epilogue
+    LOG(Message) << "Grid is finalizing now" << std::endl;
+    Grid_finalize();
+    
+    return EXIT_SUCCESS;
+}

extras/Hadrons/Makefile.am

@@ -0,0 +1,29 @@
+lib_LIBRARIES = libHadrons.a
+bin_PROGRAMS  = HadronsXmlRun HadronsXmlSchedule
+
+include modules.inc
+
+libHadrons_a_SOURCES = \
+    $(modules_cc)      \
+    Application.cc     \
+    Environment.cc     \
+    Global.cc          \
+    Module.cc
+libHadrons_adir = $(pkgincludedir)/Hadrons
+nobase_libHadrons_a_HEADERS = \
+	$(modules_hpp)            \
+	Application.hpp           \
+	Environment.hpp           \
+	Factory.hpp               \
+	GeneticScheduler.hpp      \
+	Global.hpp                \
+	Graph.hpp                 \
+	Module.hpp                \
+	Modules.hpp               \
+	ModuleFactory.hpp
+
+HadronsXmlRun_SOURCES = HadronsXmlRun.cc
+HadronsXmlRun_LDADD   = libHadrons.a -lGrid
+
+HadronsXmlSchedule_SOURCES = HadronsXmlSchedule.cc
+HadronsXmlSchedule_LDADD   = libHadrons.a -lGrid

extras/Hadrons/Module.cc

@@ -0,0 +1,71 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Module.cc
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#include <Grid/Hadrons/Module.hpp>
+
+using namespace Grid;
+using namespace QCD;
+using namespace Hadrons;
+
+/******************************************************************************
+ *                       ModuleBase implementation                            *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+ModuleBase::ModuleBase(const std::string name)
+: name_(name)
+, env_(Environment::getInstance())
+{}
+
+// access //////////////////////////////////////////////////////////////////////
+std::string ModuleBase::getName(void) const
+{
+    return name_;
+}
+
+Environment & ModuleBase::env(void) const
+{
+    return env_;
+}
+
+// get factory registration name if available
+std::string ModuleBase::getRegisteredName(void)
+{
+    HADRON_ERROR("module '" + getName() + "' has a type not registered"
+                 + " in the factory");
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+void ModuleBase::operator()(void)
+{
+    setup();
+    if (!env().isDryRun())
+    {
+        execute();
+    }
+}

extras/Hadrons/Module.hpp

@@ -0,0 +1,198 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Module.hpp
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_Module_hpp_
+#define Hadrons_Module_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Environment.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+// module registration macros
+#define MODULE_REGISTER(mod, base)\
+class mod: public base\
+{\
+public:\
+    typedef base Base;\
+    using Base::Base;\
+    virtual std::string getRegisteredName(void)\
+    {\
+        return std::string(#mod);\
+    }\
+};\
+class mod##ModuleRegistrar\
+{\
+public:\
+    mod##ModuleRegistrar(void)\
+    {\
+        ModuleFactory &modFac = ModuleFactory::getInstance();\
+        modFac.registerBuilder(#mod, [&](const std::string name)\
+                              {\
+                                  return std::unique_ptr<mod>(new mod(name));\
+                              });\
+    }\
+};\
+static mod##ModuleRegistrar mod##ModuleRegistrarInstance;
+
+#define MODULE_REGISTER_NS(mod, base, ns)\
+class mod: public base\
+{\
+public:\
+    typedef base Base;\
+    using Base::Base;\
+    virtual std::string getRegisteredName(void)\
+    {\
+        return std::string(#ns "::" #mod);\
+    }\
+};\
+class ns##mod##ModuleRegistrar\
+{\
+public:\
+    ns##mod##ModuleRegistrar(void)\
+    {\
+        ModuleFactory &modFac = ModuleFactory::getInstance();\
+        modFac.registerBuilder(#ns "::" #mod, [&](const std::string name)\
+                              {\
+                                  return std::unique_ptr<ns::mod>(new ns::mod(name));\
+                              });\
+    }\
+};\
+static ns##mod##ModuleRegistrar ns##mod##ModuleRegistrarInstance;
+
+#define ARG(...) __VA_ARGS__
+
+/******************************************************************************
+ *                            Module class                                    *
+ ******************************************************************************/
+// base class
+class ModuleBase
+{
+public:
+    // constructor
+    ModuleBase(const std::string name);
+    // destructor
+    virtual ~ModuleBase(void) = default;
+    // access
+    std::string getName(void) const;
+    Environment &env(void) const;
+    // get factory registration name if available
+    virtual std::string getRegisteredName(void);
+    // dependencies/products
+    virtual std::vector<std::string> getInput(void) = 0;
+    virtual std::vector<std::string> getOutput(void) = 0;
+    // parse parameters
+    virtual void parseParameters(XmlReader &reader, const std::string name) = 0;
+    virtual void saveParameters(XmlWriter &writer, const std::string name) = 0;
+    // setup
+    virtual void setup(void) {};
+    // execution
+    void operator()(void);
+    virtual void execute(void) = 0;
+private:
+    std::string name_;
+    Environment &env_;
+};
+
+// derived class, templating the parameter class
+template <typename P>
+class Module: public ModuleBase
+{
+public:
+    typedef P Par;
+public:
+    // constructor
+    Module(const std::string name);
+    // destructor
+    virtual ~Module(void) = default;
+    // parse parameters
+    virtual void parseParameters(XmlReader &reader, const std::string name);
+    virtual void saveParameters(XmlWriter &writer, const std::string name);
+    // parameter access
+    const P & par(void) const;
+    void      setPar(const P &par);
+private:
+    P par_;
+};
+
+// no parameter type
+class NoPar {};
+
+template <>
+class Module<NoPar>: public ModuleBase
+{
+public:
+    // constructor
+    Module(const std::string name): ModuleBase(name) {};
+    // destructor
+    virtual ~Module(void) = default;
+    // parse parameters (do nothing)
+    virtual void parseParameters(XmlReader &reader, const std::string name) {};
+    virtual void saveParameters(XmlWriter &writer, const std::string name)
+    {
+        push(writer, "options");
+        pop(writer);
+    };
+};
+
+/******************************************************************************
+ *                           Template implementation                          *
+ ******************************************************************************/
+template <typename P>
+Module<P>::Module(const std::string name)
+: ModuleBase(name)
+{}
+
+template <typename P>
+void Module<P>::parseParameters(XmlReader &reader, const std::string name)
+{
+    read(reader, name, par_);
+}
+
+template <typename P>
+void Module<P>::saveParameters(XmlWriter &writer, const std::string name)
+{
+    write(writer, name, par_);
+}
+
+template <typename P>
+const P & Module<P>::par(void) const
+{
+    return par_;
+}
+
+template <typename P>
+void Module<P>::setPar(const P &par)
+{
+    par_ = par;
+}
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Module_hpp_

extras/Hadrons/ModuleFactory.hpp

@@ -0,0 +1,49 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/ModuleFactory.hpp
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_ModuleFactory_hpp_
+#define Hadrons_ModuleFactory_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Factory.hpp>
+#include <Grid/Hadrons/Module.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                            ModuleFactory                                   *
+ ******************************************************************************/
+class ModuleFactory: public Factory<ModuleBase>
+{
+    SINGLETON_DEFCTOR(ModuleFactory)
+};
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_ModuleFactory_hpp_

extras/Hadrons/Modules.hpp

@@ -0,0 +1,40 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules.hpp
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MAction/DWF.hpp>
+#include <Grid/Hadrons/Modules/MAction/Wilson.hpp>
+#include <Grid/Hadrons/Modules/MContraction/Baryon.hpp>
+#include <Grid/Hadrons/Modules/MContraction/Meson.hpp>
+#include <Grid/Hadrons/Modules/MGauge/Load.hpp>
+#include <Grid/Hadrons/Modules/MGauge/Random.hpp>
+#include <Grid/Hadrons/Modules/MGauge/Unit.hpp>
+#include <Grid/Hadrons/Modules/MSolver/RBPrecCG.hpp>
+#include <Grid/Hadrons/Modules/MSource/Point.hpp>
+#include <Grid/Hadrons/Modules/MSource/SeqGamma.hpp>
+#include <Grid/Hadrons/Modules/MSource/Z2.hpp>
+#include <Grid/Hadrons/Modules/Quark.hpp>

extras/Hadrons/Modules/MAction/DWF.hpp

@@ -0,0 +1,134 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MAction/DWF.hpp
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_DWF_hpp_
+#define Hadrons_DWF_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                     Domain wall quark action                               *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MAction)
+
+class DWFPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(DWFPar,
+                                    std::string, gauge,
+                                    unsigned int, Ls,
+                                    double      , mass,
+                                    double      , M5);
+};
+
+template <typename FImpl>
+class TDWF: public Module<DWFPar>
+{
+public:
+    TYPE_ALIASES(FImpl,);
+public:
+    // constructor
+    TDWF(const std::string name);
+    // destructor
+    virtual ~TDWF(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(DWF, TDWF<FIMPL>, MAction);
+
+/******************************************************************************
+ *                        DWF template implementation                         *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TDWF<FImpl>::TDWF(const std::string name)
+: Module<DWFPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TDWF<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().gauge};
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TDWF<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TDWF<FImpl>::setup(void)
+{
+    unsigned int size;
+    
+    size = 2*env().template lattice4dSize<typename FImpl::DoubledGaugeField>();
+    env().registerObject(getName(), size, par().Ls);
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TDWF<FImpl>::execute(void)
+{
+    LOG(Message) << "Setting up domain wall fermion matrix with m= "
+                 << par().mass << ", M5= " << par().M5 << " and Ls= "
+                 << par().Ls << " using gauge field '" << par().gauge << "'"
+                 << std::endl;
+    env().createGrid(par().Ls);
+    auto &U      = *env().template getObject<LatticeGaugeField>(par().gauge);
+    auto &g4     = *env().getGrid();
+    auto &grb4   = *env().getRbGrid();
+    auto &g5     = *env().getGrid(par().Ls);
+    auto &grb5   = *env().getRbGrid(par().Ls);
+    FMat *fMatPt = new DomainWallFermion<FImpl>(U, g5, grb5, g4, grb4,
+                                                par().mass, par().M5);
+    env().setObject(getName(), fMatPt);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_DWF_hpp_

extras/Hadrons/Modules/MAction/Wilson.hpp

@@ -0,0 +1,126 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MAction/Wilson.hpp
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_Wilson_hpp_
+#define Hadrons_Wilson_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                            TWilson quark action                            *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MAction)
+
+class WilsonPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonPar,
+                                    std::string, gauge,
+                                    double     , mass);
+};
+
+template <typename FImpl>
+class TWilson: public Module<WilsonPar>
+{
+public:
+    TYPE_ALIASES(FImpl,);
+public:
+    // constructor
+    TWilson(const std::string name);
+    // destructor
+    virtual ~TWilson(void) = default;
+    // dependencies/products
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(Wilson, TWilson<FIMPL>, MAction);
+
+/******************************************************************************
+ *                     TWilson template implementation                        *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TWilson<FImpl>::TWilson(const std::string name)
+: Module<WilsonPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TWilson<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().gauge};
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TWilson<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TWilson<FImpl>::setup(void)
+{
+    unsigned int size;
+    
+    size = 2*env().template lattice4dSize<typename FImpl::DoubledGaugeField>();
+    env().registerObject(getName(), size);
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TWilson<FImpl>::execute()
+{
+    LOG(Message) << "Setting up TWilson fermion matrix with m= " << par().mass
+                 << " using gauge field '" << par().gauge << "'" << std::endl;
+    auto &U      = *env().template getObject<LatticeGaugeField>(par().gauge);
+    auto &grid   = *env().getGrid();
+    auto &gridRb = *env().getRbGrid();
+    FMat *fMatPt = new WilsonFermion<FImpl>(U, grid, gridRb, par().mass);
+    env().setObject(getName(), fMatPt);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Wilson_hpp_

extras/Hadrons/Modules/MContraction/Baryon.hpp

@@ -0,0 +1,131 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MContraction/Baryon.hpp
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_Baryon_hpp_
+#define Hadrons_Baryon_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                               Baryon                                       *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MContraction)
+
+class BaryonPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(BaryonPar,
+                                    std::string, q1,
+                                    std::string, q2,
+                                    std::string, q3,
+                                    std::string, output);
+};
+
+template <typename FImpl1, typename FImpl2, typename FImpl3>
+class TBaryon: public Module<BaryonPar>
+{
+public:
+    TYPE_ALIASES(FImpl1, 1);
+    TYPE_ALIASES(FImpl2, 2);
+    TYPE_ALIASES(FImpl3, 3);
+    class Result: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
+                                        std::vector<std::vector<std::vector<Complex>>>, corr);
+    };
+public:
+    // constructor
+    TBaryon(const std::string name);
+    // destructor
+    virtual ~TBaryon(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(Baryon, ARG(TBaryon<FIMPL, FIMPL, FIMPL>), MContraction);
+
+/******************************************************************************
+ *                         TBaryon implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2, typename FImpl3>
+TBaryon<FImpl1, FImpl2, FImpl3>::TBaryon(const std::string name)
+: Module<BaryonPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2, typename FImpl3>
+std::vector<std::string> TBaryon<FImpl1, FImpl2, FImpl3>::getInput(void)
+{
+    std::vector<std::string> input = {par().q1, par().q2, par().q3};
+    
+    return input;
+}
+
+template <typename FImpl1, typename FImpl2, typename FImpl3>
+std::vector<std::string> TBaryon<FImpl1, FImpl2, FImpl3>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2, typename FImpl3>
+void TBaryon<FImpl1, FImpl2, FImpl3>::execute(void)
+{
+    LOG(Message) << "Computing baryon contractions '" << getName() << "' using"
+                 << " quarks '" << par().q1 << "', '" << par().q2 << "', and '"
+                 << par().q3 << "'" << std::endl;
+    
+    XmlWriter             writer(par().output);
+    PropagatorField1      &q1 = *env().template getObject<PropagatorField1>(par().q1);
+    PropagatorField2      &q2 = *env().template getObject<PropagatorField2>(par().q2);
+    PropagatorField3      &q3 = *env().template getObject<PropagatorField3>(par().q2);
+    LatticeComplex        c(env().getGrid());
+    Result                result;
+    
+    // FIXME: do contractions
+    
+    write(writer, "meson", result);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Baryon_hpp_

extras/Hadrons/Modules/MContraction/Meson.hpp

@@ -0,0 +1,148 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MContraction/Meson.hpp
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_Meson_hpp_
+#define Hadrons_Meson_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                                TMeson                                       *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MContraction)
+
+class MesonPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(MesonPar,
+                                    std::string, q1,
+                                    std::string, q2,
+                                    std::string, output);
+};
+
+template <typename FImpl1, typename FImpl2>
+class TMeson: public Module<MesonPar>
+{
+public:
+    TYPE_ALIASES(FImpl1, 1);
+    TYPE_ALIASES(FImpl2, 2);
+    class Result: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
+                                        std::vector<std::vector<std::vector<Complex>>>, corr);
+    };
+public:
+    // constructor
+    TMeson(const std::string name);
+    // destructor
+    virtual ~TMeson(void) = default;
+    // dependencies/products
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(Meson, ARG(TMeson<FIMPL, FIMPL>), MContraction);
+
+/******************************************************************************
+ *                           TMeson implementation                            *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2>
+TMeson<FImpl1, FImpl2>::TMeson(const std::string name)
+: Module<MesonPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2>
+std::vector<std::string> TMeson<FImpl1, FImpl2>::getInput(void)
+{
+    std::vector<std::string> input = {par().q1, par().q2};
+    
+    return input;
+}
+
+template <typename FImpl1, typename FImpl2>
+std::vector<std::string> TMeson<FImpl1, FImpl2>::getOutput(void)
+{
+    std::vector<std::string> output = {getName()};
+    
+    return output;
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2>
+void TMeson<FImpl1, FImpl2>::execute(void)
+{
+    LOG(Message) << "Computing meson contractions '" << getName() << "' using"
+                 << " quarks '" << par().q1 << "' and '" << par().q2 << "'"
+                 << std::endl;
+    
+    XmlWriter             writer(par().output);
+    PropagatorField1      &q1 = *env().template getObject<PropagatorField1>(par().q1);
+    PropagatorField2      &q2 = *env().template getObject<PropagatorField2>(par().q2);
+    LatticeComplex        c(env().getGrid());
+    SpinMatrix            g[Ns*Ns], g5;
+    std::vector<TComplex> buf;
+    Result                result;
+    
+    g5 = makeGammaProd(Ns*Ns - 1);
+    result.corr.resize(Ns*Ns);
+    for (unsigned int i = 0; i < Ns*Ns; ++i)
+    {
+        g[i] = makeGammaProd(i);
+    }
+    for (unsigned int iSink = 0; iSink < Ns*Ns; ++iSink)
+    {
+        result.corr[iSink].resize(Ns*Ns);
+        for (unsigned int iSrc = 0; iSrc < Ns*Ns; ++iSrc)
+        {
+            c = trace(g[iSink]*q1*g[iSrc]*g5*adj(q2)*g5);
+            sliceSum(c, buf, Tp);
+            result.corr[iSink][iSrc].resize(buf.size());
+            for (unsigned int t = 0; t < buf.size(); ++t)
+            {
+                result.corr[iSink][iSrc][t] = TensorRemove(buf[t]);
+            }
+        }
+    }
+    write(writer, "meson", result);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Meson_hpp_

extras/Hadrons/Modules/MGauge/Load.cc

@@ -0,0 +1,78 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MGauge/Load.cc
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#include <Grid/Hadrons/Modules/MGauge/Load.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MGauge;
+
+/******************************************************************************
+*                           TLoad implementation                               *
+******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+TLoad::TLoad(const std::string name)
+: Module<LoadPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+std::vector<std::string> TLoad::getInput(void)
+{
+    std::vector<std::string> in;
+    
+    return in;
+}
+
+std::vector<std::string> TLoad::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+void TLoad::setup(void)
+{
+    env().registerLattice<LatticeGaugeField>(getName());
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+void TLoad::execute(void)
+{
+    NerscField  header;
+    std::string fileName = par().file + "."
+                           + std::to_string(env().getTrajectory());
+    
+    LOG(Message) << "Loading NERSC configuration from file '" << fileName
+                 << "'" << std::endl;
+    LatticeGaugeField &U = *env().createLattice<LatticeGaugeField>(getName());
+    NerscIO::readConfiguration(U, header, fileName);
+    LOG(Message) << "NERSC header:" << std::endl;
+    dump_nersc_header(header, LOG(Message));
+}

extras/Hadrons/Modules/MGauge/Load.hpp

@@ -0,0 +1,73 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MGauge/Load.hpp
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_Load_hpp_
+#define Hadrons_Load_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         Load a NERSC configuration                         *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MGauge)
+
+class LoadPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(LoadPar,
+                                    std::string, file);
+};
+
+class TLoad: public Module<LoadPar>
+{
+public:
+    // constructor
+    TLoad(const std::string name);
+    // destructor
+    virtual ~TLoad(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(Load, TLoad, MGauge);
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Load_hpp_

extras/Hadrons/Modules/MGauge/Random.cc

@@ -0,0 +1,69 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MGauge/Random.cc
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#include <Grid/Hadrons/Modules/MGauge/Random.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MGauge;
+
+/******************************************************************************
+*                           TRandom implementation                            *
+******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+TRandom::TRandom(const std::string name)
+: Module<NoPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+std::vector<std::string> TRandom::getInput(void)
+{
+    return std::vector<std::string>();
+}
+
+std::vector<std::string> TRandom::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+void TRandom::setup(void)
+{
+    env().registerLattice<LatticeGaugeField>(getName());
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+void TRandom::execute(void)
+{
+    LOG(Message) << "Generating random gauge configuration" << std::endl;
+    LatticeGaugeField &U = *env().createLattice<LatticeGaugeField>(getName());
+    SU3::HotConfiguration(*env().get4dRng(), U);
+}

extras/Hadrons/Modules/MGauge/Random.hpp

@@ -0,0 +1,66 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MGauge/Random.hpp
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_Random_hpp_
+#define Hadrons_Random_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                             Random gauge                                   *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MGauge)
+
+class TRandom: public Module<NoPar>
+{
+public:
+    // constructor
+    TRandom(const std::string name);
+    // destructor
+    virtual ~TRandom(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(Random, TRandom, MGauge);
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Random_hpp_

extras/Hadrons/Modules/MGauge/Unit.cc

@@ -0,0 +1,69 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MGauge/Unit.cc
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#include <Grid/Hadrons/Modules/MGauge/Unit.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MGauge;
+
+/******************************************************************************
+*                            TUnit implementation                             *
+******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+TUnit::TUnit(const std::string name)
+: Module<NoPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+std::vector<std::string> TUnit::getInput(void)
+{
+    return std::vector<std::string>();
+}
+
+std::vector<std::string> TUnit::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+void TUnit::setup(void)
+{
+    env().registerLattice<LatticeGaugeField>(getName());
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+void TUnit::execute(void)
+{
+    LOG(Message) << "Creating unit gauge configuration" << std::endl;
+    LatticeGaugeField &U = *env().createLattice<LatticeGaugeField>(getName());
+    SU3::ColdConfiguration(*env().get4dRng(), U);
+}

extras/Hadrons/Modules/MGauge/Unit.hpp

@@ -0,0 +1,66 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MGauge/Unit.hpp
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_Unit_hpp_
+#define Hadrons_Unit_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                              Unit gauge                                    *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MGauge)
+
+class TUnit: public Module<NoPar>
+{
+public:
+    // constructor
+    TUnit(const std::string name);
+    // destructor
+    virtual ~TUnit(void) = default;
+    // dependencies/products
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(Unit, TUnit, MGauge);
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Unit_hpp_

extras/Hadrons/Modules/MSolver/RBPrecCG.hpp

@@ -0,0 +1,132 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MSolver/RBPrecCG.hpp
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_RBPrecCG_hpp_
+#define Hadrons_RBPrecCG_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                     Schur red-black preconditioned CG                      *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MSolver)
+
+class RBPrecCGPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(RBPrecCGPar,
+                                    std::string, action,
+                                    double     , residual);
+};
+
+template <typename FImpl>
+class TRBPrecCG: public Module<RBPrecCGPar>
+{
+public:
+    TYPE_ALIASES(FImpl,);
+public:
+    // constructor
+    TRBPrecCG(const std::string name);
+    // destructor
+    virtual ~TRBPrecCG(void) = default;
+    // dependencies/products
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(RBPrecCG, TRBPrecCG<FIMPL>, MSolver);
+
+/******************************************************************************
+ *                      TRBPrecCG template implementation                     *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TRBPrecCG<FImpl>::TRBPrecCG(const std::string name)
+: Module(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TRBPrecCG<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().action};
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TRBPrecCG<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TRBPrecCG<FImpl>::setup(void)
+{
+    auto Ls = env().getObjectLs(par().action);
+    
+    env().registerObject(getName(), 0, Ls);
+    env().addOwnership(getName(), par().action);
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TRBPrecCG<FImpl>::execute(void)
+{
+    auto &mat   = *(env().template getObject<FMat>(par().action));
+    auto solver = [&mat, this](FermionField &sol, const FermionField &source)
+    {
+        ConjugateGradient<FermionField>           cg(par().residual, 10000);
+        SchurRedBlackDiagMooeeSolve<FermionField> schurSolver(cg);
+        
+        schurSolver(mat, source, sol);
+    };
+    
+    LOG(Message) << "setting up Schur red-black preconditioned CG for"
+                 << " action '" << par().action << "' with residual "
+                 << par().residual << std::endl;
+    env().setObject(getName(), new SolverFn(solver));
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_RBPrecCG_hpp_

extras/Hadrons/Modules/MSource/Point.hpp

@@ -0,0 +1,135 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MSource/Point.hpp
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_Point_hpp_
+#define Hadrons_Point_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/*
+ 
+ Point source
+ ------------
+ * src_x = delta_x,position
+ 
+ * options:
+ - position: space-separated integer sequence (e.g. "0 1 1 0")
+ 
+ */
+
+/******************************************************************************
+ *                                  TPoint                                     *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MSource)
+
+class PointPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(PointPar,
+                                    std::string, position);
+};
+
+template <typename FImpl>
+class TPoint: public Module<PointPar>
+{
+public:
+    TYPE_ALIASES(FImpl,);
+public:
+    // constructor
+    TPoint(const std::string name);
+    // destructor
+    virtual ~TPoint(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(Point, TPoint<FIMPL>, MSource);
+
+/******************************************************************************
+ *                       TPoint template implementation                       *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TPoint<FImpl>::TPoint(const std::string name)
+: Module<PointPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TPoint<FImpl>::getInput(void)
+{
+    std::vector<std::string> in;
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TPoint<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TPoint<FImpl>::setup(void)
+{
+    env().template registerLattice<PropagatorField>(getName());
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TPoint<FImpl>::execute(void)
+{
+    std::vector<int> position = strToVec<int>(par().position);
+    typename SitePropagator::scalar_object id;
+    
+    LOG(Message) << "Creating point source at position [" << par().position
+                 << "]" << std::endl;
+    PropagatorField &src = *env().template createLattice<PropagatorField>(getName());
+    id  = 1.;
+    src = zero;
+    pokeSite(id, src, position);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Point_hpp_

extras/Hadrons/Modules/MSource/SeqGamma.hpp

@@ -0,0 +1,164 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MSource/SeqGamma.hpp
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_SeqGamma_hpp_
+#define Hadrons_SeqGamma_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/*
+ 
+ Sequential source
+ -----------------------------
+ * src_x = q_x * theta(x_3 - tA) * theta(tB - x_3) * gamma * exp(i x.mom)
+ 
+ * options:
+ - q: input propagator (string)
+ - tA: begin timeslice (integer)
+ - tB: end timesilce (integer)
+ - gamma: gamma product to insert (integer)
+ - mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0.")
+ 
+ */
+
+/******************************************************************************
+ *                         SeqGamma                                 *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MSource)
+
+class SeqGammaPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(SeqGammaPar,
+                                    std::string,  q,
+                                    unsigned int, tA,
+                                    unsigned int, tB,
+                                    unsigned int, gamma,
+                                    std::string,  mom);
+};
+
+template <typename FImpl>
+class TSeqGamma: public Module<SeqGammaPar>
+{
+public:
+    TYPE_ALIASES(FImpl,);
+public:
+    // constructor
+    TSeqGamma(const std::string name);
+    // destructor
+    virtual ~TSeqGamma(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(SeqGamma, TSeqGamma<FIMPL>, MSource);
+
+/******************************************************************************
+ *                         TSeqGamma implementation                           *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TSeqGamma<FImpl>::TSeqGamma(const std::string name)
+: Module<SeqGammaPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TSeqGamma<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().q};
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TSeqGamma<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TSeqGamma<FImpl>::setup(void)
+{
+    env().template registerLattice<PropagatorField>(getName());
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TSeqGamma<FImpl>::execute(void)
+{
+    if (par().tA == par().tB)
+    {
+        LOG(Message) << "Generating gamma_" << par().gamma
+                     << " sequential source at t= " << par().tA << std::endl;
+    }
+    else
+    {
+        LOG(Message) << "Generating gamma_" << par().gamma
+                     << " sequential source for "
+                     << par().tA << " <= t <= " << par().tB << std::endl;
+    }
+    PropagatorField &src = *env().template createLattice<PropagatorField>(getName());
+    PropagatorField &q   = *env().template getObject<PropagatorField>(par().q);
+    Lattice<iScalar<vInteger>> t(env().getGrid());
+    LatticeComplex             ph(env().getGrid()), coor(env().getGrid());
+    SpinMatrix                 g;
+    std::vector<Real>          p;
+    Complex                    i(0.0,1.0);
+    
+    g  = makeGammaProd(par().gamma);
+    p  = strToVec<Real>(par().mom);
+    ph = zero;
+    for(unsigned int mu = 0; mu < Nd; mu++)
+    {
+        LatticeCoordinate(coor, mu);
+        ph = ph + p[mu]*coor;
+    }
+    ph = exp(i*ph);
+    LatticeCoordinate(t, Tp);
+    src = where((t >= par().tA) and (t <= par().tB), g*ph*q, 0.*q);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_SeqGamma_hpp_

extras/Hadrons/Modules/MSource/Z2.hpp

@@ -0,0 +1,151 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MSource/Z2.hpp
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_Z2_hpp_
+#define Hadrons_Z2_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/*
+ 
+ Z_2 stochastic source
+ -----------------------------
+ * src_x = eta_x * theta(x_3 - tA) * theta(tB - x_3)
+ 
+ the eta_x are independent uniform random numbers in {+/- 1 +/- i}
+ 
+ * options:
+ - tA: begin timeslice (integer)
+ - tB: end timesilce (integer)
+ 
+ */
+ 
+/******************************************************************************
+ *                          Z2 stochastic source                              *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MSource)
+
+class Z2Par: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(Z2Par,
+                                    unsigned int, tA,
+                                    unsigned int, tB);
+};
+
+template <typename FImpl>
+class TZ2: public Module<Z2Par>
+{
+public:
+    TYPE_ALIASES(FImpl,);
+public:
+    // constructor
+    TZ2(const std::string name);
+    // destructor
+    virtual ~TZ2(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(Z2, TZ2<FIMPL>, MSource);
+
+/******************************************************************************
+ *                       TZ2 template implementation                          *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TZ2<FImpl>::TZ2(const std::string name)
+: Module<Z2Par>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TZ2<FImpl>::getInput(void)
+{
+    std::vector<std::string> in;
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TZ2<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TZ2<FImpl>::setup(void)
+{
+    env().template registerLattice<PropagatorField>(getName());
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TZ2<FImpl>::execute(void)
+{
+    Lattice<iScalar<vInteger>> t(env().getGrid());
+    LatticeComplex             eta(env().getGrid());
+    Complex                    shift(1., 1.);
+    
+    if (par().tA == par().tB)
+    {
+        LOG(Message) << "Generating Z_2 wall source at t= " << par().tA
+        << std::endl;
+    }
+    else
+    {
+        LOG(Message) << "Generating Z_2 band for " << par().tA << " <= t <= "
+        << par().tB << std::endl;
+    }
+    PropagatorField &src = *env().template createLattice<PropagatorField>(getName());
+    LatticeCoordinate(t, Tp);
+    bernoulli(*env().get4dRng(), eta);
+    eta = (2.*eta - shift)*(1./::sqrt(2.));
+    eta = where((t >= par().tA) and (t <= par().tB), eta, 0.*eta);
+    src = 1.;
+    src = src*eta;
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Z2_hpp_

extras/Hadrons/Modules/Quark.hpp

@@ -0,0 +1,185 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/Quark.hpp
+
+Copyright (C) 2015
+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
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_Quark_hpp_
+#define Hadrons_Quark_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                               TQuark                                       *
+ ******************************************************************************/
+class QuarkPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(QuarkPar,
+                                    std::string, source,
+                                    std::string, solver);
+};
+
+template <typename FImpl>
+class TQuark: public Module<QuarkPar>
+{
+public:
+    TYPE_ALIASES(FImpl,);
+public:
+    // constructor
+    TQuark(const std::string name);
+    // destructor
+    virtual ~TQuark(void) = default;
+    // dependencies/products
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+private:
+    unsigned int Ls_;
+    SolverFn     *solver_{nullptr};
+};
+
+MODULE_REGISTER(Quark, TQuark<FIMPL>);
+
+/******************************************************************************
+ *                          TQuark implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TQuark<FImpl>::TQuark(const std::string name)
+: Module(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TQuark<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().source, par().solver};
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TQuark<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName(), getName() + "_5d"};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TQuark<FImpl>::setup(void)
+{
+    Ls_ = env().getObjectLs(par().solver);
+    env().template registerLattice<PropagatorField>(getName());
+    if (Ls_ > 1)
+    {
+        env().template registerLattice<PropagatorField>(getName() + "_5d", Ls_);
+    }
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TQuark<FImpl>::execute(void)
+{
+    LOG(Message) << "Computing quark propagator '" << getName() << "'"
+                 << std::endl;
+    
+    FermionField    source(env().getGrid(Ls_)), sol(env().getGrid(Ls_)),
+                    tmp(env().getGrid());
+    std::string     propName = (Ls_ == 1) ? getName() : (getName() + "_5d");
+    PropagatorField &prop    = *env().template createLattice<PropagatorField>(propName);
+    PropagatorField &fullSrc = *env().template getObject<PropagatorField>(par().source);
+    SolverFn        &solver  = *env().template getObject<SolverFn>(par().solver);
+    if (Ls_ > 1)
+    {
+        env().template createLattice<PropagatorField>(getName());
+    }
+    
+    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 < Nc; ++c)
+    {
+        LOG(Message) << "Inversion for spin= " << s << ", color= " << c
+        << std::endl;
+        // source conversion for 4D sources
+        if (!env().isObject5d(par().source))
+        {
+            if (Ls_ == 1)
+            {
+                PropToFerm(source, fullSrc, s, c);
+            }
+            else
+            {
+                source = zero;
+                PropToFerm(tmp, fullSrc, s, c);
+                InsertSlice(tmp, source, 0, 0);
+                InsertSlice(tmp, source, Ls_-1, 0);
+                axpby_ssp_pplus(source, 0., source, 1., source, 0, 0);
+                axpby_ssp_pminus(source, 0., source, 1., source, Ls_-1, Ls_-1);
+            }
+        }
+        // source conversion for 5D sources
+        else
+        {
+            if (Ls_ != env().getObjectLs(par().source))
+            {
+                HADRON_ERROR("Ls mismatch between quark action and source");
+            }
+            else
+            {
+                PropToFerm(source, fullSrc, s, c);
+            }
+        }
+        sol = zero;
+        solver(sol, source);
+        FermToProp(prop, sol, s, c);
+        // create 4D propagators from 5D one if necessary
+        if (Ls_ > 1)
+        {
+            PropagatorField &p4d =
+                *env().template getObject<PropagatorField>(getName());
+            
+            axpby_ssp_pminus(sol, 0., sol, 1., sol, 0, 0);
+            axpby_ssp_pplus(sol, 0., sol, 1., sol, 0, Ls_-1);
+            ExtractSlice(tmp, sol, 0, 0);
+            FermToProp(p4d, tmp, s, c);
+        }
+    }
+}
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Quark_hpp_

extras/Hadrons/Modules/templates/Module.cc.template

@@ -0,0 +1,39 @@
+#include <Grid/Hadrons/Modules/___FILEBASENAME___.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+
+/******************************************************************************
+*                  T___FILEBASENAME___ implementation                             *
+******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+T___FILEBASENAME___::T___FILEBASENAME___(const std::string name)
+: Module<___FILEBASENAME___Par>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+std::vector<std::string> T___FILEBASENAME___::getInput(void)
+{
+    std::vector<std::string> in;
+    
+    return in;
+}
+
+std::vector<std::string> T___FILEBASENAME___::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+void T___FILEBASENAME___::setup(void)
+{
+
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+void T___FILEBASENAME___::execute(void)
+{
+
+}

extras/Hadrons/Modules/templates/Module.hpp.template

@@ -0,0 +1,40 @@
+#ifndef Hadrons____FILEBASENAME____hpp_
+#define Hadrons____FILEBASENAME____hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         ___FILEBASENAME___                                 *
+ ******************************************************************************/
+class ___FILEBASENAME___Par: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(___FILEBASENAME___Par, 
+                                    unsigned int, i);
+};
+
+class T___FILEBASENAME___: public Module<___FILEBASENAME___Par>
+{
+public:
+    // constructor
+    T___FILEBASENAME___(const std::string name);
+    // destructor
+    virtual ~T___FILEBASENAME___(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER(___FILEBASENAME___, T___FILEBASENAME___);
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons____FILEBASENAME____hpp_

extras/Hadrons/Modules/templates/Module_in_NS.cc.template

@@ -0,0 +1,40 @@
+#include <Grid/Hadrons/Modules/___NAMESPACE___/___FILEBASENAME___.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace ___NAMESPACE___;
+
+/******************************************************************************
+*                  T___FILEBASENAME___ implementation                             *
+******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+T___FILEBASENAME___::T___FILEBASENAME___(const std::string name)
+: Module<___FILEBASENAME___Par>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+std::vector<std::string> T___FILEBASENAME___::getInput(void)
+{
+    std::vector<std::string> in;
+    
+    return in;
+}
+
+std::vector<std::string> T___FILEBASENAME___::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+void T___FILEBASENAME___::setup(void)
+{
+
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+void T___FILEBASENAME___::execute(void)
+{
+
+}

extras/Hadrons/Modules/templates/Module_in_NS.hpp.template

@@ -0,0 +1,44 @@
+#ifndef Hadrons____FILEBASENAME____hpp_
+#define Hadrons____FILEBASENAME____hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         ___FILEBASENAME___                                 *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(___NAMESPACE___)
+
+class ___FILEBASENAME___Par: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(___FILEBASENAME___Par,
+                                    unsigned int, i);
+};
+
+class T___FILEBASENAME___: public Module<___FILEBASENAME___Par>
+{
+public:
+    // constructor
+    T___FILEBASENAME___(const std::string name);
+    // destructor
+    virtual ~T___FILEBASENAME___(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(___FILEBASENAME___, T___FILEBASENAME___, ___NAMESPACE___);
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons____FILEBASENAME____hpp_

extras/Hadrons/Modules/templates/Module_tmp.hpp.template

@@ -0,0 +1,81 @@
+#ifndef Hadrons____FILEBASENAME____hpp_
+#define Hadrons____FILEBASENAME____hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         ___FILEBASENAME___                                 *
+ ******************************************************************************/
+class ___FILEBASENAME___Par: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(___FILEBASENAME___Par,
+                                    unsigned int, i);
+};
+
+template <typename FImpl>
+class T___FILEBASENAME___: public Module<___FILEBASENAME___Par>
+{
+public:
+    // constructor
+    T___FILEBASENAME___(const std::string name);
+    // destructor
+    virtual ~T___FILEBASENAME___(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER(___FILEBASENAME___, T___FILEBASENAME___<FIMPL>);
+
+/******************************************************************************
+ *                 T___FILEBASENAME___ implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+T___FILEBASENAME___<FImpl>::T___FILEBASENAME___(const std::string name)
+: Module<___FILEBASENAME___Par>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> T___FILEBASENAME___<FImpl>::getInput(void)
+{
+    std::vector<std::string> in;
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> T___FILEBASENAME___<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void T___FILEBASENAME___<FImpl>::setup(void)
+{
+    
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void T___FILEBASENAME___<FImpl>::execute(void)
+{
+    
+}
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons____FILEBASENAME____hpp_

extras/Hadrons/Modules/templates/Module_tmp_in_NS.hpp.template

@@ -0,0 +1,85 @@
+#ifndef Hadrons____FILEBASENAME____hpp_
+#define Hadrons____FILEBASENAME____hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         ___FILEBASENAME___                                 *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(___NAMESPACE___)
+
+class ___FILEBASENAME___Par: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(___FILEBASENAME___Par,
+                                    unsigned int, i);
+};
+
+template <typename FImpl>
+class T___FILEBASENAME___: public Module<___FILEBASENAME___Par>
+{
+public:
+    // constructor
+    T___FILEBASENAME___(const std::string name);
+    // destructor
+    virtual ~T___FILEBASENAME___(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(___FILEBASENAME___, T___FILEBASENAME___<FIMPL>, ___NAMESPACE___);
+
+/******************************************************************************
+ *                 T___FILEBASENAME___ implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+T___FILEBASENAME___<FImpl>::T___FILEBASENAME___(const std::string name)
+: Module<___FILEBASENAME___Par>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> T___FILEBASENAME___<FImpl>::getInput(void)
+{
+    std::vector<std::string> in;
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> T___FILEBASENAME___<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void T___FILEBASENAME___<FImpl>::setup(void)
+{
+    
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void T___FILEBASENAME___<FImpl>::execute(void)
+{
+    
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons____FILEBASENAME____hpp_

extras/Hadrons/add_module.sh

@@ -0,0 +1,31 @@
+#!/usr/bin/env bash
+
+if (( $# != 1 && $# != 2)); then
+    echo "usage: `basename $0` <module name> [<namespace>]" 1>&2
+    exit 1
+fi
+NAME=$1
+NS=$2
+
+if (( $# == 1 )); then
+	if [ -e "Modules/${NAME}.cc" ] || [ -e "Modules/${NAME}.hpp" ]; then
+	    echo "error: files Modules/${NAME}.* already exists" 1>&2
+	    exit 1
+	fi
+	sed "s/___FILEBASENAME___/${NAME}/g" Modules/templates/Module.cc.template > Modules/${NAME}.cc
+	sed "s/___FILEBASENAME___/${NAME}/g" Modules/templates/Module.hpp.template > Modules/${NAME}.hpp
+elif (( $# == 2 )); then
+	mkdir -p Modules/${NS}
+	if [ -e "Modules/${NS}/${NAME}.cc" ] || [ -e "Modules/${NS}/${NAME}.hpp" ]; then
+	    echo "error: files Modules/${NS}/${NAME}.* already exists" 1>&2
+	    exit 1
+	fi
+	TMPCC=".${NS}.${NAME}.tmp.cc"
+	TMPHPP=".${NS}.${NAME}.tmp.hpp"
+	sed "s/___FILEBASENAME___/${NAME}/g" Modules/templates/Module_in_NS.cc.template  > ${TMPCC}
+	sed "s/___FILEBASENAME___/${NAME}/g" Modules/templates/Module_in_NS.hpp.template > ${TMPHPP}
+	sed "s/___NAMESPACE___/${NS}/g" ${TMPCC}  > Modules/${NS}/${NAME}.cc
+	sed "s/___NAMESPACE___/${NS}/g" ${TMPHPP} > Modules/${NS}/${NAME}.hpp
+	rm -f ${TMPCC} ${TMPHPP}
+fi
+./make_module_list.sh

extras/Hadrons/add_module_template.sh

@@ -0,0 +1,28 @@
+#!/usr/bin/env bash
+
+if (( $# != 1 && $# != 2)); then
+    echo "usage: `basename $0` <module name> [<namespace>]" 1>&2
+    exit 1
+fi
+NAME=$1
+NS=$2
+
+if (( $# == 1 )); then
+	if [ -e "Modules/${NAME}.cc" ] || [ -e "Modules/${NAME}.hpp" ]; then
+	    echo "error: files Modules/${NAME}.* already exists" 1>&2
+	    exit 1
+	fi
+	sed "s/___FILEBASENAME___/${NAME}/g" Modules/templates/Module_tmp.hpp.template > Modules/${NAME}.hpp
+elif (( $# == 2 )); then
+	mkdir -p Modules/${NS}
+	if [ -e "Modules/${NS}/${NAME}.cc" ] || [ -e "Modules/${NS}/${NAME}.hpp" ]; then
+	    echo "error: files Modules/${NS}/${NAME}.* already exists" 1>&2
+	    exit 1
+	fi
+	TMPCC=".${NS}.${NAME}.tmp.cc"
+	TMPHPP=".${NS}.${NAME}.tmp.hpp"
+	sed "s/___FILEBASENAME___/${NAME}/g" Modules/templates/Module_tmp_in_NS.hpp.template > ${TMPHPP}
+	sed "s/___NAMESPACE___/${NS}/g" ${TMPHPP} > Modules/${NS}/${NAME}.hpp
+	rm -f ${TMPCC} ${TMPHPP}
+fi
+./make_module_list.sh

extras/Hadrons/make_module_list.sh

@@ -0,0 +1,12 @@
+#!/usr/bin/env bash
+
+echo 'modules_cc =\' > modules.inc
+find Modules -name '*.cc' -type f -print | sed 's/^/  /;$q;s/$/ \\/' >> modules.inc
+echo '' >> modules.inc
+echo 'modules_hpp =\' >> modules.inc
+find Modules -name '*.hpp' -type f -print | sed 's/^/  /;$q;s/$/ \\/' >> modules.inc
+echo '' >> modules.inc
+rm -f Modules.hpp
+for f in `find Modules -name '*.hpp'`; do
+	echo "#include <Grid/Hadrons/${f}>" >> Modules.hpp
+done

extras/Hadrons/modules.inc

@@ -0,0 +1,19 @@
+modules_cc =\
+  Modules/MGauge/Load.cc \
+  Modules/MGauge/Random.cc \
+  Modules/MGauge/Unit.cc
+
+modules_hpp =\
+  Modules/MAction/DWF.hpp \
+  Modules/MAction/Wilson.hpp \
+  Modules/MContraction/Baryon.hpp \
+  Modules/MContraction/Meson.hpp \
+  Modules/MGauge/Load.hpp \
+  Modules/MGauge/Random.hpp \
+  Modules/MGauge/Unit.hpp \
+  Modules/MSolver/RBPrecCG.hpp \
+  Modules/MSource/Point.hpp \
+  Modules/MSource/SeqGamma.hpp \
+  Modules/MSource/Z2.hpp \
+  Modules/Quark.hpp
+

extras/Makefile.am

@@ -0,0 +1 @@
+SUBDIRS = Hadrons

lib/Hadrons

@@ -0,0 +1 @@
+../extras/Hadrons

lib/Log.h

@@ -110,8 +110,8 @@ public:
   friend std::ostream& operator<< (std::ostream& stream, Logger& log){
 
     if ( log.active ) {
-      stream << log.background()<< log.topName << log.background()<< " : ";
+      stream << log.background()<< std::setw(10) << std::left << log.topName << log.background()<< " : ";
-      stream << log.colour() <<std::setw(14) << std::left << log.name << log.background() << " : ";
+      stream << log.colour() << std::setw(14) << std::left << log.name << log.background() << " : ";
       if ( log.timestamp ) {
 	StopWatch.Stop();
 	GridTime now = StopWatch.Elapsed();

lib/Threads.h

@@ -46,11 +46,13 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #endif
 #define PARALLEL_NESTED_LOOP2 _Pragma("omp parallel for collapse(2)")
 #define PARALLEL_REGION       _Pragma("omp parallel")
+#define PARALLEL_CRITICAL     _Pragma("omp critical")
 #else
 #define PARALLEL_FOR_LOOP
 #define PARALLEL_FOR_LOOP_INTERN
 #define PARALLEL_NESTED_LOOP2
 #define PARALLEL_REGION
+#define PARALLEL_CRITICAL
 #endif
 
 namespace Grid {

lib/qcd/action/fermion/FermionOperatorImpl.h

@@ -48,8 +48,10 @@ namespace QCD {
   //    typedef typename XXX         GaugeField;
   //    typedef typename XXX      GaugeActField;
   //    typedef typename XXX       FermionField;
+  //    typedef typename XXX    PropagatorField;
   //    typedef typename XXX  DoubledGaugeField;
   //    typedef typename XXX         SiteSpinor;
+  //    typedef typename XXX     SitePropagator;
   //    typedef typename XXX     SiteHalfSpinor;	
   //    typedef typename XXX         Compressor;	
   //
@@ -95,13 +97,15 @@ namespace QCD {
   
 #define INHERIT_FIMPL_TYPES(Impl)\
   typedef typename Impl::FermionField           FermionField;		\
+  typedef typename Impl::PropagatorField     PropagatorField;		\
   typedef typename Impl::DoubledGaugeField DoubledGaugeField;		\
   typedef typename Impl::SiteSpinor               SiteSpinor;		\
+  typedef typename Impl::SitePropagator       SitePropagator;		\
   typedef typename Impl::SiteHalfSpinor       SiteHalfSpinor;		\
   typedef typename Impl::Compressor               Compressor;		\
   typedef typename Impl::StencilImpl             StencilImpl;		\
   typedef typename Impl::ImplParams               ImplParams;	        \
-  typedef typename Impl::Coeff_t       Coeff_t;
+  typedef typename Impl::Coeff_t                     Coeff_t;           \
   
 #define INHERIT_IMPL_TYPES(Base) \
   INHERIT_GIMPL_TYPES(Base)	 \
@@ -127,14 +131,17 @@ namespace QCD {
     INHERIT_GIMPL_TYPES(Gimpl);
       
     template <typename vtype> using iImplSpinor            = iScalar<iVector<iVector<vtype, Dimension>, Ns> >;
+    template <typename vtype> using iImplPropagator        = iScalar<iMatrix<iMatrix<vtype, Dimension>, Ns> >;
     template <typename vtype> using iImplHalfSpinor        = iScalar<iVector<iVector<vtype, Dimension>, Nhs> >;
     template <typename vtype> using iImplDoubledGaugeField = iVector<iScalar<iMatrix<vtype, Dimension> >, Nds>;
     
     typedef iImplSpinor<Simd>            SiteSpinor;
+    typedef iImplPropagator<Simd>        SitePropagator;
     typedef iImplHalfSpinor<Simd>        SiteHalfSpinor;
     typedef iImplDoubledGaugeField<Simd> SiteDoubledGaugeField;
     
     typedef Lattice<SiteSpinor>            FermionField;
+    typedef Lattice<SitePropagator>        PropagatorField;
     typedef Lattice<SiteDoubledGaugeField> DoubledGaugeField;
     
     typedef WilsonCompressor<SiteHalfSpinor, SiteSpinor> Compressor;
@@ -216,14 +223,17 @@ class DomainWallVec5dImpl :  public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepres
   INHERIT_GIMPL_TYPES(Gimpl);
   
   template <typename vtype> using iImplSpinor            = iScalar<iVector<iVector<vtype, Nrepresentation>, Ns> >;
+  template <typename vtype> using iImplPropagator        = iScalar<iMatrix<iMatrix<vtype, Nrepresentation>, Ns> >;
   template <typename vtype> using iImplHalfSpinor        = iScalar<iVector<iVector<vtype, Nrepresentation>, Nhs> >;
   template <typename vtype> using iImplDoubledGaugeField = iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nds>;
   template <typename vtype> using iImplGaugeField        = iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nd>;
   template <typename vtype> using iImplGaugeLink         = iScalar<iScalar<iMatrix<vtype, Nrepresentation> > >;
   
   typedef iImplSpinor<Simd> SiteSpinor;
+  typedef iImplPropagator<Simd> SitePropagator;
   typedef iImplHalfSpinor<Simd> SiteHalfSpinor;
   typedef Lattice<SiteSpinor> FermionField;
+  typedef Lattice<SitePropagator> PropagatorField;
   
   // Make the doubled gauge field a *scalar*
   typedef iImplDoubledGaugeField<typename Simd::scalar_type>  SiteDoubledGaugeField;  // This is a scalar
@@ -315,14 +325,17 @@ class GparityWilsonImpl : public ConjugateGaugeImpl<GaugeImplTypes<S, Nrepresent
  INHERIT_GIMPL_TYPES(Gimpl);
       
  template <typename vtype> using iImplSpinor            = iVector<iVector<iVector<vtype, Nrepresentation>, Ns>, Ngp>;
+ template <typename vtype> using iImplPropagator        = iVector<iMatrix<iMatrix<vtype, Nrepresentation>, Ns>, Ngp >;
  template <typename vtype> using iImplHalfSpinor        = iVector<iVector<iVector<vtype, Nrepresentation>, Nhs>, Ngp>;
  template <typename vtype> using iImplDoubledGaugeField = iVector<iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nds>, Ngp>;
       
  typedef iImplSpinor<Simd> SiteSpinor;
+ typedef iImplPropagator<Simd> SitePropagator;
  typedef iImplHalfSpinor<Simd> SiteHalfSpinor;
  typedef iImplDoubledGaugeField<Simd> SiteDoubledGaugeField;
  
  typedef Lattice<SiteSpinor> FermionField;
+ typedef Lattice<SitePropagator> PropagatorField;
  typedef Lattice<SiteDoubledGaugeField> DoubledGaugeField;
  
  typedef WilsonCompressor<SiteHalfSpinor, SiteSpinor> Compressor;

lib/serialisation/BaseIO.h

@@ -67,9 +67,8 @@ namespace Grid {
     return os;
   }
   
-  class Serializable {};
-  
   // static polymorphism implemented using CRTP idiom
+  class Serializable;
   
   // Static abstract writer
   template <typename T>
@@ -122,6 +121,27 @@ namespace Grid {
     T *upcast;
   };
   
+  // serializable base class
+  class Serializable
+  {
+  public:
+    template <typename T>
+    static inline void write(Writer<T> &WR,const std::string &s,
+                             const Serializable &obj)
+    {}
+    
+    template <typename T>
+    static inline void read(Reader<T> &RD,const std::string &s,
+                            Serializable &obj)
+    {}
+    
+    friend inline std::ostream & operator<<(std::ostream &os,
+                                            const Serializable &obj)
+    {
+      return os;
+    }
+  };
+  
   // Generic writer interface
   template <typename T>
   inline void push(Writer<T> &w, const std::string &s)

lib/stencil/Lebesgue.cc

@@ -65,7 +65,7 @@ void LebesgueOrder::CartesianBlocking(void)
 {
   _LebesgueReorder.resize(0);
 
-  std::cout << GridLogDebug << " CartesianBlocking ";
+  //    std::cout << GridLogDebug << " CartesianBlocking ";
   //    for(int d=0;d<Block.size();d++) std::cout <<Block[d]<<" ";
   //    std::cout<<std::endl; 
 

scripts/filelist

@@ -24,10 +24,11 @@ for subdir in $dirs; do
     echo "tests: ${TESTLIST} ${SUB}" > Make.inc
     echo ${PREF}_PROGRAMS = ${TESTLIST} >> Make.inc
     echo >> Make.inc
+    HADLINK=`[ $subdir = './hadrons' ] && echo '-lHadrons '`
     for f in $TESTS; do
 	   BNAME=`basename $f .cc`
 	   echo ${BNAME}_SOURCES=$f >> Make.inc
-	   echo ${BNAME}_LDADD=-lGrid>> Make.inc
+	   echo ${BNAME}_LDADD=${HADLINK}-lGrid  >> Make.inc
 	   echo >> Make.inc
     done
     if [ $subdir != '.' ]; then

tests/hadrons/Makefile.am

@@ -0,0 +1,3 @@
+AM_LDFLAGS += -L../../extras/hadrons
+
+include Make.inc

tests/hadrons/Test_hadrons_meson_3pt.cc

@@ -0,0 +1,170 @@
+/*******************************************************************************
+ Grid physics library, www.github.com/paboyle/Grid
+ 
+ Source file: tests/hadrons/Test_hadrons_meson_3pt.cc
+ 
+ Copyright (C) 2015
+ 
+ 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 <Grid/Hadrons/Application.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+
+int main(int argc, char *argv[])
+{
+    // initialization //////////////////////////////////////////////////////////
+    Grid_init(&argc, &argv);
+    HadronsLogError.Active(GridLogError.isActive());
+    HadronsLogWarning.Active(GridLogWarning.isActive());
+    HadronsLogMessage.Active(GridLogMessage.isActive());
+    HadronsLogIterative.Active(GridLogIterative.isActive());
+    HadronsLogDebug.Active(GridLogDebug.isActive());
+    LOG(Message) << "Grid initialized" << std::endl;
+    
+    // run setup ///////////////////////////////////////////////////////////////
+    Application              application;
+    std::vector<std::string> flavour = {"l", "s", "c1", "c2", "c3"};
+    std::vector<double>      mass    = {.01, .04, .2  , .25 , .3  };
+    unsigned int             nt      = GridDefaultLatt()[Tp];
+    
+    // global parameters
+    Application::GlobalPar globalPar;
+    globalPar.trajCounter.start    = 1500;
+    globalPar.trajCounter.end      = 1520;
+    globalPar.trajCounter.step     = 20;
+    globalPar.seed                 = "1 2 3 4";
+    globalPar.genetic.maxGen       = 1000;
+    globalPar.genetic.maxCstGen    = 200;
+    globalPar.genetic.popSize      = 20;
+    globalPar.genetic.mutationRate = .1;
+    application.setPar(globalPar);
+    
+    // gauge field
+    application.createModule<MGauge::Unit>("gauge");
+    for (unsigned int i = 0; i < flavour.size(); ++i)
+    {
+        // actions
+        MAction::DWF::Par actionPar;
+        actionPar.gauge = "gauge";
+        actionPar.Ls    = 12;
+        actionPar.M5    = 1.8;
+        actionPar.mass  = mass[i];
+        application.createModule<MAction::DWF>("DWF_" + flavour[i], actionPar);
+        
+        // solvers
+        MSolver::RBPrecCG::Par solverPar;
+        solverPar.action   = "DWF_" + flavour[i];
+        solverPar.residual = 1.0e-8;
+        application.createModule<MSolver::RBPrecCG>("CG_" + flavour[i],
+                                                    solverPar);
+    }
+    for (unsigned int t = 0; t < nt; t += 1)
+    {
+        std::string                           srcName;
+        std::vector<std::string>              qName;
+        std::vector<std::vector<std::string>> seqName;
+        
+        // Z2 source
+        MSource::Z2::Par z2Par;
+        z2Par.tA = t;
+        z2Par.tB = t;
+        srcName  = "z2_" + std::to_string(t);
+        application.createModule<MSource::Z2>(srcName, z2Par);
+        for (unsigned int i = 0; i < flavour.size(); ++i)
+        {
+            // sequential sources
+            MSource::SeqGamma::Par seqPar;
+            qName.push_back("QZ2_" + flavour[i] + "_" + std::to_string(t));
+            seqPar.q   = qName[i];
+            seqPar.tA  = (t + nt/4) % nt;
+            seqPar.tB  = (t + nt/4) % nt;
+            seqPar.mom = "1. 0. 0. 0.";
+            seqName.push_back(std::vector<std::string>(Nd));
+            for (unsigned int mu = 0; mu < Nd; ++mu)
+            {
+                seqPar.gamma   = 0x1 << mu;
+                seqName[i][mu] = "G" + std::to_string(seqPar.gamma)
+                                 + "_" + std::to_string(seqPar.tA) + "-"
+                                 + qName[i];
+                application.createModule<MSource::SeqGamma>(seqName[i][mu], seqPar);
+            }
+            
+            // propagators
+            Quark::Par quarkPar;
+            quarkPar.solver = "CG_" + flavour[i];
+            quarkPar.source = srcName;
+            application.createModule<Quark>(qName[i], quarkPar);
+            for (unsigned int mu = 0; mu < Nd; ++mu)
+            {
+                quarkPar.source = seqName[i][mu];
+                seqName[i][mu]  = "Q_" + flavour[i] + "-" + seqName[i][mu];
+                application.createModule<Quark>(seqName[i][mu], quarkPar);
+            }
+        }
+        
+        // contractions
+        MContraction::Meson::Par mesPar;
+        for (unsigned int i = 0; i < flavour.size(); ++i)
+        for (unsigned int j = i; j < flavour.size(); ++j)
+        {
+            mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j];
+            mesPar.q1     = qName[i];
+            mesPar.q2     = qName[j];
+            application.createModule<MContraction::Meson>("meson_Z2_"
+                                                          + std::to_string(t)
+                                                          + "_"
+                                                          + flavour[i]
+                                                          + flavour[j],
+                                                          mesPar);
+        }
+        for (unsigned int i = 0; i < flavour.size(); ++i)
+        for (unsigned int j = 0; j < flavour.size(); ++j)
+        for (unsigned int mu = 0; mu < Nd; ++mu)
+        {
+            MContraction::Meson::Par mesPar;
+            
+            mesPar.output = "3pt/Z2_" + flavour[i] + flavour[j] + "_"
+                            + std::to_string(mu);
+            mesPar.q1     = qName[i];
+            mesPar.q2     = seqName[j][mu];
+            application.createModule<MContraction::Meson>("3pt_Z2_"
+                                                          + std::to_string(t)
+                                                          + "_"
+                                                          + flavour[i]
+                                                          + flavour[j]
+                                                          + "_"
+                                                          + std::to_string(mu),
+                                                          mesPar);
+        }
+    }
+    
+    // execution
+    application.saveParameterFile("meson3pt.xml");
+    application.run();
+    
+    // epilogue
+    LOG(Message) << "Grid is finalizing now" << std::endl;
+    Grid_finalize();
+    
+    return EXIT_SUCCESS;
+}

tests/hadrons/Test_hadrons_spectrum.cc

@@ -0,0 +1,132 @@
+/*******************************************************************************
+ Grid physics library, www.github.com/paboyle/Grid
+ 
+ Source file: tests/hadrons/Test_hadrons_spectrum.cc
+ 
+ Copyright (C) 2015
+ 
+ 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 <Grid/Hadrons/Application.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+
+int main(int argc, char *argv[])
+{
+    // initialization //////////////////////////////////////////////////////////
+    Grid_init(&argc, &argv);
+    HadronsLogError.Active(GridLogError.isActive());
+    HadronsLogWarning.Active(GridLogWarning.isActive());
+    HadronsLogMessage.Active(GridLogMessage.isActive());
+    HadronsLogIterative.Active(GridLogIterative.isActive());
+    HadronsLogDebug.Active(GridLogDebug.isActive());
+    LOG(Message) << "Grid initialized" << std::endl;
+    
+    // run setup ///////////////////////////////////////////////////////////////
+    Application              application;
+    std::vector<std::string> flavour = {"l", "s", "c1", "c2", "c3"};
+    std::vector<double>      mass    = {.01, .04, .2  , .25 , .3  };
+    
+    // global parameters
+    Application::GlobalPar globalPar;
+    globalPar.trajCounter.start = 1500;
+    globalPar.trajCounter.end   = 1520;
+    globalPar.trajCounter.step  = 20;
+    globalPar.seed              = "1 2 3 4";
+    application.setPar(globalPar);
+    // gauge field
+    application.createModule<MGauge::Unit>("gauge");
+    // sources
+    MSource::Z2::Par z2Par;
+    z2Par.tA = 0;
+    z2Par.tB = 0;
+    application.createModule<MSource::Z2>("z2", z2Par);
+    MSource::Point::Par ptPar;
+    ptPar.position = "0 0 0 0";
+    application.createModule<MSource::Point>("pt", ptPar);
+    for (unsigned int i = 0; i < flavour.size(); ++i)
+    {
+        // actions
+        MAction::DWF::Par actionPar;
+        actionPar.gauge = "gauge";
+        actionPar.Ls    = 12;
+        actionPar.M5    = 1.8;
+        actionPar.mass  = mass[i];
+        application.createModule<MAction::DWF>("DWF_" + flavour[i], actionPar);
+        
+        // solvers
+        MSolver::RBPrecCG::Par solverPar;
+        solverPar.action   = "DWF_" + flavour[i];
+        solverPar.residual = 1.0e-8;
+        application.createModule<MSolver::RBPrecCG>("CG_" + flavour[i],
+                                                    solverPar);
+        
+        // propagators
+        Quark::Par quarkPar;
+        quarkPar.solver = "CG_" + flavour[i];
+        quarkPar.source = "pt";
+        application.createModule<Quark>("Qpt_" + flavour[i], quarkPar);
+        quarkPar.source = "z2";
+        application.createModule<Quark>("QZ2_" + flavour[i], quarkPar);
+    }
+    for (unsigned int i = 0; i < flavour.size(); ++i)
+    for (unsigned int j = i; j < flavour.size(); ++j)
+    {
+        MContraction::Meson::Par mesPar;
+        
+        mesPar.output = "mesons/pt_" + flavour[i] + flavour[j];
+        mesPar.q1     = "Qpt_" + flavour[i];
+        mesPar.q2     = "Qpt_" + flavour[j];
+        application.createModule<MContraction::Meson>("meson_pt_"
+                                                      + flavour[i] + flavour[j],
+                                                      mesPar);
+        mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j];
+        mesPar.q1     = "QZ2_" + flavour[i];
+        mesPar.q2     = "QZ2_" + flavour[j];
+        application.createModule<MContraction::Meson>("meson_Z2_"
+                                                      + flavour[i] + flavour[j],
+                                                      mesPar);
+    }
+    for (unsigned int i = 0; i < flavour.size(); ++i)
+    for (unsigned int j = i; j < flavour.size(); ++j)
+    for (unsigned int k = j; k < flavour.size(); ++k)
+    {
+        MContraction::Baryon::Par barPar;
+        
+        barPar.output = "baryons/pt_" + flavour[i] + flavour[j] + flavour[k];
+        barPar.q1     = "Qpt_" + flavour[i];
+        barPar.q2     = "Qpt_" + flavour[j];
+        barPar.q3     = "Qpt_" + flavour[k];
+        application.createModule<MContraction::Baryon>(
+            "baryon_pt_" + flavour[i] + flavour[j] + flavour[k], barPar);
+    }
+    
+    // execution
+    application.saveParameterFile("spectrum.xml");
+    application.run();
+    
+    // epilogue
+    LOG(Message) << "Grid is finalizing now" << std::endl;
+    Grid_finalize();
+    
+    return EXIT_SUCCESS;
+}