In the middle of adding a general binary writer

lib/parallelIO/BinaryIO.h

@@ -539,6 +539,12 @@ class BinaryIO {
     int ieee64big = (format == std::string("IEEE64BIG"));
     int ieee64    = (format == std::string("IEEE64"));
 
+    if(!(ieee32big || ieee32 || ieee64big || ieee64)){
+      std::cout << GridLogError << "Unrecognized file format " << format << std::endl;
+      std::cout << GridLogError << "Allowed: IEEE32BIG | IEEE32 | IEEE64BIG | IEEE64" << std::endl; 
+      exit(0);
+    }
+
     int nd = grid->_ndimension;
     for(int d=0;d<nd;d++){
       assert(grid->CheckerBoarded(d) == 0);

lib/parallelIO/NerscIO.h

@@ -137,51 +137,53 @@ inline void NerscMachineCharacteristics(NerscField &header)
       const int y=1;
       const int z=2;
       for(int mu=0;mu<4;mu++){
-	cm(mu)()(2,x) = adj(cm(mu)()(0,y)*cm(mu)()(1,z)-cm(mu)()(0,z)*cm(mu)()(1,y)); //x= yz-zy
+  cm(mu)()(2,x) = adj(cm(mu)()(0,y)*cm(mu)()(1,z)-cm(mu)()(0,z)*cm(mu)()(1,y)); //x= yz-zy
-	cm(mu)()(2,y) = adj(cm(mu)()(0,z)*cm(mu)()(1,x)-cm(mu)()(0,x)*cm(mu)()(1,z)); //y= zx-xz
+  cm(mu)()(2,y) = adj(cm(mu)()(0,z)*cm(mu)()(1,x)-cm(mu)()(0,x)*cm(mu)()(1,z)); //y= zx-xz
-	cm(mu)()(2,z) = adj(cm(mu)()(0,x)*cm(mu)()(1,y)-cm(mu)()(0,y)*cm(mu)()(1,x)); //z= xy-yx
+  cm(mu)()(2,z) = adj(cm(mu)()(0,x)*cm(mu)()(1,y)-cm(mu)()(0,y)*cm(mu)()(1,x)); //z= xy-yx
       }
     }
 
     template<class fobj,class sobj>
     struct NerscSimpleMunger{
-
+      void operator()(fobj &in, sobj &out, uint32_t &csum) {
-      void operator() (fobj &in,sobj &out,uint32_t &csum){
+        for (int mu = 0; mu < Nd; mu++) {
-
+          for (int i = 0; i < Nc; i++) {
-      for(int mu=0;mu<4;mu++){
+            for (int j = 0; j < Nc; j++) {
-      for(int i=0;i<3;i++){
+              out(mu)()(i, j) = in(mu)()(i, j);
-      for(int j=0;j<3;j++){
+            }
-	out(mu)()(i,j) = in(mu)()(i,j);
+          }
-      }}}
+        }
-      NerscChecksum((uint32_t *)&in,sizeof(in),csum); 
+        NerscChecksum((uint32_t *)&in, sizeof(in), csum);
       };
     };
 
-    template<class fobj,class sobj>
+    template <class fobj, class sobj>
-    struct NerscSimpleUnmunger{
+    struct NerscSimpleUnmunger {
-      void operator() (sobj &in,fobj &out,uint32_t &csum){
+      void operator()(sobj &in, fobj &out, uint32_t &csum) {
-	for(int mu=0;mu<Nd;mu++){
+        for (int mu = 0; mu < Nd; mu++) {
-	for(int i=0;i<Nc;i++){
+          for (int i = 0; i < Nc; i++) {
-	for(int j=0;j<Nc;j++){
+            for (int j = 0; j < Nc; j++) {
-	  out(mu)()(i,j) = in(mu)()(i,j);
+              out(mu)()(i, j) = in(mu)()(i, j);
-	}}}
+            }
-	NerscChecksum((uint32_t *)&out,sizeof(out),csum); 
+          }
+        }
+        NerscChecksum((uint32_t *)&out, sizeof(out), csum);
       };
     };
- 
+
     template<class fobj,class sobj>
     struct Nersc3x2munger{
       void operator() (fobj &in,sobj &out,uint32_t &csum){
      
-	NerscChecksum((uint32_t *)&in,sizeof(in),csum); 
+  NerscChecksum((uint32_t *)&in,sizeof(in),csum); 
 
-	for(int mu=0;mu<4;mu++){
+  for(int mu=0;mu<4;mu++){
-	  for(int i=0;i<2;i++){
+    for(int i=0;i<2;i++){
-	    for(int j=0;j<3;j++){
+      for(int j=0;j<3;j++){
-	      out(mu)()(i,j) = in(mu)(i)(j);
+        out(mu)()(i,j) = in(mu)(i)(j);
-	    }}
+      }}
-	}
+  }
-	reconstruct3(out);
+  reconstruct3(out);
       }
     };
 
@@ -191,14 +193,14 @@ inline void NerscMachineCharacteristics(NerscField &header)
       void operator() (sobj &in,fobj &out,uint32_t &csum){
 
 
-	for(int mu=0;mu<4;mu++){
+  for(int mu=0;mu<4;mu++){
-	  for(int i=0;i<2;i++){
+    for(int i=0;i<2;i++){
-	    for(int j=0;j<3;j++){
+      for(int j=0;j<3;j++){
-	      out(mu)(i)(j) = in(mu)()(i,j);
+        out(mu)(i)(j) = in(mu)()(i,j);
-	    }}
+      }}
-	}
+  }
 
-	NerscChecksum((uint32_t *)&out,sizeof(out),csum); 
+  NerscChecksum((uint32_t *)&out,sizeof(out),csum); 
 
       }
     };
@@ -346,24 +348,24 @@ static inline void readConfiguration(Lattice<iLorentzColourMatrix<vsimd> > &Umu,
   if ( header.data_type == std::string("4D_SU3_GAUGE") ) {
     if ( ieee32 || ieee32big ) {
       //      csum=BinaryIO::readObjectSerial<iLorentzColourMatrix<vsimd>, LorentzColour2x3F> 
-	csum=BinaryIO::readObjectParallel<iLorentzColourMatrix<vsimd>, LorentzColour2x3F> 
+  csum=BinaryIO::readObjectParallel<iLorentzColourMatrix<vsimd>, LorentzColour2x3F> 
-	(Umu,file,Nersc3x2munger<LorentzColour2x3F,LorentzColourMatrix>(), offset,format);
+  (Umu,file,Nersc3x2munger<LorentzColour2x3F,LorentzColourMatrix>(), offset,format);
     }
     if ( ieee64 || ieee64big ) {
       //csum=BinaryIO::readObjectSerial<iLorentzColourMatrix<vsimd>, LorentzColour2x3D> 
       csum=BinaryIO::readObjectParallel<iLorentzColourMatrix<vsimd>, LorentzColour2x3D> 
-      	(Umu,file,Nersc3x2munger<LorentzColour2x3D,LorentzColourMatrix>(),offset,format);
+        (Umu,file,Nersc3x2munger<LorentzColour2x3D,LorentzColourMatrix>(),offset,format);
     }
   } else if ( header.data_type == std::string("4D_SU3_GAUGE_3x3") ) {
     if ( ieee32 || ieee32big ) {
       //csum=BinaryIO::readObjectSerial<iLorentzColourMatrix<vsimd>,LorentzColourMatrixF>
       csum=BinaryIO::readObjectParallel<iLorentzColourMatrix<vsimd>,LorentzColourMatrixF>
-	(Umu,file,NerscSimpleMunger<LorentzColourMatrixF,LorentzColourMatrix>(),offset,format);
+  (Umu,file,NerscSimpleMunger<LorentzColourMatrixF,LorentzColourMatrix>(),offset,format);
     }
     if ( ieee64 || ieee64big ) {
       //      csum=BinaryIO::readObjectSerial<iLorentzColourMatrix<vsimd>,LorentzColourMatrixD>
       csum=BinaryIO::readObjectParallel<iLorentzColourMatrix<vsimd>,LorentzColourMatrixD>
-	(Umu,file,NerscSimpleMunger<LorentzColourMatrixD,LorentzColourMatrix>(),offset,format);
+  (Umu,file,NerscSimpleMunger<LorentzColourMatrixD,LorentzColourMatrix>(),offset,format);
     }
   } else {
     assert(0);

lib/qcd/hmc/BinaryCheckpointer.h

@@ -0,0 +1,106 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/qcd/hmc/NerscCheckpointer.h
+
+    Copyright (C) 2015
+
+Author: paboyle <paboyle@ph.ed.ac.uk>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+    *************************************************************************************/
+    /*  END LEGAL */
+#ifndef NERSC_CHECKPOINTER
+#define NERSC_CHECKPOINTER
+
+#include <string>
+#include <iostream>
+#include <sstream>
+
+
+namespace Grid{
+  namespace QCD{
+    
+    template<class fobj,class sobj>
+    struct BinarySimpleUnmunger{
+      void operator() (sobj &in,fobj &out,uint32_t &csum){
+        fobj = sobj;
+        csum =0;
+    };
+
+
+    // Only for gauge fields
+    template<class Impl> 
+    class BinaryHmcCheckpointer : public HmcObservable<typename Impl::Field> {
+    private:
+      std::string configStem;
+      std::string rngStem;
+      int SaveInterval;
+    public:
+      INHERIT_FIELD_TYPES(Impl); // The Field is a Lattice object
+
+      typedef typename Field::vector_object vobj;
+      typedef typename vobj::scalar_object sobj;
+ 
+
+      BinaryHmcCheckpointer(std::string cf, std::string rn,int savemodulo) {
+        configStem  = cf;
+        rngStem     = rn;
+        SaveInterval= savemodulo;
+      };
+
+      void TrajectoryComplete(int traj, Field &U, GridSerialRNG &sRNG,
+                              GridParallelRNG &pRNG) {
+        if ((traj % SaveInterval) == 0) {
+          std::string rng;
+          {
+            std::ostringstream os;
+            os << rngStem << "." << traj;
+            rng = os.str();
+          }
+          std::string config;
+          {
+            std::ostringstream os;
+            os << configStem << "." << traj;
+            config = os.str();
+          }
+
+          int prec = getPrecision<Field>::value; //get precision of oject
+          std::string floating_point = "IEEE64BIG"; //default precision
+          if (prec == 1)
+            floating_point = "IEEE32BIG"
+
+          BinarySimpleUnmunger<sobj, sobj> munge;
+          BinaryIO::writeRNGSerial(sRNG, pRNG, rng, 0);
+          BinaryIO::writeObjectParallel<vobj, sobj>(U, config, munge, 0 ,floating_point);
+        }
+      };
+
+      void CheckpointRestore(int traj, GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG & pRNG ){
+
+  std::string rng;   { std::ostringstream os; os << rngStem     <<"."<< traj; rng = os.str(); }
+  std::string config;{ std::ostringstream os; os << configStem  <<"."<< traj; config = os.str();}
+
+  NerscField header;
+  NerscIO::readRNGState(sRNG,pRNG,header,rng);
+  NerscIO::readConfiguration(U,header,config);
+      };
+
+    };
+}}
+#endif

lib/qcd/hmc/integrators/Integrator.h

@@ -41,8 +41,6 @@ directory
 #ifndef INTEGRATOR_INCLUDED
 #define INTEGRATOR_INCLUDED
 
-// class Observer;
-
 #include <memory>
 
 namespace Grid {
@@ -77,24 +75,20 @@ struct IntegratorParameters {
 template <class FieldImplementation, class SmearingPolicy, class RepresentationPolicy>
 class Integrator {
  protected:
-  typedef IntegratorParameters ParameterType;
   typedef typename FieldImplementation::Field MomentaField;  //for readability
   typedef typename FieldImplementation::Field Field;
 
+  int levels;  // number of integration levels
+  double t_U;  // Track time passing on each level and for U and for P
+  std::vector<double> t_P;  
+
+  MomentaField P;
+  SmearingPolicy& Smearer;
+  RepresentationPolicy Representations;
   IntegratorParameters Params;
 
   const ActionSet<Field, RepresentationPolicy> as;
 
-  int levels;  //
-  double t_U;  // Track time passing on each level and for U and for P
-  std::vector<double> t_P;  //
-
-  MomentaField P;
-
-  SmearingPolicy& Smearer;
-
-  RepresentationPolicy Representations;
-
   void update_P(Field& U, int level, double ep) {
     t_P[level] += ep;
     update_P(P, U, level, ep);