Options to use mersenne twister OR ranlux48 via --enable-rng flag at configure time.

Can save and restore RNG state via new (serial) I/O routines in a NERSC header style file.
Store a Parallel (one per site) and a single serial RNG file.

configure

@@ -749,6 +749,7 @@ enable_openmp
 enable_simd
 enable_precision
 enable_comms
+enable_rng
 enable_chroma
 '
       ac_precious_vars='build_alias
@@ -1395,6 +1396,8 @@ Optional Features:
   --enable-precision=single|double
                           Select default word size of Real
   --enable-comms=none|mpi Select communications
+  --enable-rng=ranlux48|mt19937
+                          Select Random Number Generator to be used
   --enable-chroma         Expect chroma compiled under c++11
 
 Some influential environment variables:
@@ -6387,6 +6390,10 @@ done
 #Please install or provide the correct path to your installation
 #Info at: http://www.mpfr.org/)])
 
+#
+# SIMD instructions selection
+#
+
 # Check whether --enable-simd was given.
 if test "${enable_simd+set}" = set; then :
   enableval=$enable_simd; \
@@ -6519,6 +6526,10 @@ $as_echo "#define GRID_DEFAULT_PRECISION_DOUBLE 1" >>confdefs.h
      ;;
 esac
 
+#
+# Comms selection
+#
+
 # Check whether --enable-comms was given.
 if test "${enable_comms+set}" = set; then :
   enableval=$enable_comms; ac_COMMS=${enable_comms}
@@ -6562,6 +6573,35 @@ else
 fi
 
 
+#
+# RNG selection
+#
+# Check whether --enable-rng was given.
+if test "${enable_rng+set}" = set; then :
+  enableval=$enable_rng; \
+	ac_RNG=${enable_rng}
+else
+  ac_RNG=ranlux48
+fi
+
+case ${ac_RNG} in
+     ranlux48)
+
+$as_echo "#define RNG_RANLUX 1" >>confdefs.h
+
+     ;;
+     mt19937)
+
+$as_echo "#define RNG_MT19937 1" >>confdefs.h
+
+     ;;
+     *)
+     as_fn_error $? "${ac_RNG} unsupported --enable-rng option" "$LINENO" 5;
+     ;;
+esac
+#
+# Chroma regression tests
+#
 # Check whether --enable-chroma was given.
 if test "${enable_chroma+set}" = set; then :
   enableval=$enable_chroma; ac_CHROMA=yes

configure.ac

@@ -66,6 +66,10 @@ AC_CHECK_FUNCS([gettimeofday])
 #Please install or provide the correct path to your installation
 #Info at: http://www.mpfr.org/)])
 
+#
+# SIMD instructions selection
+#
+
 AC_ARG_ENABLE([simd],[AC_HELP_STRING([--enable-simd=SSE4|AVX|AVXFMA4|AVX2|AVX512|IMCI],\
 	[Select instructions to be SSE4.0, AVX 1.0, AVX 2.0+FMA, AVX 512, IMCI])],\
 	[ac_SIMD=${enable_simd}],[ac_SIMD=AVX2])
@@ -160,6 +164,10 @@ case ${ac_PRECISION} in
      ;;
 esac
 
+#
+# Comms selection
+#
+
 AC_ARG_ENABLE([comms],[AC_HELP_STRING([--enable-comms=none|mpi],[Select communications])],[ac_COMMS=${enable_comms}],[ac_COMMS=none])
 
 case ${ac_COMMS} in
@@ -179,6 +187,26 @@ esac
 AM_CONDITIONAL(BUILD_COMMS_MPI,[ test "X${ac_COMMS}X" == "XmpiX" ])
 AM_CONDITIONAL(BUILD_COMMS_NONE,[ test "X${ac_COMMS}X" == "XnoneX" ])
 
+#
+# RNG selection
+#
+AC_ARG_ENABLE([rng],[AC_HELP_STRING([--enable-rng=ranlux48|mt19937],\
+	[Select Random Number Generator to be used])],\
+	[ac_RNG=${enable_rng}],[ac_RNG=ranlux48])
+case ${ac_RNG} in
+     ranlux48)
+     AC_DEFINE([RNG_RANLUX],[1],[RNG_RANLUX] )
+     ;;
+     mt19937)
+     AC_DEFINE([RNG_MT19937],[1],[RNG_MT19937] )
+     ;;
+     *)
+     AC_MSG_ERROR([${ac_RNG} unsupported --enable-rng option]); 
+     ;;
+esac
+#
+# Chroma regression tests
+#
 AC_ARG_ENABLE([chroma],[AC_HELP_STRING([--enable-chroma],[Expect chroma compiled under c++11 ])],ac_CHROMA=yes,ac_CHROMA=no)
 
 case ${ac_CHROMA} in

lib/Config.h.in

@@ -143,6 +143,12 @@
 /* Define to the version of this package. */
 #undef PACKAGE_VERSION
 
+/* RNG_MT19937 */
+#undef RNG_MT19937
+
+/* RNG_RANLUX */
+#undef RNG_RANLUX
+
 /* SSE4 Intrinsics */
 #undef SSE4
 

lib/lattice/Lattice_rng.h

@@ -81,11 +81,39 @@ namespace Grid {
     int _seeded;
     // One generator per site.
     // Uniform and Gaussian distributions from these generators.
-    std::vector<std::ranlux48>             _generators;
+#ifdef RNG_RANLUX
+    typedef uint64_t      RngStateType;
+    typedef std::ranlux48 RngEngine;
+    static const int RngStateCount = 15;
+#else
+    typedef std::mt19937 RngEngine;
+    typedef uint32_t     RngStateType;
+    static const int     RngStateCount = std::mt19937::state_size;
+#endif
+    std::vector<RngEngine>             _generators;
     std::uniform_real_distribution<double> _uniform;
     std::normal_distribution<double>       _gaussian;
 
-
+    void GetState(std::vector<RngStateType> & saved,int gen) {
+      saved.resize(RngStateCount);
+      std::stringstream ss;
+      //      std::cout << _generators[gen]<<std::endl;
+      ss<<_generators[gen];
+      ss.seekg(0,ss.beg);
+      for(int i=0;i<RngStateCount;i++){
+	ss>>saved[i];
+      }
+    }
+    void SetState(std::vector<RngStateType> & saved,int gen){
+      assert(saved.size()==RngStateCount);
+      std::stringstream ss;
+      for(int i=0;i<RngStateCount;i++){
+	ss<< saved[i]<<" ";
+      }
+      ss.seekg(0,ss.beg);
+      ss>>_generators[gen];
+      //      std::cout << _generators[gen]<<std::endl;
+    }
   };
 
   class GridSerialRNG : public GridRNGbase {
@@ -98,7 +126,7 @@ namespace Grid {
     {
       typename source::result_type init = src();
       CartesianCommunicator::BroadcastWorld(0,(void *)&init,sizeof(init));
-      _generators[0] = std::ranlux48(init);
+      _generators[0] = RngEngine(init);
       _seeded=1;
     }    
 
@@ -185,8 +213,6 @@ namespace Grid {
 
   class GridParallelRNG : public GridRNGbase {
   public:
-    // One generator per site.
-    std::vector<std::ranlux48>             _generators;
     
     // Uniform and Gaussian distributions from these generators.
     std::uniform_real_distribution<double> _uniform;
@@ -220,7 +246,7 @@ namespace Grid {
       int gsites = _grid->_gsites;
 
       typename source::result_type init = src();
-      std::ranlux48 pseeder(init);
+      RngEngine pseeder(init);
       std::uniform_int_distribution<uint64_t> ui;
 
       for(int gidx=0;gidx<gsites;gidx++){
@@ -241,7 +267,7 @@ namespace Grid {
 	if( rank == _grid->ThisRank() ){
 	  fixedSeed ssrc(site_seeds);
 	  typename source::result_type sinit = ssrc();
-	  _generators[l_idx] = std::ranlux48(sinit);
+	  _generators[l_idx] = RngEngine(sinit);
 	}
       }
       _seeded=1;

lib/parallelIO/BinaryIO.h

@@ -232,6 +232,110 @@ class BinaryIO {
     return csum;
   }
 
+  static inline uint32_t writeRNGSerial(GridSerialRNG &serial,GridParallelRNG &parallel,std::string file,int offset)
+  {
+    typedef typename GridSerialRNG::RngStateType RngStateType;
+    const int RngStateCount = GridSerialRNG::RngStateCount;
+
+    GridBase *grid = parallel._grid;
+    int gsites = grid->_gsites;
+
+    //////////////////////////////////////////////////
+    // Serialise through node zero
+    //////////////////////////////////////////////////
+    std::cout<< GridLogMessage<< "Serial RNG write I/O "<< file<<std::endl;
+
+    std::ofstream fout;
+    if ( grid->IsBoss() ) {
+      fout.open(file,std::ios::binary|std::ios::out|std::ios::in);
+      fout.seekp(offset);
+    }
+    
+    uint32_t csum=0;
+    std::vector<RngStateType> saved(RngStateCount);
+    int bytes = sizeof(RngStateType)*saved.size();
+    std::vector<int> gcoor;
+
+    for(int gidx=0;gidx<gsites;gidx++){
+
+      int rank,o_idx,i_idx;
+      grid->GlobalIndexToGlobalCoor(gidx,gcoor);
+      grid->GlobalCoorToRankIndex(rank,o_idx,i_idx,gcoor);
+      int l_idx=parallel.generator_idx(o_idx,i_idx);
+
+      if( rank == grid->ThisRank() ){
+	//	std::cout << "rank" << rank<<" Getting state for index "<<l_idx<<std::endl;
+	parallel.GetState(saved,l_idx);
+      }
+
+      grid->Broadcast(rank,(void *)&saved[0],bytes);
+
+      if ( grid->IsBoss() ) {
+	Uint32Checksum((uint32_t *)&saved[0],bytes,csum);
+	fout.write((char *)&saved[0],bytes);
+      }
+
+    }
+
+    if ( grid->IsBoss() ) {
+      serial.GetState(saved,0);
+      Uint32Checksum((uint32_t *)&saved[0],bytes,csum);
+      fout.write((char *)&saved[0],bytes);
+    }
+
+    return csum;
+  }
+  static inline uint32_t readRNGSerial(GridSerialRNG &serial,GridParallelRNG &parallel,std::string file,int offset)
+  {
+    typedef typename GridSerialRNG::RngStateType RngStateType;
+    const int RngStateCount = GridSerialRNG::RngStateCount;
+
+    GridBase *grid = parallel._grid;
+    int gsites = grid->_gsites;
+
+    //////////////////////////////////////////////////
+    // Serialise through node zero
+    //////////////////////////////////////////////////
+    std::cout<< GridLogMessage<< "Serial RNG read I/O "<< file<<std::endl;
+
+    std::ifstream fin(file,std::ios::binary|std::ios::in);
+    fin.seekg(offset);
+    
+    uint32_t csum=0;
+    std::vector<RngStateType> saved(RngStateCount);
+    int bytes = sizeof(RngStateType)*saved.size();
+    std::vector<int> gcoor;
+
+    for(int gidx=0;gidx<gsites;gidx++){
+
+      int rank,o_idx,i_idx;
+      grid->GlobalIndexToGlobalCoor(gidx,gcoor);
+      grid->GlobalCoorToRankIndex(rank,o_idx,i_idx,gcoor);
+      int l_idx=parallel.generator_idx(o_idx,i_idx);
+
+      if ( grid->IsBoss() ) {
+	fin.read((char *)&saved[0],bytes);
+	Uint32Checksum((uint32_t *)&saved[0],bytes,csum);
+      }
+
+      grid->Broadcast(0,(void *)&saved[0],bytes);
+
+      if( rank == grid->ThisRank() ){
+	parallel.SetState(saved,l_idx);
+      }
+
+    }
+
+    if ( grid->IsBoss() ) {
+      fin.read((char *)&saved[0],bytes);
+      serial.SetState(saved,0);
+      Uint32Checksum((uint32_t *)&saved[0],bytes,csum);
+    }
+
+    return csum;
+  }
+
+
   template<class vobj,class fobj,class munger>
   static inline uint32_t readObjectParallel(Lattice<vobj> &Umu,std::string file,munger munge,int offset,const std::string &format)
   {
@@ -284,6 +388,7 @@ class BinaryIO {
     {
       uint32_t tmp = IOnode;
       grid->GlobalSum(tmp);
+      std::cout<< std::dec ;
       std::cout<< GridLogMessage<< "Parallel read I/O to "<< file << " with " <<tmp<< " IOnodes for subslice ";
       for(int d=0;d<grid->_ndimension;d++){
 	std::cout<< range[d];

lib/parallelIO/NerscIO.h

@@ -181,9 +181,12 @@ inline void NerscMachineCharacteristics(NerscField &header)
 class NerscIO : public BinaryIO { 
  public:
 
+  static inline void truncate(std::string file){
+    std::ofstream fout(file,std::ios::out);
+  }
   static inline unsigned int writeHeader(NerscField &field,std::string file)
   {
-    std::ofstream fout(file,std::ios::out);
+    std::ofstream fout(file,std::ios::out|std::ios::in);
   
     fout.seekp(0,std::ios::beg);
     fout << "BEGIN_HEADER"      << std::endl;
@@ -201,8 +204,7 @@ class NerscIO : public BinaryIO {
       fout << "BOUNDARY_"<<i+1<<" = " << field.boundary[i] << std::endl;
     }
 
-    fout << "CHECKSUM = "<< std::hex << std::setw(10) << field.checksum << std::endl;
+    fout << "CHECKSUM = "<< std::hex << std::setw(10) << field.checksum << std::dec<<std::endl;
-    fout << std::dec;
 
     fout << "ENSEMBLE_ID = "     << field.ensemble_id      << std::endl;
     fout << "ENSEMBLE_LABEL = "  << field.ensemble_label   << std::endl;
@@ -287,9 +289,6 @@ static inline int readHeader(std::string file,GridBase *grid,  NerscField &field
   return field.data_start;
 }
 
-
-
-
 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 // Now the meat: the object readers
 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -377,6 +376,8 @@ static inline void writeConfiguration(Lattice<iLorentzColourMatrix<vsimd> > &Umu
   uint32_t csum;
   int offset;
   
+  truncate(file);
+
   if ( two_row ) { 
 
     header.floating_point = std::string("IEEE64BIG");
@@ -392,8 +393,7 @@ static inline void writeConfiguration(Lattice<iLorentzColourMatrix<vsimd> > &Umu
 
     
     std::cout << GridLogMessage << " TESTING PARALLEL WRITE offsets " << offset1 << " "<< offset << std::endl;
-    std::cout << GridLogMessage <<std::hex<< " TESTING PARALLEL WRITE csums   " << csum1 << " "<< csum << std::endl;
+    std::cout << GridLogMessage << " TESTING PARALLEL WRITE csums   " << csum1 << " "<<std::hex<< csum << std::dec<< std::endl;
-    std::cout << std::dec;
 
     assert(offset1==offset);  
     assert(csum1==csum);  
@@ -410,7 +410,81 @@ static inline void writeConfiguration(Lattice<iLorentzColourMatrix<vsimd> > &Umu
   std::cout<<GridLogMessage <<"Written NERSC Configuration "<<file<< " checksum "<<std::hex<<csum<< std::dec<<" plaq "<< header.plaquette <<std::endl;
 
  }
+
+
+    ///////////////////////////////
+    // RNG state
+    ///////////////////////////////
+static inline void writeRNGState(GridSerialRNG &serial,GridParallelRNG &parallel,std::string file)
+{
+  typedef typename GridParallelRNG::RngStateType RngStateType;
+
+  // Following should become arguments
+  NerscField header;
+  header.sequence_number = 1;
+  header.ensemble_id     = "UKQCD";
+  header.ensemble_label  = "DWF";
+
+  GridBase *grid = parallel._grid;
+
+  NerscGrid(grid,header);
+  header.link_trace=0.0;
+  header.plaquette=0.0;
+  NerscMachineCharacteristics(header);
+
+  uint32_t csum;
+  int offset;
+  
+#ifdef RNG_RANLUX
+    header.floating_point = std::string("UINT64");
+    header.data_type      = std::string("RANLUX48");
+#else
+    header.floating_point = std::string("UINT32");
+    header.data_type      = std::string("MT19937");
+#endif
+
+  truncate(file);
+  offset = writeHeader(header,file);
+  csum=BinaryIO::writeRNGSerial(serial,parallel,file,offset);
+  header.checksum = csum;
+  offset = writeHeader(header,file);
+
+  std::cout<<GridLogMessage <<"Written NERSC RNG STATE "<<file<< " checksum "<<std::hex<<csum<<std::dec<<std::endl;
+
+ }
+    
+static inline void readRNGState(GridSerialRNG &serial,GridParallelRNG & parallel,NerscField& header,std::string file)
+{
+  typedef typename GridParallelRNG::RngStateType RngStateType;
+
+  GridBase *grid = parallel._grid;
+
+  int offset = readHeader(file,grid,header);
+
+  NerscField clone(header);
+
+  std::string format(header.floating_point);
+  std::string data_type(header.data_type);
+
+#ifdef RNG_RANLUX
+  assert(format == std::string("UINT64"));
+  assert(data_type == std::string("RANLUX48"));
+#else
+  assert(format == std::string("UINT32"));
+  assert(data_type == std::string("MT19937"));
+#endif
+
+  // depending on datatype, set up munger;
+  // munger is a function of <floating point, Real, data_type>
+  uint32_t csum=BinaryIO::readRNGSerial(serial,parallel,file,offset);
+
+  assert(csum == header.checksum );
+
+  std::cout<<GridLogMessage <<"Read NERSC RNG file "<<file<< " format "<< data_type <<std::endl;
+}
+
 };
 
+
 }}
 #endif

tests/Test_main.cc

@@ -56,19 +56,25 @@ int main (int argc, char ** argv)
     GridRedBlackCartesian rbFine(latt_size,simd_layout,mpi_layout);
     GridParallelRNG       FineRNG(&Fine);
     GridSerialRNG       SerialRNG;
+    GridSerialRNG       SerialRNG1;
 
     FineRNG.SeedRandomDevice();
+    SerialRNG.SeedRandomDevice();
 
     std::cout <<"SerialRNG" << SerialRNG._generators[0] <<std::endl;
-    std::stringstream output(std::stringstream::out|std::stringstream::binary);
-    output <<SerialRNG._generators[0]<<std::endl;
-    std::cout << output.str();
 
-    {
+    std::vector<typename GridSerialRNG::RngStateType> saved;
-      std::ofstream of("rngstate",std::ios::out|std::ios::binary);
+    SerialRNG.GetState(saved,0);
-      of << SerialRNG._generators[0];
+    SerialRNG1.SetState(saved,0);
+
+    RealD dd1,dd2;
+
+    std::cout << "Testing RNG state save restore"<<std::endl;
+    for(int i=0;i<10;i++){
+      random(SerialRNG,dd1);
+      random(SerialRNG1,dd2);
+      std::cout << "i "<<i<<" "<<dd1<< " " <<dd2<<std::endl;
     }
-
     LatticeColourMatrix Foo(&Fine);
     LatticeColourMatrix Bar(&Fine);
 

tests/Test_nersc_io.cc

@@ -20,6 +20,30 @@ int main (int argc, char ** argv)
   GridCartesian     Fine(latt_size,simd_layout,mpi_layout);
   GridCartesian     Coarse(clatt_size,simd_layout,mpi_layout);
 
+  GridParallelRNG   pRNGa(&Fine);
+  GridParallelRNG   pRNGb(&Fine);
+  GridSerialRNG     sRNGa;
+  GridSerialRNG     sRNGb;
+
+  pRNGa.SeedRandomDevice();
+  sRNGa.SeedRandomDevice();
+  
+  std::string rfile("./ckpoint_rng.4000");
+  NerscIO::writeRNGState(sRNGa,pRNGa,rfile);
+  NerscField rngheader;
+  NerscIO::readRNGState (sRNGb,pRNGb,rngheader,rfile);
+
+  LatticeComplex tmpa(&Fine); random(pRNGa,tmpa);
+  LatticeComplex tmpb(&Fine); random(pRNGb,tmpb);
+  tmpa = tmpa - tmpb;
+  std::cout << " difference between restored randoms and orig "<<norm2( tmpa ) <<" / "<< norm2(tmpb)<<std::endl;
+
+  ComplexD a,b;
+
+  random(sRNGa,a);
+  random(sRNGb,b);
+  std::cout << " serial RNG numbers "<<a<<" "<<b<<std::endl;
+
   LatticeGaugeField Umu(&Fine);
   LatticeGaugeField Umu_diff(&Fine);
   LatticeGaugeField Umu_saved(&Fine);