Rearranging files in hmc

lib/Make.inc

@@ -1,4 +1,4 @@
 
-HFILES=./Cshift.h ./simd/Grid_avx.h ./simd/Grid_vector_types.h ./simd/Grid_sse4.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_vector_unops.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./Tensors.h ./Algorithms.h ./communicator/Communicator_base.h ./lattice/Lattice_rng.h ./lattice/Lattice_reduction.h ./lattice/Lattice_transfer.h ./lattice/Lattice_unary.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_comparison.h ./lattice/Lattice_overload.h ./lattice/Lattice_reality.h ./lattice/Lattice_local.h ./lattice/Lattice_conformable.h ./lattice/Lattice_where.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_ET.h ./lattice/Lattice_transpose.h ./lattice/Lattice_trace.h ./Stencil.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_exp.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_class.h ./tensors/Tensor_logical.h ./tensors/Tensor_transpose.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_reality.h ./tensors/Tensor_trace.h ./tensors/Tensor_index.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_outer.h ./tensors/Tensor_inner.h ./tensors/Tensor_traits.h ./tensors/Tensor_Ta.h ./tensors/Tensor_unary.h ./tensors/Tensor_determinant.h ./tensors/Tensor_arith.h ./tensors/Tensor_extract_merge.h ./Communicator.h ./Cartesian.h ./parallelIO/NerscIO.h ./qcd/QCD.h ./qcd/hmc/integrators/Integrator_base.h ./qcd/hmc/HMC.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/LinalgUtils.h ./qcd/utils/CovariantCshift.h ./qcd/utils/WilsonLoops.h ./qcd/action/ActionBase.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/spin/TwoSpinor.h ./qcd/spin/Dirac.h ./cshift/Cshift_common.h ./cshift/Cshift_none.h ./cshift/Cshift_mpi.h ./Simd.h ./GridConfig.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_red_black.h ./cartesian/Cartesian_full.h ./AlignedAllocator.h ./Lattice.h ./Old/Tensor_poke.h ./Old/Tensor_peek.h ./Threads.h ./Grid.h ./algorithms/Preconditioner.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/AdefGeneric.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/Zolotarev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Remez.h ./algorithms/LinearOperator.h ./algorithms/SparseMatrix.h ./algorithms/CoarsenedMatrix.h ./stencil/Lebesgue.h
+HFILES=./Cshift.h ./simd/Grid_avx.h ./simd/Grid_vector_types.h ./simd/Grid_sse4.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_vector_unops.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./Tensors.h ./Algorithms.h ./communicator/Communicator_base.h ./lattice/Lattice_rng.h ./lattice/Lattice_reduction.h ./lattice/Lattice_transfer.h ./lattice/Lattice_unary.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_comparison.h ./lattice/Lattice_overload.h ./lattice/Lattice_reality.h ./lattice/Lattice_local.h ./lattice/Lattice_conformable.h ./lattice/Lattice_where.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_ET.h ./lattice/Lattice_transpose.h ./lattice/Lattice_trace.h ./Stencil.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_exp.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_class.h ./tensors/Tensor_logical.h ./tensors/Tensor_transpose.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_reality.h ./tensors/Tensor_trace.h ./tensors/Tensor_index.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_outer.h ./tensors/Tensor_inner.h ./tensors/Tensor_traits.h ./tensors/Tensor_Ta.h ./tensors/Tensor_unary.h ./tensors/Tensor_determinant.h ./tensors/Tensor_arith.h ./tensors/Tensor_extract_merge.h ./Communicator.h ./Cartesian.h ./parallelIO/NerscIO.h ./qcd/QCD.h ./qcd/hmc/integrators/Integrator.h ./qcd/hmc/integrators/Integrator_algorithm.h ./qcd/hmc/HMC.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/LinalgUtils.h ./qcd/utils/CovariantCshift.h ./qcd/utils/WilsonLoops.h ./qcd/action/ActionBase.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/spin/TwoSpinor.h ./qcd/spin/Dirac.h ./cshift/Cshift_common.h ./cshift/Cshift_none.h ./cshift/Cshift_mpi.h ./Simd.h ./GridConfig.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_red_black.h ./cartesian/Cartesian_full.h ./AlignedAllocator.h ./Lattice.h ./Old/Tensor_poke.h ./Old/Tensor_peek.h ./Threads.h ./Grid.h ./algorithms/Preconditioner.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/AdefGeneric.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/Zolotarev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Remez.h ./algorithms/LinearOperator.h ./algorithms/SparseMatrix.h ./algorithms/CoarsenedMatrix.h ./stencil/Lebesgue.h
 
-CCFILES=./qcd/hmc/integrators/Integrator_base.cc ./qcd/hmc/HMC.cc ./qcd/utils/SpaceTimeGrid.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/spin/Dirac.cc ./GridInit.cc ./algorithms/approx/MultiShiftFunction.cc ./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc
+CCFILES=./qcd/hmc/integrators/Integrator.cc ./qcd/hmc/HMC.cc ./qcd/utils/SpaceTimeGrid.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/spin/Dirac.cc ./GridInit.cc ./algorithms/approx/MultiShiftFunction.cc ./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc

lib/qcd/QCD.h

@@ -413,7 +413,8 @@ namespace QCD {
 #include <qcd/utils/WilsonLoops.h>
 #include <qcd/utils/SUn.h>
 #include <qcd/action/Actions.h>
-#include <qcd/hmc/integrators/Integrator_base.h>
+#include <qcd/hmc/integrators/Integrator.h>
+#include <qcd/hmc/integrators/Integrator_algorithm.h>
 #include <qcd/hmc/HMC.h>
 
 

lib/qcd/hmc/HMC.h

@@ -25,11 +25,11 @@ namespace Grid{
       HMCparameters();
     };
     
-    template <class IntegType> 
+    template <class Algorithm> 
     class HybridMonteCarlo{
       const HMCparameters Params;
       GridSerialRNG sRNG;
-      Integrator<IntegType>& MD;
+      Integrator<Algorithm>& MD;
       
       bool metropolis_test(const RealD DeltaH){
 	RealD rn_test;
@@ -51,23 +51,21 @@ namespace Grid{
 
       RealD evolve_step(LatticeLorentzColourMatrix& U){
 
-	MD.init(U);     // set U and initialize P and phi's 
+	MD.init(U); // set U and initialize P and phi's 
-	RealD H0 = MD.S(U);     // current state            
+	RealD H0 = MD.S(U); // initial state action  
 	std::cout<<"Total H before = "<< H0 << "\n";
       
-	MD.integrate(U,0);
+	MD.integrate(U);
       
-	RealD H1 = MD.S(U);     // updated state            
+	RealD H1 = MD.S(U); // updated state action            
 	std::cout<<"Total H after = "<< H1 << "\n";
       
 	return (H1-H0);
       }
-
-      
       
     public:
     HybridMonteCarlo(HMCparameters Pms, 
-		     Integrator<IntegType>& MolDyn):
+		     Integrator<Algorithm>& MolDyn):
       Params(Pms),MD(MolDyn){
 	//FIXME
 
@@ -87,7 +85,6 @@ namespace Grid{
 	
 	  DeltaH = evolve_step(Uin);
 	  std::cout<< " dH = "<< DeltaH << "\n";
-	
 	}
 
 	// Actual updates (evolve a copy Ucopy then copy back eventually)

lib/qcd/hmc/integrators/Integrator.h

@@ -0,0 +1,148 @@
+//--------------------------------------------------------------------
+/*! @file Integrator.h
+ * @brief Declaration of classes for the Molecular Dynamics integrator
+ *
+ * @author Guido Cossu
+ */
+//--------------------------------------------------------------------
+
+#ifndef INTEGRATOR_INCLUDED
+#define INTEGRATOR_INCLUDED
+
+class Observer;
+
+#include <memory>
+
+namespace Grid{
+  namespace QCD{
+
+    typedef Action<LatticeLorentzColourMatrix>*  ActPtr; // now force the same colours as the rest of the code
+    typedef std::vector<ActPtr> ActionLevel;
+    typedef std::vector<ActionLevel> ActionSet;
+    typedef std::vector<Observer*> ObserverList;
+    
+    struct IntegratorParameters{
+      int Nexp;
+      int MDsteps;  // number of outer steps
+      RealD trajL;  // trajectory length 
+      RealD stepsize;
+
+      IntegratorParameters(int Nexp_,
+			   int MDsteps_, 
+			   RealD trajL_):
+      Nexp(Nexp_),MDsteps(MDsteps_),trajL(trajL_),stepsize(trajL/MDsteps){};
+    };
+
+
+    namespace MDutils{
+      void generate_momenta(LatticeLorentzColourMatrix&,GridParallelRNG&);
+      void generate_momenta_su3(LatticeLorentzColourMatrix&,GridParallelRNG&);
+    }
+
+    /*! @brief Class for Molecular Dynamics management */   
+    template< class IntegratorAlgorithm >
+    class Integrator{
+    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
+     
+      IntegratorAlgorithm TheIntegrator;
+
+      void register_observers();
+      void notify_observers();
+
+      void update_P(LatticeLorentzColourMatrix&U, int level,double ep){
+	for(int a=0; a<as[level].size(); ++a){
+	  LatticeLorentzColourMatrix force(U._grid);
+	  as[level].at(a)->deriv(U,force);
+	  *P -= force*ep;
+	}
+      }
+
+
+      void update_U(LatticeLorentzColourMatrix&U, double ep){
+	//rewrite exponential to deal automatically  with the lorentz index?
+	LatticeColourMatrix Umu(U._grid);
+	LatticeColourMatrix Pmu(U._grid);
+	for (int mu = 0; mu < Nd; mu++){
+	  Umu=peekLorentz(U, mu);
+	  Pmu=peekLorentz(*P, mu);
+	  Umu = expMat(Pmu, ep, Params.Nexp)*Umu;
+	  pokeLorentz(U, Umu, mu);
+	}
+
+      }
+      
+
+      
+      friend void IntegratorAlgorithm::step (LatticeLorentzColourMatrix& U, 
+					     int level, std::vector<int>& clock,
+					     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());
+	pRNG.SeedRandomDevice();
+      };
+      
+      ~Integrator(){}
+
+
+      //Initialization of momenta and actions
+      void init(LatticeLorentzColourMatrix& U){
+	std::cout<< "Integrator init\n";
+
+	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);
+	  }
+	}
+      }
+
+      
+      // Calculate action
+      RealD S(LatticeLorentzColourMatrix& U){
+	LatticeComplex Hloc(U._grid);
+	Hloc = zero;
+	// Momenta
+	for (int mu=0; mu <Nd; mu++){
+	  LatticeColourMatrix Pmu = peekLorentz(*P, mu);
+	  Hloc -= trace(Pmu*Pmu);
+	}
+	Complex Hsum = sum(Hloc);
+	
+	RealD H = Hsum.real();
+
+	std::cout << "H_p = "<< H << "\n";
+
+	// 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);
+	
+	return H;
+      }
+
+      void integrate(LatticeLorentzColourMatrix& U){
+	std::vector<int> clock;
+	clock.resize(as.size(),0);
+	for(int step=0; step< Params.MDsteps; ++step)   // MD step
+	  TheIntegrator.step(U,0,clock, (this));
+      }
+    };
+    
+
+
+
+
+
+    
+  }
+}
+#endif//INTEGRATOR_INCLUDED

lib/qcd/hmc/integrators/Integrator_algorithm.h

@@ -0,0 +1,152 @@
+//--------------------------------------------------------------------
+/*! @file Integrator_algorithm.h
+ * @brief Declaration of classes for the Molecular Dynamics algorithms
+ *
+ * @author Guido Cossu
+ */
+//--------------------------------------------------------------------
+
+#ifndef INTEGRATOR_ALG_INCLUDED
+#define INTEGRATOR_ALG_INCLUDED
+
+
+namespace Grid{
+  namespace QCD{
+
+
+    class MinimumNorm2{
+      const double lambda = 0.1931833275037836;
+    public:
+      void step (LatticeLorentzColourMatrix& U, 
+		 int level, std::vector<int>& clock,
+		 Integrator<MinimumNorm2>* Integ){
+	// level  : current level
+	// fl     : final level
+	// eps    : current step size
+
+	int fl = Integ->as.size() -1;
+	double eps = Integ->Params.stepsize;
+	
+	for(int l=0; l<=level; ++l) eps/= 2.0*Integ->Nrel[l];
+	
+	int fin = Integ->Nrel[0];
+	for(int l=1; l<=level; ++l) fin*= 2.0*Integ->Nrel[l];
+	fin = 3*Integ->Params.MDsteps*fin -1;
+	
+	
+	for(int e=0; e<Integ->Nrel[level]; ++e){
+	  
+	  if(clock[level] == 0){    // initial half step 
+	    Integ->update_P(U,level,lambda*eps);
+	    ++clock[level];
+	    for(int l=0; l<level;++l) std::cout<<"   ";
+	    std::cout<<"P "<< clock[level] <<std::endl;
+	  }
+	  
+	  if(level == fl){          // lowest level 
+	    Integ->update_U(U,0.5*eps);
+	    
+	    for(int l=0; l<level;++l) std::cout<<"   ";
+	    std::cout<<"U "<< (clock[level]+1) <<std::endl;
+	  }else{                 // recursive function call 
+	    step(U,level+1,clock, Integ);
+	  }
+	  
+	  Integ->update_P(U,level,(1.0-2.0*lambda)*eps);
+	  ++clock[level];
+	  for(int l=0; l<level;++l) std::cout<<"   ";
+	  std::cout<<"P "<< (clock[level]) <<std::endl;
+	  
+	  if(level == fl){          // lowest level 
+	    Integ->update_U(U,0.5*eps);
+	    
+	    for(int l=0; l<level;++l) std::cout<<"   ";
+	    std::cout<<"U "<< (clock[level]+1) <<std::endl;
+	  }else{                 // recursive function call 
+	    step(U,level+1,clock, Integ);
+	  }    
+	  
+	  
+	  if(clock[level] == fin){  // final half step
+	    Integ->update_P(U,level,lambda*eps);
+	    
+	    ++clock[level];
+	    for(int l=0; l<level;++l) std::cout<<"   ";
+	    std::cout<<"P "<< clock[level] <<std::endl;
+	  }else{                  // bulk step
+	    Integ->update_P(U,level,lambda*2.0*eps);
+	    
+	    clock[level]+=2;
+	    for(int l=0; l<level;++l) std::cout<<"   ";
+	    std::cout<<"P "<< clock[level] <<std::endl;
+	  }
+	}
+	
+	
+	
+      }
+      
+    };
+    
+    class LeapFrog{
+    public:
+      void step (LatticeLorentzColourMatrix& U, 
+		 int level, std::vector<int>& clock,
+		 Integrator<LeapFrog>* Integ){
+	// level  : current level
+	// fl     : final level
+	// eps    : current step size
+	
+	int fl = Integ->as.size() -1;
+	double eps = Integ->Params.stepsize;
+	
+	// Get current level step size
+	for(int l=0; l<=level; ++l) eps/= Integ->Nrel[l];
+	
+	int fin = 1;
+	for(int l=0; l<=level; ++l) fin*= Integ->Nrel[l];
+	fin = 2*Integ->Params.MDsteps*fin - 1;
+	
+	for(int e=0; e<Integ->Nrel[level]; ++e){
+	  
+	  if(clock[level] == 0){    // initial half step
+	    Integ->update_P(U, level,eps/2.0);
+	    ++clock[level];
+	    for(int l=0; l<level;++l) std::cout<<"   ";
+	    std::cout<<"P "<< 0.5*clock[level] <<std::endl;
+	  }
+	  if(level == fl){          // lowest level
+	    Integ->update_U(U, eps);
+	    for(int l=0; l<level;++l) std::cout<<"   ";
+	    std::cout<<"U "<< 0.5*(clock[level]+1) <<std::endl;
+	  }else{                 // recursive function call
+	    step(U, level+1,clock, Integ);
+	  }
+	  if(clock[level] == fin){  // final half step
+	    Integ->update_P(U, level,eps/2.0);
+	    
+	    ++clock[level];
+	    for(int l=0; l<level;++l) std::cout<<"   ";
+	    std::cout<<"P "<< 0.5*clock[level] <<std::endl;
+	  }else{                  // bulk step
+	    Integ->update_P(U, level,eps);
+	    
+	    clock[level]+=2;
+	    for(int l=0; l<level;++l) std::cout<<"   ";
+	    std::cout<<"P "<< 0.5*clock[level] <<std::endl;
+	  }
+	}
+
+
+
+
+      }
+    };
+
+
+  }
+}
+
+
+
+#endif//INTEGRATOR_INCLUDED

lib/qcd/hmc/integrators/Integrator_base.h

@@ -1,280 +0,0 @@
-//--------------------------------------------------------------------
-/*! @file Integrator_base.h
- * @brief Declaration of classes for the abstract Molecular Dynamics integrator
- *
- * @author Guido Cossu
- */
-//--------------------------------------------------------------------
-
-#ifndef INTEGRATOR_INCLUDED
-#define INTEGRATOR_INCLUDED
-
-#include <memory>
-
-class Observer;
-
-
-/*! @brief Abstract base class for Molecular Dynamics management */
-
-namespace Grid{
-  namespace QCD{
-
-    typedef Action<LatticeLorentzColourMatrix>*  ActPtr; // now force the same size as the rest of the code
-    typedef std::vector<ActPtr> ActionLevel;
-    typedef std::vector<ActionLevel> ActionSet;
-    typedef std::vector<Observer*> ObserverList;
-    
-    struct IntegratorParameters{
-      int Nexp;
-      int MDsteps;  // number of outer steps
-      RealD trajL;  // trajectory length 
-      RealD stepsize;
-
-      IntegratorParameters(int Nexp_,
-			   int MDsteps_, 
-			   RealD trajL_):
-	Nexp(Nexp_),MDsteps(MDsteps_),trajL(trajL_),stepsize(trajL/MDsteps){};
-    };
-
-
-    namespace MDutils{
-      void generate_momenta(LatticeLorentzColourMatrix&,GridParallelRNG&);
-      void generate_momenta_su3(LatticeLorentzColourMatrix&,GridParallelRNG&);
-    }
-
-    
-    template< class IntegratorPolicy >
-    class Integrator{
-    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
-     
-      IntegratorPolicy TheIntegrator;
-
-      void register_observers();
-      void notify_observers();
-
-      void update_P(LatticeLorentzColourMatrix&U, int level,double ep){
-	for(int a=0; a<as[level].size(); ++a){
-	  LatticeLorentzColourMatrix force(U._grid);
-	  as[level].at(a)->deriv(U,force);
-	  *P -= force*ep;
-	}
-      }
-
-
-      void update_U(LatticeLorentzColourMatrix&U, double ep){
-	//rewrite exponential to deal automatically  with the lorentz index?
-	LatticeColourMatrix Umu(U._grid);
-	LatticeColourMatrix Pmu(U._grid);
-	for (int mu = 0; mu < Nd; mu++){
-	  Umu=peekLorentz(U, mu);
-	  Pmu=peekLorentz(*P, mu);
-	  Umu = expMat(Pmu, ep, Params.Nexp)*Umu;
-	  pokeLorentz(U, Umu, mu);
-	}
-
-      }
-      
-
-      
-      friend void IntegratorPolicy::step (LatticeLorentzColourMatrix& U, 
-				     int level, std::vector<int>& clock,
-				     Integrator<IntegratorPolicy>* 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());
-	pRNG.SeedRandomDevice();
-      };
-      
-      ~Integrator(){}
-
-
-      //Initialization of momenta and actions
-      void init(LatticeLorentzColourMatrix& U){
-	std::cout<< "Integrator init\n";
-
-	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);
-	  }
-	}
-      }
-
-      
-      // Calculate action
-      RealD S(LatticeLorentzColourMatrix& U){
-	LatticeComplex Hloc(U._grid);
-	Hloc = zero;
-	// Momenta
-	for (int mu=0; mu <Nd; mu++){
-	  LatticeColourMatrix Pmu = peekLorentz(*P, mu);
-	  Hloc -= trace(Pmu*Pmu);
-	}
-	Complex Hsum = sum(Hloc);
-	
-	RealD H = Hsum.real();
-
-	std::cout << "H_p = "<< H << "\n";
-
-	// 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);
-	
-	return H;
-      }
-      
-      
-
-      void integrate(LatticeLorentzColourMatrix& U, int level){
-	std::vector<int> clock;
-	clock.resize(as.size(),0);
-	for(int step=0; step< Params.MDsteps; ++step)   // MD step
-	  TheIntegrator.step(U,0,clock, (this));
-      }
-    };
-    
-
-    class MinimumNorm2{
-      const double lambda = 0.1931833275037836;
-    public:
-      void step (LatticeLorentzColourMatrix& U, 
-		 int level, std::vector<int>& clock,
-		 Integrator<MinimumNorm2>* Integ){
-	// level  : current level
-	// fl  : final level
-	// eps : current step size
-
-	int fl = Integ->as.size() -1;
-	double eps = Integ->Params.stepsize;
-	
-	for(int l=0; l<=level; ++l) eps/= 2.0*Integ->Nrel[l];
-	
-	int fin = Integ->Nrel[0];
-	for(int l=1; l<=level; ++l) fin*= 2.0*Integ->Nrel[l];
-	fin = 3*Integ->Params.MDsteps*fin -1;
-	
-	
-	for(int e=0; e<Integ->Nrel[level]; ++e){
-	  
-	  if(clock[level] == 0){    // initial half step 
-	    Integ->update_P(U,level,lambda*eps);
-	    ++clock[level];
-	    for(int l=0; l<level;++l) std::cout<<"   ";
-	    std::cout<<"P "<< clock[level] <<std::endl;
-	  }
-	  
-	  if(level == fl){          // lowest level 
-	    Integ->update_U(U,0.5*eps);
-	    
-	    for(int l=0; l<level;++l) std::cout<<"   ";
-	    std::cout<<"U "<< (clock[level]+1) <<std::endl;
-	  }else{                 // recursive function call 
-	    step(U,level+1,clock, Integ);
-	  }
-	  
-	  Integ->update_P(U,level,(1.0-2.0*lambda)*eps);
-	  ++clock[level];
-	  for(int l=0; l<level;++l) std::cout<<"   ";
-	  std::cout<<"P "<< (clock[level]) <<std::endl;
-	  
-	  if(level == fl){          // lowest level 
-	    Integ->update_U(U,0.5*eps);
-	    
-	    for(int l=0; l<level;++l) std::cout<<"   ";
-	    std::cout<<"U "<< (clock[level]+1) <<std::endl;
-	  }else{                 // recursive function call 
-	    step(U,level+1,clock, Integ);
-	  }    
-	  
-	  
-	  if(clock[level] == fin){  // final half step
-	    Integ->update_P(U,level,lambda*eps);
-	    
-	    ++clock[level];
-	    for(int l=0; l<level;++l) std::cout<<"   ";
-	    std::cout<<"P "<< clock[level] <<std::endl;
-	  }else{                  // bulk step
-	    Integ->update_P(U,level,lambda*2.0*eps);
-	    
-	    clock[level]+=2;
-	    for(int l=0; l<level;++l) std::cout<<"   ";
-	    std::cout<<"P "<< clock[level] <<std::endl;
-	  }
-	}
-	
-	
-	
-      }
-      
-    };
-    
-    class LeapFrog{
-    public:
-      void step (LatticeLorentzColourMatrix& U, 
-		 int level, std::vector<int>& clock,
-		 Integrator<LeapFrog>* Integ){
-	// fl  : final level
-	// eps : current step size
-	
-	int fl = Integ->as.size() -1;
-	double eps = Integ->Params.stepsize;
-	
-	// Get current level step size
-	for(int l=0; l<=level; ++l) eps/= Integ->Nrel[l];
-	
-	int fin = 1;
-	for(int l=0; l<=level; ++l) fin*= Integ->Nrel[l];
-	fin = 2*Integ->Params.MDsteps*fin - 1;
-	
-	for(int e=0; e<Integ->Nrel[level]; ++e){
-	  
-	  if(clock[level] == 0){    // initial half step
-	    Integ->update_P(U, level,eps/2.0);
-	    ++clock[level];
-	    for(int l=0; l<level;++l) std::cout<<"   ";
-	    std::cout<<"P "<< 0.5*clock[level] <<std::endl;
-	  }
-	  if(level == fl){          // lowest level
-	    Integ->update_U(U, eps);
-	    for(int l=0; l<level;++l) std::cout<<"   ";
-	    std::cout<<"U "<< 0.5*(clock[level]+1) <<std::endl;
-	  }else{                 // recursive function call
-	    step(U, level+1,clock, Integ);
-	  }
-	  if(clock[level] == fin){  // final half step
-	    Integ->update_P(U, level,eps/2.0);
-	    
-	    ++clock[level];
-	    for(int l=0; l<level;++l) std::cout<<"   ";
-	    std::cout<<"P "<< 0.5*clock[level] <<std::endl;
-	  }else{                  // bulk step
-	    Integ->update_P(U, level,eps);
-	    
-	    clock[level]+=2;
-	    for(int l=0; l<level;++l) std::cout<<"   ";
-	    std::cout<<"P "<< 0.5*clock[level] <<std::endl;
-	  }
-	}
-
-
-
-
-      }
-    };
-
-
-
-
-    
-  }
-}
-#endif//INTEGRATOR_INCLUDED

tests/Test_hmc_WilsonGauge.cc

@@ -40,7 +40,7 @@ int main (int argc, char ** argv)
 
   // Create integrator
   typedef MinimumNorm2  IntegratorAlgorithm;// change here to change the algorithm
-  IntegratorParameters MDpar(12,30,1.0);
+  IntegratorParameters MDpar(12,5,1.0);
   std::vector<int> rel ={1};
   Integrator<IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet,rel);