Merge branch 'develop' of https://github.com/paboyle/Grid into feature/Lanczos

configure.ac

@@ -13,6 +13,10 @@ m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])])
 ################ Get git info
 #AC_REVISION([m4_esyscmd_s([./scripts/configure.commit])])
 
+################ Set flags
+# do not move!
+CXXFLAGS="-O3 $CXXFLAGS"
+
 ############### Checks for programs
 AC_PROG_CXX
 AC_PROG_RANLIB
@@ -27,7 +31,6 @@ AX_GXX_VERSION
 AC_DEFINE_UNQUOTED([GXX_VERSION],["$GXX_VERSION"],
       [version of g++ that will compile the code])
 
-CXXFLAGS="-O3 $CXXFLAGS"
 
 
 ############### Checks for typedefs, structures, and compiler characteristics
@@ -55,6 +58,10 @@ AC_CHECK_HEADERS(numaif.h)
 AC_CHECK_DECLS([ntohll],[], [], [[#include <arpa/inet.h>]])
 AC_CHECK_DECLS([be64toh],[], [], [[#include <arpa/inet.h>]])
 
+############## Standard libraries
+AC_CHECK_LIB([m],[cos])
+AC_CHECK_LIB([stdc++],[abort])
+
 ############### GMP and MPFR
 AC_ARG_WITH([gmp],
     [AS_HELP_STRING([--with-gmp=prefix],

lib/cartesian/Cartesian_base.h

@@ -185,17 +185,18 @@ public:
     ////////////////////////////////////////////////////////////////
 
     void show_decomposition(){
-      std::cout << GridLogMessage << "Full Dimensions    : " << _fdimensions << std::endl;
+      std::cout << GridLogMessage << "\tFull Dimensions    : " << _fdimensions << std::endl;
-      std::cout << GridLogMessage << "Global Dimensions  : " << _gdimensions << std::endl;
+      std::cout << GridLogMessage << "\tSIMD layout        : " << _simd_layout << std::endl;
-      std::cout << GridLogMessage << "Local Dimensions   : " << _ldimensions << std::endl;
+      std::cout << GridLogMessage << "\tGlobal Dimensions  : " << _gdimensions << std::endl;
-      std::cout << GridLogMessage << "Reduced Dimensions : " << _rdimensions << std::endl;
+      std::cout << GridLogMessage << "\tLocal Dimensions   : " << _ldimensions << std::endl;
-      std::cout << GridLogMessage << "Outer strides      : " << _ostride << std::endl;
+      std::cout << GridLogMessage << "\tReduced Dimensions : " << _rdimensions << std::endl;
-      std::cout << GridLogMessage << "Inner strides      : " << _istride << std::endl;
+      std::cout << GridLogMessage << "\tOuter strides      : " << _ostride << std::endl;
-      std::cout << GridLogMessage << "iSites             : " << _isites << std::endl;
+      std::cout << GridLogMessage << "\tInner strides      : " << _istride << std::endl;
-      std::cout << GridLogMessage << "oSites             : " << _osites << std::endl;
+      std::cout << GridLogMessage << "\tiSites             : " << _isites << std::endl;
-      std::cout << GridLogMessage << "lSites             : " << lSites() << std::endl;        
+      std::cout << GridLogMessage << "\toSites             : " << _osites << std::endl;
-      std::cout << GridLogMessage << "gSites             : " << gSites() << std::endl;
+      std::cout << GridLogMessage << "\tlSites             : " << lSites() << std::endl;        
-      std::cout << GridLogMessage << "Nd                 : " << _ndimension << std::endl;             
+      std::cout << GridLogMessage << "\tgSites             : " << gSites() << std::endl;
+      std::cout << GridLogMessage << "\tNd                 : " << _ndimension << std::endl;             
     } 
 
     ////////////////////////////////////////////////////////////////

lib/parallelIO/BinaryIO.h

@@ -99,13 +99,14 @@ class BinaryIO {
     NerscChecksum(grid,scalardata,nersc_csum);
   }
 
-  template<class fobj> static inline void NerscChecksum(GridBase *grid,std::vector<fobj> &fbuf,uint32_t &nersc_csum)
+  template <class fobj>
+  static inline void NerscChecksum(GridBase *grid, std::vector<fobj> &fbuf, uint32_t &nersc_csum)
   {
     const uint64_t size32 = sizeof(fobj) / sizeof(uint32_t);
 
-
     uint64_t lsites = grid->lSites();
-    if (fbuf.size()==1) {
+    if (fbuf.size() == 1)
+    {
       lsites = 1;
     }
 
@@ -114,9 +115,11 @@ class BinaryIO {
       uint32_t nersc_csum_thr = 0;
 
       #pragma omp for
-      for(uint64_t local_site=0;local_site<lsites;local_site++){
+      for (uint64_t local_site = 0; local_site < lsites; local_site++)
+      {
         uint32_t *site_buf = (uint32_t *)&fbuf[local_site];
-	for(uint64_t j=0;j<size32;j++){
+        for (uint64_t j = 0; j < size32; j++)
+        {
           nersc_csum_thr = nersc_csum_thr + site_buf[j];
         }
       }
@@ -127,6 +130,7 @@ class BinaryIO {
       }
     }
   }
+
   template<class fobj> static inline void ScidacChecksum(GridBase *grid,std::vector<fobj> &fbuf,uint32_t &scidac_csuma,uint32_t &scidac_csumb)
   {
     const uint64_t size32 = sizeof(fobj)/sizeof(uint32_t);
@@ -367,12 +371,16 @@ class BinaryIO {
                   << iodata.size() * sizeof(fobj) << " bytes" << std::endl;
         std::ifstream fin;
         fin.open(file, std::ios::binary | std::ios::in);
-	if ( control & BINARYIO_MASTER_APPEND )  {
+        if (control & BINARYIO_MASTER_APPEND)
+        {
           fin.seekg(-sizeof(fobj), fin.end);
-	} else { 
+        }
+        else
+        {
           fin.seekg(offset + myrank * lsites * sizeof(fobj));
         }
-	fin.read((char *)&iodata[0],iodata.size()*sizeof(fobj));assert( fin.fail()==0);
+        fin.read((char *)&iodata[0], iodata.size() * sizeof(fobj));
+        assert(fin.fail() == 0);
         fin.close();
       }
       timer.Stop();
@@ -406,9 +414,29 @@ class BinaryIO {
       if ( (control & BINARYIO_LEXICOGRAPHIC) && (nrank > 1) ) {
 #ifdef USE_MPI_IO
         std::cout << GridLogMessage << "MPI write I/O " << file << std::endl;
-	ierr=MPI_File_open(grid->communicator,(char *) file.c_str(), MPI_MODE_RDWR|MPI_MODE_CREATE,MPI_INFO_NULL, &fh); assert(ierr==0);
+        ierr = MPI_File_open(grid->communicator, (char *)file.c_str(), MPI_MODE_RDWR | MPI_MODE_CREATE, MPI_INFO_NULL, &fh);
-	ierr=MPI_File_set_view(fh, disp, mpiObject, fileArray, "native", MPI_INFO_NULL);                        assert(ierr==0);
+        std::cout << GridLogMessage << "Checking for errors" << std::endl;
-	ierr=MPI_File_write_all(fh, &iodata[0], 1, localArray, &status);                                        assert(ierr==0);
+        if (ierr != MPI_SUCCESS)
+        {
+          char error_string[BUFSIZ];
+          int length_of_error_string, error_class;
+
+          MPI_Error_class(ierr, &error_class);
+          MPI_Error_string(error_class, error_string, &length_of_error_string);
+          fprintf(stderr, "%3d: %s\n", myrank, error_string);
+          MPI_Error_string(ierr, error_string, &length_of_error_string);
+          fprintf(stderr, "%3d: %s\n", myrank, error_string);
+          MPI_Abort(MPI_COMM_WORLD, 1); //assert(ierr == 0);
+        }
+
+        std::cout << GridLogDebug << "MPI read I/O set view " << file << std::endl;
+        ierr = MPI_File_set_view(fh, disp, mpiObject, fileArray, "native", MPI_INFO_NULL);
+        assert(ierr == 0);
+
+        std::cout << GridLogDebug << "MPI read I/O write all " << file << std::endl;
+        ierr = MPI_File_write_all(fh, &iodata[0], 1, localArray, &status);
+        assert(ierr == 0);
+
         MPI_File_close(&fh);
         MPI_Type_free(&fileArray);
         MPI_Type_free(&localArray);
@@ -416,15 +444,43 @@ class BinaryIO {
 	assert(0);
 #endif
       } else { 
-	std::ofstream fout; fout.open(file,std::ios::binary|std::ios::out|std::ios::in);
+        
+	std::ofstream fout; 
+  fout.exceptions ( std::fstream::failbit | std::fstream::badbit );
+  try {
+    fout.open(file,std::ios::binary|std::ios::out|std::ios::in);
+  } catch (const std::fstream::failure& exc) {
+    std::cout << GridLogError << "Error in opening the file " << file << " for output" <<std::endl;
+    std::cout << GridLogError << "Exception description: " << exc.what() << std::endl;
+    std::cout << GridLogError << "Probable cause: wrong path, inaccessible location "<< std::endl;
+    #ifdef USE_MPI_IO
+    MPI_Abort(MPI_COMM_WORLD,1);
+    #else
+    exit(1);
+    #endif
+  }
 	std::cout << GridLogMessage<< "C++ write I/O "<< file<<" : "
 		        << iodata.size()*sizeof(fobj)<<" bytes"<<std::endl;
+	
   if ( control & BINARYIO_MASTER_APPEND )  {
 	  fout.seekp(0,fout.end);
 	} else {
 	  fout.seekp(offset+myrank*lsites*sizeof(fobj));
 	}
-	fout.write((char *)&iodata[0],iodata.size()*sizeof(fobj));assert( fout.fail()==0);
+  
+  try {
+  	fout.write((char *)&iodata[0],iodata.size()*sizeof(fobj));//assert( fout.fail()==0);
+  }
+  catch (const std::fstream::failure& exc) {
+    std::cout << "Exception in writing file " << file << std::endl;
+    std::cout << GridLogError << "Exception description: "<< exc.what() << std::endl;
+    #ifdef USE_MPI_IO
+    MPI_Abort(MPI_COMM_WORLD,1);
+    #else
+    exit(1);
+    #endif
+  }
+
 	fout.close();
   }
   timer.Stop();
@@ -442,12 +498,15 @@ class BinaryIO {
     //////////////////////////////////////////////////////////////////////////////
     // Safety check
     //////////////////////////////////////////////////////////////////////////////
+    // if the data size is 1 we do not want to sum over the MPI ranks
+    if (iodata.size() != 1){
       grid->Barrier();
       grid->GlobalSum(nersc_csum);
       grid->GlobalXOR(scidac_csuma);
       grid->GlobalXOR(scidac_csumb);
       grid->Barrier();
     }
+  }
 
   /////////////////////////////////////////////////////////////////////////////
   // Read a Lattice of object
@@ -546,9 +605,9 @@ class BinaryIO {
     int gsites = grid->gSites();
     int lsites = grid->lSites();
 
-    uint32_t nersc_csum_tmp;
+    uint32_t nersc_csum_tmp   = 0;
-    uint32_t scidac_csuma_tmp;
+    uint32_t scidac_csuma_tmp = 0;
-    uint32_t scidac_csumb_tmp;
+    uint32_t scidac_csumb_tmp = 0;
 
     GridStopWatch timer;
 

lib/qcd/hmc/HMCResourceManager.h

@@ -165,7 +165,7 @@ class HMCResourceManager {
   // Grids
   //////////////////////////////////////////////////////////////
 
-  void AddGrid(std::string s, GridModule& M) {
+  void AddGrid(const std::string s, GridModule& M) {
     // Check for name clashes
     auto search = Grids.find(s);
     if (search != Grids.end()) {
@@ -174,14 +174,24 @@ class HMCResourceManager {
       exit(1);
     }
     Grids[s] = std::move(M);
+    std::cout << GridLogMessage << "::::::::::::::::::::::::::::::::::::::::" <<std::endl;
+    std::cout << GridLogMessage << "HMCResourceManager:" << std::endl;
+    std::cout << GridLogMessage << "Created grid set with name '" << s << "' and decomposition for the full cartesian " << std::endl;
+    Grids[s].show_full_decomposition();
+    std::cout << GridLogMessage << "::::::::::::::::::::::::::::::::::::::::" <<std::endl;
   }
 
   // Add a named grid set, 4d shortcut
-  void AddFourDimGrid(std::string s) {
+  void AddFourDimGrid(const std::string s) {
     GridFourDimModule<vComplex> Mod;
     AddGrid(s, Mod);
   }
 
+  // Add a named grid set, 4d shortcut + tweak simd lanes
+  void AddFourDimGrid(const std::string s, const std::vector<int> simd_decomposition) {
+    GridFourDimModule<vComplex> Mod(simd_decomposition);
+    AddGrid(s, Mod);
+  }
 
 
   GridCartesian* GetCartesian(std::string s = "") {

lib/qcd/hmc/HMC_GridModules.h

@@ -43,10 +43,11 @@ public:
   std::string, lattice,
   std::string, mpi);
 
-  std::vector<int> getLattice(){return strToVec<int>(lattice);}
+  std::vector<int> getLattice() const {return strToVec<int>(lattice);}
-  std::vector<int> getMpi()    {return strToVec<int>(mpi);}
+  std::vector<int> getMpi()     const {return strToVec<int>(mpi);}
 
-  void check(){
+
+  void check() const {
     if (getLattice().size() != getMpi().size() ) {
       std::cout << GridLogError
                 << "Error in GridModuleParameters: lattice and mpi dimensions "
@@ -84,6 +85,8 @@ class GridModule {
 
   void set_full(GridCartesian* grid) { grid_.reset(grid); }
   void set_rb(GridRedBlackCartesian* rbgrid) { rbgrid_.reset(rbgrid); }
+  void show_full_decomposition(){ grid_->show_decomposition(); }
+  void show_rb_decomposition(){ rbgrid_->show_decomposition(); }
 
  protected:
   std::unique_ptr<GridCartesian> grid_;
@@ -96,27 +99,68 @@ class GridModule {
 ////////////////////////////////////
 // Note: the space time grid should be out of the QCD namespace
 template <class vector_type>
-class GridFourDimModule : public GridModule {
+class GridFourDimModule : public GridModule
+{
 public:
-  GridFourDimModule() {
+  GridFourDimModule()
+  {
     using namespace QCD;
     set_full(SpaceTimeGrid::makeFourDimGrid(
-        GridDefaultLatt(), GridDefaultSimd(4, vector_type::Nsimd()),
+        GridDefaultLatt(), 
+        GridDefaultSimd(4, vector_type::Nsimd()),
         GridDefaultMpi()));
     set_rb(SpaceTimeGrid::makeFourDimRedBlackGrid(grid_.get()));
   }
 
-  GridFourDimModule(GridModuleParameters Params) {
+  GridFourDimModule(const std::vector<int> tweak_simd)
+  {
+    using namespace QCD;
+    if (tweak_simd.size() != 4)
+    {
+      std::cout << GridLogError
+                << "Error in GridFourDimModule: SIMD size different from 4" 
+                << std::endl;
+      exit(1);
+    }
+
+    // Checks that the product agrees with the expectation
+    int simd_sum = 1;
+    for (auto &n : tweak_simd)
+      simd_sum *= n;
+    std::cout << GridLogDebug << "TweakSIMD: " << tweak_simd << "  Sum: " << simd_sum << std::endl;
+
+    if (simd_sum == vector_type::Nsimd())
+    {
+      set_full(SpaceTimeGrid::makeFourDimGrid(
+          GridDefaultLatt(), 
+          tweak_simd, 
+          GridDefaultMpi()));
+      set_rb(SpaceTimeGrid::makeFourDimRedBlackGrid(grid_.get()));
+    }
+    else
+    {
+      std::cout << GridLogError 
+                << "Error in GridFourDimModule: SIMD lanes must sum to " 
+                << vector_type::Nsimd() 
+                << std::endl;
+    }
+  }
+
+  GridFourDimModule(const GridModuleParameters Params)
+  {
     using namespace QCD;
-    Params.check();
     std::vector<int> lattice_v = Params.getLattice();
     std::vector<int> mpi_v = Params.getMpi();
-    if (lattice_v.size() == 4) {
+    if (lattice_v.size() == 4)
+    {
       set_full(SpaceTimeGrid::makeFourDimGrid(
-          lattice_v, GridDefaultSimd(4, vector_type::Nsimd()),
+          lattice_v, 
+          GridDefaultSimd(4, vector_type::Nsimd()),
           mpi_v));
       set_rb(SpaceTimeGrid::makeFourDimRedBlackGrid(grid_.get()));
-    } else {
+    }
+    else
+    {
       std::cout << GridLogError
                 << "Error in GridFourDimModule: lattice dimension different from 4"
                 << std::endl;

lib/qcd/modules/ObservableModules.h

@@ -84,8 +84,6 @@ class PlaquetteMod: public ObservableModule<PlaquetteLogger<Impl>, NoParameters>
   typedef ObservableModule<PlaquetteLogger<Impl>, NoParameters> ObsBase;
   using ObsBase::ObsBase; // for constructors
 
-
-
   // acquire resource
   virtual void initialize(){
     this->ObservablePtr.reset(new PlaquetteLogger<Impl>());
@@ -94,23 +92,22 @@ class PlaquetteMod: public ObservableModule<PlaquetteLogger<Impl>, NoParameters>
   PlaquetteMod(): ObsBase(NoParameters()){}
 };
 
+
 template < class Impl >
-class TopologicalChargeMod: public ObservableModule<TopologicalCharge<Impl>, NoParameters>{
+class TopologicalChargeMod: public ObservableModule<TopologicalCharge<Impl>, TopologyObsParameters>{
-  typedef ObservableModule<TopologicalCharge<Impl>, NoParameters> ObsBase;
+  typedef ObservableModule<TopologicalCharge<Impl>, TopologyObsParameters> ObsBase;
   using ObsBase::ObsBase; // for constructors
 
-
-
   // acquire resource
   virtual void initialize(){
-    this->ObservablePtr.reset(new TopologicalCharge<Impl>());
+    this->ObservablePtr.reset(new TopologicalCharge<Impl>(this->Par_));
   }
   public:
-  TopologicalChargeMod(): ObsBase(NoParameters()){}
+  TopologicalChargeMod(TopologyObsParameters Par): ObsBase(Par){}
+  TopologicalChargeMod(): ObsBase(){}
 };
 
 
-
 }// QCD temporarily here
 
 

lib/qcd/observables/topological_charge.h

@@ -33,9 +33,45 @@ directory
 namespace Grid {
 namespace QCD {
 
+struct TopologySmearingParameters : Serializable {
+    GRID_SERIALIZABLE_CLASS_MEMBERS(TopologySmearingParameters,
+    int, steps,
+    float, step_size,
+    int, meas_interval,
+    float, maxTau);
+
+    TopologySmearingParameters(int s = 0, float ss = 0.0f, int mi = 0, float mT = 0.0f):
+        steps(s), step_size(ss), meas_interval(mi), maxTau(mT){}
+
+    template < class ReaderClass >
+    TopologySmearingParameters(Reader<ReaderClass>& Reader){
+        read(Reader, "Smearing", *this);  
+    }  
+};
+
+
+
+struct TopologyObsParameters : Serializable {
+    GRID_SERIALIZABLE_CLASS_MEMBERS(TopologyObsParameters,
+      int, interval,
+      bool, do_smearing,
+      TopologySmearingParameters, Smearing);  
+
+    TopologyObsParameters(int interval = 1, bool smearing = false):
+        interval(interval), Smearing(smearing){}
+
+    template <class ReaderClass >
+      TopologyObsParameters(Reader<ReaderClass>& Reader){
+        read(Reader, "TopologyMeasurement", *this);
+  }
+};
+
+
 // this is only defined for a gauge theory
 template <class Impl>
 class TopologicalCharge : public HmcObservable<typename Impl::Field> {
+    TopologyObsParameters Pars;
+
  public:
     // here forces the Impl to be of gauge fields
     // if not the compiler will complain
@@ -44,21 +80,40 @@ class TopologicalCharge : public HmcObservable<typename Impl::Field> {
     // necessary for HmcObservable compatibility
     typedef typename Impl::Field Field;
 
+    TopologicalCharge(int interval = 1, bool do_smearing = false):
+        Pars(interval, do_smearing){}
+    
+    TopologicalCharge(TopologyObsParameters P):Pars(P){
+        std::cout << GridLogDebug << "Creating TopologicalCharge " << std::endl;
+    }
+
     void TrajectoryComplete(int traj,
                             Field &U,
                             GridSerialRNG &sRNG,
                             GridParallelRNG &pRNG) {
 
-    Real q = WilsonLoops<Impl>::TopologicalCharge(U);
+    if (traj%Pars.interval == 0){
-
+        // Smearing
+        Field Usmear = U;
         int def_prec = std::cout.precision();
         
+        if (Pars.do_smearing){
+            // using wilson flow by default here
+            WilsonFlow<PeriodicGimplR> WF(Pars.Smearing.steps, Pars.Smearing.step_size, Pars.Smearing.meas_interval);
+            WF.smear_adaptive(Usmear, U, Pars.Smearing.maxTau);
+            Real T0   = WF.energyDensityPlaquette(Usmear);
+            std::cout << GridLogMessage << std::setprecision(std::numeric_limits<Real>::digits10 + 1)
+                      << "T0                : [ " << traj << " ] "<< T0 << std::endl;
+        }
+
+        Real q    = WilsonLoops<Impl>::TopologicalCharge(Usmear);
         std::cout << GridLogMessage
             << std::setprecision(std::numeric_limits<Real>::digits10 + 1)
             << "Topological Charge: [ " << traj << " ] "<< q << std::endl;
 
         std::cout.precision(def_prec);
         }
+    }
 
 };
 }

lib/qcd/utils/GaugeFix.h

@@ -26,8 +26,8 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     /*  END LEGAL */
 //#include <Grid/Grid.h>
 
-using namespace Grid;
+namespace Grid {
-using namespace Grid::QCD;
+namespace QCD {
 
 template <class Gimpl> 
 class FourierAcceleratedGaugeFixer  : public Gimpl {
@@ -186,3 +186,5 @@ class FourierAcceleratedGaugeFixer  : public Gimpl {
   }  
 };
 
+}
+}

lib/tensors/Tensor_arith_mul.h

@@ -98,7 +98,9 @@ template<class rtype,class vtype,class mtype,int N>
 strong_inline void mult(iVector<rtype,N> * __restrict__ ret,
                  const iVector<vtype,N> * __restrict__ rhs,
                  const iScalar<mtype> * __restrict__ lhs){
-    mult(ret,lhs,rhs);
+    for(int c1=0;c1<N;c1++){
+        mult(&ret->_internal[c1],&rhs->_internal[c1],&lhs->_internal);
+    }                 
 }
     
 

lib/util/Init.cc

@@ -377,7 +377,7 @@ void Grid_init(int *argc,char ***argv)
   std::cout << GridLogDebug << "Requesting "<< CartesianCommunicator::MAX_MPI_SHM_BYTES <<" byte stencil comms buffers "<<std::endl;
 
   if( GridCmdOptionExists(*argv,*argv+*argc,"--decomposition") ){
-    std::cout<<GridLogMessage<<"Grid Decomposition\n";
+    std::cout<<GridLogMessage<<"Grid Default Decomposition patterns\n";
     std::cout<<GridLogMessage<<"\tOpenMP threads : "<<GridThread::GetThreads()<<std::endl;
     std::cout<<GridLogMessage<<"\tMPI tasks      : "<<GridCmdVectorIntToString(GridDefaultMpi())<<std::endl;
     std::cout<<GridLogMessage<<"\tvRealF         : "<<sizeof(vRealF)*8    <<"bits ; " <<GridCmdVectorIntToString(GridDefaultSimd(4,vRealF::Nsimd()))<<std::endl;

tests/hmc/Test_hmc_EOMobiusRatio.cc

@@ -40,12 +40,6 @@ namespace Grid{
 				    double, StoppingCondition,
 				    int, MaxCGIterations,
 				    bool, ApplySmearing);
-
-    //template <class ReaderClass >
-    //FermionParameters(Reader<ReaderClass>& Reader){
-    //  read(Reader, "Mobius", *this);
-    //}
-
   };
 
   
@@ -113,9 +107,17 @@ int main(int argc, char **argv) {
   bool ApplySmearing = MyParams.Mobius.ApplySmearing;
   
   
+  // Use this if you want to tweak the default decomposition
+  // commented out as very architecture speficic
+  
+  //std::vector<int> simd_lanes({2,2,1,1});
+
+  // Grid from the command line arguments --grid and --mpi
+  // drop the simd_lanes argument to fall back to the default decomposition for the SIMD lanes
+  
+  //TheHMC.Resources.AddFourDimGrid("gauge", simd_lanes); // tweak the SIMD lanes
+  TheHMC.Resources.AddFourDimGrid("gauge"); // use default simd lanes decomposition
   
-  // Grid from the command line
-  TheHMC.Resources.AddFourDimGrid("gauge");
   // Possibile to create the module by hand 
   // hardcoding parameters or using a Reader
 

tests/hmc/Test_hmc_WilsonGauge.cc

@@ -66,7 +66,14 @@ int main(int argc, char **argv) {
   typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
   typedef TopologicalChargeMod<HMCWrapper::ImplPolicy> QObs;
   TheHMC.Resources.AddObservable<PlaqObs>();
-  TheHMC.Resources.AddObservable<QObs>();
+  TopologyObsParameters TopParams;
+  TopParams.interval = 5;
+  TopParams.do_smearing = true;
+  TopParams.Smearing.steps = 200;
+  TopParams.Smearing.step_size = 0.01;
+  TopParams.Smearing.meas_interval = 50;
+  TopParams.Smearing.maxTau = 2.0; 
+  TheHMC.Resources.AddObservable<QObs>(TopParams);
   //////////////////////////////////////////////
 
   /////////////////////////////////////////////////////////////