Global edit on HMC sector -- making GaugeField a template parameter and

preparing to pass integrator, smearing, bc's as policy classes to hmc.

Propose to unify "integrator" and integrator algorithm in a base/derived
way to override step. Want to read through ForceGradient to ensure
that abstraction covers the force gradient case.

lib/Make.inc

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

lib/parallelIO/NerscIO.h

@@ -66,11 +66,11 @@ inline void NerscGrid(GridBase *grid,NerscField &header)
     header.boundary[d] = std::string("PERIODIC");
   }
 }
-template<class GaugeMatrix,class GaugeField>
+template<class GaugeField>
 inline void NerscStatistics(GaugeField & data,NerscField &header)
 {
-  header.link_trace=Grid::QCD::WilsonLoops<GaugeMatrix,GaugeField>::linkTrace(data);
+  header.link_trace=Grid::QCD::WilsonLoops<GaugeField>::linkTrace(data);
-  header.plaquette =Grid::QCD::WilsonLoops<GaugeMatrix,GaugeField>::avgPlaquette(data);
+  header.plaquette =Grid::QCD::WilsonLoops<GaugeField>::avgPlaquette(data);
 }
 
 inline void NerscMachineCharacteristics(NerscField &header)
@@ -298,7 +298,6 @@ template<class vsimd>
 static inline void readConfiguration(Lattice<iLorentzColourMatrix<vsimd> > &Umu,NerscField& header,std::string file)
 {
   typedef Lattice<iLorentzColourMatrix<vsimd> > GaugeField;
-  typedef Lattice<iColourMatrix<vsimd> > GaugeMatrix;
 
   GridBase *grid = Umu._grid;
   int offset = readHeader(file,Umu._grid,header);
@@ -341,7 +340,7 @@ static inline void readConfiguration(Lattice<iLorentzColourMatrix<vsimd> > &Umu,
     assert(0);
   }
 
-  NerscStatistics<GaugeMatrix,GaugeField>(Umu,clone);
+  NerscStatistics<GaugeField>(Umu,clone);
 
   assert(fabs(clone.plaquette -header.plaquette ) < 1.0e-5 );
   assert(fabs(clone.link_trace-header.link_trace) < 1.0e-6 );
@@ -354,7 +353,7 @@ template<class vsimd>
 static inline void writeConfiguration(Lattice<iLorentzColourMatrix<vsimd> > &Umu,std::string file, int two_row,int bits32)
 {
   typedef Lattice<iLorentzColourMatrix<vsimd> > GaugeField;
-  typedef Lattice<iColourMatrix<vsimd> > GaugeMatrix;
+
   typedef iLorentzColourMatrix<vsimd> vobj;
   typedef typename vobj::scalar_object sobj;
 
@@ -372,7 +371,7 @@ static inline void writeConfiguration(Lattice<iLorentzColourMatrix<vsimd> > &Umu
   GridBase *grid = Umu._grid;
 
   NerscGrid(grid,header);
-  NerscStatistics<GaugeMatrix,GaugeField>(Umu,header);
+  NerscStatistics<GaugeField>(Umu,header);
   NerscMachineCharacteristics(header);
 
   uint32_t csum;

lib/qcd/QCD.h

@@ -62,7 +62,6 @@ namespace QCD {
     template<typename vtype> using iGparitySpinColourVector       = iVector<iVector<iVector<vtype, Nc>, Nhs>, Ngp >;
     template<typename vtype> using iGparityHalfSpinColourVector   = iVector<iVector<iVector<vtype, Nc>, Nhs>, Ngp >;
 
-
     // Spin matrix
     typedef iSpinMatrix<Complex  >          SpinMatrix;
     typedef iSpinMatrix<ComplexF >          SpinMatrixF;
@@ -154,7 +153,7 @@ namespace QCD {
     typedef iHalfSpinColourVector<vComplexD> vHalfSpinColourVectorD;
     
     // singlets
-    typedef iSinglet<Complex >         TComplex;    // FIXME This is painful. Tensor singlet complex type.
+    typedef iSinglet<Complex >         TComplex;     // FIXME This is painful. Tensor singlet complex type.
     typedef iSinglet<ComplexF>         TComplexF;    // FIXME This is painful. Tensor singlet complex type.
     typedef iSinglet<ComplexD>         TComplexD;    // FIXME This is painful. Tensor singlet complex type.
 
@@ -162,7 +161,7 @@ namespace QCD {
     typedef iSinglet<vComplexF>        vTComplexF;   // what if we don't know the tensor structure
     typedef iSinglet<vComplexD>        vTComplexD;   // what if we don't know the tensor structure
 
-    typedef iSinglet<Real >            TReal;       // Shouldn't need these; can I make it work without?
+    typedef iSinglet<Real >            TReal;        // Shouldn't need these; can I make it work without?
     typedef iSinglet<RealF>            TRealF;       // Shouldn't need these; can I make it work without?
     typedef iSinglet<RealD>            TRealD;       // Shouldn't need these; can I make it work without?
 
@@ -251,6 +250,8 @@ namespace QCD {
     typedef LatticeDoubleStoredColourMatrixF       LatticeDoubledGaugeFieldF;
     typedef LatticeDoubleStoredColourMatrixD       LatticeDoubledGaugeFieldD;
 
+    template<class GF> using LorentzScalar = Lattice<iScalar<typename GF::vector_object::element> >;
+
     // Uhgg... typing this hurt  ;)
     // (my keyboard got burning hot when I typed this, must be the anti-Fermion)
     typedef Lattice<vColourVector>          LatticeStaggeredFermion;    

lib/qcd/action/Actions.h

@@ -16,10 +16,6 @@
 #include <qcd/action/ActionBase.h>
 #include <qcd/action/ActionParams.h>
 
-////////////////////////////////////////////
-// Gauge Actions
-////////////////////////////////////////////
-#include <qcd/action/gauge/WilsonGaugeAction.h>
 
 
 ////////////////////////////////////////////
@@ -31,6 +27,17 @@
 #include <qcd/action/fermion/WilsonKernels.h>        //used by all wilson type fermions
 
 
+////////////////////////////////////////////
+// Gauge Actions
+////////////////////////////////////////////
+#include <qcd/action/gauge/WilsonGaugeAction.h>
+namespace Grid {
+namespace QCD {
+typedef WilsonGaugeAction<LatticeGaugeField>     WilsonGaugeActionR;
+typedef WilsonGaugeAction<LatticeGaugeFieldF>    WilsonGaugeActionF;
+typedef WilsonGaugeAction<LatticeGaugeFieldD>    WilsonGaugeActionD;
+}}
+
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 // Explicit explicit template instantiation is still required in the .cc files
 //

lib/qcd/action/gauge/WilsonGaugeAction.h

@@ -7,8 +7,12 @@ namespace Grid{
     ////////////////////////////////////////////////////////////////////////
     // Wilson Gauge Action .. should I template the Nc etc..
     ////////////////////////////////////////////////////////////////////////
-    template<class GaugeField, class MatrixField>
+    template<class GaugeField>
-      class WilsonGaugeAction : public Action<GaugeField> {
+    class WilsonGaugeAction : public Action<GaugeField> {
+    public:
+
+      typedef LorentzScalar<GaugeField> GaugeLinkField;
+
     private:
       RealD beta;
     public:
@@ -17,23 +21,26 @@ namespace Grid{
       virtual void init(const GaugeField &U, GridParallelRNG& pRNG) {};
       
       virtual RealD S(const GaugeField &U) {
-	RealD plaq = WilsonLoops<MatrixField,GaugeField>::avgPlaquette(U);
+	RealD plaq = WilsonLoops<GaugeField>::avgPlaquette(U);
 	std::cout<<GridLogMessage << "Plaq : "<<plaq << "\n";
 	RealD vol = U._grid->gSites();
 	RealD action=beta*(1.0 -plaq)*(Nd*(Nd-1.0))*vol*0.5;
 	std::cout << GridLogMessage << "WilsonGauge action "<<action<<std::endl;
 	return action;
       };
+
       virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
 	//not optimal implementation FIXME
 	//extend Ta to include Lorentz indexes
 	RealD factor = 0.5*beta/RealD(Nc);
-	MatrixField dSdU_mu(U._grid);
+	GaugeLinkField Umu(U._grid);
-	MatrixField Umu(U._grid);
+	GaugeLinkField dSdU_mu(U._grid);
 	for (int mu=0; mu < Nd; mu++){
+
 	  Umu = PeekIndex<LorentzIndex>(U,mu);
+
 	  // Staple in direction mu
-	  WilsonLoops<MatrixField,GaugeField>::Staple(dSdU_mu,U,mu);
+	  WilsonLoops<GaugeField>::Staple(dSdU_mu,U,mu);
 	  dSdU_mu = Ta(Umu*adj(dSdU_mu))*factor;
 	  pokeLorentz(dSdU, dSdU_mu, mu);
 	}

lib/qcd/hmc/HMC.h

@@ -15,6 +15,7 @@
 namespace Grid{
   namespace QCD{
     
+
     struct HMCparameters{
       Integer Nsweeps; /* @brief Number of sweeps in this run */
       Integer TotalSweeps; /* @brief If provided, the total number of sweeps */
@@ -26,14 +27,17 @@ namespace Grid{
       HMCparameters();
     };
     
-    template <class Algorithm> 
+    //    template <class GaugeField, class Integrator, class Smearer, class Boundary> 
-    class HybridMonteCarlo{
+    template <class GaugeField, class Algorithm>
+    class HybridMonteCarlo {
+    private:
 
       const HMCparameters Params;
-
+      
-      GridSerialRNG sRNG; // Fixme: need a RNG management strategy.
+      GridSerialRNG   &sRNG;   // Fixme: need a RNG management strategy.
-
+      GridParallelRNG &pRNG; // Fixme: need a RNG management strategy.
-      Integrator<Algorithm>& MD;
+      typedef Integrator<GaugeField,Algorithm> IntegratorType;
+      IntegratorType &TheIntegrator;
 
       /////////////////////////////////////////////////////////
       // Metropolis step
@@ -63,16 +67,20 @@ namespace Grid{
       /////////////////////////////////////////////////////////
       // Evolution
       /////////////////////////////////////////////////////////
-      RealD evolve_step(LatticeGaugeField& U){
+      RealD evolve_step(GaugeField& U){
-	MD.init(U); // set U and initialize P and phi's 
+
+	TheIntegrator.init(U); // set U and initialize P and phi's 
+
+	RealD H0 = TheIntegrator.S(U); // initial state action  
 
-	RealD H0 = MD.S(U); // initial state action  
 	std::cout<<GridLogMessage<<"Total H before = "<< H0 << "\n";
 
-	MD.integrate(U);
+	TheIntegrator.integrate(U);
       
-	RealD H1 = MD.S(U); // updated state action            
+	RealD H1 = TheIntegrator.S(U); // updated state action            
+
 	std::cout<<GridLogMessage<<"Total H after = "<< H1 << "\n";
+
 	return (H1-H0);
       }
       
@@ -81,16 +89,18 @@ namespace Grid{
       /////////////////////////////////////////
       // Constructor
       /////////////////////////////////////////
-      HybridMonteCarlo(HMCparameters Pms,  Integrator<Algorithm>& MolDyn): Params(Pms),MD(MolDyn) {
+     HybridMonteCarlo(HMCparameters Pms,  IntegratorType &_Int, GridSerialRNG &_sRNG, GridParallelRNG &_pRNG ) :
-
+        Params(Pms), 
-	//FIXME...  initialize RNGs also with seed ; RNG management strategy
+	TheIntegrator(_Int), 
-	sRNG.SeedRandomDevice();
+	sRNG(_sRNG),
-
+	pRNG(_pRNG)
+      {
       }
       ~HybridMonteCarlo(){};
 
 
-      void evolve(LatticeGaugeField& Uin){
+      void evolve(GaugeField& Uin){
+
 	Real DeltaH;
 	
 	// Thermalizations
@@ -102,7 +112,7 @@ namespace Grid{
 	}
 
 	// Actual updates (evolve a copy Ucopy then copy back eventually)
-	LatticeGaugeField Ucopy(Uin._grid);
+	GaugeField Ucopy(Uin._grid);
 	for(int iter=Params.StartingConfig; iter < Params.Nsweeps+Params.StartingConfig; ++iter){
 	  std::cout<<GridLogMessage << "-- # Sweep = "<< iter <<  "\n";
 	  

lib/qcd/hmc/integrators/Integrator.cc

@@ -1,28 +0,0 @@
-/*! 
-  @file Integrator_base.cc
-  @brief utilities for MD including funcs to generate initial HMC momentum
- */
-
-#include <Grid.h>
-
-namespace Grid{
-  namespace QCD{
-
-    void MDutils::generate_momenta(LatticeLorentzColourMatrix& P,GridParallelRNG& pRNG){
-      // for future support of different groups
-      MDutils::generate_momenta_su3(P, pRNG);
-    }
-
-    void MDutils::generate_momenta_su3(LatticeLorentzColourMatrix& P,GridParallelRNG& pRNG){
-      LatticeColourMatrix Pmu(P._grid);
-      Pmu = zero;
-      for(int mu=0;mu<Nd;mu++){
-	SU3::GaussianLieAlgebraMatrix(pRNG, Pmu);
-	PokeIndex<LorentzIndex>(P, Pmu, mu);
-      }
-      
-    }
-
-
-  }
-}

lib/qcd/hmc/integrators/Integrator.h

@@ -10,27 +10,15 @@
 #ifndef INTEGRATOR_INCLUDED
 #define INTEGRATOR_INCLUDED
 
-class Observer;
+//class Observer;
 
 #include <memory>
 
 namespace Grid{
   namespace QCD{
 
-    typedef Action<LatticeGaugeField>*  ActPtr; // now force the same colours as the rest of the code
-    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;
-    
     struct IntegratorParameters{
+
       int Nexp;
       int MDsteps;  // number of outer steps
       RealD trajL;  // trajectory length 
@@ -39,94 +27,133 @@ namespace Grid{
       IntegratorParameters(int Nexp_,
 			   int MDsteps_, 
 			   RealD trajL_):
-      Nexp(Nexp_),MDsteps(MDsteps_),trajL(trajL_),stepsize(trajL/MDsteps){};
+        Nexp(Nexp_),
+	MDsteps(MDsteps_),
+	trajL(trajL_),
+	stepsize(trajL/MDsteps)
+        {
+	  // empty body constructor
+	};
+
     };
 
-    namespace MDutils{
+    // Should match any legal (SU(n)) gauge field
-      void generate_momenta(LatticeGaugeField&,GridParallelRNG&);
+    // Need to use this template to match Ncol to pass to SU<N> class
-      void generate_momenta_su3(LatticeGaugeField&,GridParallelRNG&);
+    template<int Ncol,class vec> void generate_momenta(Lattice< iVector< iScalar< iMatrix<vec,Ncol> >, Nd> > & P,GridParallelRNG& pRNG){
+      typedef Lattice< iScalar< iScalar< iMatrix<vec,Ncol> > > > GaugeLinkField;
+      GaugeLinkField Pmu(P._grid);
+      Pmu = zero;
+      for(int mu=0;mu<Nd;mu++){
+	SU<Ncol>::GaussianLieAlgebraMatrix(pRNG, Pmu);
+	PokeIndex<LorentzIndex>(P, Pmu, mu);
+      }
     }
 
+    template<class GaugeField> struct ActionLevel{
+      public:
+      
+	typedef Action<GaugeField>*  ActPtr; // now force the same colours as the rest of the code
+
+	int multiplier;
+
+	std::vector<ActPtr> actions;
+
+        ActionLevel(int mul = 1) : multiplier(mul) {
+	  assert (mul > 0);
+	};
+
+	void push_back(ActPtr ptr){
+	  actions.push_back(ptr);
+	}
+    };
+
+    template<class GaugeField> using ActionSet = std::vector<ActionLevel< GaugeField > >;
+
     /*! @brief Class for Molecular Dynamics management */   
-    template< class IntegratorAlgorithm >
+    template<class GaugeField, class Algorithm >
-    class Integrator{
+    class Integrator : public Algorithm {
+
     private:
 
-      int levels;
-      std::vector<double> t_P;
-      double t_U;
-
       IntegratorParameters Params;
-      const ActionSet as;
+
-      std::unique_ptr<LatticeGaugeField> P;
+      GridParallelRNG pRNG;        // Store this somewhere more sensible and pass as reference
-      GridParallelRNG pRNG;
+      const ActionSet<GaugeField> as;
+
+      int levels;              //
+      double t_U;              // Track time passing on each level and for U and for P
+      std::vector<double> t_P; //
+
+      GaugeField P;  // is a pointer really necessary?
 
       //ObserverList observers; // not yet
-     
+      //      typedef std::vector<Observer*> ObserverList;
-      IntegratorAlgorithm TheIntegrator;
+      //      void register_observers();
+      //      void notify_observers();
 
-      void register_observers();
 
-      void notify_observers();
+      void update_P(GaugeField&U, int level,double ep){
-
-      void update_P(LatticeGaugeField&U, int level,double ep){
 	t_P[level]+=ep;
 	for(int a=0; a<as[level].actions.size(); ++a){
-	  LatticeGaugeField force(U._grid);
+	  GaugeField force(U._grid);
 	  as[level].actions.at(a)->deriv(U,force);
-	  *P -= force*ep;
+	  P = P - force*ep;
 	}
+
 	std::cout<<GridLogMessage;
 	for(int l=0; l<level;++l) std::cout<<"   ";	    
 	std::cout<<"["<<level<<"] P " << " dt "<< ep <<" : t_P "<< t_P[level] <<std::endl;
+
       }
 
-      void update_U(LatticeGaugeField&U, double ep){
+      void update_U(GaugeField&U, double ep){
 	//rewrite exponential to deal automatically  with the lorentz index?
-	LatticeColourMatrix Umu(U._grid);
+	//	GaugeLinkField Umu(U._grid);
-	LatticeColourMatrix Pmu(U._grid);
+	//	GaugeLinkField Pmu(U._grid);
 	for (int mu = 0; mu < Nd; mu++){
-	  Umu=PeekIndex<LorentzIndex>(U, mu);
+	  auto Umu=PeekIndex<LorentzIndex>(U, mu);
-	  Pmu=PeekIndex<LorentzIndex>(*P, mu);
+	  auto Pmu=PeekIndex<LorentzIndex>(P, mu);
 	  Umu = expMat(Pmu, ep, Params.Nexp)*Umu;
 	  PokeIndex<LorentzIndex>(U, Umu, mu);
 	}
+
 	t_U+=ep;
 	int fl = levels-1;
 	std::cout<<GridLogMessage<<"   ";
 	for(int l=0; l<fl;++l) std::cout<<"   ";	    
 	std::cout<<"["<<fl<<"] U " << " dt "<< ep <<" : t_U "<< t_U <<std::endl;
+
       }
       
-      friend void IntegratorAlgorithm::step (LatticeGaugeField& U, 
+      friend void Algorithm::step (GaugeField& U, 
-					     int level, std::vector<int>& clock,
+				   int level, 
-					     Integrator<IntegratorAlgorithm>* Integ);
+				   std::vector<int>& clock,
+				   Integrator<GaugeField,Algorithm>* Integ);
+
     public:
 
       Integrator(GridBase* grid, 
 		 IntegratorParameters Par,
-		 ActionSet& Aset):
+		 ActionSet<GaugeField> & Aset):
           Params(Par),
     	  as(Aset),
-	  P(new LatticeGaugeField(grid)),
+	  P(grid),
 	  pRNG(grid),
 	  levels(Aset.size())
       {
-	
+	std::vector<int> seeds({1,2,3,4,5}); // Fixme; Pass it the RNG as a ref
-	std::vector<int> seeds({1,2,3,4,5});
 	pRNG.SeedFixedIntegers(seeds);
 
 	t_P.resize(levels,0.0);
 	t_U=0.0;
-
       };
       
       ~Integrator(){}
 
       //Initialization of momenta and actions
-      void init(LatticeGaugeField& U){
+      void init(GaugeField& U){
 	std::cout<<GridLogMessage<< "Integrator init\n";
-	MDutils::generate_momenta(*P,pRNG);
+	generate_momenta(P,pRNG);
 	for(int level=0; level< as.size(); ++level){
 	  for(int actionID=0; actionID<as[level].actions.size(); ++actionID){
 	    as[level].actions.at(actionID)->init(U, pRNG);
@@ -135,12 +162,12 @@ namespace Grid{
       }
 
       // Calculate action
-      RealD S(LatticeGaugeField& U){
+      RealD S(GaugeField& U){
-	LatticeComplex Hloc(U._grid);
+
-	Hloc = zero;
+	LatticeComplex Hloc(U._grid);	Hloc = zero;
 	// Momenta
 	for (int mu=0; mu <Nd; mu++){
-	  LatticeColourMatrix Pmu = peekLorentz(*P, mu);
+	  auto Pmu = PeekIndex<LorentzIndex>(P, mu);
 	  Hloc -= trace(Pmu*Pmu);
 	}
 	Complex Hsum = sum(Hloc);
@@ -162,16 +189,23 @@ namespace Grid{
 	return H;
       }
 
-      void integrate(LatticeGaugeField& U){
+      void integrate(GaugeField& U){
+
 	std::vector<int> clock;
+
 	clock.resize(as.size(),0);
+
+	// All the clock stuff is removed if we pass first, last to the step down the way
 	for(int step=0; step< Params.MDsteps; ++step){   // MD step
-	  TheIntegrator.step(U,0,clock, (this));
+	  Algorithm::step(U,0,clock, (this));
 	}
+
+	// Check the clocks all match
 	for(int level=0; level<as.size(); ++level){
-	  assert(fabs(t_U - t_P[level])<1.0e-3);
+	  assert(fabs(t_U - t_P[level])<1.0e-4);
 	  std::cout<<GridLogMessage<<" times["<<level<<"]= "<<t_P[level]<< " " << t_U <<std::endl;
 	}	
+
       }
     };
     

lib/qcd/hmc/integrators/Integrator_algorithm.h

@@ -12,111 +12,13 @@
 namespace Grid{
   namespace QCD{
 
-    /*
+   /* PAB:
-Chroma: Recursive min norm
+    *
-00132     Real dtau = traj_length / Real(n_steps);
+    * Recursive leapfrog; explanation of nested stepping
-00133     Real lambda_dt = dtau*lambda;
+    *
-00134     Real dtauby2 = dtau / Real(2);
-00135     Real one_minus_2lambda_dt = (Real(1)-Real(2)*lambda)*dtau;
-00136     Real two_lambda_dt = lambda_dt*Real(2);
-00137 
-00138     // Its sts so:
-00139     expSdt(s, lambda_dt); 
-00140     for(int i=0; i < n_steps-1; i++) {  // N-1 full steps
-00141       // Roll the exp(lambda_dt T) here and start
-00142       // Next iter into one
-00143       subIntegrator(s, dtauby2);         <--- either leapU or next integrator
-00144       expSdt(s, one_minus_2lambda_dt);  
-00145       subIntegrator(s, dtauby2);         <--- either leapU or next integrator
-00146       expSdt(s, two_lambda_dt); 
-00147     }
-00148     // Last step, can't roll the first and last exp(lambda_dt T) 
-00149     // together.
-00150     subIntegrator(s, dtauby2);
-00151     expSdt(s, one_minus_2lambda_dt);
-00152     subIntegrator(s, dtauby2);
-00153     expSdt(s, lambda_dt);
-
-    * 
-    */
-
-
-    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->as[l].multiplier;
-	
-	// which is final half step
-	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;
-
-	int multiplier = Integ->as[level].multiplier;
-	for(int e=0; e<multiplier; ++e){       // steps per step
-
-	  if(clock[level] == 0){    // initial half step 
-	    Integ->update_P(U,level,lambda*eps);   ++clock[level];
-	  }
-	  
-	  if(level == fl){          // lowest level 
-	    Integ->update_U(U,0.5*eps);
-	  }else{                 // recursive function call 
-	    step(U,level+1,clock, Integ);
-	  }
-	  
-	  Integ->update_P(U,level,(1.0-2.0*lambda)*eps); ++clock[level];
-	  
-	  if(level == fl){          // lowest level 
-	    Integ->update_U(U,0.5*eps);
-	  }else{                 // recursive function call 
-	    step(U,level+1,clock, Integ);
-	  }    
-	  
-	  // Handle the final half step
-	  std::cout << GridLogMessage << " P[" <<level<<"] clock = "<<clock[level]<<"/"<<fin<<std::endl;
-	  int mm = (clock[level]==fin) ? 1 : 2;
-	  Integ->update_P(U,level,lambda*eps*mm); clock[level]+=mm;
-
-	}
-		
-      }
-      
-    };
-
-
-
-    /*
-
-Chroma: Recursive leapfrog
-
-
-00124     Real dtau = traj_length / n_steps;
-00125     Real dtauby2 = dtau/2;
-00126 
-00127     // Its sts so:
-00128     expSdt(s, dtauby2);  // First half step
-00129     for(int i=0; i < n_steps-1; i++) {  // N-1 full steps
-00130       subIntegrator(s, dtau);  <--- either leapU or next integrator
-00131       expSdt(s, dtau);
-00132     }
-00133     subIntegrator(s, dtau);     // Last Full Step
-00134     expSdt(s, dtauby2);  // Last Half Step
-
     * Nested 1:4; units in dt for top level integrator
-    * CHROMA                           GUIDO
+    *
+    * CHROMA                           IroIro
     *   0        1                      0              
     *  P 1/2                           P 1/2
     *          P 1/16                                  P1/16
@@ -129,7 +31,7 @@ Chroma: Recursive leapfrog
     *                 U 1/8                                 U1/8
     *          P 1/16                                  P1/8
     *  P 1                             P 1
-    *          P 1/16                    * skipped
+    *          P 1/16                    * skipped --- avoids revaluating force
     *                 U 1/8                                 U1/8
     *          P 1/8                                   P1/8
     *                 U 1/8                                 U1/8
@@ -159,13 +61,16 @@ Chroma: Recursive leapfrog
     *                 U 1/8                                 U1/8
     *          P 1/16                                  P1/16
     *  P 1/2                            P 1/2
-    * Total
     */    
-    class LeapFrog{
+
+    template<class GaugeField> class LeapFrog {
     public:
-      void step (LatticeLorentzColourMatrix& U, 
+
+      typedef LeapFrog<GaugeField> Algorithm;
+
+      void step (GaugeField& U, 
 		 int level, std::vector<int>& clock,
-		 Integrator<LeapFrog>* Integ){
+		 Integrator<GaugeField,Algorithm> * Integ){
 
 	// level  : current level
 	// fl     : final level
@@ -186,13 +91,7 @@ Chroma: Recursive leapfrog
 
 	  int first_step,last_step;
 
-	  if ( level==0 ) {
+	  first_step = (clock[level]==0);
-	    first_step = (clock[level]==0);
-	    last_step  = (clock[level]==fin);
-	  } else {
-	    first_step = (e==0);
-	    last_step  = (e==multiplier-1);
-	  }
 
 	  if(first_step){    // initial half step
 	    Integ->update_P(U, level,eps/2.0);	    ++clock[level];
@@ -204,8 +103,7 @@ Chroma: Recursive leapfrog
 	    step(U, level+1,clock, Integ);
 	  }
 
-	  // Handle the final half step
+	  last_step  = (clock[level]==fin);
-	  std::cout << GridLogMessage << " P[" <<level<<"] clock = "<<clock[level]<<"/"<<fin<<std::endl;
 	  int mm = last_step ? 1 : 2;
 	  Integ->update_P(U, level,mm*eps/2.0);	    
 	  clock[level]+=mm;
@@ -214,6 +112,65 @@ Chroma: Recursive leapfrog
       }
     };
 
+    template<class GaugeField> class MinimumNorm2 {
+    public:
+      typedef MinimumNorm2<GaugeField> Algorithm;
+
+    private:
+      const double lambda = 0.1931833275037836;
+    public:
+
+      void step (GaugeField& U, 
+		 int level, std::vector<int>& clock,
+		 Integrator<GaugeField,Algorithm>* 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->as[l].multiplier;
+	
+	// which is final half step
+	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;
+
+	int multiplier = Integ->as[level].multiplier;
+	for(int e=0; e<multiplier; ++e){       // steps per step
+
+	  int first_step,last_step;
+
+	  first_step = (clock[level]==0);
+
+	  if(first_step){    // initial half step 
+	    Integ->update_P(U,level,lambda*eps);   ++clock[level];
+	  }
+	  
+	  if(level == fl){          // lowest level 
+	    Integ->update_U(U,0.5*eps);
+	  }else{                 // recursive function call 
+	    step(U,level+1,clock, Integ);
+	  }
+	  
+	  Integ->update_P(U,level,(1.0-2.0*lambda)*eps); ++clock[level];
+	  
+	  if(level == fl){          // lowest level 
+	    Integ->update_U(U,0.5*eps);
+	  }else{                 // recursive function call 
+	    step(U,level+1,clock, Integ);
+	  }    
+	  
+	  last_step  = (clock[level]==fin);
+	  int mm = (last_step) ? 1 : 2;
+	  Integ->update_P(U,level,lambda*eps*mm); clock[level]+=mm;
+
+	}
+      }
+    };
 
   }
 }

lib/qcd/utils/WilsonLoops.h

@@ -4,9 +4,12 @@ namespace Grid {
 namespace QCD {
 
 // Common wilson loop observables
-template<class GaugeMat,class GaugeLorentz>
+template<class GaugeLorentz>
 class WilsonLoops {
 public:
+
+  typedef LorentzScalar<GaugeLorentz> GaugeMat;
+
   //////////////////////////////////////////////////
   // directed plaquette oriented in mu,nu plane
   //////////////////////////////////////////////////
@@ -141,10 +144,10 @@ void siteRectangle(GaugeMat &plaq,const std::vector<GaugeMat> &U, const int mu,
 };
 
 
- typedef WilsonLoops<LatticeColourMatrix,LatticeGaugeField> ColourWilsonLoops;
+ typedef WilsonLoops<LatticeGaugeField> ColourWilsonLoops;
- typedef WilsonLoops<LatticeColourMatrix,LatticeGaugeField> U1WilsonLoops;
+ typedef WilsonLoops<LatticeGaugeField> U1WilsonLoops;
- typedef WilsonLoops<LatticeColourMatrix,LatticeGaugeField> SU2WilsonLoops;
+ typedef WilsonLoops<LatticeGaugeField> SU2WilsonLoops;
- typedef WilsonLoops<LatticeColourMatrix,LatticeGaugeField> SU3WilsonLoops;
+ typedef WilsonLoops<LatticeGaugeField> SU3WilsonLoops;
 
 }}
 

lib/tensors/Tensor_class.h

@@ -23,6 +23,7 @@ template<class vtype> class iScalar
 public:
   vtype _internal;
 
+  typedef vtype element;
   typedef typename GridTypeMapper<vtype>::scalar_type scalar_type;
   typedef typename GridTypeMapper<vtype>::vector_type vector_type;
   typedef typename GridTypeMapper<vtype>::tensor_reduced tensor_reduced_v;
@@ -124,6 +125,7 @@ template<class vtype,int N> class iVector
 public:
   vtype _internal[N];
 
+  typedef vtype element;
   typedef typename GridTypeMapper<vtype>::scalar_type scalar_type;
   typedef typename GridTypeMapper<vtype>::vector_type vector_type;
   typedef typename GridTypeMapper<vtype>::tensor_reduced tensor_reduced_v;
@@ -219,6 +221,7 @@ template<class vtype,int N> class iMatrix
 public:
   vtype _internal[N][N];
 
+  typedef vtype element;
   typedef typename GridTypeMapper<vtype>::scalar_type scalar_type;
   typedef typename GridTypeMapper<vtype>::vector_type vector_type;
   typedef typename GridTypeMapper<vtype>::tensor_reduced tensor_reduced_v;

tests/Test_hmc_EODWFRatio.cc

@@ -17,7 +17,9 @@ int main (int argc, char ** argv)
   GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
   GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
 
+  GridSerialRNG    sRNG;
   GridParallelRNG  pRNG(UGrid);
+  sRNG.SeedRandomDevice();
   pRNG.SeedRandomDevice();
 
   LatticeLorentzColourMatrix     U(UGrid);
@@ -25,7 +27,7 @@ int main (int argc, char ** argv)
   SU3::HotConfiguration(pRNG, U);
 
   // simplify template declaration? Strip the lorentz from the second template
-  WilsonGaugeAction<LatticeLorentzColourMatrix, LatticeColourMatrix> Waction(5.6);
+  WilsonGaugeActionR Waction(5.6);
 
   Real mass=0.04;
   Real pv  =1.0;
@@ -38,21 +40,21 @@ int main (int argc, char ** argv)
   TwoFlavourEvenOddRatioPseudoFermionAction<WilsonImplR> Nf2(NumOp, DenOp,CG,CG);
   
   //Collect actions
-  ActionLevel Level1;
+  ActionLevel<LatticeGaugeField> Level1;
   Level1.push_back(&Nf2);
   Level1.push_back(&Waction);
-  ActionSet FullSet;
+  ActionSet<LatticeGaugeField> FullSet;
   FullSet.push_back(Level1);
 
   // Create integrator
   //  typedef LeapFrog  IntegratorAlgorithm;// change here to change the algorithm
-  typedef MinimumNorm2  IntegratorAlgorithm;// change here to change the algorithm
+  typedef MinimumNorm2<LatticeGaugeField>  IntegratorAlgorithm;// change here to change the algorithm
   IntegratorParameters MDpar(12,20,1.0);
-  Integrator<IntegratorAlgorithm> MDynamics(UGrid,MDpar, FullSet);
+  Integrator<LatticeGaugeField,IntegratorAlgorithm> MDynamics(UGrid,MDpar, FullSet);
 
   // Create HMC
   HMCparameters HMCpar;
-  HybridMonteCarlo<IntegratorAlgorithm>  HMC(HMCpar, MDynamics);
+  HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm>  HMC(HMCpar, MDynamics,sRNG,pRNG);
 
   // Run it
   HMC.evolve(U);

tests/Test_hmc_EOWilsonFermionGauge.cc

@@ -15,15 +15,21 @@ int main (int argc, char ** argv)
   
   GridCartesian            Fine(latt_size,simd_layout,mpi_layout);
   GridRedBlackCartesian  RBFine(latt_size,simd_layout,mpi_layout);
-  GridParallelRNG  pRNG(&Fine);
+
   std::vector<int> seeds({6,7,8,80});
+  GridParallelRNG  pRNG(&Fine);
   pRNG.SeedFixedIntegers(seeds);
+
+  std::vector<int> seedsS({1,2,3,4});
+  GridSerialRNG    sRNG;
+  sRNG.SeedFixedIntegers(seedsS);
+
   LatticeLorentzColourMatrix     U(&Fine);
 
   SU3::HotConfiguration(pRNG, U);
 
   // simplify template declaration? Strip the lorentz from the second template
-  WilsonGaugeAction<LatticeLorentzColourMatrix, LatticeColourMatrix> Waction(5.6);
+  WilsonGaugeActionR Waction(5.6);
 
   Real mass=-0.77;
   WilsonFermionR FermOp(U,Fine,RBFine,mass);
@@ -33,25 +39,25 @@ int main (int argc, char ** argv)
   TwoFlavourEvenOddPseudoFermionAction<WilsonImplR> WilsonNf2(FermOp,CG,CG);
   
   //Collect actions
-  ActionLevel Level1(1);
+  ActionLevel<LatticeGaugeField> Level1(1);
-  ActionLevel Level2(4);
+  ActionLevel<LatticeGaugeField> Level2(4);
   Level1.push_back(&WilsonNf2);
   Level2.push_back(&Waction);
 
-  ActionSet FullSet;
+  ActionSet<LatticeGaugeField> FullSet;
   FullSet.push_back(Level1);
   FullSet.push_back(Level2);
 
   // Create integrator
   //  typedef LeapFrog  IntegratorAlgorithm;// change here to change the algorithm
   //  IntegratorParameters MDpar(12,40,1.0);
-  typedef MinimumNorm2  IntegratorAlgorithm;// change here to change the algorithm
+  typedef MinimumNorm2<LatticeGaugeField>  IntegratorAlgorithm;// change here to change the algorithm
   IntegratorParameters MDpar(12,10,1.0);
-  Integrator<IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet);
+  Integrator<LatticeGaugeField,IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet);
 
   // Create HMC
   HMCparameters HMCpar;
-  HybridMonteCarlo<IntegratorAlgorithm>  HMC(HMCpar, MDynamics);
+  HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm>  HMC(HMCpar, MDynamics,sRNG,pRNG);
 
   HMC.evolve(U);
 

tests/Test_hmc_EOWilsonRatio.cc

@@ -15,14 +15,21 @@ int main (int argc, char ** argv)
   
   GridCartesian            Fine(latt_size,simd_layout,mpi_layout);
   GridRedBlackCartesian  RBFine(latt_size,simd_layout,mpi_layout);
+
+  std::vector<int> seeds({6,7,8,80});
   GridParallelRNG  pRNG(&Fine);
-  pRNG.SeedRandomDevice();
+  pRNG.SeedFixedIntegers(seeds);
+
+  std::vector<int> seedsS({1,2,3,4});
+  GridSerialRNG    sRNG;
+  sRNG.SeedFixedIntegers(seedsS);
+
   LatticeLorentzColourMatrix     U(&Fine);
 
   SU3::HotConfiguration(pRNG, U);
 
   // simplify template declaration? Strip the lorentz from the second template
-  WilsonGaugeAction<LatticeLorentzColourMatrix, LatticeColourMatrix> Waction(5.6);
+  WilsonGaugeActionR Waction(5.6);
 
   Real mass=-0.77;
   Real pv  =0.0;
@@ -34,22 +41,22 @@ int main (int argc, char ** argv)
   TwoFlavourEvenOddRatioPseudoFermionAction<WilsonImplR> WilsonNf2(NumOp, DenOp,CG,CG);
   
   //Collect actions
-  ActionLevel Level1;
+  ActionLevel<LatticeGaugeField> Level1;
   Level1.push_back(&WilsonNf2);
   Level1.push_back(&Waction);
-  ActionSet FullSet;
+  ActionSet<LatticeGaugeField> FullSet;
   FullSet.push_back(Level1);
 
 
   // Create integrator
-  typedef MinimumNorm2  IntegratorAlgorithm;// change here to change the algorithm
+  typedef MinimumNorm2<LatticeGaugeField>  IntegratorAlgorithm;// change here to change the algorithm
   //  typedef LeapFrog  IntegratorAlgorithm;// change here to change the algorithm
   IntegratorParameters MDpar(12,20,1.0);
-  Integrator<IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet);
+  Integrator<LatticeGaugeField,IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet);
 
   // Create HMC
   HMCparameters HMCpar;
-  HybridMonteCarlo<IntegratorAlgorithm>  HMC(HMCpar, MDynamics);
+  HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm>  HMC(HMCpar, MDynamics,sRNG,pRNG);
 
   HMC.evolve(U);
 

tests/Test_hmc_WilsonFermionGauge.cc

@@ -15,16 +15,21 @@ int main (int argc, char ** argv)
   
   GridCartesian            Fine(latt_size,simd_layout,mpi_layout);
   GridRedBlackCartesian  RBFine(latt_size,simd_layout,mpi_layout);
-  GridParallelRNG  pRNG(&Fine);
 
-  std::vector<int> seeds({5,6,7,8});
+  std::vector<int> seeds({6,7,8,80});
+  GridParallelRNG  pRNG(&Fine);
   pRNG.SeedFixedIntegers(seeds);
+
+  std::vector<int> seedsS({1,2,3,4});
+  GridSerialRNG    sRNG;
+  sRNG.SeedFixedIntegers(seedsS);
+
   LatticeLorentzColourMatrix     U(&Fine);
 
   SU3::HotConfiguration(pRNG, U);
 
   // simplify template declaration? Strip the lorentz from the second template
-  WilsonGaugeAction<LatticeLorentzColourMatrix, LatticeColourMatrix> Waction(5.6);
+  WilsonGaugeActionR Waction(5.6);
 
   Real mass=-0.77;
   WilsonFermionR FermOp(U,Fine,RBFine,mass);
@@ -34,25 +39,25 @@ int main (int argc, char ** argv)
   TwoFlavourPseudoFermionAction<WilsonImplR> WilsonNf2(FermOp,CG,CG);
   
   //Collect actions
-  ActionLevel Level1(1);
+  ActionLevel<LatticeGaugeField> Level1(1);
-  ActionLevel Level2(4);
+  ActionLevel<LatticeGaugeField> Level2(4);
   Level1.push_back(&WilsonNf2);
   Level1.push_back(&Waction);
   //  Level1.push_back(&Waction);
 
-  ActionSet FullSet;
+  ActionSet<LatticeGaugeField> FullSet;
   FullSet.push_back(Level1);
   //  FullSet.push_back(Level2);
 
   // Create integrator
-  typedef MinimumNorm2  IntegratorAlgorithm;// change here to change the algorithm
+  typedef MinimumNorm2<LatticeGaugeField>  IntegratorAlgorithm;// change here to change the algorithm
   //  typedef LeapFrog  IntegratorAlgorithm;// change here to change the algorithm
   IntegratorParameters MDpar(12,20,1.0);
-  Integrator<IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet);
+  Integrator<LatticeGaugeField,IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet);
 
   // Create HMC
   HMCparameters HMCpar;
-  HybridMonteCarlo<IntegratorAlgorithm>  HMC(HMCpar, MDynamics);
+  HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm>  HMC(HMCpar, MDynamics,sRNG,pRNG);
 
   HMC.evolve(U);
 

tests/Test_hmc_WilsonGauge.cc

@@ -22,31 +22,37 @@ int main (int argc, char ** argv)
   double volume = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
   
   GridCartesian           Fine(latt_size,simd_layout,mpi_layout);
+
+
+  std::vector<int> seeds({6,7,8,80});
   GridParallelRNG  pRNG(&Fine);
-  std::vector<int> seeds({1,2,3,4,5,6,7,8});
   pRNG.SeedFixedIntegers(seeds);
+
+  std::vector<int> seedsS({1,2,3,4});
+  GridSerialRNG    sRNG;
+  sRNG.SeedFixedIntegers(seedsS);
+
   LatticeGaugeField U(&Fine);
 
   SU3::HotConfiguration(pRNG, U);
- 
 
   // simplify template declaration? Strip the lorentz from the second template
-  WilsonGaugeAction<LatticeGaugeField, LatticeColourMatrix> Waction(6.0);
+  WilsonGaugeActionR Waction(6.0);
 
   //Collect actions
-  ActionLevel Level1;
+  ActionLevel<LatticeGaugeField> Level1;
   Level1.push_back(&Waction);
-  ActionSet FullSet;
+  ActionSet<LatticeGaugeField> FullSet;
   FullSet.push_back(Level1);
 
   // Create integrator
-  typedef MinimumNorm2  IntegratorAlgorithm;// change here to modify the algorithm
+  typedef MinimumNorm2<LatticeGaugeField>  IntegratorAlgorithm;// change here to modify the algorithm
   IntegratorParameters MDpar(12,20,1.0);
-  Integrator<IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet);
+  Integrator<LatticeGaugeField,IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet);
 
   // Create HMC
   HMCparameters HMCpar;
-  HybridMonteCarlo<IntegratorAlgorithm>  HMC(HMCpar, MDynamics);
+  HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm>  HMC(HMCpar, MDynamics, sRNG, pRNG);
 
   HMC.evolve(U);
 

tests/Test_hmc_WilsonRatio.cc

@@ -15,14 +15,21 @@ int main (int argc, char ** argv)
   
   GridCartesian            Fine(latt_size,simd_layout,mpi_layout);
   GridRedBlackCartesian  RBFine(latt_size,simd_layout,mpi_layout);
+
+  std::vector<int> seeds({6,7,8,80});
   GridParallelRNG  pRNG(&Fine);
-  pRNG.SeedRandomDevice();
+  pRNG.SeedFixedIntegers(seeds);
+
+  std::vector<int> seedsS({1,2,3,4});
+  GridSerialRNG    sRNG;
+  sRNG.SeedFixedIntegers(seedsS);
+
   LatticeLorentzColourMatrix     U(&Fine);
 
   SU3::HotConfiguration(pRNG, U);
 
   // simplify template declaration? Strip the lorentz from the second template
-  WilsonGaugeAction<LatticeLorentzColourMatrix, LatticeColourMatrix> Waction(5.6);
+  WilsonGaugeActionR Waction(5.6);
 
   Real mass=-0.77;
   Real pv  =0.0;
@@ -34,22 +41,22 @@ int main (int argc, char ** argv)
   TwoFlavourRatioPseudoFermionAction<WilsonImplR> WilsonNf2(NumOp,FermOp,CG,CG);
   
   //Collect actions
-  ActionLevel Level1;
+  ActionLevel<LatticeGaugeField> Level1;
   Level1.push_back(&WilsonNf2);
   Level1.push_back(&Waction);
-  ActionSet FullSet;
+  ActionSet<LatticeGaugeField> FullSet;
   FullSet.push_back(Level1);
 
 
   // Create integrator
-  typedef MinimumNorm2  IntegratorAlgorithm;// change here to change the algorithm
+  typedef MinimumNorm2<LatticeGaugeField>  IntegratorAlgorithm;// change here to change the algorithm
   //  typedef LeapFrog  IntegratorAlgorithm;// change here to change the algorithm
   IntegratorParameters MDpar(12,20,1.0);
-  Integrator<IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet);
+  Integrator<LatticeGaugeField,IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet);
 
   // Create HMC
   HMCparameters HMCpar;
-  HybridMonteCarlo<IntegratorAlgorithm>  HMC(HMCpar, MDynamics);
+  HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm>  HMC(HMCpar, MDynamics, sRNG, pRNG);
 
   HMC.evolve(U);
 

tests/Test_rhmc_EOWilson1p1.cc

@@ -15,14 +15,22 @@ int main (int argc, char ** argv)
   
   GridCartesian            Fine(latt_size,simd_layout,mpi_layout);
   GridRedBlackCartesian  RBFine(latt_size,simd_layout,mpi_layout);
+
+
+  std::vector<int> seeds({6,7,8,80});
   GridParallelRNG  pRNG(&Fine);
-  pRNG.SeedRandomDevice();
+  pRNG.SeedFixedIntegers(seeds);
+
+  std::vector<int> seedsS({1,2,3,4});
+  GridSerialRNG    sRNG;
+  sRNG.SeedFixedIntegers(seedsS);
+
   LatticeLorentzColourMatrix     U(&Fine);
 
   SU3::HotConfiguration(pRNG, U);
 
   // simplify template declaration? Strip the lorentz from the second template
-  WilsonGaugeAction<LatticeLorentzColourMatrix, LatticeColourMatrix> Waction(5.6);
+  WilsonGaugeActionR Waction(5.6);
 
   Real mass=-0.77;
   WilsonFermionR FermOp(U,Fine,RBFine,mass);
@@ -33,23 +41,23 @@ int main (int argc, char ** argv)
   OneFlavourEvenOddRationalPseudoFermionAction<WilsonImplR> WilsonNf1b(FermOp,Params);
   
   //Collect actions
-  ActionLevel Level1;
+  ActionLevel<LatticeGaugeField> Level1;
   Level1.push_back(&WilsonNf1a);
   Level1.push_back(&WilsonNf1b);
   Level1.push_back(&Waction);
 
-  ActionSet FullSet;
+  ActionSet<LatticeGaugeField> FullSet;
   FullSet.push_back(Level1);
 
   // Create integrator
-  typedef MinimumNorm2  IntegratorAlgorithm;// change here to change the algorithm
+  typedef MinimumNorm2<LatticeGaugeField>  IntegratorAlgorithm;// change here to change the algorithm
 
   IntegratorParameters MDpar(12,20,1.0);
-  Integrator<IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet);
+  Integrator<LatticeGaugeField,IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet);
 
   // Create HMC
   HMCparameters HMCpar;
-  HybridMonteCarlo<IntegratorAlgorithm>  HMC(HMCpar, MDynamics);
+  HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm>  HMC(HMCpar, MDynamics, sRNG, pRNG);
   HMC.evolve(U);
 
 }

tests/Test_rhmc_EOWilsonRatio.cc

@@ -15,14 +15,22 @@ int main (int argc, char ** argv)
   
   GridCartesian            Fine(latt_size,simd_layout,mpi_layout);
   GridRedBlackCartesian  RBFine(latt_size,simd_layout,mpi_layout);
+
+
+  std::vector<int> seeds({6,7,8,80});
   GridParallelRNG  pRNG(&Fine);
-  pRNG.SeedRandomDevice();
+  pRNG.SeedFixedIntegers(seeds);
+
+  std::vector<int> seedsS({1,2,3,4});
+  GridSerialRNG    sRNG;
+  sRNG.SeedFixedIntegers(seedsS);
+
   LatticeLorentzColourMatrix     U(&Fine);
 
   SU3::HotConfiguration(pRNG, U);
 
   // simplify template declaration? Strip the lorentz from the second template
-  WilsonGaugeAction<LatticeLorentzColourMatrix, LatticeColourMatrix> Waction(5.6);
+  WilsonGaugeActionR Waction(5.6);
 
   Real mass=-0.77;
   Real pv  =0.0;
@@ -35,23 +43,23 @@ int main (int argc, char ** argv)
   OneFlavourEvenOddRatioRationalPseudoFermionAction<WilsonImplR> WilsonNf1b(NumOp,DenOp,Params);
   
   //Collect actions
-  ActionLevel Level1;
+  ActionLevel<LatticeGaugeField> Level1;
   Level1.push_back(&WilsonNf1a);
   Level1.push_back(&WilsonNf1b);
   Level1.push_back(&Waction);
-  ActionSet FullSet;
+  ActionSet<LatticeGaugeField> FullSet;
   FullSet.push_back(Level1);
 
 
   // Create integrator
-  typedef MinimumNorm2  IntegratorAlgorithm;// change here to change the algorithm
+  typedef MinimumNorm2<LatticeGaugeField>  IntegratorAlgorithm;// change here to change the algorithm
   //  typedef LeapFrog  IntegratorAlgorithm;// change here to change the algorithm
   IntegratorParameters MDpar(12,20,1.0);
-  Integrator<IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet);
+  Integrator<LatticeGaugeField,IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet);
 
   // Create HMC
   HMCparameters HMCpar;
-  HybridMonteCarlo<IntegratorAlgorithm>  HMC(HMCpar, MDynamics);
+  HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm>  HMC(HMCpar, MDynamics, sRNG, pRNG);
 
   HMC.evolve(U);
 

tests/Test_rhmc_Wilson1p1.cc

@@ -15,14 +15,22 @@ int main (int argc, char ** argv)
   
   GridCartesian            Fine(latt_size,simd_layout,mpi_layout);
   GridRedBlackCartesian  RBFine(latt_size,simd_layout,mpi_layout);
+
+
+  std::vector<int> seeds({6,7,8,80});
   GridParallelRNG  pRNG(&Fine);
-  pRNG.SeedRandomDevice();
+  pRNG.SeedFixedIntegers(seeds);
+
+  std::vector<int> seedsS({1,2,3,4});
+  GridSerialRNG    sRNG;
+  sRNG.SeedFixedIntegers(seedsS);
+
   LatticeLorentzColourMatrix     U(&Fine);
 
   SU3::HotConfiguration(pRNG, U);
 
   // simplify template declaration? Strip the lorentz from the second template
-  WilsonGaugeAction<LatticeLorentzColourMatrix, LatticeColourMatrix> Waction(5.6);
+  WilsonGaugeActionR Waction(5.6);
 
   Real mass=-0.77;
   WilsonFermionR FermOp(U,Fine,RBFine,mass);
@@ -33,23 +41,23 @@ int main (int argc, char ** argv)
   OneFlavourRationalPseudoFermionAction<WilsonImplR> WilsonNf1b(FermOp,Params);
   
   //Collect actions
-  ActionLevel Level1;
+  ActionLevel<LatticeGaugeField> Level1;
   Level1.push_back(&WilsonNf1a);
   Level1.push_back(&WilsonNf1b);
   Level1.push_back(&Waction);
 
-  ActionSet FullSet;
+  ActionSet<LatticeGaugeField> FullSet;
   FullSet.push_back(Level1);
 
   // Create integrator
-  typedef MinimumNorm2  IntegratorAlgorithm;// change here to change the algorithm
+  typedef MinimumNorm2<LatticeGaugeField>  IntegratorAlgorithm;// change here to change the algorithm
 
   IntegratorParameters MDpar(12,20,1.0);
-  Integrator<IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet);
+  Integrator<LatticeGaugeField,IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet);
 
   // Create HMC
   HMCparameters HMCpar;
-  HybridMonteCarlo<IntegratorAlgorithm>  HMC(HMCpar, MDynamics);
+  HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm>  HMC(HMCpar, MDynamics, sRNG, pRNG);
 
   HMC.evolve(U);
 

tests/Test_rhmc_WilsonRatio.cc

@@ -15,14 +15,22 @@ int main (int argc, char ** argv)
   
   GridCartesian            Fine(latt_size,simd_layout,mpi_layout);
   GridRedBlackCartesian  RBFine(latt_size,simd_layout,mpi_layout);
+
+
+  std::vector<int> seeds({6,7,8,80});
   GridParallelRNG  pRNG(&Fine);
-  pRNG.SeedRandomDevice();
+  pRNG.SeedFixedIntegers(seeds);
+
+  std::vector<int> seedsS({1,2,3,4});
+  GridSerialRNG    sRNG;
+  sRNG.SeedFixedIntegers(seedsS);
+
   LatticeLorentzColourMatrix     U(&Fine);
 
   SU3::HotConfiguration(pRNG, U);
 
   // simplify template declaration? Strip the lorentz from the second template
-  WilsonGaugeAction<LatticeLorentzColourMatrix, LatticeColourMatrix> Waction(5.6);
+  WilsonGaugeActionR Waction(5.6);
 
   Real mass=-0.77;
   Real pv  =0.0;
@@ -35,23 +43,23 @@ int main (int argc, char ** argv)
   OneFlavourRatioRationalPseudoFermionAction<WilsonImplR> WilsonNf1b(NumOp,DenOp,Params);
   
   //Collect actions
-  ActionLevel Level1;
+  ActionLevel<LatticeGaugeField> Level1;
   Level1.push_back(&WilsonNf1a);
   Level1.push_back(&WilsonNf1b);
   Level1.push_back(&Waction);
-  ActionSet FullSet;
+  ActionSet<LatticeGaugeField> FullSet;
   FullSet.push_back(Level1);
 
 
   // Create integrator
-  typedef MinimumNorm2  IntegratorAlgorithm;// change here to change the algorithm
+  typedef MinimumNorm2<LatticeGaugeField>  IntegratorAlgorithm;// change here to change the algorithm
   //  typedef LeapFrog  IntegratorAlgorithm;// change here to change the algorithm
   IntegratorParameters MDpar(12,20,1.0);
-  Integrator<IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet);
+  Integrator<LatticeGaugeField,IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet);
 
   // Create HMC
   HMCparameters HMCpar;
-  HybridMonteCarlo<IntegratorAlgorithm>  HMC(HMCpar, MDynamics);
+  HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm>  HMC(HMCpar, MDynamics, sRNG, pRNG);
 
   HMC.evolve(U);