From 52148463418c22fd114dee45c84b11e303010897 Mon Sep 17 00:00:00 2001
From: Guido Cossu <guido.cossu@ed.ac.uk>
Date: Thu, 22 Dec 2016 12:41:56 +0000
Subject: [PATCH] Adding a resource manager

---
 lib/Grid.h                               |   2 +
 lib/qcd/action/ActionBase.h              |  49 ++---
 lib/qcd/action/gauge/WilsonGaugeAction.h |  17 +-
 lib/qcd/hmc/GenericHMCrunner.h           | 257 +++++++++++------------
 lib/qcd/hmc/HMC.h                        |   6 +-
 lib/qcd/hmc/HMCModules.h                 |  96 +++++++++
 lib/qcd/hmc/HMCResourceManager.h         | 110 ++++++++++
 lib/qcd/hmc/integrators/Integrator.h     |  19 +-
 tests/hmc/Test_hmc_WilsonGauge_Binary.cc | 188 ++++++-----------
 tests/solver/Test_dwf_cg_prec_LsVec.cc   |   2 +-
 10 files changed, 435 insertions(+), 311 deletions(-)
 create mode 100644 lib/qcd/hmc/HMCModules.h
 create mode 100644 lib/qcd/hmc/HMCResourceManager.h

diff --git a/lib/Grid.h b/lib/Grid.h
index 78b5a6bd..7eaf8511 100644
--- a/lib/Grid.h
+++ b/lib/Grid.h
@@ -87,6 +87,8 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #include <Grid/qcd/hmc/NerscCheckpointer.h>
 #include <Grid/qcd/hmc/BinaryCheckpointer.h>
 #include <Grid/qcd/hmc/ILDGCheckpointer.h>
+#include <Grid/qcd/hmc/HMCModules.h>
+#include <Grid/qcd/hmc/HMCResourceManager.h>
 #include <Grid/qcd/hmc/HmcRunner.h>
 #include <Grid/qcd/hmc/GenericHMCrunner.h>
 
diff --git a/lib/qcd/action/ActionBase.h b/lib/qcd/action/ActionBase.h
index 371ccd24..6ff4bcb1 100644
--- a/lib/qcd/action/ActionBase.h
+++ b/lib/qcd/action/ActionBase.h
@@ -36,14 +36,13 @@ template <class GaugeField>
 class Action {
  public:
   bool is_smeared = false;
-  // Boundary conditions? // Heatbath?
-  virtual void refresh(const GaugeField& U,
-                       GridParallelRNG& pRNG) = 0;  // refresh pseudofermions
-  virtual RealD S(const GaugeField& U) = 0;         // evaluate the action
-  virtual void deriv(const GaugeField& U,
-                     GaugeField& dSdU) = 0;  // evaluate the action derivative
-  virtual std::string action_name() = 0; // return the action name
-  virtual ~Action(){};
+  // Heatbath?
+  virtual void refresh(const GaugeField& U, const GridParallelRNG& pRNG) = 0; // refresh pseudofermions
+  virtual RealD S(const GaugeField& U) = 0;                       // evaluate the action
+  virtual void deriv(const GaugeField& U, GaugeField& dSdU) = 0;  // evaluate the action derivative
+  virtual std::string action_name()    = 0;                       // return the action name
+  virtual std::string LogParameters()  = 0;  											// prints action parameters
+  virtual ~Action(){}
 };
 
 // Indexing of tuple types
@@ -60,32 +59,10 @@ struct Index<T, std::tuple<U, Types...>> {
   static const std::size_t value = 1 + Index<T, std::tuple<Types...>>::value;
 };
 
-/*
-template <class GaugeField>
-struct ActionLevel {
- public:
-  typedef Action<GaugeField>*
-      ActPtr;  // now force the same colours as the rest of the code
-
-  //Add supported representations here
-
-
-  unsigned int multiplier;
-
-  std::vector<ActPtr> actions;
-
-  ActionLevel(unsigned int mul = 1) : actions(0), multiplier(mul) {
-    assert(mul >= 1);
-  };
-
-  void push_back(ActPtr ptr) { actions.push_back(ptr); }
-};
-*/
-
 template <class Field, class Repr = NoHirep >
 struct ActionLevel {
  public:
-  unsigned int multiplier; 
+  unsigned int multiplier;
 
   // Fundamental repr actions separated because of the smearing
   typedef Action<Field>* ActPtr;
@@ -98,15 +75,13 @@ struct ActionLevel {
 
   std::vector<ActPtr>& actions;
 
-  ActionLevel(unsigned int mul = 1) : actions(std::get<0>(actions_hirep)), multiplier(mul) {
+  explicit ActionLevel(unsigned int mul = 1) : actions(std::get<0>(actions_hirep)), multiplier(mul) {
     // initialize the hirep vectors to zero.
-    //apply(this->resize, actions_hirep, 0); //need a working resize
+    // apply(this->resize, actions_hirep, 0); //need a working resize
     assert(mul >= 1);
-  };
-
-  //void push_back(ActPtr ptr) { actions.push_back(ptr); }
-
+  }
 
+  // void push_back(ActPtr ptr) { actions.push_back(ptr); }
 
   template < class GenField >
   void push_back(Action<GenField>* ptr) {
diff --git a/lib/qcd/action/gauge/WilsonGaugeAction.h b/lib/qcd/action/gauge/WilsonGaugeAction.h
index aa61abf2..fd281e20 100644
--- a/lib/qcd/action/gauge/WilsonGaugeAction.h
+++ b/lib/qcd/action/gauge/WilsonGaugeAction.h
@@ -43,18 +43,22 @@ class WilsonGaugeAction : public Action<typename Gimpl::GaugeField> {
  public:
   INHERIT_GIMPL_TYPES(Gimpl);
 
-  //      typedef LorentzScalar<GaugeField> GaugeLinkField;
-
  private:
   RealD beta;
 
  public:
-  WilsonGaugeAction(RealD b) : beta(b){};
+  explicit WilsonGaugeAction(RealD b) : beta(b){}
 
-  virtual std::string action_name(){return "WilsonGaugeAction";}
+  virtual std::string action_name() {return "WilsonGaugeAction";}
+
+  virtual std::string LogParameters(){
+  	std::stringstream sstream;
+  	sstream << GridLogMessage << "[WilsonGaugeAction] Beta: " << beta << std::endl;
+  	return sstream.str();
+  }
 
   virtual void refresh(const GaugeField &U,
-                       GridParallelRNG &pRNG){};  // noop as no pseudoferms
+                       const GridParallelRNG &pRNG){};  // noop as no pseudoferms
 
   virtual RealD S(const GaugeField &U) {
     RealD plaq = WilsonLoops<Gimpl>::avgPlaquette(U);
@@ -80,8 +84,7 @@ class WilsonGaugeAction : public Action<typename Gimpl::GaugeField> {
 
       PokeIndex<LorentzIndex>(dSdU, dSdU_mu, mu);
     }
-
-  };
+  }
 };
 
 
diff --git a/lib/qcd/hmc/GenericHMCrunner.h b/lib/qcd/hmc/GenericHMCrunner.h
index 2530559f..1d328f8e 100644
--- a/lib/qcd/hmc/GenericHMCrunner.h
+++ b/lib/qcd/hmc/GenericHMCrunner.h
@@ -30,176 +30,167 @@ with this program; if not, write to the Free Software Foundation, Inc.,
 #ifndef GRID_GENERIC_HMC_RUNNER
 #define GRID_GENERIC_HMC_RUNNER
 
+
+#include <unordered_map>
+
 namespace Grid {
 namespace QCD {
 
-// Virtual Class for HMC specific for gauge theories
-// implement a specific theory by defining the BuildTheAction
-template <class Implementation, class RepresentationsPolicy = NoHirep>
+template <class Implementation, 
+          template < typename, typename, typename > class Integrator, 
+          class RepresentationsPolicy = NoHirep >
 class BinaryHmcRunnerTemplate {
 public:
-    INHERIT_FIELD_TYPES(Implementation);
-    typedef Implementation ImplPolicy;
+ INHERIT_FIELD_TYPES(Implementation);
+	    typedef Implementation ImplPolicy; // visible from outside
+      template < typename S = NoSmearing<Implementation> > 
+      using IntegratorType = Integrator<Implementation,S,RepresentationsPolicy>;
 
-    enum StartType_t { ColdStart,
-                       HotStart,
-                       TepidStart,
-                       CheckpointStart };
+      enum StartType_t 
+      { 
+        ColdStart,
+        HotStart,
+        TepidStart,
+        CheckpointStart,
+        FilenameStart 
+      };
 
-    ActionSet<Field, RepresentationsPolicy> TheAction;
+      struct HMCPayload
+      {
+        StartType_t StartType;
+        HMCparameters Parameters;
 
-    // A vector of HmcObservable
-    // that can be injected from outside
-    std::vector<HmcObservable<typename Implementation::Field> *>
-        ObservablesList;
+        HMCPayload() { StartType = HotStart; }
+      };
 
-    IntegratorParameters MDparameters;
+     // These can be rationalised, some private
+	  HMCPayload Payload; // Parameters 
+	  HMCResourceManager Resources;
+	  IntegratorParameters MDparameters;
 
-    GridCartesian *        UGrid;
-    GridRedBlackCartesian *UrbGrid;
+	  ActionSet<Field, RepresentationsPolicy> TheAction;
+
+    // A vector of HmcObservable that can be injected from outside
+	  std::vector<HmcObservable<typename Implementation::Field> *> ObservablesList;
+
+	  //GridCartesian *        UGrid;
 
     // These two are unnecessary, eliminate
-    GridCartesian *        FGrid;
-    GridRedBlackCartesian *FrbGrid;
+	  // GridRedBlackCartesian *UrbGrid;
+	  // GridCartesian *        FGrid;
+	  // GridRedBlackCartesian *FrbGrid;
 
-    std::vector<int> SerialSeed;
-    std::vector<int> ParallelSeed;
+	  void ReadCommandLine(int argc, char ** argv) {
+	  	std::string arg;
 
-    void RNGSeeds(std::vector<int> S, std::vector<int> P) {
-      SerialSeed   = S;
-      ParallelSeed = P;
-    }
+	  	if (GridCmdOptionExists(argv, argv + argc, "--StartType")) {
+	  		arg = GridCmdOptionPayload(argv, argv + argc, "--StartType");
+	  		if (arg == "HotStart") {
+	  			Payload.StartType = HotStart;
+	  		} else if (arg == "ColdStart") {
+	  			Payload.StartType = ColdStart;
+	  		} else if (arg == "TepidStart") {
+	  			Payload.StartType = TepidStart;
+	  		} else if (arg == "CheckpointStart") {
+	  			Payload.StartType = CheckpointStart;
+	  		} else {
+	  			std::cout << GridLogError << "Unrecognized option in --StartType\n";
+	  			std::cout
+	  			<< GridLogError
+	  			<< "Valid [HotStart, ColdStart, TepidStart, CheckpointStart]\n";
+	  			assert(0);
+	  		}
+	  	}
 
-    virtual void BuildTheAction(int argc, char **argv) = 0;  // necessary?
+	  	if (GridCmdOptionExists(argv, argv + argc, "--StartTrajectory")) {
+	  		arg = GridCmdOptionPayload(argv, argv + argc, "--StartTrajectory");
+	  		std::vector<int> ivec(0);
+	  		GridCmdOptionIntVector(arg, ivec);
+	  		Payload.Parameters.StartTrajectory = ivec[0];
+	  	}
 
-    // A couple of wrapper classes
-    template <class IOCheckpointer>
-    void Run(int argc, char **argv, IOCheckpointer &Checkpoint)  {
-      NoSmearing<Implementation> S;
-      Runner(argc, argv, Checkpoint, S);
-    }
+	  	if (GridCmdOptionExists(argv, argv + argc, "--Trajectories")) {
+	  		arg = GridCmdOptionPayload(argv, argv + argc, "--Trajectories");
+	  		std::vector<int> ivec(0);
+	  		GridCmdOptionIntVector(arg, ivec);
+	  		Payload.Parameters.Trajectories = ivec[0];
+	  	}
 
-    template <class IOCheckpointer, class SmearingPolicy>
-    void Run(int argc, char **argv, IOCheckpointer &CP, SmearingPolicy &S) {
-      Runner(argc, argv, CP, S);
-    }
-    //////////////////////////////
+	  	if (GridCmdOptionExists(argv, argv + argc, "--Thermalizations")) {
+	  		arg = GridCmdOptionPayload(argv, argv + argc, "--Thermalizations");
+	  		std::vector<int> ivec(0);
+	  		GridCmdOptionIntVector(arg, ivec);
+	  		Payload.Parameters.NoMetropolisUntil = ivec[0];
+	  	}
 
-    
+	  }
 
+  // A couple of wrapper functions
+  template <class IOCheckpointer> void Run(IOCheckpointer &CP)  {
+    NoSmearing<Implementation> S;
+    Runner(CP, S);
+  }
 
-    template <class SmearingPolicy, class IOCheckpointer>
-    void Runner(int             argc,
-                char **         argv,
-                IOCheckpointer &Checkpoint,
-                SmearingPolicy &Smearing) {
-      StartType_t StartType = HotStart;
+  template <class IOCheckpointer, class SmearingPolicy> void Run(IOCheckpointer &CP, SmearingPolicy &S) {
+    Runner(CP, S);
+  }
 
-      std::string arg;
+  //////////////////////////////////////////////////////////////////
 
-      if (GridCmdOptionExists(argv, argv + argc, "--StartType")) {
-        arg = GridCmdOptionPayload(argv, argv + argc, "--StartType");
-        if (arg == "HotStart") {
-          StartType = HotStart;
-        } else if (arg == "ColdStart") {
-          StartType = ColdStart;
-        } else if (arg == "TepidStart") {
-          StartType = TepidStart;
-        } else if (arg == "CheckpointStart") {
-          StartType = CheckpointStart;
-        } else {
-          std::cout << GridLogError << "Unrecognized option in --StartType\n";
-          std::cout
-              << GridLogError
-              << "Valid [HotStart, ColdStart, TepidStart, CheckpointStart]\n";
-          assert(0);
-        }
-      }
+private:
+  template <class SmearingPolicy, class IOCheckpointer>
+    void Runner(IOCheckpointer &Checkpoint,	SmearingPolicy &Smearing) {
+  auto UGrid = Resources.GetCartesian();
+  Resources.AddRNGs();
+  Field U(UGrid);
 
-      int StartTraj = 0;
-      if (GridCmdOptionExists(argv, argv + argc, "--StartTrajectory")) {
-        arg = GridCmdOptionPayload(argv, argv + argc, "--StartTrajectory");
-        std::vector<int> ivec(0);
-        GridCmdOptionIntVector(arg, ivec);
-        StartTraj = ivec[0];
-      }
+  typedef IntegratorType<SmearingPolicy> TheIntegrator;
+  TheIntegrator MDynamics(UGrid, MDparameters, TheAction, Smearing);
 
-      int NumTraj = 1;
-      if (GridCmdOptionExists(argv, argv + argc, "--Trajectories")) {
-        arg = GridCmdOptionPayload(argv, argv + argc, "--Trajectories");
-        std::vector<int> ivec(0);
-        GridCmdOptionIntVector(arg, ivec);
-        NumTraj = ivec[0];
-      }
-
-      int NumThermalizations = 10;
-      if (GridCmdOptionExists(argv, argv + argc, "--Thermalizations")) {
-        arg = GridCmdOptionPayload(argv, argv + argc, "--Thermalizations");
-        std::vector<int> ivec(0);
-        GridCmdOptionIntVector(arg, ivec);
-        NumThermalizations = ivec[0];
-      }
-
-      GridSerialRNG   sRNG;
-      GridParallelRNG pRNG(UGrid);
-      Field           U(UGrid);
-
-
-      typedef MinimumNorm2<Implementation, SmearingPolicy, RepresentationsPolicy> IntegratorType;  // change here to change the algorithm
-      IntegratorType MDynamics(UGrid, MDparameters, TheAction, Smearing);
-
-      HMCparameters HMCpar;
-      HMCpar.StartTrajectory   = StartTraj;
-      HMCpar.Trajectories      = NumTraj;
-      HMCpar.NoMetropolisUntil = NumThermalizations;
-
-      if (StartType == HotStart) {
+  if (Payload.StartType == HotStart) {
         // Hot start
-        HMCpar.MetropolisTest = true;
-        sRNG.SeedFixedIntegers(SerialSeed);
-        pRNG.SeedFixedIntegers(ParallelSeed);
-        Implementation::HotConfiguration(pRNG, U);
-      } else if (StartType == ColdStart) {
+   Payload.Parameters.MetropolisTest = true;
+   Resources.SeedFixedIntegers();
+   Implementation::HotConfiguration(Resources.GetParallelRNG(), U);
+ } else if (Payload.StartType == ColdStart) {
         // Cold start
-        HMCpar.MetropolisTest = true;
-        sRNG.SeedFixedIntegers(SerialSeed);
-        pRNG.SeedFixedIntegers(ParallelSeed);
-        Implementation::ColdConfiguration(pRNG, U);
-      } else if (StartType == TepidStart) {
+   Payload.Parameters.MetropolisTest = true;
+   Resources.SeedFixedIntegers();
+   Implementation::ColdConfiguration(Resources.GetParallelRNG(), U);
+ } else if (Payload.StartType == TepidStart) {
         // Tepid start
-        HMCpar.MetropolisTest = true;
-        sRNG.SeedFixedIntegers(SerialSeed);
-        pRNG.SeedFixedIntegers(ParallelSeed);
-        Implementation::TepidConfiguration(pRNG, U);
-      } else if (StartType == CheckpointStart) {
-        HMCpar.MetropolisTest = true;
+   Payload.Parameters.MetropolisTest = true;
+   Resources.SeedFixedIntegers();
+   Implementation::TepidConfiguration(Resources.GetParallelRNG(), U);
+ } else if (Payload.StartType == CheckpointStart) {
+   Payload.Parameters.MetropolisTest = true;
         // CheckpointRestart
-        Checkpoint.CheckpointRestore(StartTraj, U, sRNG, pRNG);
-      }
+   Checkpoint.CheckpointRestore(Payload.Parameters.StartTrajectory, U, Resources.GetSerialRNG(), Resources.GetParallelRNG());
+ }
 
-      Smearing.set_Field(U);
+ Smearing.set_Field(U);
 
-      HybridMonteCarlo<IntegratorType> HMC(HMCpar, MDynamics, sRNG, pRNG, U);
+ HybridMonteCarlo<TheIntegrator> HMC(Payload.Parameters, MDynamics, Resources.GetSerialRNG(), Resources.GetParallelRNG(), U);
 
-      for (int obs = 0; obs < ObservablesList.size(); obs++)
-        HMC.AddObservable(ObservablesList[obs]);
+ for (int obs = 0; obs < ObservablesList.size(); obs++)
+   HMC.AddObservable(ObservablesList[obs]);
 
       // Run it
-      HMC.evolve();
-    }
+ HMC.evolve();
+}
 };
 
-  // These are for gauge fields
-  typedef BinaryHmcRunnerTemplate<PeriodicGimplR> BinaryHmcRunner;
-  typedef BinaryHmcRunnerTemplate<PeriodicGimplF> BinaryHmcRunnerF;
-  typedef BinaryHmcRunnerTemplate<PeriodicGimplD> BinaryHmcRunnerD;
+// These are for gauge fields, default integrator MinimumNorm2
+template <template <typename, typename, typename> class Integrator >  using BinaryHmcRunner = BinaryHmcRunnerTemplate<PeriodicGimplR, Integrator > ;
+template <template <typename, typename, typename> class Integrator >  using BinaryHmcRunnerF = BinaryHmcRunnerTemplate<PeriodicGimplF, Integrator > ;
+template <template <typename, typename, typename> class Integrator >  using BinaryHmcRunnerD = BinaryHmcRunnerTemplate<PeriodicGimplD, Integrator > ;
 
-  template <class RepresentationsPolicy>
-  using BinaryHmcRunnerTemplateHirep = BinaryHmcRunnerTemplate<PeriodicGimplR, RepresentationsPolicy>;
+template <class RepresentationsPolicy, template <typename, typename, typename> class Integrator >
+using BinaryHmcRunnerTemplateHirep = BinaryHmcRunnerTemplate<PeriodicGimplR, Integrator, RepresentationsPolicy>;
 
-  typedef BinaryHmcRunnerTemplate<ScalarImplR, ScalarFields>
-      ScalarBinaryHmcRunner;
+typedef BinaryHmcRunnerTemplate<ScalarImplR, MinimumNorm2, ScalarFields> ScalarBinaryHmcRunner;
 
 }  // namespace QCD
 }  // namespace Grid
-#endif 
+
+#endif  // GRID_GENERIC_HMC_RUNNER
diff --git a/lib/qcd/hmc/HMC.h b/lib/qcd/hmc/HMC.h
index e39faae0..6253ec77 100644
--- a/lib/qcd/hmc/HMC.h
+++ b/lib/qcd/hmc/HMC.h
@@ -55,7 +55,7 @@ struct HMCparameters {
     MetropolisTest = true;
     NoMetropolisUntil = 10;
     StartTrajectory = 0;
-    Trajectories = 200;
+    Trajectories = 10;
     /////////////////////////////////
   }
 
@@ -206,10 +206,10 @@ class HybridMonteCarlo {
 
     Params.print_parameters();
     TheIntegrator.print_parameters();
+    TheIntegrator.print_actions();
 
     // Actual updates (evolve a copy Ucopy then copy back eventually)
-    for (int traj = Params.StartTrajectory;
-         traj < Params.Trajectories + Params.StartTrajectory; ++traj) {
+    for (int traj = Params.StartTrajectory; traj < Params.Trajectories + Params.StartTrajectory; ++traj) {
       std::cout << GridLogMessage << "-- # Trajectory = " << traj << "\n";
       Ucopy = Ucur;
 
diff --git a/lib/qcd/hmc/HMCModules.h b/lib/qcd/hmc/HMCModules.h
new file mode 100644
index 00000000..84f586f6
--- /dev/null
+++ b/lib/qcd/hmc/HMCModules.h
@@ -0,0 +1,96 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/hmc/GenericHmcRunner.h
+
+Copyright (C) 2015
+Copyright (C) 2016
+
+Author: Guido Cossu <guido.cossu@ed.ac.uk>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
+#ifndef GRID_HMC_MODULES
+#define GRID_HMC_MODULES
+
+namespace Grid {
+namespace QCD {
+
+///////////////////////////////////////////////////
+// Modules
+class GridModule {
+ public:
+  GridCartesian* get_full() { return grid_.get(); }
+  GridRedBlackCartesian* get_rb() { return rbgrid_.get(); }
+
+  void set_full(GridCartesian* grid) { grid_.reset(grid); }
+  void set_rb(GridRedBlackCartesian* rbgrid) { rbgrid_.reset(rbgrid); }
+
+ protected:
+  std::unique_ptr<GridCartesian> grid_;
+  std::unique_ptr<GridRedBlackCartesian> rbgrid_;
+};
+
+// helpers
+class GridFourDimModule : public GridModule {
+ public:
+  GridFourDimModule() {
+    set_full(SpaceTimeGrid::makeFourDimGrid(
+        GridDefaultLatt(), GridDefaultSimd(4, vComplex::Nsimd()),
+        GridDefaultMpi()));
+    set_rb(SpaceTimeGrid::makeFourDimRedBlackGrid(grid_.get()));
+  }
+
+};
+
+class RNGModule{
+   // Random number generators
+   GridSerialRNG sRNG_;
+   std::unique_ptr<GridParallelRNG> pRNG_;
+   std::vector<int> SerialSeed_;
+   std::vector<int> ParallelSeed_;
+
+public:
+   void set_pRNG(GridParallelRNG* pRNG){
+      pRNG_.reset(pRNG);
+   }
+
+   void set_RNGSeeds(const std::vector<int> S, const std::vector<int> P) {
+    SerialSeed_   = S;
+    ParallelSeed_ = P;
+  }
+
+  GridSerialRNG& get_sRNG(){return sRNG_;}
+  GridParallelRNG& get_pRNG(){return *pRNG_.get();}
+  void seed(){
+    sRNG_.SeedFixedIntegers(SerialSeed_);
+    pRNG_->SeedFixedIntegers(ParallelSeed_);
+  }
+};
+
+
+
+
+
+
+}  // namespace QCD
+}  // namespace Grid
+
+#endif  // GRID_HMC_MODULES
\ No newline at end of file
diff --git a/lib/qcd/hmc/HMCResourceManager.h b/lib/qcd/hmc/HMCResourceManager.h
new file mode 100644
index 00000000..5cfd354d
--- /dev/null
+++ b/lib/qcd/hmc/HMCResourceManager.h
@@ -0,0 +1,110 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/hmc/GenericHmcRunner.h
+
+Copyright (C) 2015
+
+Author: paboyle <paboyle@ph.ed.ac.uk>
+Author: Guido Cossu <guido.cossu@ed.ac.uk>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+  See the full license in the file "LICENSE" in the top level distribution
+  directory
+  *************************************************************************************/
+/*  END LEGAL */
+#ifndef HMC_RESOURCE_MANAGER_H
+#define HMC_RESOURCE_MANAGER_H
+
+
+#include <unordered_map>
+
+namespace Grid {
+namespace QCD {
+
+  // HMC Resource manager
+  class HMCResourceManager {
+	// Storage for grid pairs (std + red-black)
+   std::unordered_map<std::string, GridModule> Grids;
+   RNGModule RNGs; 
+
+   bool have_RNG;
+
+ public:	
+   HMCResourceManager():have_RNG(false){}
+   void AddGrid(std::string s, GridModule& M){
+    // Check for name clashes
+    auto search = Grids.find(s);
+    if(search != Grids.end()) {
+      std::cout << GridLogError << "Grid with name \"" << search->first << "\" already present. Terminating\n" ;
+      exit(1);        
+    }
+    Grids[s] = std::move(M);
+  }
+
+  // Add a named grid set
+  void AddFourDimGrid(std::string s) {
+    GridFourDimModule Mod;
+    AddGrid(s,Mod);
+  }
+
+  GridCartesian* GetCartesian(std::string s="") {
+    if (s.empty()) s = Grids.begin()->first;
+    std::cout << GridLogDebug << "Getting cartesian grid from: "<< s << std::endl;
+    return Grids[s].get_full();
+  }
+
+  GridRedBlackCartesian* GetRBCartesian(std::string s="") {
+    if (s.empty()) s = Grids.begin()->first;
+    std::cout << GridLogDebug << "Getting rb-cartesian grid from: "<< s << std::endl;
+    return Grids[s].get_rb();
+  }
+
+  void AddRNGs(std::string s="") {
+		// Couple the RNGs to the GridModule tagged by s
+    // default is the first grid set
+    assert(Grids.size()>0 && !have_RNG );
+    if (s.empty()) s = Grids.begin()->first;
+    std::cout << GridLogDebug << "Adding RNG to grid: "<< s << std::endl;
+    RNGs.set_pRNG(new GridParallelRNG(GetCartesian(s)));
+    //pRNG.reset(new GridParallelRNG(GetCartesian(s)));
+    have_RNG = true;
+  }
+
+  void AddRNGSeeds(const std::vector<int> S, const std::vector<int> P) {
+    RNGs.set_RNGSeeds(S,P);
+    //SerialSeed   = S;
+    //ParallelSeed = P;
+  }
+
+  GridSerialRNG& GetSerialRNG() {return RNGs.get_sRNG();}
+  GridParallelRNG& GetParallelRNG() {
+    assert(have_RNG);
+    return RNGs.get_pRNG();
+  }
+  void SeedFixedIntegers() {
+    assert(have_RNG);
+    RNGs.seed();
+    //sRNG.SeedFixedIntegers(SerialSeed);
+    //pRNG->SeedFixedIntegers(ParallelSeed);
+  }
+};
+
+}
+}
+
+#endif // HMC_RESOURCE_MANAGER_H
\ No newline at end of file
diff --git a/lib/qcd/hmc/integrators/Integrator.h b/lib/qcd/hmc/integrators/Integrator.h
index 48ec746c..5d72a05f 100644
--- a/lib/qcd/hmc/integrators/Integrator.h
+++ b/lib/qcd/hmc/integrators/Integrator.h
@@ -57,7 +57,7 @@ struct IntegratorParameters {
 
 
   void print_parameters() {
-        std::cout << GridLogMessage << "[Integrator] Trajectory length  : " << trajL << std::endl;
+    std::cout << GridLogMessage << "[Integrator] Trajectory length  : " << trajL << std::endl;
     std::cout << GridLogMessage << "[Integrator] Number of MD steps : " << MDsteps << std::endl;
     std::cout << GridLogMessage << "[Integrator] Step size          : " << stepsize << std::endl;
   }
@@ -173,6 +173,20 @@ class Integrator {
     Params.print_parameters();
   }
 
+  void print_actions(){
+  	std::cout << GridLogMessage << ":::::::::::::::::::::::::::::::::::::::::" << std::endl;
+  	std::cout << GridLogMessage << "[Integrator] Action summary: "<<std::endl;
+  	for (int level = 0; level < as.size(); ++level) {
+  		std::cout << GridLogMessage << "[Integrator] ---- Level: "<< level << std::endl;
+  		for (int actionID = 0; actionID < as[level].actions.size(); ++actionID) {
+  			std::cout << GridLogMessage << "["<< as[level].actions.at(actionID)->action_name() << "] ID: " << actionID << std::endl;
+  			std::cout << as[level].actions.at(actionID)->LogParameters();
+  		}
+  	}
+  	std::cout << GridLogMessage << ":::::::::::::::::::::::::::::::::::::::::"<< std::endl;
+
+  }
+
   // to be used by the actionlevel class to iterate
   // over the representations
   struct _refresh {
@@ -189,8 +203,7 @@ class Integrator {
 
   // Initialization of momenta and actions
   void refresh(Field& U, GridParallelRNG& pRNG) {
-    //assert(P._grid == U._grid);
-    P.reset(U._grid);
+    assert(P._grid == U._grid);
     std::cout << GridLogIntegrator << "Integrator refresh\n";
     FieldImplementation::generate_momenta(P, pRNG);
  
diff --git a/tests/hmc/Test_hmc_WilsonGauge_Binary.cc b/tests/hmc/Test_hmc_WilsonGauge_Binary.cc
index d5228ce8..e2d05ee5 100644
--- a/tests/hmc/Test_hmc_WilsonGauge_Binary.cc
+++ b/tests/hmc/Test_hmc_WilsonGauge_Binary.cc
@@ -35,140 +35,52 @@ using namespace Grid;
 using namespace Grid::QCD;
 
 namespace Grid {
-namespace QCD {
+  namespace QCD {
 
-//Change here the type of reader
-typedef Grid::TextReader InputFileReader; 
+    //Change here the type of reader
+    typedef Grid::TextReader InputFileReader; 
 
 
-class HMCRunnerParameters : Serializable {
- public:
-  GRID_SERIALIZABLE_CLASS_MEMBERS(HMCRunnerParameters,
-                                  double, beta,
-                                  int, MDsteps,
-                                  double, TrajectorLength,
-                                  int, SaveInterval,
-                                  std::string, format,
-                                  std::string, conf_prefix,
-                                  std::string, rng_prefix,
-                                  std::string, serial_seeds,
-                                  std::string, parallel_seeds,
-                                  );
+    class HMCRunnerParameters : Serializable {
+    public:
+      GRID_SERIALIZABLE_CLASS_MEMBERS(HMCRunnerParameters,
+        double, beta,
+        int, MDsteps,
+        double, TrajectorLength,
+        int, SaveInterval,
+        std::string, format,
+        std::string, conf_prefix,
+        std::string, rng_prefix,
+        std::string, serial_seeds,
+        std::string, parallel_seeds,
+        );
 
-  HMCRunnerParameters() {}
-};
-
-
-// Derive from the BinaryHmcRunner (templated for gauge fields)
-class HmcRunner : public BinaryHmcRunner {
- public:
-    HMCRunnerParameters HMCPar;
-    void BuildTheAction(int argc, char **argv){}
-};
-/*
-
-    // eliminate arcg and argv from here
-    void BuildTheAction(int argc, char **argv)
-
-  {
-    // Typedefs to simplify notation
-    typedef WilsonGaugeActionR GaugeAction;
-    typedef WilsonImplR ImplPolicy;
-    typedef WilsonFermionR FermionAction;
-    typedef typename FermionAction::FermionField FermionField;
-
-    // this can be simplified too. MakeDefaultGrid(Nd)
-    UGrid = SpaceTimeGrid::makeFourDimGrid(
-        GridDefaultLatt(), 
-        GridDefaultSimd(Nd, vComplex::Nsimd()), 
-        GridDefaultMpi());
-    
-
-    // Gauge action
-    std::cout << GridLogMessage << "Beta: " << HMCPar.beta << std::endl;
-    GaugeAction Waction(HMCPar.beta);
-
-    // Collect actions
-    ActionLevel<Field> Level1(1);
-    Level1.push_back(&Waction);
-    TheAction.push_back(Level1);
-
-    // Add observables
-    // options for checkpointers
-    // this can be moved outside the BuildTheAction
-    //BinaryHmcCheckpointer
-    //ILDGHmcCheckpointer
-    //NerscHmcCheckpointer
-    NerscHmcCheckpointer<BinaryHmcRunner::ImplPolicy> Checkpoint(
-        HMCPar.conf_prefix, HMCPar.rng_prefix, HMCPar.SaveInterval, HMCPar.format);
-    // Can implement also a specific function in the hmcrunner
-    // AddCheckpoint (...) that takes the same parameters + a string/tag
-    // defining the type of the checkpointer
-    // with tags can be implemented by overloading and no ifs
-    // Then force all checkpoint to have few common functions
-    // return an object that is then passed to the Run function
-
-    PlaquetteLogger<BinaryHmcRunner::ImplPolicy> PlaqLog(
-        std::string("Plaquette"));
-    ObservablesList.push_back(&PlaqLog);
-    ObservablesList.push_back(&Checkpoint);
-
-    // This must run from here so that the grids are defined
-    Run(argc, argv, Checkpoint);  // no smearing
-  };
-};
-*/
-}
+      HMCRunnerParameters() {}
+    };
+  }
 }
 
 int main(int argc, char **argv) {
   Grid_init(&argc, &argv);
+   // Typedefs to simplify notation
+  typedef BinaryHmcRunner<MinimumNorm2> HMCWrapper;// Uses the default minimum norm
+  typedef WilsonGaugeActionR GaugeAction;
 
   int threads = GridThread::GetThreads();
-  std::cout << GridLogMessage << "Grid is setup to use " << threads
-            << " threads" << std::endl;
+  std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
 
-  HmcRunner TheHMC;
+  //////////////////////////////////////////////////////////////
+  // Input file section 
   // make input file name general
+  HMCRunnerParameters HMCPar;
   InputFileReader Reader("input.wilson_gauge.params");
-  read(Reader, "HMC", TheHMC.HMCPar);
-
-  std::cout << GridLogMessage << TheHMC.HMCPar << std::endl;
+  read(Reader, "HMC", HMCPar);
+  std::cout << GridLogMessage << HMCPar << std::endl;
 
   // Seeds for the random number generators
   // generalise
-  std::vector<int> SerSeed = strToVec<int>(TheHMC.HMCPar.serial_seeds);
-  std::vector<int> ParSeed = strToVec<int>(TheHMC.HMCPar.parallel_seeds);
-
-  TheHMC.RNGSeeds(SerSeed, ParSeed);
-
-  TheHMC.MDparameters.set(TheHMC.HMCPar.MDsteps, TheHMC.HMCPar.TrajectorLength);
-
-  //TheHMC.BuildTheAction(argc, argv);
-
-
-
-     // Typedefs to simplify notation
-    typedef WilsonGaugeActionR GaugeAction;
-    typedef WilsonImplR ImplPolicy;
-    typedef WilsonFermionR FermionAction;
-    typedef typename FermionAction::FermionField FermionField;
-
-    // this can be simplified too. MakeDefaultGrid(Nd)
-    TheHMC.UGrid = SpaceTimeGrid::makeFourDimGrid(
-        GridDefaultLatt(), 
-        GridDefaultSimd(Nd, vComplex::Nsimd()), 
-        GridDefaultMpi());
-    
-
-    // Gauge action
-    std::cout << GridLogMessage << "Beta: " << TheHMC.HMCPar.beta << std::endl;
-    GaugeAction Waction(TheHMC.HMCPar.beta);
-
-    // Collect actions
-    ActionLevel<BinaryHmcRunner::Field> Level1(1);
-    Level1.push_back(&Waction);
-    TheHMC.TheAction.push_back(Level1);
+  std::vector<int> SerSeed = strToVec<int>(HMCPar.serial_seeds);
+  std::vector<int> ParSeed = strToVec<int>(HMCPar.parallel_seeds);
 
     // Add observables
     // options for checkpointers
@@ -176,8 +88,8 @@ int main(int argc, char **argv) {
     //BinaryHmcCheckpointer
     //ILDGHmcCheckpointer
     //NerscHmcCheckpointer
-    NerscHmcCheckpointer<BinaryHmcRunner::ImplPolicy> Checkpoint(
-        TheHMC.HMCPar.conf_prefix, TheHMC.HMCPar.rng_prefix, TheHMC.HMCPar.SaveInterval, TheHMC.HMCPar.format);
+  BinaryHmcCheckpointer<HMCWrapper::ImplPolicy> Checkpoint(HMCPar.conf_prefix, HMCPar.rng_prefix, 
+                                                           HMCPar.SaveInterval, HMCPar.format);
     // Can implement also a specific function in the hmcrunner
     // AddCheckpoint (...) that takes the same parameters + a string/tag
     // defining the type of the checkpointer
@@ -185,16 +97,38 @@ int main(int argc, char **argv) {
     // Then force all checkpoint to have few common functions
     // return an object that is then passed to the Run function
 
-    PlaquetteLogger<BinaryHmcRunner::ImplPolicy> PlaqLog(
-        std::string("Plaquette"));
-    TheHMC.ObservablesList.push_back(&PlaqLog);
-    TheHMC.ObservablesList.push_back(&Checkpoint);
-
-    // This must run from here so that the grids are defined
-    TheHMC.Run(argc, argv, Checkpoint);  // no smearing
+  /////////////////////////////////////////////////////////////
+  HMCWrapper TheHMC;
+  TheHMC.Resources.AddFourDimGrid("gauge");
 
+  /////////////////////////////////////////////////////////////
+  // Collect actions, here use more encapsulation
+
+  // Gauge action
+  std::cout << GridLogMessage << "Beta: " << HMCPar.beta << std::endl;
+  GaugeAction Waction(HMCPar.beta);
+
+  ActionLevel<HMCWrapper::Field> Level1(1);
+  Level1.push_back(&Waction);
+  TheHMC.TheAction.push_back(Level1);
+  /////////////////////////////////////////////////////////////
+
+  // Construct observables 
+  PlaquetteLogger<HMCWrapper::ImplPolicy> PlaqLog("Plaquette");
+  TheHMC.ObservablesList.push_back(&PlaqLog);
+  TheHMC.ObservablesList.push_back(&Checkpoint);
+  //////////////////////////////////////////////
+
+  // Fill resources
+  TheHMC.Resources.AddRNGSeeds(SerSeed, ParSeed);
+  TheHMC.MDparameters.set(HMCPar.MDsteps, HMCPar.TrajectorLength);
+
+  // eventually smearing here
+  ////////////////////////////////////////////////////////////////
 
 
+  TheHMC.ReadCommandLine(argc, argv);
+  TheHMC.Run(Checkpoint);  // no smearing
 
   Grid_finalize();
 }
diff --git a/tests/solver/Test_dwf_cg_prec_LsVec.cc b/tests/solver/Test_dwf_cg_prec_LsVec.cc
index dfa0a194..9fedfb58 100644
--- a/tests/solver/Test_dwf_cg_prec_LsVec.cc
+++ b/tests/solver/Test_dwf_cg_prec_LsVec.cc
@@ -110,4 +110,4 @@ int main(int argc, char** argv) {
   Ddwf.Report();
 
   Grid_finalize();
-}  
+}