From bb580ae077493b656e8ef3e5caf08016d6edb987 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Sat, 7 May 2016 13:19:38 -0700 Subject: [PATCH] Hadrons: significant overhaul of the object registration system, previous version didn't allow dry runs --- programs/Hadrons/AWilson.cc | 12 +- programs/Hadrons/AWilson.hpp | 2 + programs/Hadrons/Application.cc | 38 ++++-- programs/Hadrons/Environment.cc | 204 +++++++++++++++++++------------ programs/Hadrons/Environment.hpp | 113 ++++++++--------- programs/Hadrons/GLoad.cc | 10 +- programs/Hadrons/GLoad.hpp | 7 +- programs/Hadrons/GRandom.cc | 10 +- programs/Hadrons/GRandom.hpp | 6 +- programs/Hadrons/GUnit.cc | 10 +- programs/Hadrons/GUnit.hpp | 6 +- programs/Hadrons/Graph.hpp | 95 ++++++++++++-- programs/Hadrons/MQuark.cc | 46 +++---- programs/Hadrons/MQuark.hpp | 3 - programs/Hadrons/Module.cc | 1 - programs/Hadrons/Module.hpp | 2 - programs/Hadrons/SolRBPrecCG.cc | 4 +- programs/Hadrons/SrcPoint.cc | 14 +-- programs/Hadrons/SrcPoint.hpp | 7 +- programs/Hadrons/SrcZ2.cc | 16 +-- programs/Hadrons/SrcZ2.hpp | 7 +- 21 files changed, 364 insertions(+), 249 deletions(-) diff --git a/programs/Hadrons/AWilson.cc b/programs/Hadrons/AWilson.cc index 0be43ee7..5a75a257 100644 --- a/programs/Hadrons/AWilson.cc +++ b/programs/Hadrons/AWilson.cc @@ -59,6 +59,15 @@ std::vector AWilson::getOutput(void) return out; } +// setup /////////////////////////////////////////////////////////////////////// +void AWilson::setup(void) +{ + unsigned int size; + + size = 3*env().lattice4dSize(); + 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(); - env().addFermionMatrix(getName(), fMatPt, size); + env().addFermionMatrix(getName(), fMatPt); } diff --git a/programs/Hadrons/AWilson.hpp b/programs/Hadrons/AWilson.hpp index df2d6f41..1eaf0a36 100644 --- a/programs/Hadrons/AWilson.hpp +++ b/programs/Hadrons/AWilson.hpp @@ -56,6 +56,8 @@ public: // dependencies/products virtual std::vector getInput(void); virtual std::vector getOutput(void); + // setup + virtual void setup(void); // execution virtual void execute(void); private: diff --git a/programs/Hadrons/Application.cc b/programs/Hadrons/Application.cc index 7fb69ab1..aa7e66bc 100644 --- a/programs/Hadrons/Application.cc +++ b/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 moduleGraph; + // memory peak and invers memory peak functions + auto memPeak = [this](const std::vector &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 &program) + { + return 1./memPeak(program); + }; + + Graph moduleGraph; LOG(Message) << "Scheduling computation..." << std::endl; @@ -161,7 +185,8 @@ void Application::schedule(void) // } // env_.dryRun(false); std::vector 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 &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 &program) (*module_[program[i]])(); size = env_.getTotalSize(); // print used memory after execution - LOG(Message) << "Allocated objects: " << sizeString(size*locVol_) - << " (" << sizeString(size) << "/site)" << std::endl; + LOG(Message) << "Allocated objects: " << MEM_MSG(size) << std::endl; if (size > memPeak) { memPeak = size; @@ -237,7 +262,7 @@ unsigned int Application::execute(const std::vector &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 &program) } // print used memory after garbage collection size = env_.getTotalSize(); - LOG(Message) << "Allocated objects: " << sizeString(size*locVol_) - << " (" << sizeString(size) << "/site)" << std::endl; + LOG(Message) << "Allocated objects: " << MEM_MSG(size) << std::endl; } return memPeak; diff --git a/programs/Hadrons/Environment.cc b/programs/Hadrons/Environment.cc index 0f870ba2..ee0b842d 100644 --- a/programs/Hadrons/Environment.cc +++ b/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, - const unsigned int size) +void Environment::addFermionMatrix(const std::string name, FMat *fMat) { - fMat_[name].reset(fMat); - addSize(name, size); + if (hasObject(name)) + { + 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, - const std::string actionName) +void Environment::addSolver(const std::string name, Solver s) { - solver_[name] = s; - solverAction_[name] = actionName; + if (hasObject(name)) + { + 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()); -} - - -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; - } + return (hasObject(name) and (lattice_.find(name) != lattice_.end())); } // 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; -} diff --git a/programs/Hadrons/Environment.hpp b/programs/Hadrons/Environment.hpp index f5265181..7f55308a 100644 --- a/programs/Hadrons/Environment.hpp +++ b/programs/Hadrons/Environment.hpp @@ -47,6 +47,11 @@ public: typedef std::unique_ptr RngPt; typedef std::unique_ptr FMatPt; typedef std::unique_ptr 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, - const unsigned int size); + void addFermionMatrix(const std::string name, FMat *mat); 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, - const std::string actionName); + void addSolver(const std::string name, Solver s); 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 - unsigned int lattice4dSize(void) const; - template - void create(const std::string name, - const unsigned int Ls = 1); + T * create(const std::string name); template 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; - unsigned int getLatticeLs(const std::string name) const; + void freeLattice(const std::string name); + template + 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 + 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}; - unsigned int traj_; - GridPt grid4d_; - std::map grid5d_; - GridRbPt gridRb4d_; - std::map gridRb5d_; - RngPt rng4d_; - std::map fMat_; - std::map solver_; - std::map solverAction_; - std::map lattice_; - std::map objectSize_; + bool dryRun_{false}; + unsigned int traj_; + GridPt grid4d_; + std::map grid5d_; + GridRbPt gridRb4d_; + std::map gridRb5d_; + RngPt rng4d_; + std::map object_; + std::map lattice_; + std::map fMat_; + std::map solver_; + std::map solverAction_; std::map> owners_; std::map> properties_; }; @@ -123,41 +135,13 @@ unsigned int Environment::lattice4dSize(void) const } template -void Environment::create(const std::string name, const unsigned int Ls) +T * Environment::create(const std::string name) { - GridCartesian *g4 = getGrid(); - GridCartesian *g; + GridCartesian *g = getGrid(getObjectLs(name)); - if (hasLattice(name)) - { - HADRON_ERROR("object '" + name + "' already exists"); - } - 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()*Ls); + lattice_[name].reset(new T(g)); + + return dynamic_cast(lattice_[name].get()); } template @@ -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 +void Environment::registerLattice(const std::string name, const unsigned int Ls) +{ + createGrid(Ls); + registerObject(name, Ls*lattice4dSize()); +} + END_HADRONS_NAMESPACE #endif // Hadrons_Environment_hpp_ diff --git a/programs/Hadrons/GLoad.cc b/programs/Hadrons/GLoad.cc index 55715fbb..cff9f3ed 100644 --- a/programs/Hadrons/GLoad.cc +++ b/programs/Hadrons/GLoad.cc @@ -59,11 +59,10 @@ std::vector GLoad::getOutput(void) return out; } -// allocation ////////////////////////////////////////////////////////////////// -void GLoad::allocate(void) +// setup /////////////////////////////////////////////////////////////////////// +void GLoad::setup(void) { - env().create(getName()); - gauge_ = env().get(getName()); + env().registerLattice(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(getName()); + NerscIO::readConfiguration(U, header, fileName); LOG(Message) << "NERSC header:" << std::endl; dump_nersc_header(header, LOG(Message)); } diff --git a/programs/Hadrons/GLoad.hpp b/programs/Hadrons/GLoad.hpp index 4527c33c..cb076bb5 100644 --- a/programs/Hadrons/GLoad.hpp +++ b/programs/Hadrons/GLoad.hpp @@ -55,13 +55,12 @@ public: // dependency relation virtual std::vector getInput(void); virtual std::vector getOutput(void); - // allocation - virtual void allocate(void); + // setup + virtual void setup(void); // execution virtual void execute(void); private: - Par par_; - LatticeGaugeField *gauge_ = nullptr; + Par par_; }; MODULE_REGISTER(GLoad); diff --git a/programs/Hadrons/GRandom.cc b/programs/Hadrons/GRandom.cc index 53007045..7f81624c 100644 --- a/programs/Hadrons/GRandom.cc +++ b/programs/Hadrons/GRandom.cc @@ -51,16 +51,16 @@ std::vector GRandom::getOutput(void) return out; } -// allocation ////////////////////////////////////////////////////////////////// -void GRandom::allocate(void) +// setup /////////////////////////////////////////////////////////////////////// +void GRandom::setup(void) { - env().create(getName()); - gauge_ = env().get(getName()); + env().registerLattice(getName()); } // execution /////////////////////////////////////////////////////////////////// void GRandom::execute(void) { LOG(Message) << "Generating random gauge configuration" << std::endl; - SU3::HotConfiguration(*env().get4dRng(), *gauge_); + LatticeGaugeField &U = *env().create(getName()); + SU3::HotConfiguration(*env().get4dRng(), U); } diff --git a/programs/Hadrons/GRandom.hpp b/programs/Hadrons/GRandom.hpp index 12d5be5f..5d0f4197 100644 --- a/programs/Hadrons/GRandom.hpp +++ b/programs/Hadrons/GRandom.hpp @@ -47,12 +47,10 @@ public: // dependency relation virtual std::vector getInput(void); virtual std::vector getOutput(void); - // allocation - virtual void allocate(void); + // setup + virtual void setup(void); // execution virtual void execute(void); -private: - LatticeGaugeField *gauge_ = nullptr; }; MODULE_REGISTER(GRandom); diff --git a/programs/Hadrons/GUnit.cc b/programs/Hadrons/GUnit.cc index 5a86aa00..75bdc597 100644 --- a/programs/Hadrons/GUnit.cc +++ b/programs/Hadrons/GUnit.cc @@ -51,16 +51,16 @@ std::vector GUnit::getOutput(void) return out; } -// allocation ////////////////////////////////////////////////////////////////// -void GUnit::allocate(void) +// setup /////////////////////////////////////////////////////////////////////// +void GUnit::setup(void) { - env().create(getName()); - gauge_ = env().get(getName()); + env().registerLattice(getName()); } // execution /////////////////////////////////////////////////////////////////// void GUnit::execute(void) { LOG(Message) << "Creating unit gauge configuration" << std::endl; - SU3::ColdConfiguration(*env().get4dRng(), *gauge_); + LatticeGaugeField &U = *env().create(getName()); + SU3::ColdConfiguration(*env().get4dRng(), U); } diff --git a/programs/Hadrons/GUnit.hpp b/programs/Hadrons/GUnit.hpp index 4dfc268d..3c72bd07 100644 --- a/programs/Hadrons/GUnit.hpp +++ b/programs/Hadrons/GUnit.hpp @@ -47,12 +47,10 @@ public: // dependencies/products virtual std::vector getInput(void); virtual std::vector getOutput(void); - // allocation - virtual void allocate(void); + // setup + virtual void setup(void); // execution virtual void execute(void); -private: - LatticeGaugeField *gauge_ = nullptr; }; MODULE_REGISTER(GUnit); diff --git a/programs/Hadrons/Graph.hpp b/programs/Hadrons/Graph.hpp index e43f98d1..6022dd3a 100644 --- a/programs/Hadrons/Graph.hpp +++ b/programs/Hadrons/Graph.hpp @@ -52,7 +52,7 @@ public: typedef std::pair Edge; public: // constructor - Graph(void) = default; + Graph(void); // destructor virtual ~Graph(void) = default; // access @@ -71,7 +71,7 @@ public: std::vector getParents(const T &value) const; std::vector getRoots(void) const; std::vector> getConnectedComponents(void) const; - std::vector topoSort(void); + std::vector topoSort(const bool randomize = false); std::vector> allTopoSort(void); // I/O friend std::ostream & operator<<(std::ostream &out, const Graph &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 isMarked_; - std::set edgeSet_; + std::map isMarked_; + std::set edgeSet_; + std::mt19937 gen_; }; // build depedency matrix from topological sorts @@ -121,6 +124,15 @@ makeDependencyMatrix(const std::vector> &topSort); * in the worst case E = V^2 */ +// constructor ///////////////////////////////////////////////////////////////// +template +Graph::Graph(void) +{ + std::random_device rd; + + gen_.seed(rd()); +} + // access ////////////////////////////////////////////////////////////////////// // complexity: log(V) template @@ -297,6 +309,37 @@ const T * Graph::getFirstMarked(const bool isMarked) const } } +// complexity: O(log(V)) +template +const T * Graph::getRandomMarked(const bool isMarked) +{ + auto pred = [&isMarked](const std::pair &v) + { + return (v.second == isMarked); + }; + std::uniform_int_distribution 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 const T * Graph::getFirstUnmarked(void) const @@ -304,6 +347,13 @@ const T * Graph::getFirstUnmarked(void) const return getFirstMarked(false); } +// complexity: O(log(V)) +template +const T * Graph::getRandomUnmarked(void) +{ + return getRandomMarked(false); +} + // prune marked/unmarked vertices ////////////////////////////////////////////// // complexity: O(V^2*log(V)) template @@ -462,10 +512,10 @@ std::vector> Graph::getConnectedComponents(void) const return res; } -// topological sort using Tarjan's algorithm +// topological sort using a directed DFS algorithm // complexity: O(V*log(V)) template -std::vector Graph::topoSort(void) +std::vector Graph::topoSort(const bool randomize) { std::stack buf; std::vector res; @@ -479,16 +529,20 @@ std::vector Graph::topoSort(void) { HADRON_ERROR("cannot topologically sort a cyclic graph"); } - if (!this->isMarked(v)) + if (!isMarked(v)) { - std::vector child = this->getChildren(v); - + std::vector 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 Graph::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 Graph::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 std::vector> Graph::allTopoSort(void) { diff --git a/programs/Hadrons/MQuark.cc b/programs/Hadrons/MQuark.cc index d27a545e..c6af53f4 100644 --- a/programs/Hadrons/MQuark.cc +++ b/programs/Hadrons/MQuark.cc @@ -63,27 +63,11 @@ std::vector MQuark::getOutput(void) // setup /////////////////////////////////////////////////////////////////////// void MQuark::setup(void) { - auto dim = env().getFermionMatrix(par_.solver)->Grid()->GlobalDimensions(); - - if (dim.size() == Nd) - { - Ls_ = 1; - } - else - { - Ls_ = dim[0]; - } -} - -// allocation ////////////////////////////////////////////////////////////////// -void MQuark::allocate(void) -{ - env().create(getName()); - quark_ = env().get(getName()); + Ls_ = env().getObjectLs(env().getSolverAction(par_.solver)); + env().registerLattice(getName()); if (Ls_ > 1) { - env().create(getName() + "_5d", Ls_); - quark5d_ = env().get(getName() + "_5d"); + env().registerLattice(getName() + "_5d", Ls_); } } @@ -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(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(*quark_, sol, s, c); - } - else - { - HADRON_ERROR("MQuark not implemented with 5D actions"); - } + FermToProp(prop, sol, s, c); + } + // create 4D propagators from 5D one if necessary + if (Ls_ > 1) + { + HADRON_ERROR("MQuark not implemented with 5D actions"); } } diff --git a/programs/Hadrons/MQuark.hpp b/programs/Hadrons/MQuark.hpp index aaa3864e..14ced177 100644 --- a/programs/Hadrons/MQuark.hpp +++ b/programs/Hadrons/MQuark.hpp @@ -58,14 +58,11 @@ public: virtual std::vector 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}; }; diff --git a/programs/Hadrons/Module.cc b/programs/Hadrons/Module.cc index 3233c210..a2dedb52 100644 --- a/programs/Hadrons/Module.cc +++ b/programs/Hadrons/Module.cc @@ -55,7 +55,6 @@ Environment & Module::env(void) const void Module::operator()(void) { setup(); - allocate(); if (!env().isDryRun()) { execute(); diff --git a/programs/Hadrons/Module.hpp b/programs/Hadrons/Module.hpp index 605aec06..2b4bfabe 100644 --- a/programs/Hadrons/Module.hpp +++ b/programs/Hadrons/Module.hpp @@ -69,8 +69,6 @@ public: virtual std::vector getOutput(void) = 0; // setup virtual void setup(void) {}; - // allocation - virtual void allocate(void) {}; // execution void operator()(void); virtual void execute(void) = 0; diff --git a/programs/Hadrons/SolRBPrecCG.cc b/programs/Hadrons/SolRBPrecCG.cc index 8487991e..fc1d00aa 100644 --- a/programs/Hadrons/SolRBPrecCG.cc +++ b/programs/Hadrons/SolRBPrecCG.cc @@ -63,7 +63,9 @@ std::vector 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); } diff --git a/programs/Hadrons/SrcPoint.cc b/programs/Hadrons/SrcPoint.cc index 40046cd0..559e56d3 100644 --- a/programs/Hadrons/SrcPoint.cc +++ b/programs/Hadrons/SrcPoint.cc @@ -59,11 +59,10 @@ std::vector SrcPoint::getOutput(void) return out; } -// allocation ////////////////////////////////////////////////////////////////// -void SrcPoint::allocate(void) +// setup /////////////////////////////////////////////////////////////////////// +void SrcPoint::setup(void) { - env().create(getName()); - src_ = env().get(getName()); + env().registerLattice(getName()); } // execution /////////////////////////////////////////////////////////////////// @@ -74,7 +73,8 @@ void SrcPoint::execute(void) LOG(Message) << "Creating point source at position [" << par_.position << "]" << std::endl; - id = 1.; - *src_ = zero; - pokeSite(id, *src_, position); + LatticePropagator &src = *env().create(getName()); + id = 1.; + src = zero; + pokeSite(id, src, position); } diff --git a/programs/Hadrons/SrcPoint.hpp b/programs/Hadrons/SrcPoint.hpp index 75cd9aa6..8cfa89c9 100644 --- a/programs/Hadrons/SrcPoint.hpp +++ b/programs/Hadrons/SrcPoint.hpp @@ -66,13 +66,12 @@ public: // dependency relation virtual std::vector getInput(void); virtual std::vector getOutput(void); - // allocation - virtual void allocate(void); + // setup + virtual void setup(void); // execution virtual void execute(void); private: - Par par_; - LatticePropagator *src_{nullptr}; + Par par_; }; MODULE_REGISTER(SrcPoint); diff --git a/programs/Hadrons/SrcZ2.cc b/programs/Hadrons/SrcZ2.cc index 42d5dbb2..0ed6227d 100644 --- a/programs/Hadrons/SrcZ2.cc +++ b/programs/Hadrons/SrcZ2.cc @@ -59,11 +59,10 @@ std::vector SrcZ2::getOutput(void) return out; } -// allocation ////////////////////////////////////////////////////////////////// -void SrcZ2::allocate(void) +// setup /////////////////////////////////////////////////////////////////////// +void SrcZ2::setup(void) { - env().create(getName()); - src_ = env().get(getName()); + env().registerLattice(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(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; + eta = (2.*eta - shift)*(1./::sqrt(2.)); + eta = where((t >= par_.tA) and (t <= par_.tB), eta, 0.*eta); + src = 1.; + src = src*eta; } diff --git a/programs/Hadrons/SrcZ2.hpp b/programs/Hadrons/SrcZ2.hpp index dc4c35c6..a1458693 100644 --- a/programs/Hadrons/SrcZ2.hpp +++ b/programs/Hadrons/SrcZ2.hpp @@ -68,13 +68,12 @@ public: // dependency relation virtual std::vector getInput(void); virtual std::vector getOutput(void); - // allocation - virtual void allocate(void); + // setup + virtual void setup(void); // execution virtual void execute(void); private: - Par par_; - LatticePropagator *src_; + Par par_; }; MODULE_REGISTER(SrcZ2);