Merge branch 'develop' into feature/CG-reliable-update

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/parallelIO/BinaryIO.h

@@ -98,35 +98,39 @@ 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);
+    const uint64_t size32 = sizeof(fobj) / sizeof(uint32_t);
 
-
+    uint64_t lsites = grid->lSites();
-    uint64_t lsites              =grid->lSites();
+    if (fbuf.size() == 1)
-    if (fbuf.size()==1) {
+    {
-      lsites=1;
+      lsites = 1;
     }
 
-#pragma omp parallel
+    #pragma omp parallel
-    { 
+    {
-      uint32_t nersc_csum_thr=0;
+      uint32_t nersc_csum_thr = 0;
 
-#pragma omp for
+      #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++){
+        uint32_t *site_buf = (uint32_t *)&fbuf[local_site];
-	  nersc_csum_thr=nersc_csum_thr+site_buf[j];
+        for (uint64_t j = 0; j < size32; j++)
-	}
+        {
+          nersc_csum_thr = nersc_csum_thr + site_buf[j];
+        }
       }
 
-#pragma omp critical
+      #pragma omp critical
       {
-	nersc_csum  += nersc_csum_thr;
+        nersc_csum += nersc_csum_thr;
       }
     }
   }
+
   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);
@@ -266,7 +270,7 @@ class BinaryIO {
     grid->Barrier();
     GridStopWatch timer; 
     GridStopWatch bstimer;
-
+    
     nersc_csum=0;
     scidac_csuma=0;
     scidac_csumb=0;
@@ -362,18 +366,22 @@ class BinaryIO {
 #else 
 	assert(0);
 #endif
-      } else { 
+      } else {
-	std::cout<< GridLogMessage<< "C++ read I/O "<< file<<" : "
+        std::cout << GridLogMessage << "C++ read I/O " << file << " : "
-		 << iodata.size()*sizeof(fobj)<<" bytes"<<std::endl;
+                  << iodata.size() * sizeof(fobj) << " bytes" << std::endl;
-	std::ifstream fin;
+        std::ifstream fin;
-	fin.open(file,std::ios::binary|std::ios::in);
+        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 { 
+          fin.seekg(-sizeof(fobj), fin.end);
-	  fin.seekg(offset+myrank*lsites*sizeof(fobj));
+        }
-	}
+        else
-	fin.read((char *)&iodata[0],iodata.size()*sizeof(fobj));assert( fin.fail()==0);
+        {
-	fin.close();
+          fin.seekg(offset + myrank * lsites * sizeof(fobj));
+        }
+        fin.read((char *)&iodata[0], iodata.size() * sizeof(fobj));
+        assert(fin.fail() == 0);
+        fin.close();
       }
       timer.Stop();
 
@@ -416,19 +424,47 @@ class BinaryIO {
 	assert(0);
 #endif
       } else { 
-	std::ofstream fout; fout.open(file,std::ios::binary|std::ios::out|std::ios::in);
+        
-	std::cout<< GridLogMessage<< "C++ write I/O "<< file<<" : "
+	std::ofstream fout; 
-		 << iodata.size()*sizeof(fobj)<<" bytes"<<std::endl;
+  fout.exceptions ( std::fstream::failbit | std::fstream::badbit );
-	if ( control & BINARYIO_MASTER_APPEND )  {
+  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();
+  timer.Stop();
-    }
+  }
 
     std::cout<<GridLogMessage<<"IOobject: ";
     if ( control & BINARYIO_READ) std::cout << " read  ";
@@ -442,11 +478,14 @@ class BinaryIO {
     //////////////////////////////////////////////////////////////////////////////
     // Safety check
     //////////////////////////////////////////////////////////////////////////////
-    grid->Barrier();
+    // if the data size is 1 we do not want to sum over the MPI ranks
-    grid->GlobalSum(nersc_csum);
+    if (iodata.size() != 1){
-    grid->GlobalXOR(scidac_csuma);
+      grid->Barrier();
-    grid->GlobalXOR(scidac_csumb);
+      grid->GlobalSum(nersc_csum);
-    grid->Barrier();
+      grid->GlobalXOR(scidac_csuma);
+      grid->GlobalXOR(scidac_csumb);
+      grid->Barrier();
+    }
   }
 
   /////////////////////////////////////////////////////////////////////////////
@@ -546,9 +585,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/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,20 +80,39 @@ 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;
+        }
 
-    int def_prec = std::cout.precision();
+        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 << GridLogMessage
+        std::cout.precision(def_prec);
-        << std::setprecision(std::numeric_limits<Real>::digits10 + 1)
+        }
-        << "Topological Charge: [ " << traj << " ] "<< q << std::endl;
-
-    std::cout.precision(def_prec);
     }
 
 };

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);
+    }                 
 }
     
 

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);
   //////////////////////////////////////////////
 
   /////////////////////////////////////////////////////////////