Merge branch 'hadrons' into feature/hadrons

Makefile.am

@@ -1,5 +1,5 @@
 # additional include paths necessary to compile the C++ library
 AM_CXXFLAGS = -I$(top_srcdir)/
-SUBDIRS = lib tests benchmarks
+SUBDIRS = lib tests benchmarks programs
 
 filelist: $(SUBDIRS)

configure.ac

@@ -287,6 +287,8 @@ AC_CONFIG_FILES(lib/Makefile)
 AC_CONFIG_FILES(tests/Makefile)
 AC_CONFIG_FILES(tests/qdpxx/Makefile)
 AC_CONFIG_FILES(benchmarks/Makefile)
+AC_CONFIG_FILES(programs/Makefile)
+AC_CONFIG_FILES(programs/Hadrons/Makefile)
 AC_OUTPUT
 
 

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)

programs/Hadrons/Application.cc

@@ -0,0 +1,195 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/Application.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 <Hadrons/Application.hpp>
+#include <Hadrons/GeneticScheduler.hpp>
+
+using namespace Grid;
+using namespace QCD;
+using namespace Hadrons;
+
+#define BIG_SEP "==============="
+#define SEP     "---------------"
+
+/******************************************************************************
+ *                       Application implementation                           *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+Application::Application(const std::string parameterFileName)
+: parameterFileName_(parameterFileName)
+, 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;
+}
+
+// destructor //////////////////////////////////////////////////////////////////
+Application::~Application(void)
+{}
+
+// execute /////////////////////////////////////////////////////////////////////
+void Application::run(void)
+{
+    parseParameterFile();
+    schedule();
+    configLoop();
+}
+
+// parse parameter file ////////////////////////////////////////////////////////
+class ObjectId: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(ObjectId,
+                                    std::string, name,
+                                    std::string, type);
+};
+
+void Application::parseParameterFile(void)
+{
+    XmlReader reader(parameterFileName_);
+    ObjectId  id;
+    
+    LOG(Message) << "Reading '" << parameterFileName_ << "'..." << std::endl;
+    read(reader, "parameters", 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);
+    env_.setSeed(strToVec<int>(par_.seed));
+    env_.printContent();
+}
+
+// schedule computation ////////////////////////////////////////////////////////
+#define MEM_MSG(size)\
+sizeString((size)*locVol_) << " (" << sizeString(size)  << "/site)"
+
+void Application::schedule(void)
+{
+    // memory peak function
+    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;
+    };
+    
+    // constrained topological sort using a genetic algorithm
+    LOG(Message) << "Scheduling computation..." << std::endl;
+    constexpr unsigned int maxGen = 200, maxCstGen = 50;
+    unsigned int           k = 0, gen, prevPeak, nCstPeak = 0;
+    auto                   graph = env_.makeModuleGraph();
+    auto                   con = graph.getConnectedComponents();
+    std::random_device     rd;
+    GeneticScheduler<unsigned int>::Parameters par;
+
+    par.popSize      = 20;
+    par.mutationRate = .1;
+    par.seed         = rd();
+    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
+        {
+            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 < maxGen) and (nCstPeak < maxCstGen));
+        auto &t = scheduler.getMinSchedule();
+        LOG(Message) << "Program " << i + 1 << " (memory peak: "
+                     << MEM_MSG(scheduler.getMinValue()) << "):" << std::endl;
+        for (unsigned int j = 0; j < t.size(); ++j)
+        {
+            program_.push_back(t[j]);
+            LOG(Message) << std::setw(4) << std::right << k + 1 << ": "
+                         << env_.getModuleName(program_[k]) << std::endl;
+            k++;
+        }
+    }
+}
+
+// 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;
+}

programs/Hadrons/Application.hpp

@@ -0,0 +1,85 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/Application.hpp
+
+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.
+*******************************************************************************/
+
+#ifndef Hadrons_Application_hpp_
+#define Hadrons_Application_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Environment.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+class TrajRange: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(TrajRange,
+                                    unsigned int, start,
+                                    unsigned int, end,
+                                    unsigned int, step);
+};
+
+class GlobalPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(GlobalPar,
+                                    TrajRange,   trajCounter,
+                                    std::string, seed);
+};
+
+/******************************************************************************
+ *                         Main program manager                               *
+ ******************************************************************************/
+class Application
+{
+public:
+
+public:
+    // constructor
+    Application(const std::string parameterFileName);
+    // destructor
+    virtual ~Application(void);
+    // execute
+    void run(void);
+private:
+    // parse parameter file
+    void parseParameterFile(void);
+    // schedule computation
+    void schedule(void);
+    // loop on configurations
+    void configLoop(void);
+private:
+    long unsigned int                               locVol_;
+    std::string                                     parameterFileName_;
+    GlobalPar                                       par_;
+    Environment                                     &env_;
+    std::vector<unsigned int>                       program_;
+};
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Application_hpp_

programs/Hadrons/Environment.cc

@@ -0,0 +1,684 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/Environment.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 <Hadrons/Environment.hpp>
+#include <Hadrons/Module.hpp>
+#include <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::createModule(const std::string name, const std::string type,
+                               XmlReader &reader)
+{
+    if (!hasModule(name))
+    {
+        auto                      &factory = ModuleFactory::getInstance();
+        std::vector<unsigned int> inputAddress;
+        ModuleInfo                m;
+
+        m.data = factory.create(type, name);
+        m.type = typeIdPt(*m.data.get());
+        m.name = name;
+        m.data->parseParameters(reader, "options");
+        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));
+        moduleAddress_[name] = module_.size() - 1;
+        for (auto &out: output)
+        {
+            if (!hasObject(out))
+            {
+                addObject(out , module_.size() - 1);
+            }
+            else
+            {
+                if (object_[objectAddress_[out]].module < 0)
+                {
+                    object_[objectAddress_[out]].module = module_.size() - 1;
+                }
+                else
+                {
+                    HADRON_ERROR("object '" + out
+                        + "' is already produced by module '"
+                        + module_[object_[getObjectAddress(out)].module].name
+                        + "' (while creating module '" + name + "')");
+                }
+            }
+        }
+    }
+    else
+    {
+        HADRON_ERROR("module '" + name + "' already exists");
+    }
+}
+
+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 module_[address].type->name();
+    }
+    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)"
+
+unsigned int Environment::executeProgram(const std::vector<unsigned int> &p)
+{
+    unsigned int                        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;
+}
+
+unsigned int 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] = 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 object_[address].type->name();
+    }
+    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));
+}
+
+unsigned int 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));
+    }
+}
+
+unsigned int 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);
+}
+
+long unsigned int Environment::getTotalSize(void) const
+{
+    long unsigned int 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())
+        {
+            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) << std::right << i << ": "
+                     << getModuleName(i) << " ("
+                     << getModuleType(i) << ")" << std::endl;
+    }
+    LOG(Message) << "Objects: " << std::endl;
+    for (unsigned int i = 0; i < object_.size(); ++i)
+    {
+        LOG(Message) << std::setw(4) << std::right << i << ": "
+                     << getObjectName(i) << std::endl;
+    }
+}

programs/Hadrons/Environment.hpp

@@ -0,0 +1,351 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/Environment.hpp
+
+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.
+*******************************************************************************/
+
+#ifndef Hadrons_Environment_hpp_
+#define Hadrons_Environment_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Graph.hpp>
+
+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 std::unique_ptr<ModuleBase>                 ModPt;
+    typedef std::unique_ptr<GridCartesian>              GridPt;
+    typedef std::unique_ptr<GridRedBlackCartesian>      GridRbPt;
+    typedef FermionOperator<WilsonImplR>                FMat;
+    typedef std::unique_ptr<FMat>                       FMatPt;
+    typedef std::function<void(LatticeFermion &,
+                               const LatticeFermion &)> Solver;
+    typedef std::unique_ptr<GridParallelRNG>            RngPt;
+    typedef std::unique_ptr<LatticeBase>                LatticePt;
+private:
+    struct ModuleInfo
+    {
+        const std::type_info        *type{nullptr};
+        std::string                 name;
+        std::unique_ptr<ModuleBase> data{nullptr};
+        std::vector<unsigned int>   input;
+    };
+    struct ObjInfo
+    {
+        unsigned int            size{0}, 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                    createModule(const std::string name,
+                                         const std::string type,
+                                         XmlReader &reader);
+    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;
+    unsigned int            executeProgram(const std::vector<unsigned int> &p);
+    unsigned int            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;
+    unsigned int            getObjectSize(const unsigned int address) const;
+    unsigned int            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;
+    long unsigned int       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_;
+    // fermion matrix store
+    std::map<unsigned int, FMatPt>         fMat_;
+    // solver store & solver/action map
+    std::map<unsigned int, Solver>         solver_;
+    std::map<std::string, std::string>     solverAction_;
+    // object store
+    std::vector<ObjInfo>                   object_;
+    std::map<std::string, unsigned int>    objectAddress_;
+};
+
+/******************************************************************************
+ *                        template implementation                             *
+ ******************************************************************************/
+template <typename T>
+Holder<T>::Holder(T *pt)
+: objPt_(pt)
+{}
+
+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);
+}
+
+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_

programs/Hadrons/Factory.hpp

@@ -0,0 +1,104 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/Factory.hpp
+
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution 
+directory.
+*******************************************************************************/
+
+#ifndef Hadrons_Factory_hpp_
+#define Hadrons_Factory_hpp_
+
+#include <Hadrons/Global.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                        abstract factory class                              *
+ ******************************************************************************/
+template <typename T>
+class Factory
+{
+public:
+    typedef std::function<std::unique_ptr<T>(const std::string)> Func;
+public:
+    // constructor
+    Factory(void) = default;
+    // destructor
+    virtual ~Factory(void) = default;
+    // registration
+    void registerBuilder(const std::string type, const Func &f);
+    // get 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_

programs/Hadrons/GeneticScheduler.hpp

@@ -0,0 +1,236 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/GeneticScheduler.hpp
+
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution 
+directory.
+*******************************************************************************/
+
+#ifndef Hadrons_GeneticScheduler_hpp_
+#define Hadrons_GeneticScheduler_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Graph.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                   Scheduler based on a genetic algorithm                   *
+ ******************************************************************************/
+template <typename T>
+class GeneticScheduler
+{
+public:
+    typedef std::function<int(const std::vector<T> &)> 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 std::vector<T> & getMinSchedule(void);
+    int                    getMinValue(void);
+    // breed a new generation
+    void nextGeneration(void);
+    // print population
+    friend std::ostream & operator<<(std::ostream &out,
+                                     const GeneticScheduler<T> &s)
+    {
+        for (auto &p: s.population_)
+        {
+            out << p.second << ": " << p.first << std::endl;
+        }
+        
+        return out;
+    }
+private:
+    // randomly initialize population
+    void initPopulation(void);
+    // genetic operators
+    const std::vector<T> & selection(void);
+    void                   crossover(const std::vector<T> &c1,
+                                     const std::vector<T> &c2);
+    void                   mutation(std::vector<T> &c);
+private:
+    Graph<T>                           &graph_;
+    const ObjFunc                      &func_;
+    const Parameters                   par_;
+    std::multimap<int, std::vector<T>> 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 std::vector<T> & 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)
+{
+    std::uniform_real_distribution<double> dis(0., 1.);
+    
+    // random initialization of the population if necessary
+    if (population_.size() != par_.popSize)
+    {
+        initPopulation();
+    }
+
+    // mating
+    for (unsigned int i = 0; i < par_.popSize/2; ++i)
+    {
+        auto &p1 = selection(), &p2 = selection();
+        crossover(p1, p2);
+    }
+    
+    // random mutations
+    auto buf = population_;
+    population_.clear();
+    for (auto &c: buf)
+    {
+        if (dis(gen_) < par_.mutationRate)
+        {
+            mutation(c.second);
+        }
+        population_.emplace(func_(c.second), c.second);
+    }
+    
+    // grim reaper
+    auto it = population_.begin();
+    
+    std::advance(it, par_.popSize);
+    population_.erase(it, population_.end());
+}
+
+// randomly initialize population //////////////////////////////////////////////
+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);
+    }
+}
+
+// genetic operators ///////////////////////////////////////////////////////////
+template <typename T>
+const std::vector<T> & GeneticScheduler<T>::selection(void)
+{
+    std::vector<double> prob;
+    
+    for (auto &c: population_)
+    {
+        prob.push_back(1./c.first);
+    }
+    std::discrete_distribution<unsigned int> dis(prob.begin(), prob.end());
+    auto rIt = population_.begin();
+    std::advance(rIt, dis(gen_));
+    
+    return rIt->second;
+}
+
+template <typename T>
+void GeneticScheduler<T>::crossover(const std::vector<T> &p1,
+                                    const std::vector<T> &p2)
+{
+    std::uniform_int_distribution<unsigned int> dis(0, p1.size() - 1);
+    unsigned int                                cut = dis(gen_);
+    std::vector<T>                              c1, c2, buf;
+    
+    auto cross = [&buf, cut](std::vector<T> &c, const std::vector<T> &p1,
+                             const std::vector<T> &p2)
+    {
+        buf = p2;
+        for (unsigned int i = 0; i < cut; ++i)
+        {
+            c.push_back(p1[i]);
+            buf.erase(std::find(buf.begin(), buf.end(), p1[i]));
+        }
+        for (unsigned int i = 0; i < buf.size(); ++i)
+        {
+            c.push_back(buf[i]);
+        }
+    };
+    
+    cross(c1, p1, p2);
+    cross(c2, p2, p1);
+    population_.emplace(func_(c1), c1);
+    population_.emplace(func_(c2), c2);
+}
+
+template <typename T>
+void GeneticScheduler<T>::mutation(std::vector<T> &c)
+{
+    std::uniform_int_distribution<unsigned int> dis(0, c.size() - 1);
+    unsigned int                                cut = dis(gen_);
+    Graph<T>                                    g = graph_;
+    std::vector<T>                              buf;
+    
+    for (unsigned int i = cut; i < c.size(); ++i)
+    {
+        g.removeVertex(c[i]);
+    }
+    if (g.size() > 0)
+    {
+        buf = g.topoSort(gen_);
+    }
+    for (unsigned int i = cut; i < c.size(); ++i)
+    {
+        buf.push_back(c[i]);
+    }
+    c = buf;
+}
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_GeneticScheduler_hpp_

programs/Hadrons/Global.cc

@@ -0,0 +1,65 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/Global.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 <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);
+}

programs/Hadrons/Global.hpp

@@ -0,0 +1,128 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/Global.hpp
+
+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.
+*******************************************************************************/
+
+#ifndef Hadrons_Global_hpp_
+#define Hadrons_Global_hpp_
+
+#include <set>
+#include <stack>
+#include <cxxabi.h>
+#include <Grid.h>
+
+#define BEGIN_HADRONS_NAMESPACE \
+namespace Grid {\
+using namespace QCD;\
+namespace Hadrons {\
+using Grid::operator<<;
+#define END_HADRONS_NAMESPACE }}
+
+/* the 'using Grid::operator<<;' statement prevents a very nasty compilation
+ * error with GCC (clang compiles fine without it).
+ */
+
+BEGIN_HADRONS_NAMESPACE
+
+class HadronsLogger: public Logger
+{
+public:
+    HadronsLogger(int on, std::string nm): Logger("Hadrons", on, nm,
+                                                  Logger::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);
+
+template <typename T>
+std::string typeName(const T &x)
+{
+    std::string name(typeid(x).name());
+
+    return name;
+}
+
+template <typename T>
+std::string typeName(void)
+{
+    std::string name(typeid(T).name());
+
+    return name;
+}
+
+template <typename T>
+const std::type_info * typeIdPt(const T &x)
+{
+    return &typeid(x);
+}
+
+template <typename T>
+const std::type_info * typeName(void)
+{
+    return &typeid(T);
+}
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Global_hpp_

programs/Hadrons/Graph.hpp

@@ -0,0 +1,758 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/Graph.hpp
+
+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.
+*******************************************************************************/
+
+#ifndef Hadrons_Graph_hpp_
+#define Hadrons_Graph_hpp_
+
+#include <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_

programs/Hadrons/Hadrons.cc

@@ -0,0 +1,66 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/Hadrons.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 <Hadrons/Application.hpp>
+
+using namespace Grid;
+using namespace QCD;
+using namespace Hadrons;
+
+int main(int argc, char *argv[])
+{
+    // parse command line
+    std::string parameterFileName;
+    
+    if (argc < 2)
+    {
+        std::cerr << "usage: " << argv[0] << " <parameter file> [Grid options]";
+        std::cerr << std::endl;
+        std::exit(EXIT_FAILURE);
+    }
+    parameterFileName = argv[1];
+    
+    // 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.run();
+    
+    // epilogue
+    LOG(Message) << "Grid is finalizing now" << std::endl;
+    Grid_finalize();
+    
+    return EXIT_SUCCESS;
+}

programs/Hadrons/Make.inc

@@ -0,0 +1,4 @@
+bin_PROGRAMS = Hadrons
+
+Hadrons_SOURCES = Hadrons.cc
+Hadrons_LDADD   = -lGrid

programs/Hadrons/Makefile.am

@@ -0,0 +1,16 @@
+AM_CXXFLAGS = -I$(top_srcdir)/programs -I../$(top_srcdir)/programs -I$(top_srcdir)/lib
+AM_LDFLAGS  = -L$(top_builddir)/lib
+
+bin_PROGRAMS = Hadrons
+
+include modules.inc
+
+Hadrons_SOURCES =   \
+    $(modules)      \
+    Application.cc  \
+    Environment.cc  \
+    Global.cc       \
+    Hadrons.cc      \
+    Module.cc
+
+Hadrons_LDADD   = -lGrid

programs/Hadrons/Module.cc

@@ -0,0 +1,62 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/Module.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 <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_;
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+void ModuleBase::operator()(void)
+{
+    setup();
+    if (!env().isDryRun())
+    {
+        execute();
+    }
+}

programs/Hadrons/Module.cc.template

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

programs/Hadrons/Module.hpp

@@ -0,0 +1,148 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/Module.hpp
+
+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.
+*******************************************************************************/
+
+#ifndef Hadrons_Module_hpp_
+#define Hadrons_Module_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Environment.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+// module registration macro
+#define MODULE_REGISTER(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;
+
+/******************************************************************************
+ *                            Module class                                    *
+ ******************************************************************************/
+// base class
+class ModuleBase
+{
+public:
+    // convenient type shortcuts
+    typedef Environment::FMat   FMat;
+    typedef Environment::Solver Solver;
+public:
+    // constructor
+    ModuleBase(const std::string name);
+    // destructor
+    virtual ~ModuleBase(void) = default;
+    // access
+    std::string getName(void) const;
+    Environment &env(void) const;
+    // 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;
+    // 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);
+    // 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) {};
+};
+
+/******************************************************************************
+ *                           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_);read(reader, 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_

programs/Hadrons/Module.hpp.template

@@ -0,0 +1,40 @@
+#ifndef Hadrons____FILEBASENAME____hpp_
+#define Hadrons____FILEBASENAME____hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         ___FILEBASENAME___                                 *
+ ******************************************************************************/
+class ___FILEBASENAME___Par: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(___FILEBASENAME___Par, 
+                                    unsigned int, i);
+};
+
+class ___FILEBASENAME___: public Module<___FILEBASENAME___Par>
+{
+public:
+    // constructor
+    ___FILEBASENAME___(const std::string name);
+    // destructor
+    virtual ~___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___);
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons____FILEBASENAME____hpp_

programs/Hadrons/ModuleFactory.hpp

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

programs/Hadrons/Modules/ADWF.cc

@@ -0,0 +1,81 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/ADWF.cc
+
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution 
+directory.
+*******************************************************************************/
+
+#include <Hadrons/Modules/ADWF.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+
+/******************************************************************************
+*                          ADWF implementation                                *
+******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+ADWF::ADWF(const std::string name)
+: Module<ADWFPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+std::vector<std::string> ADWF::getInput(void)
+{
+    std::vector<std::string> in = {par().gauge};
+    
+    return in;
+}
+
+std::vector<std::string> ADWF::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+void ADWF::setup(void)
+{
+    unsigned int size;
+    
+    size = 3*env().lattice4dSize<WilsonFermionR::DoubledGaugeField>();
+    env().registerObject(getName(), size, par().Ls);
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+void ADWF::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().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 DomainWallFermionR(U, g5, grb5, g4, grb4, par().mass,
+                                          par().M5);
+    env().setObject(getName(), fMatPt);
+}

programs/Hadrons/Modules/ADWF.hpp

@@ -0,0 +1,70 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/ADWF.hpp
+
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution 
+directory.
+*******************************************************************************/
+
+#ifndef Hadrons_ADWF_hpp_
+#define Hadrons_ADWF_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                     Domain wall quark action                               *
+ ******************************************************************************/
+class ADWFPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(ADWFPar, 
+                                    std::string, gauge,
+                                    unsigned int, Ls,
+                                    double      , mass,
+                                    double      , M5);
+};
+
+class ADWF: public Module<ADWFPar>
+{
+public:
+    // constructor
+    ADWF(const std::string name);
+    // destructor
+    virtual ~ADWF(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(ADWF);
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_ADWF_hpp_

programs/Hadrons/Modules/AWilson.cc

@@ -0,0 +1,75 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/AWilson.cc
+
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution 
+directory.
+*******************************************************************************/
+
+#include <Hadrons/Modules/AWilson.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+
+/******************************************************************************
+*                         AWilson implementation                              *
+******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+AWilson::AWilson(const std::string name)
+: Module<AWilsonPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+std::vector<std::string> AWilson::getInput(void)
+{
+    std::vector<std::string> in = {par().gauge};
+    
+    return in;
+}
+
+std::vector<std::string> AWilson::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+void AWilson::setup(void)
+{
+    unsigned int size;
+    
+    size = 3*env().lattice4dSize<WilsonFermionR::DoubledGaugeField>();
+    env().registerObject(getName(), size);
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+void AWilson::execute()
+{
+    LOG(Message) << "Setting up Wilson fermion matrix with m= " << par().mass
+                 << " using gauge field '" << par().gauge << "'" << std::endl;
+    auto &U      = *env().getObject<LatticeGaugeField>(par().gauge);
+    auto &grid   = *env().getGrid();
+    auto &gridRb = *env().getRbGrid();
+    FMat *fMatPt = new WilsonFermionR(U, grid, gridRb, par().mass);
+    env().setObject(getName(), fMatPt);
+}

programs/Hadrons/Modules/AWilson.hpp

@@ -0,0 +1,68 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/AWilson.hpp
+
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution 
+directory.
+*******************************************************************************/
+
+#ifndef Hadrons_AWilson_hpp_
+#define Hadrons_AWilson_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                            Wilson quark action                             *
+ ******************************************************************************/
+class AWilsonPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(AWilsonPar,
+                                    std::string, gauge,
+                                    double     , mass);
+};
+
+class AWilson: public Module<AWilsonPar>
+{
+public:
+    // constructor
+    AWilson(const std::string name);
+    // destructor
+    virtual ~AWilson(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(AWilson);
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_AWilson_hpp_

programs/Hadrons/Modules/CMeson.cc

@@ -0,0 +1,94 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/CMeson.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 <Hadrons/Modules/CMeson.hpp>
+
+using namespace Grid;
+using namespace QCD;
+using namespace Hadrons;
+
+/******************************************************************************
+ *                          CMeson implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+CMeson::CMeson(const std::string name)
+: Module<CMesonPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+std::vector<std::string> CMeson::getInput(void)
+{
+    std::vector<std::string> input = {par().q1, par().q2};
+    
+    return input;
+}
+
+std::vector<std::string> CMeson::getOutput(void)
+{
+    std::vector<std::string> output = {getName()};
+    
+    return output;
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+void CMeson::execute(void)
+{
+    LOG(Message) << "Computing meson contraction '" << getName() << "' using"
+                 << " quarks '" << par().q1 << "' and '" << par().q2 << "'"
+                 << std::endl;
+    
+    XmlWriter             writer(par().output);
+    LatticePropagator     &q1 = *env().getObject<LatticePropagator>(par().q1);
+    LatticePropagator     &q2 = *env().getObject<LatticePropagator>(par().q2);
+    LatticeComplex        c(env().getGrid());
+    SpinMatrix            g[Ns*Ns], g5;
+    std::vector<TComplex> buf;
+    Result                result;
+    unsigned int          nt = env().getGrid()->GlobalDimensions()[Tp];
+    
+    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);
+}

programs/Hadrons/Modules/CMeson.hpp

@@ -0,0 +1,74 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/CMeson.hpp
+
+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.
+*******************************************************************************/
+
+#ifndef Hadrons_CMeson_hpp_
+#define Hadrons_CMeson_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                               CMeson                                       *
+ ******************************************************************************/
+class CMesonPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(CMesonPar,
+                                    std::string, q1,
+                                    std::string, q2,
+                                    std::string, output);
+};
+
+class CMeson: public Module<CMesonPar>
+{
+public:
+    class Result: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
+                        std::vector<std::vector<std::vector<Complex>>>, corr);
+    };
+public:
+    // constructor
+    CMeson(const std::string name);
+    // destructor
+    virtual ~CMeson(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(CMeson);
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_CMeson_hpp_

programs/Hadrons/Modules/GLoad.cc

@@ -0,0 +1,75 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/GLoad.cc
+
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution 
+directory.
+*******************************************************************************/
+
+#include <Hadrons/Modules/GLoad.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+
+/******************************************************************************
+*                          GLoad implementation                               *
+******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+GLoad::GLoad(const std::string name)
+: Module<GLoadPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+std::vector<std::string> GLoad::getInput(void)
+{
+    std::vector<std::string> in;
+    
+    return in;
+}
+
+std::vector<std::string> GLoad::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+void GLoad::setup(void)
+{
+    env().registerLattice<LatticeGaugeField>(getName());
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+void GLoad::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));
+}

programs/Hadrons/Modules/GLoad.hpp

@@ -0,0 +1,67 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/GLoad.hpp
+
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution 
+directory.
+*******************************************************************************/
+
+#ifndef Hadrons_GLoad_hpp_
+#define Hadrons_GLoad_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         Load a NERSC configuration                         *
+ ******************************************************************************/
+class GLoadPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(GLoadPar,
+                                    std::string, file);
+};
+
+class GLoad: public Module<GLoadPar>
+{
+public:
+    // constructor
+    GLoad(const std::string name);
+    // destructor
+    virtual ~GLoad(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(GLoad);
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_GLoad_hpp_

programs/Hadrons/Modules/GRandom.cc

@@ -0,0 +1,66 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/GRandom.cc
+
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution 
+directory.
+*******************************************************************************/
+
+#include <Hadrons/Modules/GRandom.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+
+/******************************************************************************
+*                          GRandom implementation                             *
+******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+GRandom::GRandom(const std::string name)
+: Module<NoPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+std::vector<std::string> GRandom::getInput(void)
+{
+    return std::vector<std::string>();
+}
+
+std::vector<std::string> GRandom::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+void GRandom::setup(void)
+{
+    env().registerLattice<LatticeGaugeField>(getName());
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+void GRandom::execute(void)
+{
+    LOG(Message) << "Generating random gauge configuration" << std::endl;
+    LatticeGaugeField &U = *env().createLattice<LatticeGaugeField>(getName());
+    SU3::HotConfiguration(*env().get4dRng(), U);
+}

programs/Hadrons/Modules/GRandom.hpp

@@ -0,0 +1,60 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/GRandom.hpp
+
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution 
+directory.
+*******************************************************************************/
+
+#ifndef Hadrons_GRandom_hpp_
+#define Hadrons_GRandom_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                             Random gauge                                   *
+ ******************************************************************************/
+class GRandom: public Module<NoPar>
+{
+public:
+    // constructor
+    GRandom(const std::string name);
+    // destructor
+    virtual ~GRandom(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(GRandom);
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_GRandom_hpp_

programs/Hadrons/Modules/GUnit.cc

@@ -0,0 +1,66 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/GUnit.cc
+
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution 
+directory.
+*******************************************************************************/
+
+#include <Hadrons/Modules/GUnit.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+
+/******************************************************************************
+*                            GUnit implementation                             *
+******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+GUnit::GUnit(const std::string name)
+: Module<NoPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+std::vector<std::string> GUnit::getInput(void)
+{
+    return std::vector<std::string>();
+}
+
+std::vector<std::string> GUnit::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+void GUnit::setup(void)
+{
+    env().registerLattice<LatticeGaugeField>(getName());
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+void GUnit::execute(void)
+{
+    LOG(Message) << "Creating unit gauge configuration" << std::endl;
+    LatticeGaugeField &U = *env().createLattice<LatticeGaugeField>(getName());
+    SU3::ColdConfiguration(*env().get4dRng(), U);
+}

programs/Hadrons/Modules/GUnit.hpp

@@ -0,0 +1,60 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/GUnit.hpp
+
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution 
+directory.
+*******************************************************************************/
+
+#ifndef Hadrons_GUnit_hpp_
+#define Hadrons_GUnit_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                              Unit gauge                                    *
+ ******************************************************************************/
+class GUnit: public Module<NoPar>
+{
+public:
+    // constructor
+    GUnit(const std::string name);
+    // destructor
+    virtual ~GUnit(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(GUnit);
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_GUnit_hpp_

programs/Hadrons/Modules/MQuark.cc

@@ -0,0 +1,135 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/MQuark.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 <Hadrons/Modules/MQuark.hpp>
+
+using namespace Grid;
+using namespace QCD;
+using namespace Hadrons;
+
+/******************************************************************************
+ *                          MQuark implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+MQuark::MQuark(const std::string name)
+: Module(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+std::vector<std::string> MQuark::getInput(void)
+{
+    std::vector<std::string> in = {par().source, par().solver};
+    
+    return in;
+}
+
+std::vector<std::string> MQuark::getOutput(void)
+{
+    std::vector<std::string> out = {getName(), getName() + "_5d"};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+void MQuark::setup(void)
+{
+    Ls_ = env().getObjectLs(par().solver);
+    env().registerLattice<LatticePropagator>(getName());
+    if (Ls_ > 1)
+    {
+        env().registerLattice<LatticePropagator>(getName() + "_5d", Ls_);
+    }
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+void MQuark::execute(void)
+{
+    LatticeFermion source(env().getGrid(Ls_)), sol(env().getGrid(Ls_)),
+                   tmp(env().getGrid());
+    std::string    propName = (Ls_ == 1) ? getName() : (getName() + "_5d");
+    
+    LOG(Message) << "Computing quark propagator '" << getName() << "'"
+                 << std::endl;
+    LatticePropagator   &prop    = *env().createLattice<LatticePropagator>(propName);
+    LatticePropagator   &fullSrc = *env().getObject<LatticePropagator>(par().source);
+    Environment::Solver &solver  = *env().getObject<Environment::Solver>(par().solver);
+    if (Ls_ > 1)
+    {
+        env().createLattice<LatticePropagator>(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)
+        {
+            LatticePropagator &p4d = *env().getObject<LatticePropagator>(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);
+        }
+    }
+}

programs/Hadrons/Modules/MQuark.hpp

@@ -0,0 +1,71 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/MQuark.hpp
+
+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.
+*******************************************************************************/
+
+#ifndef Hadrons_MQuark_hpp_
+#define Hadrons_MQuark_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                               MQuark                                       *
+ ******************************************************************************/
+class MQuarkPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(MQuarkPar,
+                                    std::string, source,
+                                    std::string, solver);
+};
+
+class MQuark: public Module<MQuarkPar>
+{
+public:
+    // constructor
+    MQuark(const std::string name);
+    // destructor
+    virtual ~MQuark(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_;
+    Environment::Solver *solver_{nullptr};
+};
+
+MODULE_REGISTER(MQuark);
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MQuark_hpp_

programs/Hadrons/Modules/SolRBPrecCG.cc

@@ -0,0 +1,83 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/SolRBPrecCG.cc
+
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution 
+directory.
+*******************************************************************************/
+
+#include <Hadrons/Modules/SolRBPrecCG.hpp>
+
+using namespace Grid;
+using namespace QCD;
+using namespace Hadrons;
+
+/******************************************************************************
+*                       SolRBPrecCG implementation                            *
+******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+SolRBPrecCG::SolRBPrecCG(const std::string name)
+: Module(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+std::vector<std::string> SolRBPrecCG::getInput(void)
+{
+    std::vector<std::string> in = {par().action};
+    
+    return in;
+}
+
+std::vector<std::string> SolRBPrecCG::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+void SolRBPrecCG::setup(void)
+{
+    auto Ls = env().getObjectLs(par().action);
+    
+    env().registerObject(getName(), 0, Ls);
+    env().addOwnership(getName(), par().action);
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+void SolRBPrecCG::execute(void)
+{
+    auto &mat   = *(env().getObject<Environment::FMat>(par().action));
+    auto solver = [&mat, this](LatticeFermion &sol,
+                               const LatticeFermion &source)
+    {
+        ConjugateGradient<LatticeFermion>           cg(par().residual, 10000);
+        SchurRedBlackDiagMooeeSolve<LatticeFermion> 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 Environment::Solver(solver));
+}

programs/Hadrons/Modules/SolRBPrecCG.hpp

@@ -0,0 +1,68 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/SolRBPrecCG.hpp
+
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution 
+directory.
+*******************************************************************************/
+
+#ifndef Hadrons_SolRBPrecCG_hpp_
+#define Hadrons_SolRBPrecCG_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                     Schur red-black preconditioned CG                      *
+ ******************************************************************************/
+class SolRBPrecCGPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(SolRBPrecCGPar,
+                                    std::string, action,
+                                    double     , residual);
+};
+
+class SolRBPrecCG: public Module<SolRBPrecCGPar>
+{
+public:
+    // constructor
+    SolRBPrecCG(const std::string name);
+    // destructor
+    virtual ~SolRBPrecCG(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(SolRBPrecCG);
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_SolRBPrecCG_hpp_

programs/Hadrons/Modules/SrcPoint.cc

@@ -0,0 +1,74 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/SrcPoint.cc
+
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution 
+directory.
+*******************************************************************************/
+
+#include <Hadrons/Modules/SrcPoint.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+
+/******************************************************************************
+*                         SrcPoint implementation                             *
+******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+SrcPoint::SrcPoint(const std::string name)
+: Module<SrcPointPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+std::vector<std::string> SrcPoint::getInput(void)
+{
+    std::vector<std::string> in;
+    
+    return in;
+}
+
+std::vector<std::string> SrcPoint::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+void SrcPoint::setup(void)
+{
+    env().registerLattice<LatticePropagator>(getName());
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+void SrcPoint::execute(void)
+{
+    std::vector<int> position = strToVec<int>(par().position);
+    SpinColourMatrix id;
+    
+    LOG(Message) << "Creating point source at position [" << par().position
+                 << "]" << std::endl;
+    LatticePropagator &src = *env().createLattice<LatticePropagator>(getName());
+    id  = 1.;
+    src = zero;
+    pokeSite(id, src, position);
+}

programs/Hadrons/Modules/SrcPoint.hpp

@@ -0,0 +1,78 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/SrcPoint.hpp
+
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution 
+directory.
+*******************************************************************************/
+
+#ifndef Hadrons_SrcPoint_hpp_
+#define Hadrons_SrcPoint_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/*
+ 
+ Point source
+ ------------
+ * src_x = delta_x,o
+ 
+ * options: o
+ - position: space-separated integer sequence (e.g. "0 1 1 0")
+ 
+ */
+
+/******************************************************************************
+ *                                SrcPoint                                    *
+ ******************************************************************************/
+class SrcPointPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(SrcPointPar,
+                                    std::string, position);
+};
+
+class SrcPoint: public Module<SrcPointPar>
+{
+public:
+    // constructor
+    SrcPoint(const std::string name);
+    // destructor
+    virtual ~SrcPoint(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(SrcPoint);
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_SrcPoint_hpp_

programs/Hadrons/Modules/SrcZ2.cc

@@ -0,0 +1,87 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/SrcZ2.cc
+
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution 
+directory.
+*******************************************************************************/
+
+#include <Hadrons/Modules/SrcZ2.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+
+/******************************************************************************
+*                           SrcZ2 implementation                              *
+******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+SrcZ2::SrcZ2(const std::string name)
+: Module<SrcZ2Par>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+std::vector<std::string> SrcZ2::getInput(void)
+{
+    std::vector<std::string> in;
+    
+    return in;
+}
+
+std::vector<std::string> SrcZ2::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+void SrcZ2::setup(void)
+{
+    env().registerLattice<LatticePropagator>(getName());
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+void SrcZ2::execute(void)
+{
+    Lattice<iScalar<vInteger>> t(env().getGrid());
+    LatticeComplex             eta(env().getGrid());
+    LatticeFermion             phi(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;
+    }
+    LatticePropagator &src = *env().createLattice<LatticePropagator>(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;
+}

programs/Hadrons/Modules/SrcZ2.hpp

@@ -0,0 +1,80 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/SrcZ2.hpp
+
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution 
+directory.
+*******************************************************************************/
+
+#ifndef Hadrons_SrcZ2_hpp_
+#define Hadrons_SrcZ2_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/*
+ 
+ Z_2 stochastic source
+ -----------------------------
+ * src_x = eta_x * theta(x_3 - ta) * theta(tb - x_3)
+ 
+ * options:
+ - tA: begin timeslice (integer)
+ - tB: end timesilce (integer)
+ 
+ */
+ 
+/******************************************************************************
+ *                               SrcZ2                                        *
+ ******************************************************************************/
+class SrcZ2Par: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(SrcZ2Par,
+                                    unsigned int, tA,
+                                    unsigned int, tB);
+};
+
+class SrcZ2: public Module<SrcZ2Par>
+{
+public:
+    // constructor
+    SrcZ2(const std::string name);
+    // destructor
+    virtual ~SrcZ2(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(SrcZ2);
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_SrcZ2_hpp_

programs/Hadrons/add_module.sh

@@ -0,0 +1,15 @@
+#!/usr/bin/env bash
+
+if (( $# != 1 )); then
+    echo "usage: `basename $0` <module name>" 1>&2
+    exit 1
+fi
+NAME=$1
+
+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" Module.cc.template > Modules/${NAME}.cc
+sed "s/___FILEBASENAME___/${NAME}/g" Module.hpp.template > Modules/${NAME}.hpp
+./make_module_list.sh

programs/Hadrons/make_module_list.sh

@@ -0,0 +1,4 @@
+#!/usr/bin/env bash
+
+echo 'modules =\' > modules.inc
+find Modules -name '*.cc' -type f -print | sed 's/^/  /;$q;s/$/ \\/' >> modules.inc

programs/Hadrons/modules.inc

@@ -0,0 +1,11 @@
+modules =\
+  Modules/ADWF.cc \
+  Modules/AWilson.cc \
+  Modules/CMeson.cc \
+  Modules/GLoad.cc \
+  Modules/GRandom.cc \
+  Modules/GUnit.cc \
+  Modules/MQuark.cc \
+  Modules/SolRBPrecCG.cc \
+  Modules/SrcPoint.cc \
+  Modules/SrcZ2.cc

programs/Makefile.am

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

scripts/copyright

@@ -5,13 +5,12 @@ while (( "$#" )); do
 echo $1
 
 cat > message  <<EOF
-    /*************************************************************************************
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
 
-    Grid physics library, www.github.com/paboyle/Grid 
+Source file: $1
 
-    Source file: $1
+Copyright (C) 2015
-
-    Copyright (C) 2015
 
 EOF
 
@@ -19,23 +18,24 @@ git log $1 | grep Author | sort -u >> message
 
 cat >> message <<EOF
 
-    This program is free software; you can redistribute it and/or modify
+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
+it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
+the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
+(at your option) any later version.
 
-    This program is distributed in the hope that it will be useful,
+This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
+but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
+GNU General Public License for more details.
 
-    You should have received a copy of the GNU General Public License along
+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.,
+with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution 
+directory.
+*******************************************************************************/
 
-    See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
-    /*  END LEGAL */
 EOF
 
 cat message > tmp.fil