From 1e9317e5cf6c46f43390861e12fa011c63e0b2b4 Mon Sep 17 00:00:00 2001
From: neo <cossu@post.kek.jp>
Date: Mon, 6 Jul 2015 18:32:20 +0900
Subject: [PATCH] Simplifying HMC syntax for the final user

---
 lib/GridConfig.h                              |  4 +--
 lib/qcd/hmc/integrators/Integrator.h          | 28 +++++++++++--------
 .../hmc/integrators/Integrator_algorithm.h    | 14 +++++-----
 tests/Test_hmc_WilsonGauge.cc                 |  5 ++--
 4 files changed, 28 insertions(+), 23 deletions(-)

diff --git a/lib/GridConfig.h b/lib/GridConfig.h
index f3a1be6a..45d07e55 100644
--- a/lib/GridConfig.h
+++ b/lib/GridConfig.h
@@ -20,10 +20,10 @@
 #define GRID_COMMS_NONE 1
 
 /* GRID_DEFAULT_PRECISION is DOUBLE */
-/* #undef GRID_DEFAULT_PRECISION_DOUBLE */
+#define GRID_DEFAULT_PRECISION_DOUBLE 1
 
 /* GRID_DEFAULT_PRECISION is SINGLE */
-#define GRID_DEFAULT_PRECISION_SINGLE 1
+/* #undef GRID_DEFAULT_PRECISION_SINGLE */
 
 /* Support Altivec instructions */
 /* #undef HAVE_ALTIVEC */
diff --git a/lib/qcd/hmc/integrators/Integrator.h b/lib/qcd/hmc/integrators/Integrator.h
index 5eb319b2..522991ac 100644
--- a/lib/qcd/hmc/integrators/Integrator.h
+++ b/lib/qcd/hmc/integrators/Integrator.h
@@ -17,7 +17,15 @@ namespace Grid{
   namespace QCD{
 
     typedef Action<LatticeLorentzColourMatrix>*  ActPtr; // now force the same colours as the rest of the code
-    typedef std::vector<ActPtr> ActionLevel;
+    struct ActionLevel{
+      int multiplier;
+    public:
+      std::vector<ActPtr> actions;
+      explicit ActionLevel(int mul = 1):multiplier(mul){assert (mul > 0);};
+      void push_back(ActPtr ptr){
+	actions.push_back(ptr);
+      }
+    };
     typedef std::vector<ActionLevel> ActionSet;
     typedef std::vector<Observer*> ObserverList;
     
@@ -45,7 +53,6 @@ namespace Grid{
     private:
       IntegratorParameters Params;
       const ActionSet as;
-      const std::vector<int> Nrel; //relative step size per level
       std::unique_ptr<LatticeLorentzColourMatrix> P;
       GridParallelRNG pRNG;
       //ObserverList observers; // not yet
@@ -56,9 +63,9 @@ namespace Grid{
       void notify_observers();
 
       void update_P(LatticeLorentzColourMatrix&U, int level,double ep){
-	for(int a=0; a<as[level].size(); ++a){
+	for(int a=0; a<as[level].actions.size(); ++a){
 	  LatticeLorentzColourMatrix force(U._grid);
-	  as[level].at(a)->deriv(U,force);
+	  as[level].actions.at(a)->deriv(U,force);
 	  *P -= force*ep;
 	}
       }
@@ -84,9 +91,8 @@ namespace Grid{
 					     Integrator<IntegratorAlgorithm>* Integ);
     public:
     Integrator(GridBase* grid, IntegratorParameters Par,
-		 ActionSet& Aset, std::vector<int> Nrel_):
-      Params(Par),as(Aset),Nrel(Nrel_),P(new LatticeLorentzColourMatrix(grid)),pRNG(grid){
-	assert(as.size() == Nrel.size());
+		 ActionSet& Aset):
+      Params(Par),as(Aset),P(new LatticeLorentzColourMatrix(grid)),pRNG(grid){
 	pRNG.SeedRandomDevice();
       };
       
@@ -99,8 +105,8 @@ namespace Grid{
 
 	MDutils::generate_momenta(*P,pRNG);
 	for(int level=0; level< as.size(); ++level){
-	  for(int actionID=0; actionID<as.at(level).size(); ++actionID){
-	    as[level].at(actionID)->init(U, pRNG);
+	  for(int actionID=0; actionID<as[level].actions.size(); ++actionID){
+	    as[level].actions.at(actionID)->init(U, pRNG);
 	  }
 	}
       }
@@ -123,8 +129,8 @@ namespace Grid{
 
 	// Actions
 	for(int level=0; level<as.size(); ++level)
-	  for(int actionID=0; actionID<as.at(level).size(); ++actionID)
-	    H += as[level].at(actionID)->S(U);
+	  for(int actionID=0; actionID<as[level].actions.size(); ++actionID)
+	    H += as[level].actions.at(actionID)->S(U);
 	
 	return H;
       }
diff --git a/lib/qcd/hmc/integrators/Integrator_algorithm.h b/lib/qcd/hmc/integrators/Integrator_algorithm.h
index 0dce2978..c192b404 100644
--- a/lib/qcd/hmc/integrators/Integrator_algorithm.h
+++ b/lib/qcd/hmc/integrators/Integrator_algorithm.h
@@ -27,14 +27,14 @@ namespace Grid{
 	int fl = Integ->as.size() -1;
 	double eps = Integ->Params.stepsize;
 	
-	for(int l=0; l<=level; ++l) eps/= 2.0*Integ->Nrel[l];
+	for(int l=0; l<=level; ++l) eps/= 2.0*Integ->as[l].multiplier;
 	
-	int fin = Integ->Nrel[0];
-	for(int l=1; l<=level; ++l) fin*= 2.0*Integ->Nrel[l];
+	int fin = Integ->as[0].multiplier;
+	for(int l=1; l<=level; ++l) fin*= 2.0*Integ->as[l].multiplier;
 	fin = 3*Integ->Params.MDsteps*fin -1;
 	
 	
-	for(int e=0; e<Integ->Nrel[level]; ++e){
+	for(int e=0; e<Integ->as[level].multiplier; ++e){
 	  
 	  if(clock[level] == 0){    // initial half step 
 	    Integ->update_P(U,level,lambda*eps);
@@ -101,13 +101,13 @@ namespace Grid{
 	double eps = Integ->Params.stepsize;
 	
 	// Get current level step size
-	for(int l=0; l<=level; ++l) eps/= Integ->Nrel[l];
+	for(int l=0; l<=level; ++l) eps/= Integ->as[l].multiplier;
 	
 	int fin = 1;
-	for(int l=0; l<=level; ++l) fin*= Integ->Nrel[l];
+	for(int l=0; l<=level; ++l) fin*= Integ->as[l].multiplier;
 	fin = 2*Integ->Params.MDsteps*fin - 1;
 	
-	for(int e=0; e<Integ->Nrel[level]; ++e){
+	for(int e=0; e<Integ->as[level].multiplier; ++e){
 	  
 	  if(clock[level] == 0){    // initial half step
 	    Integ->update_P(U, level,eps/2.0);
diff --git a/tests/Test_hmc_WilsonGauge.cc b/tests/Test_hmc_WilsonGauge.cc
index 9b9bd388..d39ec9b9 100644
--- a/tests/Test_hmc_WilsonGauge.cc
+++ b/tests/Test_hmc_WilsonGauge.cc
@@ -39,10 +39,9 @@ int main (int argc, char ** argv)
   FullSet.push_back(Level1);
 
   // Create integrator
-  typedef MinimumNorm2  IntegratorAlgorithm;// change here to change the algorithm
+  typedef MinimumNorm2  IntegratorAlgorithm;// change here to modify the algorithm
   IntegratorParameters MDpar(12,5,1.0);
-  std::vector<int> rel ={1};
-  Integrator<IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet,rel);
+  Integrator<IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet);
 
   // Create HMC
   HMCparameters HMCpar;