Hadrons: significant overhaul of the object registration system, previous version didn't allow dry runs

programs/Hadrons/AWilson.cc

@@ -59,6 +59,15 @@ std::vector<std::string> AWilson::getOutput(void)
     return out;
 }
 
+// setup ///////////////////////////////////////////////////////////////////////
+void AWilson::setup(void)
+{
+    unsigned int size;
+    
+    size = 3*env().lattice4dSize<WilsonFermionR::DoubledGaugeField>();
+    env().registerObject(getName(), size);
+}
+
 // execution ///////////////////////////////////////////////////////////////////
 void AWilson::execute()
 {
@@ -70,6 +79,5 @@ void AWilson::execute()
     
     LOG(Message) << "Setting up Wilson fermion matrix with m= " << par_.mass
                  << " using gauge field '" << par_.gauge << "'" << std::endl;
-    size = 3*env().lattice4dSize<WilsonFermionR::DoubledGaugeField>();
+    env().addFermionMatrix(getName(), fMatPt);
-    env().addFermionMatrix(getName(), fMatPt, size);
 }

programs/Hadrons/AWilson.hpp

@@ -56,6 +56,8 @@ public:
     // 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:

programs/Hadrons/Application.cc

@@ -107,9 +107,33 @@ void Application::parseParameterFile(void)
 }
 
 // schedule computation ////////////////////////////////////////////////////////
+#define MEM_MSG(size)\
+sizeString((size)*locVol_) << " (" << sizeString(size)  << "/site)"
+
 void Application::schedule(void)
 {
-    Graph<std::string> moduleGraph;
+    // memory peak and invers memory peak functions
+    auto memPeak = [this](const std::vector<std::string> &program)
+    {
+        unsigned int memPeak;
+        bool         msg;
+        
+        msg = HadronsLogMessage.isActive();
+        HadronsLogMessage.Active(false);
+        env_.dryRun(true);
+        memPeak = execute(program);
+        env_.dryRun(false);
+        env_.freeAll();
+        HadronsLogMessage.Active(true);
+        
+        return memPeak;
+    };
+    auto invMemPeak = [&memPeak](const std::vector<std::string> &program)
+    {
+        return 1./memPeak(program);
+    };
+    
+    Graph<std::string>            moduleGraph;
     
     LOG(Message) << "Scheduling computation..." << std::endl;
     
@@ -161,7 +185,8 @@ void Application::schedule(void)
 //        }
 //        env_.dryRun(false);
         std::vector<std::string> t = con[i].topoSort();
-        LOG(Message) << "Program " << i + 1 << ":" << std::endl;
+        LOG(Message) << "Program " << i + 1 << " (memory peak: "
+                     << MEM_MSG(memPeak(t)) << "):" << std::endl;
         for (unsigned int j = 0; j < t.size(); ++j)
         {
             program_.push_back(t[j]);
@@ -210,6 +235,7 @@ unsigned int Application::execute(const std::vector<std::string> &program)
             freeProg[program.rend() - it - 1].insert(n.first);
         }
     }
+    
     // program execution
     for (unsigned int i = 0; i < program.size(); ++i)
     {
@@ -220,8 +246,7 @@ unsigned int Application::execute(const std::vector<std::string> &program)
         (*module_[program[i]])();
         size = env_.getTotalSize();
         // print used memory after execution
-        LOG(Message) << "Allocated objects: " << sizeString(size*locVol_)
+        LOG(Message) << "Allocated objects: " << MEM_MSG(size) << std::endl;
-                     << " (" << sizeString(size)  << "/site)" << std::endl;
         if (size > memPeak)
         {
             memPeak = size;
@@ -237,7 +262,7 @@ unsigned int Application::execute(const std::vector<std::string> &program)
                 // continue garbage collection while there are still
                 // objects without owners
                 continueCollect = continueCollect or !env_.hasOwners(n);
-                if(env_.free(n))
+                if(env_.freeObject(n))
                 {
                     // if an object has been freed, remove it from
                     // the garbage collection schedule
@@ -256,8 +281,7 @@ unsigned int Application::execute(const std::vector<std::string> &program)
         }
         // print used memory after garbage collection
         size = env_.getTotalSize();
-        LOG(Message) << "Allocated objects: " << sizeString(size*locVol_)
+        LOG(Message) << "Allocated objects: " << MEM_MSG(size) << std::endl;
-                     << " (" << sizeString(size)  << "/site)" << std::endl;
     }
     
     return memPeak;

programs/Hadrons/Environment.cc

@@ -67,6 +67,14 @@ unsigned int Environment::getTrajectory(void) const
 }
 
 // grids ///////////////////////////////////////////////////////////////////////
+void Environment::createGrid(const unsigned int Ls)
+{
+    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
@@ -82,7 +90,7 @@ GridCartesian * Environment::getGrid(const unsigned int Ls) const
     }
     catch(std::out_of_range &)
     {
-        HADRON_ERROR("no 5D grid with Ls= " << Ls);
+        HADRON_ERROR("no grid with Ls= " << Ls);
     }
 }
 
@@ -106,11 +114,16 @@ GridRedBlackCartesian * Environment::getRbGrid(const unsigned int Ls) const
 }
 
 // fermion actions /////////////////////////////////////////////////////////////
-void Environment::addFermionMatrix(const std::string name, FMat *fMat,
+void Environment::addFermionMatrix(const std::string name, FMat *fMat)
-                                   const unsigned int size)
 {
-    fMat_[name].reset(fMat);
+    if (hasObject(name))
-    addSize(name, size);
+    {
+        fMat_[name].reset(fMat);
+    }
+    else
+    {
+        HADRON_ERROR("no object named '" << name << "'");
+    }
 }
 
 Environment::FMat * Environment::getFermionMatrix(const std::string name) const
@@ -138,41 +151,66 @@ void Environment::freeFermionMatrix(const std::string name)
     {
         LOG(Message) << "freeing fermion matrix '" << name << "'" << std::endl;
         fMat_.erase(name);
-        objectSize_.erase(name);
+        object_.erase(name);
     }
     else
     {
-        HADRON_ERROR("trying to free undefined fermion matrix '" + name + "'");
+        HADRON_ERROR("trying to free unknown fermion matrix '" + name + "'");
     }
 }
 
 bool Environment::hasFermionMatrix(const std::string name) const
 {
-    return (fMat_.find(name) != fMat_.end());
+    return (hasObject(name) and (fMat_.find(name) != fMat_.end()));
 }
 
 // solvers /////////////////////////////////////////////////////////////////////
-void Environment::addSolver(const std::string name, Solver s,
+void Environment::addSolver(const std::string name, Solver s)
-                            const std::string actionName)
 {
-    solver_[name]       = s;
+    if (hasObject(name))
-    solverAction_[name] = actionName;
+    {
+        solver_[name] = s;
+    }
+    else
+    {
+        HADRON_ERROR("no object named '" << name << "'");
+    }
 }
 
 bool Environment::hasSolver(const std::string name) const
 {
-    return (solver_.find(name) != solver_.end());
+    return (hasObject(name) and (solver_.find(name) != solver_.end()));
+}
+
+void Environment::setSolverAction(const std::string name,
+                                  const std::string actionName)
+{
+    if (hasObject(name))
+    {
+        solverAction_[name] = actionName;
+    }
+    else
+    {
+        HADRON_ERROR("no object named '" << name << "'");
+    }
 }
 
 std::string Environment::getSolverAction(const std::string name) const
 {
-    if (hasSolver(name))
+    if (hasObject(name))
     {
-        return solverAction_.at(name);
+        try
+        {
+            return solverAction_.at(name);
+        }
+        catch (std::out_of_range &)
+        {
+            HADRON_ERROR("no action registered for solver '" << name << "'")
+        }
     }
     else
     {
-        HADRON_ERROR("no solver with name '" << name << "'");
+        HADRON_ERROR("no object with name '" << name << "'");
     }
 }
 
@@ -207,47 +245,81 @@ void Environment::freeLattice(const std::string name)
     {
         LOG(Message) << "freeing lattice '" << name << "'" << std::endl;
         lattice_.erase(name);
-        objectSize_.erase(name);
+        object_.erase(name);
     }
     else
     {
-        HADRON_ERROR("trying to free undefined lattice '" + name + "'");
+        HADRON_ERROR("trying to free unknown lattice '" + name + "'");
     }
 }
 
 bool Environment::hasLattice(const std::string name) const
 {
-    return (lattice_.find(name) != lattice_.end());
+    return (hasObject(name) and (lattice_.find(name) != lattice_.end()));
-}
-
-
-bool Environment::isLattice5d(const std::string name) const
-{
-    if (hasLattice(name))
-    {
-        return (lattice_.at(name)->_grid->GlobalDimensions().size() == Nd + 1);
-    }
-    else
-    {
-        HADRON_ERROR("object '" + name + "' undefined");
-        
-        return false;
-    }
-}
-
-unsigned int Environment::getLatticeLs(const std::string name) const
-{
-    if (isLattice5d(name))
-    {
-        return lattice_.at(name)->_grid->GlobalDimensions()[0];
-    }
-    else
-    {
-        return 1;
-    }
 }
 
 // general memory management ///////////////////////////////////////////////////
+bool Environment::hasObject(const std::string name) const
+{
+    return (object_.find(name) != object_.end());
+}
+
+void Environment::registerObject(const std::string name,
+                                 const unsigned int size, const unsigned int Ls)
+{
+    if (!hasObject(name))
+    {
+        ObjInfo info{size, Ls};
+        
+        object_[name] = info;
+    }
+    else
+    {
+        HADRON_ERROR("object '" + name + "' already exists");
+    }
+}
+
+unsigned int Environment::getObjectSize(const std::string name) const
+{
+    if (hasObject(name))
+    {
+        return object_.at(name).size;
+    }
+    else
+    {
+        HADRON_ERROR("no object named '" + name + "'");
+    }
+}
+
+long unsigned int Environment::getTotalSize(void) const
+{
+    long unsigned int size = 0;
+    
+    for (auto &s: object_)
+    {
+        size += s.second.size;
+    }
+    
+    return size;
+}
+
+unsigned int Environment::getObjectLs(const std::string name) const
+{
+    if (hasObject(name))
+    {
+        return object_.at(name).Ls;
+    }
+    else
+    {
+        HADRON_ERROR("no object named '" + name + "'");
+    }
+}
+
+bool Environment::isObject5d(const std::string name) const
+{
+    return (getObjectLs(name) > 1);
+}
+
 void Environment::addOwnership(const std::string owner,
                                const std::string property)
 {
@@ -267,7 +339,7 @@ bool Environment::hasOwners(const std::string name) const
     }
 }
 
-bool Environment::free(const std::string name)
+bool Environment::freeObject(const std::string name)
 {
     if (!hasOwners(name))
     {
@@ -284,6 +356,10 @@ bool Environment::free(const std::string name)
         {
             freeFermionMatrix(name);
         }
+        else if (hasObject(name))
+        {
+            object_.erase(name);
+        }
         
         return true;
     }
@@ -295,41 +371,11 @@ bool Environment::free(const std::string name)
 
 void Environment::freeAll(void)
 {
+    object_.clear();
     lattice_.clear();
     fMat_.clear();
     solver_.clear();
-    objectSize_.clear();
+    solverAction_.clear();
     owners_.clear();
     properties_.clear();
 }
-
-unsigned int Environment::getSize(const std::string name) const
-{
-    if (hasLattice(name) or hasFermionMatrix(name))
-    {
-        return objectSize_.at(name);
-    }
-    else
-    {
-        HADRON_ERROR("object '" + name + "' undefined");
-        
-        return 0;
-    }
-}
-
-long unsigned int Environment::getTotalSize(void) const
-{
-    long unsigned int size = 0;
-    
-    for (auto &s: objectSize_)
-    {
-        size += s.second;
-    }
-    
-    return size;
-}
-
-void Environment::addSize(const std::string name, const unsigned int size)
-{
-    objectSize_[name] = size;
-}

programs/Hadrons/Environment.hpp

@@ -47,6 +47,11 @@ public:
     typedef std::unique_ptr<GridParallelRNG>            RngPt;
     typedef std::unique_ptr<FMat>                       FMatPt;
     typedef std::unique_ptr<LatticeBase>                LatticePt;
+private:
+    struct ObjInfo
+    {
+        unsigned int size, Ls;
+    };
 public:
     // dry run
     void                    dryRun(const bool isDry);
@@ -55,18 +60,19 @@ public:
     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;
     // fermion actions
-    void                    addFermionMatrix(const std::string name, FMat *mat,
+    void                    addFermionMatrix(const std::string name, FMat *mat);
-                                             const unsigned int size);
     FMat *                  getFermionMatrix(const std::string name) const;
     void                    freeFermionMatrix(const std::string name);
     bool                    hasFermionMatrix(const std::string name) const;
     // solvers
-    void                    addSolver(const std::string name, Solver s,
+    void                    addSolver(const std::string name, Solver s);
-                                      const std::string actionName);
     bool                    hasSolver(const std::string name) const;
+    void                    setSolverAction(const std::string name,
+                                            const std::string actionName);
     std::string             getSolverAction(const std::string name) const;
     void                    callSolver(const std::string name,
                                        LatticeFermion &sol,
@@ -76,39 +82,45 @@ public:
     GridParallelRNG *       get4dRng(void) const;
     // lattice store
     template <typename T>
-    unsigned int            lattice4dSize(void) const;
+    T *                     create(const std::string name);
-    template <typename T>
-    void                    create(const std::string name,
-                                   const unsigned int Ls = 1);
     template <typename T>
     T *                     get(const std::string name) const;
-    void                    freeLattice(const std::string name);
     bool                    hasLattice(const std::string name) const;
-    bool                    isLattice5d(const std::string name) const;
+    void                    freeLattice(const std::string name);
-    unsigned int            getLatticeLs(const std::string name) const;
+    template <typename T>
+    unsigned int            lattice4dSize(void) const;
     // general memory management
+    bool                    hasObject(const std::string name) const;
+    void                    registerObject(const std::string name,
+                                           const unsigned int size,
+                                           const unsigned int Ls = 1);
+    template <typename T>
+    void                    registerLattice(const std::string name,
+                                            const unsigned int Ls = 1);
+    unsigned int            getObjectSize(const std::string name) const;
+    long unsigned int       getTotalSize(void) const;
+    unsigned int            getObjectLs(const std::string name) const;
+    bool                    isObject5d(const std::string name) const;
     void                    addOwnership(const std::string owner,
                                          const std::string property);
     bool                    hasOwners(const std::string name) const;
-    bool                    free(const std::string name);
+    bool                    freeObject(const std::string name);
     void                    freeAll(void);
-    unsigned int            getSize(const std::string name) const;
-    long unsigned int       getTotalSize(void) const;
 private:
-    void addSize(const std::string name, const unsigned int size);
+    
 private:
-    bool                                dryRun_{false};
+    bool                                         dryRun_{false};
-    unsigned int                        traj_;
+    unsigned int                                 traj_;
-    GridPt                              grid4d_;
+    GridPt                                       grid4d_;
-    std::map<unsigned int, GridPt>      grid5d_;
+    std::map<unsigned int, GridPt>               grid5d_;
-    GridRbPt                            gridRb4d_;
+    GridRbPt                                     gridRb4d_;
-    std::map<unsigned int, GridRbPt>    gridRb5d_;
+    std::map<unsigned int, GridRbPt>             gridRb5d_;
-    RngPt                               rng4d_;
+    RngPt                                        rng4d_;
-    std::map<std::string, FMatPt>       fMat_;
+    std::map<std::string, ObjInfo>               object_;
-    std::map<std::string, Solver>       solver_;
+    std::map<std::string, LatticePt>             lattice_;
-    std::map<std::string, std::string>  solverAction_;
+    std::map<std::string, FMatPt>                fMat_;
-    std::map<std::string, LatticePt>    lattice_;
+    std::map<std::string, Solver>                solver_;
-    std::map<std::string, unsigned int> objectSize_;
+    std::map<std::string, std::string>           solverAction_;
     std::map<std::string, std::set<std::string>> owners_;
     std::map<std::string, std::set<std::string>> properties_;
 };
@@ -123,41 +135,13 @@ unsigned int Environment::lattice4dSize(void) const
 }
 
 template <typename T>
-void Environment::create(const std::string name, const unsigned int Ls)
+T * Environment::create(const std::string name)
 {
-    GridCartesian *g4 = getGrid();
+    GridCartesian *g = getGrid(getObjectLs(name));
-    GridCartesian *g;
     
-    if (hasLattice(name))
+    lattice_[name].reset(new T(g));
-    {
+
-        HADRON_ERROR("object '" + name + "' already exists");
+    return dynamic_cast<T *>(lattice_[name].get());
-    }
-    if (Ls > 1)
-    {
-        try
-        {
-            g = grid5d_.at(Ls).get();
-        }
-        catch(std::out_of_range &)
-        {
-            grid5d_[Ls].reset(SpaceTimeGrid::makeFiveDimGrid(Ls, g4));
-            gridRb5d_[Ls].reset(SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, g4));
-            g = grid5d_[Ls].get();
-        }
-    }
-    else
-    {
-        g = g4;
-    }
-    if (!isDryRun())
-    {
-        lattice_[name].reset(new T(g));
-    }
-    else
-    {
-        lattice_[name].reset(nullptr);
-    }
-    addSize(name, lattice4dSize<T>()*Ls);
 }
 
 template <typename T>
@@ -178,12 +162,19 @@ T * Environment::get(const std::string name) const
     }
     else
     {
-        HADRON_ERROR("object '" + name + "' undefined");
+        HADRON_ERROR("no lattice name '" + name + "'");
         
         return nullptr;
     }
 }
 
+template <typename T>
+void Environment::registerLattice(const std::string name, const unsigned int Ls)
+{
+    createGrid(Ls);
+    registerObject(name, Ls*lattice4dSize<T>());
+}
+
 END_HADRONS_NAMESPACE
 
 #endif // Hadrons_Environment_hpp_

programs/Hadrons/GLoad.cc

@@ -59,11 +59,10 @@ std::vector<std::string> GLoad::getOutput(void)
     return out;
 }
 
-// allocation //////////////////////////////////////////////////////////////////
+// setup ///////////////////////////////////////////////////////////////////////
-void GLoad::allocate(void)
+void GLoad::setup(void)
 {
-    env().create<LatticeGaugeField>(getName());
+    env().registerLattice<LatticeGaugeField>(getName());
-    gauge_ = env().get<LatticeGaugeField>(getName());
 }
 
 // execution ///////////////////////////////////////////////////////////////////
@@ -75,7 +74,8 @@ void GLoad::execute(void)
     
     LOG(Message) << "Loading NERSC configuration from file '" << fileName
                  << "'" << std::endl;
-    NerscIO::readConfiguration(*gauge_, header, fileName);
+    LatticeGaugeField &U = *env().create<LatticeGaugeField>(getName());
+    NerscIO::readConfiguration(U, header, fileName);
     LOG(Message) << "NERSC header:" << std::endl;
     dump_nersc_header(header, LOG(Message));
 }

programs/Hadrons/GLoad.hpp

@@ -55,13 +55,12 @@ public:
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
-    // allocation
+    // setup
-    virtual void allocate(void);
+    virtual void setup(void);
     // execution
     virtual void execute(void);
 private:
-    Par               par_;
+    Par par_;
-    LatticeGaugeField *gauge_ = nullptr;
 };
 
 MODULE_REGISTER(GLoad);

programs/Hadrons/GRandom.cc

@@ -51,16 +51,16 @@ std::vector<std::string> GRandom::getOutput(void)
     return out;
 }
 
-// allocation //////////////////////////////////////////////////////////////////
+// setup ///////////////////////////////////////////////////////////////////////
-void GRandom::allocate(void)
+void GRandom::setup(void)
 {
-    env().create<LatticeGaugeField>(getName());
+    env().registerLattice<LatticeGaugeField>(getName());
-    gauge_ = env().get<LatticeGaugeField>(getName());
 }
 
 // execution ///////////////////////////////////////////////////////////////////
 void GRandom::execute(void)
 {
     LOG(Message) << "Generating random gauge configuration" << std::endl;
-    SU3::HotConfiguration(*env().get4dRng(), *gauge_);
+    LatticeGaugeField &U = *env().create<LatticeGaugeField>(getName());
+    SU3::HotConfiguration(*env().get4dRng(), U);
 }

programs/Hadrons/GRandom.hpp

@@ -47,12 +47,10 @@ public:
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
-    // allocation
+    // setup
-    virtual void allocate(void);
+    virtual void setup(void);
     // execution
     virtual void execute(void);
-private:
-    LatticeGaugeField *gauge_ = nullptr;
 };
 
 MODULE_REGISTER(GRandom);

programs/Hadrons/GUnit.cc

@@ -51,16 +51,16 @@ std::vector<std::string> GUnit::getOutput(void)
     return out;
 }
 
-// allocation //////////////////////////////////////////////////////////////////
+// setup ///////////////////////////////////////////////////////////////////////
-void GUnit::allocate(void)
+void GUnit::setup(void)
 {
-    env().create<LatticeGaugeField>(getName());
+    env().registerLattice<LatticeGaugeField>(getName());
-    gauge_ = env().get<LatticeGaugeField>(getName());
 }
 
 // execution ///////////////////////////////////////////////////////////////////
 void GUnit::execute(void)
 {
     LOG(Message) << "Creating unit gauge configuration" << std::endl;
-    SU3::ColdConfiguration(*env().get4dRng(), *gauge_);
+    LatticeGaugeField &U = *env().create<LatticeGaugeField>(getName());
+    SU3::ColdConfiguration(*env().get4dRng(), U);
 }

programs/Hadrons/GUnit.hpp

@@ -47,12 +47,10 @@ public:
     // dependencies/products
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
-    // allocation
+    // setup
-    virtual void allocate(void);
+    virtual void setup(void);
     // execution
     virtual void execute(void);
-private:
-    LatticeGaugeField *gauge_ = nullptr;
 };
 
 MODULE_REGISTER(GUnit);

programs/Hadrons/Graph.hpp

@@ -52,7 +52,7 @@ public:
     typedef std::pair<T, T> Edge;
 public:
     // constructor
-    Graph(void) = default;
+    Graph(void);
     // destructor
     virtual ~Graph(void) = default;
     // access
@@ -71,7 +71,7 @@ public:
     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);
+    std::vector<T>              topoSort(const bool randomize = false);
     std::vector<std::vector<T>> allTopoSort(void);
     // I/O
     friend std::ostream & operator<<(std::ostream &out, const Graph<T> &g)
@@ -97,7 +97,9 @@ private:
     void      unmarkAll(void);
     bool      isMarked(const T &value) const;
     const T * getFirstMarked(const bool isMarked = true) const;
+    const T * getRandomMarked(const bool isMarked = true);
     const T * getFirstUnmarked(void) const;
+    const T * getRandomUnmarked(void);
     // prune marked/unmarked vertices
     void removeMarked(const bool isMarked = true);
     void removeUnmarked(void);
@@ -105,8 +107,9 @@ private:
     void depthFirstSearch(void);
     void depthFirstSearch(const T &root);
 private:
-    std::map<T, bool> isMarked_;
+    std::map<T, bool>  isMarked_;
-    std::set<Edge>    edgeSet_;
+    std::set<Edge>     edgeSet_;
+    std::mt19937       gen_;
 };
 
 // build depedency matrix from topological sorts
@@ -121,6 +124,15 @@ makeDependencyMatrix(const std::vector<std::vector<T>> &topSort);
  * in the worst case E = V^2
  */
 
+// constructor /////////////////////////////////////////////////////////////////
+template <typename T>
+Graph<T>::Graph(void)
+{
+    std::random_device rd;
+    
+    gen_.seed(rd());
+}
+
 // access //////////////////////////////////////////////////////////////////////
 // complexity: log(V)
 template <typename T>
@@ -297,6 +309,37 @@ const T * Graph<T>::getFirstMarked(const bool isMarked) const
     }
 }
 
+// complexity: O(log(V))
+template <typename T>
+const T * Graph<T>::getRandomMarked(const bool isMarked)
+{
+    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
@@ -304,6 +347,13 @@ const T * Graph<T>::getFirstUnmarked(void) const
     return getFirstMarked(false);
 }
 
+// complexity: O(log(V))
+template <typename T>
+const T * Graph<T>::getRandomUnmarked(void)
+{
+    return getRandomMarked(false);
+}
+
 // prune marked/unmarked vertices //////////////////////////////////////////////
 // complexity: O(V^2*log(V))
 template <typename T>
@@ -462,10 +512,10 @@ std::vector<Graph<T>> Graph<T>::getConnectedComponents(void) const
     return res;
 }
 
-// topological sort using Tarjan's algorithm
+// topological sort using a directed DFS algorithm
 // complexity: O(V*log(V))
 template <typename T>
-std::vector<T> Graph<T>::topoSort(void)
+std::vector<T> Graph<T>::topoSort(const bool randomize)
 {
     std::stack<T>     buf;
     std::vector<T>    res;
@@ -479,16 +529,20 @@ std::vector<T> Graph<T>::topoSort(void)
         {
             HADRON_ERROR("cannot topologically sort a cyclic graph");
         }
-        if (!this->isMarked(v))
+        if (!isMarked(v))
         {
-            std::vector<T> child = this->getChildren(v);
+            std::vector<T> child = getChildren(v);
-            
+
             tmpMarked[v] = true;
+            if (randomize)
+            {
+                std::shuffle(child.begin(), child.end(), gen_);
+            }
             for (auto &c: child)
             {
                 visit(c);
             }
-            this->mark(v);
+            mark(v);
             tmpMarked[v] = false;
             buf.push(v);
         }
@@ -501,11 +555,26 @@ std::vector<T> Graph<T>::topoSort(void)
     }
     
     // loop on unmarked vertices
-    vPt = getFirstUnmarked();
+    unmarkAll();
+    if (randomize)
+    {
+        vPt = getRandomUnmarked();
+    }
+    else
+    {
+        vPt = getFirstUnmarked();
+    }
     while (vPt)
     {
         visit(*vPt);
-        vPt = getFirstUnmarked();
+        if (randomize)
+        {
+            vPt = getRandomUnmarked();
+        }
+        else
+        {
+            vPt = getFirstUnmarked();
+        }
     }
     unmarkAll();
     
@@ -522,7 +591,7 @@ std::vector<T> Graph<T>::topoSort(void)
 // 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))
+// complexity: O(V*log(V)) (from the paper, but really ?)
 template <typename T>
 std::vector<std::vector<T>> Graph<T>::allTopoSort(void)
 {

programs/Hadrons/MQuark.cc

@@ -63,27 +63,11 @@ std::vector<std::string> MQuark::getOutput(void)
 // setup ///////////////////////////////////////////////////////////////////////
 void MQuark::setup(void)
 {
-    auto dim = env().getFermionMatrix(par_.solver)->Grid()->GlobalDimensions();
+    Ls_ = env().getObjectLs(env().getSolverAction(par_.solver));
-    
+    env().registerLattice<LatticePropagator>(getName());
-    if (dim.size() == Nd)
-    {
-        Ls_ = 1;
-    }
-    else
-    {
-        Ls_ = dim[0];
-    }
-}
-
-// allocation //////////////////////////////////////////////////////////////////
-void MQuark::allocate(void)
-{
-    env().create<LatticePropagator>(getName());
-    quark_ = env().get<LatticePropagator>(getName());
     if (Ls_ > 1)
     {
-        env().create<LatticePropagator>(getName() + "_5d", Ls_);
+        env().registerLattice<LatticePropagator>(getName() + "_5d", Ls_);
-        quark5d_ = env().get<LatticePropagator>(getName() + "_5d");
     }
 }
 
@@ -92,10 +76,13 @@ void MQuark::execute(void)
 {
     LatticePropagator *fullSource;
     LatticeFermion    source(env().getGrid(Ls_)), sol(env().getGrid(Ls_));
+    std::string       propName = (Ls_ == 1) ? getName() : (getName() + "_5d");
     
     LOG(Message) << "Computing quark propagator '" << getName() << "'"
                  << std::endl;
-    if (!env().isLattice5d(par_.source))
+    LatticePropagator &prop = *env().create<LatticePropagator>(propName);
+    // source conversion for 4D sources
+    if (!env().isObject5d(par_.source))
     {
         if (Ls_ == 1)
         {
@@ -106,15 +93,16 @@ void MQuark::execute(void)
             HADRON_ERROR("MQuark not implemented with 5D actions");
         }
     }
+    // source conversion for 5D sources
     else
     {
         if (Ls_ == 1)
         {
             HADRON_ERROR("MQuark not implemented with 5D actions");
         }
-        else if (Ls_ != env().getLatticeLs(par_.source))
+        else if (Ls_ != env().getObjectLs(par_.source))
         {
-            HADRON_ERROR("MQuark not implemented with 5D actions");
+            HADRON_ERROR("Ls mismatch between quark action and source");
         }
         else
         {
@@ -129,13 +117,11 @@ void MQuark::execute(void)
         PropToFerm(source, *fullSource, s, c);
         sol = zero;
         env().callSolver(par_.solver, sol, source);
-        if (Ls_ == 1)
+        FermToProp(prop, sol, s, c);
-        {
+    }
-            FermToProp(*quark_, sol, s, c);
+    // create 4D propagators from 5D one if necessary
-        }
+    if (Ls_ > 1)
-        else
+    {
-        {
+        HADRON_ERROR("MQuark not implemented with 5D actions");
-            HADRON_ERROR("MQuark not implemented with 5D actions");
-        }
     }
 }

programs/Hadrons/MQuark.hpp

@@ -58,14 +58,11 @@ public:
     virtual std::vector<std::string> getOutput(void);
     // setup
     virtual void setup(void);
-    // allocation
-    virtual void allocate(void);
     // execution
     virtual void execute(void);
 private:
     Par                 par_;
     unsigned int        Ls_;
-    LatticePropagator   *source_{nullptr}, *quark_{nullptr}, *quark5d_{nullptr};
     Environment::Solver *solver_{nullptr};
 };
 

programs/Hadrons/Module.cc

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

programs/Hadrons/Module.hpp

@@ -69,8 +69,6 @@ public:
     virtual std::vector<std::string> getOutput(void) = 0;
     // setup
     virtual void setup(void) {};
-    // allocation
-    virtual void allocate(void) {};
     // execution
     void operator()(void);
     virtual void execute(void) = 0;

programs/Hadrons/SolRBPrecCG.cc

@@ -63,7 +63,9 @@ std::vector<std::string> SolRBPrecCG::getOutput(void)
 // setup ///////////////////////////////////////////////////////////////////////
 void SolRBPrecCG::setup(void)
 {
+    env().registerObject(getName(), 0);
     env().addOwnership(getName(), par_.action);
+    env().setSolverAction(getName(), par_.action);
 }
 
 // execution ///////////////////////////////////////////////////////////////////
@@ -82,5 +84,5 @@ void SolRBPrecCG::execute(void)
     LOG(Message) << "setting up Schur red-black preconditioned CG for"
                  << " action '" << par_.action << "' with residual "
                  << par_.residual << std::endl;
-    env().addSolver(getName(), solver, par_.action);
+    env().addSolver(getName(), solver);
 }

programs/Hadrons/SrcPoint.cc

@@ -59,11 +59,10 @@ std::vector<std::string> SrcPoint::getOutput(void)
     return out;
 }
 
-// allocation //////////////////////////////////////////////////////////////////
+// setup ///////////////////////////////////////////////////////////////////////
-void SrcPoint::allocate(void)
+void SrcPoint::setup(void)
 {
-    env().create<LatticePropagator>(getName());
+    env().registerLattice<LatticePropagator>(getName());
-    src_ = env().get<LatticePropagator>(getName());
 }
 
 // execution ///////////////////////////////////////////////////////////////////
@@ -74,7 +73,8 @@ void SrcPoint::execute(void)
     
     LOG(Message) << "Creating point source at position [" << par_.position
                  << "]" << std::endl;
-    id    = 1.;
+    LatticePropagator &src = *env().create<LatticePropagator>(getName());
-    *src_ = zero;
+    id  = 1.;
-    pokeSite(id, *src_, position);
+    src = zero;
+    pokeSite(id, src, position);
 }

programs/Hadrons/SrcPoint.hpp

@@ -66,13 +66,12 @@ public:
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
-    // allocation
+    // setup
-    virtual void allocate(void);
+    virtual void setup(void);
     // execution
     virtual void execute(void);
 private:
-    Par               par_;
+    Par par_;
-    LatticePropagator *src_{nullptr};
 };
 
 MODULE_REGISTER(SrcPoint);

programs/Hadrons/SrcZ2.cc

@@ -59,11 +59,10 @@ std::vector<std::string> SrcZ2::getOutput(void)
     return out;
 }
 
-// allocation //////////////////////////////////////////////////////////////////
+// setup ///////////////////////////////////////////////////////////////////////
-void SrcZ2::allocate(void)
+void SrcZ2::setup(void)
 {
-    env().create<LatticePropagator>(getName());
+    env().registerLattice<LatticePropagator>(getName());
-    src_ = env().get<LatticePropagator>(getName());
 }
 
 // execution ///////////////////////////////////////////////////////////////////
@@ -84,10 +83,11 @@ void SrcZ2::execute(void)
         LOG(Message) << "Generating Z_2 band for " << par_.tA << " <= t <= "
                      << par_.tB << std::endl;
     }
+    LatticePropagator &src = *env().create<LatticePropagator>(getName());
     LatticeCoordinate(t, Tp);
     bernoulli(*env().get4dRng(), eta);
-    eta   = (2.*eta - shift)*(1./::sqrt(2.));
+    eta = (2.*eta - shift)*(1./::sqrt(2.));
-    eta   = where((t >= par_.tA) and (t <= par_.tB), eta, 0.*eta);
+    eta = where((t >= par_.tA) and (t <= par_.tB), eta, 0.*eta);
-    *src_ = 1.;
+    src = 1.;
-    *src_ = (*src_)*eta;
+    src = src*eta;
 }

programs/Hadrons/SrcZ2.hpp

@@ -68,13 +68,12 @@ public:
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
-    // allocation
+    // setup
-    virtual void allocate(void);
+    virtual void setup(void);
     // execution
     virtual void execute(void);
 private:
-    Par               par_;
+    Par par_;
-    LatticePropagator *src_;
 };
 
 MODULE_REGISTER(SrcZ2);