Hadrons: general lattice store & a lot of code cleaning

programs/Hadrons/AWilson.cc

@@ -62,12 +62,13 @@ std::vector<std::string> AWilson::getOutput(void)
 // execution ///////////////////////////////////////////////////////////////////
 void AWilson::execute(Environment &env)
 {
-    auto &U      = *env.getGauge(par_.gauge);
+    auto &U      = *env.get<LatticeGaugeField>(par_.gauge);
     auto &grid   = *env.getGrid();
     auto &gridRb = *env.getRbGrid();
     
-    LOG(Message) << "setting up Wilson fermion matrix with m= " << par_.mass
+    LOG(Message) << "Setting up Wilson fermion matrix with m= " << par_.mass
                  << " using gauge field '" << par_.gauge << "'" << std::endl;
     env.addFermionMatrix(getName(),
                          new WilsonFermionR(U, grid, gridRb, par_.mass));
+    LOG(Message) << sizeof(*env.getFermionMatrix(getName())) << std::endl;
 }

programs/Hadrons/Application.cc

@@ -47,6 +47,16 @@ Application::Application(const std::string parameterFileName)
     {
         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 //////////////////////////////////////////////////////////////////
@@ -169,8 +179,8 @@ void Application::configLoop(void)
     
     for (unsigned int t = range.start; t < range.end; t += range.step)
     {
-        LOG(Message) << "Starting measurement for trajectory " << t
+        LOG(Message) << "========== Starting measurement for trajectory " << t
-                     << std::endl;
+                     << " ==========" << std::endl;
         env_.setTrajectory(t);
         execute(program_);
         env_.freeAll();
@@ -179,7 +189,7 @@ void Application::configLoop(void)
 
 unsigned int Application::execute(const std::vector<std::string> &program)
 {
-    unsigned int                          memPeak = 0;
+    unsigned int                          memPeak = 0, size;
     std::vector<std::vector<std::string>> freeProg;
     
     freeProg.resize(program.size());
@@ -200,19 +210,53 @@ unsigned int Application::execute(const std::vector<std::string> &program)
     }
     for (unsigned int i = 0; i < program.size(); ++i)
     {
-        LOG(Message) << "Measurement step " << i+1 << "/" << program.size()
+        LOG(Message) << "---------- Measurement step " << i+1 << "/"
-                     << " (module '" << program[i] << "')" << std::endl;
+                     << program.size() << " (module '" << program[i] << "')"
+                     << " ----------" << std::endl;
         (*module_[program[i]])(env_);
-        LOG(Message) << "allocated propagators: " << env_.nProp() << std::endl;
+        size = env_.getTotalSize();
-        if (env_.nProp() > memPeak)
+        LOG(Message) << "Allocated objects: " << sizeString(size*locVol_)
+                     << " (" << sizeString(size)  << "/site)" << std::endl;
+        if (size > memPeak)
         {
-            memPeak = env_.nProp();
+            memPeak = size;
         }
+        LOG(Message) << "Garbage collection..." << std::endl;
         for (auto &n: freeProg[i])
         {
             env_.free(n);
         }
+        size = env_.getTotalSize();
+        LOG(Message) << "Allocated objects: " << sizeString(size*locVol_)
+                     << " (" << sizeString(size)  << "/site)" << std::endl;
     }
     
     return memPeak;
 }
+
+// pretty size formatting //////////////////////////////////////////////////////
+std::string Application::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/Application.hpp

@@ -85,7 +85,10 @@ private:
     // program execution
     void         configLoop(void);
     unsigned int execute(const std::vector<std::string> &program);
+    // pretty size formatting
+    std::string sizeString(long unsigned int bytes);
 private:
+    long unsigned int                               locVol_;
     std::string                                     parameterFileName_;
     GlobalPar                                       par_;
     Environment                                     &env_;

programs/Hadrons/CMeson.cc

@@ -63,13 +63,13 @@ std::vector<std::string> CMeson::getOutput(void)
 // execution ///////////////////////////////////////////////////////////////////
 void CMeson::execute(Environment &env)
 {
-    LOG(Message) << "computing meson contraction '" << getName() << "' using"
+    LOG(Message) << "Computing meson contraction '" << getName() << "' using"
                  << " quarks '" << par_.q1 << " and '" << par_.q2 << "'"
                  << std::endl;
     
     XmlWriter             writer(par_.output);
-    LatticePropagator     &q1 = *env.getProp(par_.q1);
+    LatticePropagator     &q1 = *env.get<LatticePropagator>(par_.q1);
-    LatticePropagator     &q2 = *env.getProp(par_.q2);
+    LatticePropagator     &q2 = *env.get<LatticePropagator>(par_.q2);
     LatticeComplex        c(env.getGrid());
     SpinMatrix            g[Ns*Ns], g5;
     std::vector<TComplex> buf;

programs/Hadrons/Environment.cc

@@ -151,173 +151,6 @@ void Environment::callSolver(const std::string name, LatticeFermion &sol,
     }
 }
 
-// quark propagators ///////////////////////////////////////////////////////////
-void Environment::createProp(const std::string name, const unsigned int Ls)
-{
-    GridCartesian *g4 = getGrid();
-
-    if (propExists(name))
-    {
-        HADRON_ERROR("propagator '" + name + "' already exists");
-    }
-    if (Ls > 1)
-    {
-        GridCartesian *g;
-        
-        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();
-        }
-        if (!isDryRun())
-        {
-            prop_[name].reset(new LatticePropagator(g));
-        }
-        else
-        {
-            prop_[name].reset(nullptr);
-        }
-        propSize_[name] = Ls;
-    }
-    else
-    {
-        if (!isDryRun())
-        {
-            prop_[name].reset(new LatticePropagator(g4));
-        }
-        else
-        {
-            prop_[name].reset(nullptr);
-        }
-        propSize_[name] = 1;
-    }
-}
-
-void Environment::freeProp(const std::string name)
-{
-    if (propExists(name))
-    {
-        prop_.erase(name);
-        propSize_.erase(name);
-    }
-    else
-    {
-        HADRON_ERROR("trying to free unknown propagator '" + name + "'");
-    }
-}
-
-bool Environment::isProp5d(const std::string name) const
-{
-    if (propExists(name))
-    {
-        return (getProp(name)->_grid->GlobalDimensions().size() == Nd + 1);
-    }
-    else
-    {
-        HADRON_ERROR("propagator '" + name + "' unknown");
-        
-        return false;
-    }
-}
-
-unsigned int Environment::getPropLs(const std::string name) const
-{
-    if (propExists(name))
-    {
-        if (isProp5d(name))
-        {
-            return getProp(name)->_grid->GlobalDimensions()[0];
-        }
-        else
-        {
-            return 1;
-        }
-    }
-    else
-    {
-        HADRON_ERROR("propagator '" + name + "' unknown");
-        
-        return 0;
-    }
-}
-
-LatticePropagator * Environment::getProp(const std::string name) const
-{
-    if (propExists(name))
-    {
-        return prop_.at(name).get();
-    }
-    else
-    {
-        HADRON_ERROR("propagator '" + name + "' unknown");
-        
-        return nullptr;
-    }
-}
-
-bool Environment::propExists(const std::string name) const
-{
-    return (prop_.find(name) != prop_.end());
-}
-
-unsigned int Environment::nProp(void) const
-{
-    unsigned int size = 0;
-    
-    for (auto &s: propSize_)
-    {
-        size += s.second;
-    }
-    
-    return size;
-}
-
-// gauge configuration /////////////////////////////////////////////////////////
-void Environment::createGauge(const std::string name)
-{
-    if (gaugeExists(name))
-    {
-        HADRON_ERROR("gauge field '" + name + "' already exists");
-    }
-    gauge_[name].reset(new LatticeGaugeField(getGrid()));
-}
-
-void Environment::freeGauge(const std::string name)
-{
-    if (gaugeExists(name))
-    {
-        gauge_.erase(name);
-    }
-    else
-    {
-        HADRON_ERROR("trying to free unknown gauge field '" + name + "'");
-    }
-}
-
-LatticeGaugeField * Environment::getGauge(const std::string name) const
-{
-    if (gaugeExists(name))
-    {
-        return gauge_.at(name).get();
-    }
-    else
-    {
-        HADRON_ERROR("gauge field '" + name + "' unknown");
-        
-        return nullptr;
-    }
-}
-
-bool Environment::gaugeExists(const std::string name) const
-{
-    return (gauge_.find(name) != gauge_.end());
-}
-
 // random number generator /////////////////////////////////////////////////////
 void Environment::setSeed(const std::vector<int> &seed)
 {
@@ -329,24 +162,95 @@ GridParallelRNG * Environment::get4dRng(void) const
     return rng4d_.get();
 }
 
-// general free ////////////////////////////////////////////////////////////////
+// data store //////////////////////////////////////////////////////////////////
+void Environment::freeLattice(const std::string name)
+{
+    if (hasLattice(name))
+    {
+        LOG(Message) << "freeing lattice '" << name << "'" << std::endl;
+        lattice_.erase(name);
+        objectSize_.erase(name);
+    }
+    else
+    {
+        HADRON_ERROR("trying to free undefined 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;
+    }
+}
+
+// general memory management ///////////////////////////////////////////////////
 void Environment::free(const std::string name)
 {
-    if (propExists(name))
+    if (hasLattice(name))
     {
-        LOG(Message) << "freeing propagator '" << name << "'" << std::endl;
+        freeLattice(name);
-        freeProp(name);
-    }
-    else if (gaugeExists(name))
-    {
-        LOG(Message) << "freeing gauge field '" << name << "'" << std::endl;
-        freeGauge(name);
     }
 }
 
 void Environment::freeAll(void)
 {
-    prop_.clear();
+    lattice_.clear();
-    propSize_.clear();
+    objectSize_.clear();
-    gauge_.clear();
+}
+
+void Environment::addSize(const std::string name, const unsigned int size)
+{
+    objectSize_[name] = size;
+}
+
+unsigned int Environment::getSize(const std::string name) const
+{
+    if (hasLattice(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;
 }

programs/Hadrons/Environment.hpp

@@ -46,8 +46,7 @@ public:
     typedef std::unique_ptr<GridRedBlackCartesian>      GridRbPt;
     typedef std::unique_ptr<GridParallelRNG>            RngPt;
     typedef std::unique_ptr<FMat>                       FMatPt;
-    typedef std::unique_ptr<LatticePropagator>          PropPt;
+    typedef std::unique_ptr<LatticeBase>                LatticePt;
-    typedef std::unique_ptr<LatticeGaugeField>          GaugePt;
 public:
     // dry run
     void                    dryRun(const bool isDry);
@@ -67,26 +66,27 @@ public:
     void                    callSolver(const std::string name,
                                        LatticeFermion &sol,
                                        const LatticeFermion &src) const;
-    // quark propagators
-    void                    createProp(const std::string name,
-                                       const unsigned int Ls = 1);
-    void                    freeProp(const std::string name);
-    bool                    isProp5d(const std::string name) const;
-    unsigned int            getPropLs(const std::string name) const;
-    LatticePropagator *     getProp(const std::string name) const;
-    bool                    propExists(const std::string name) const;
-    unsigned int            nProp(void) const;
-    // gauge configurations
-    void                    createGauge(const std::string name);
-    void                    freeGauge(const std::string name);
-    LatticeGaugeField *     getGauge(const std::string name) const;
-    bool                    gaugeExists(const std::string name) const;
     // random number generator
     void                    setSeed(const std::vector<int> &seed);
     GridParallelRNG *       get4dRng(void) const;
-    // general free
+    // lattice store
+    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;
+    unsigned int            getLatticeLs(const std::string name) const;
+    // general memory management
     void                    free(const std::string name);
     void                    freeAll(void);
+    void                    addSize(const std::string name,
+                                    const unsigned int size);
+    unsigned int            getSize(const std::string name) const;
+    long unsigned int       getTotalSize(void) const;
 private:
     bool                                dryRun_{false};
     unsigned int                        traj_;
@@ -98,11 +98,75 @@ private:
     std::map<std::string, FMatPt>       fMat_;
     std::map<std::string, Solver>       solver_;
     std::map<std::string, std::string>  solverAction_;
-    std::map<std::string, PropPt>       prop_;
+    std::map<std::string, LatticePt>    lattice_;
-    std::map<std::string, unsigned int> propSize_;
+    std::map<std::string, unsigned int> objectSize_;
-    std::map<std::string, GaugePt>      gauge_;
 };
 
+/******************************************************************************
+ *                        template implementation                             *
+ ******************************************************************************/
+template <typename T>
+void Environment::create(const std::string name, const unsigned int Ls)
+{
+    GridCartesian *g4 = getGrid();
+    GridCartesian *g;
+    
+    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);
+    }
+    objectSize_[name] = sizeof(typename T::vector_object)/g->Nsimd()*Ls;
+}
+
+template <typename T>
+T * Environment::get(const std::string name) const
+{
+    if (hasLattice(name))
+    {
+        try
+        {
+            return dynamic_cast<T *>(lattice_.at(name).get());
+        }
+        catch (std::bad_cast &)
+        {
+            HADRON_ERROR("object '" + name + "' does not have type "
+                         + typeid(T *).name() + "(object type: "
+                         + typeid(lattice_.at(name).get()).name() + ")");
+        }
+    }
+    else
+    {
+        HADRON_ERROR("object '" + name + "' undefined");
+        
+        return nullptr;
+    }
+}
+
 END_HADRONS_NAMESPACE
 
 #endif // Hadrons_Environment_hpp_

programs/Hadrons/GLoad.cc

@@ -62,8 +62,8 @@ std::vector<std::string> GLoad::getOutput(void)
 // allocation //////////////////////////////////////////////////////////////////
 void GLoad::allocate(Environment &env)
 {
-    env.createGauge(getName());
+    env.create<LatticeGaugeField>(getName());
-    gauge_ = env.getGauge(getName());
+    gauge_ = env.get<LatticeGaugeField>(getName());
 }
 
 // execution ///////////////////////////////////////////////////////////////////
@@ -73,7 +73,9 @@ void GLoad::execute(Environment &env)
     std::string fileName = par_.file + "."
                            + std::to_string(env.getTrajectory());
     
-    LOG(Message) << "loading NERSC configuration from file '" << fileName
+    LOG(Message) << "Loading NERSC configuration from file '" << fileName
                  << "'" << std::endl;
     NerscIO::readConfiguration(*gauge_, header, fileName);
+    LOG(Message) << "NERSC header:" << std::endl;
+    dump_nersc_header(header, LOG(Message));
 }

programs/Hadrons/GRandom.cc

@@ -54,13 +54,13 @@ std::vector<std::string> GRandom::getOutput(void)
 // allocation //////////////////////////////////////////////////////////////////
 void GRandom::allocate(Environment &env)
 {
-    env.createGauge(getName());
+    env.create<LatticeGaugeField>(getName());
-    gauge_ = env.getGauge(getName());
+    gauge_ = env.get<LatticeGaugeField>(getName());
 }
 
 // execution ///////////////////////////////////////////////////////////////////
 void GRandom::execute(Environment &env)
 {
-    LOG(Message) << "generating random gauge configuration" << std::endl;
+    LOG(Message) << "Generating random gauge configuration" << std::endl;
     SU3::HotConfiguration(*env.get4dRng(), *gauge_);
 }

programs/Hadrons/GUnit.cc

@@ -54,13 +54,13 @@ std::vector<std::string> GUnit::getOutput(void)
 // allocation //////////////////////////////////////////////////////////////////
 void GUnit::allocate(Environment &env)
 {
-    env.createGauge(getName());
+    env.create<LatticeGaugeField>(getName());
-    gauge_ = env.getGauge(getName());
+    gauge_ = env.get<LatticeGaugeField>(getName());
 }
 
 // execution ///////////////////////////////////////////////////////////////////
 void GUnit::execute(Environment &env)
 {
-    LOG(Message) << "creating unit gauge configuration" << std::endl;
+    LOG(Message) << "Creating unit gauge configuration" << std::endl;
     SU3::ColdConfiguration(*env.get4dRng(), *gauge_);
 }

programs/Hadrons/MQuark.cc

@@ -78,12 +78,12 @@ void MQuark::setup(Environment &env)
 // allocation //////////////////////////////////////////////////////////////////
 void MQuark::allocate(Environment &env)
 {
-    env.createProp(getName());
+    env.create<LatticePropagator>(getName());
-    quark_ = env.getProp(getName());
+    quark_ = env.get<LatticePropagator>(getName());
     if (Ls_ > 1)
     {
-        env.createProp(getName() + "_5d", Ls_);
+        env.create<LatticePropagator>(getName() + "_5d", Ls_);
-        quark5d_ = env.getProp(getName() + "_5d");
+        quark5d_ = env.get<LatticePropagator>(getName() + "_5d");
     }
 }
 
@@ -93,13 +93,13 @@ void MQuark::execute(Environment &env)
     LatticePropagator *fullSource;
     LatticeFermion    source(env.getGrid(Ls_)), sol(env.getGrid(Ls_));
     
-    LOG(Message) << "computing quark propagator '" << getName() << "'"
+    LOG(Message) << "Computing quark propagator '" << getName() << "'"
                  << std::endl;
-    if (!env.isProp5d(par_.source))
+    if (!env.isLattice5d(par_.source))
     {
         if (Ls_ == 1)
         {
-            fullSource = env.getProp(par_.source);
+            fullSource = env.get<LatticePropagator>(par_.source);
         }
         else
         {
@@ -112,17 +112,17 @@ void MQuark::execute(Environment &env)
         {
             HADRON_ERROR("MQuark not implemented with 5D actions");
         }
-        else if (Ls_ != env.getPropLs(par_.source))
+        else if (Ls_ != env.getLatticeLs(par_.source))
         {
             HADRON_ERROR("MQuark not implemented with 5D actions");
         }
         else
         {
-            fullSource = env.getProp(par_.source);
+            fullSource = env.get<LatticePropagator>(par_.source);
         }
     }
     
-    LOG(Message) << "inverting using solver '" << par_.solver
+    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)

programs/Hadrons/SrcPoint.cc

@@ -62,8 +62,8 @@ std::vector<std::string> SrcPoint::getOutput(void)
 // allocation //////////////////////////////////////////////////////////////////
 void SrcPoint::allocate(Environment &env)
 {
-    env.createProp(getName());
+    env.create<LatticePropagator>(getName());
-    src_ = env.getProp(getName());
+    src_ = env.get<LatticePropagator>(getName());
 }
 
 // execution ///////////////////////////////////////////////////////////////////
@@ -72,7 +72,7 @@ void SrcPoint::execute(Environment &env)
     std::vector<int> position = strToVec<int>(par_.position);
     SpinColourMatrix id;
     
-    LOG(Message) << "creating point source at position [" << par_.position
+    LOG(Message) << "Creating point source at position [" << par_.position
                  << "]" << std::endl;
     id    = 1.;
     *src_ = zero;

programs/Hadrons/SrcZ2.cc

@@ -62,8 +62,8 @@ std::vector<std::string> SrcZ2::getOutput(void)
 // allocation //////////////////////////////////////////////////////////////////
 void SrcZ2::allocate(Environment &env)
 {
-    env.createProp(getName());
+    env.create<LatticePropagator>(getName());
-    src_ = env.getProp(getName());
+    src_ = env.get<LatticePropagator>(getName());
 }
 
 // execution ///////////////////////////////////////////////////////////////////
@@ -76,12 +76,12 @@ void SrcZ2::execute(Environment &env)
     
     if (par_.tA == par_.tB)
     {
-        LOG(Message) << "generating Z_2 wall source at t= " << par_.tA
+        LOG(Message) << "Generating Z_2 wall source at t= " << par_.tA
                      << std::endl;
     }
     else
     {
-        LOG(Message) << "generating Z_2 band for " << par_.tA << " <= t <= "
+        LOG(Message) << "Generating Z_2 band for " << par_.tA << " <= t <= "
                      << par_.tB << std::endl;
     }
     LatticeCoordinate(t, Tp);