From ae4de947989d1c9299b7dbeb8c1a570f745a84d7 Mon Sep 17 00:00:00 2001
From: paboyle <paboyle@ph.ed.ac.uk>
Date: Sun, 18 Jun 2017 00:11:23 +0100
Subject: [PATCH] SciDAC I/O support

---
 lib/parallelIO/MetaData.h | 124 ++++++++++++++++++++++++++++++++++----
 lib/parallelIO/NerscIO.h  |   4 +-
 2 files changed, 114 insertions(+), 14 deletions(-)

diff --git a/lib/parallelIO/MetaData.h b/lib/parallelIO/MetaData.h
index 1bad07f2..6d45d0a5 100644
--- a/lib/parallelIO/MetaData.h
+++ b/lib/parallelIO/MetaData.h
@@ -38,9 +38,24 @@
 
 namespace Grid {
 
-    ////////////////////////////////////////////////////////////////////////////////
-    // header specification/interpretation
-    ////////////////////////////////////////////////////////////////////////////////
+  ///////////////////////////////////////////////////////
+  // Precision mapping
+  ///////////////////////////////////////////////////////
+  template<class vobj> static std::string getFormatString (void)
+  {
+    std::string format;
+    typedef typename getPrecision<vobj>::real_scalar_type stype;
+    if ( sizeof(stype) == sizeof(float) ) {
+      format = std::string("IEEE32BIG");
+    }
+    if ( sizeof(stype) == sizeof(double) ) {
+      format = std::string("IEEE64BIG");
+    }
+    return format;
+  }
+  ////////////////////////////////////////////////////////////////////////////////
+  // header specification/interpretation
+  ////////////////////////////////////////////////////////////////////////////////
     class FieldMetaData : Serializable {
     public:
 
@@ -66,8 +81,15 @@ namespace Grid {
 				      std::string, creation_date,
 				      std::string, archive_date,
 				      std::string, floating_point);
+      FieldMetaData(void) { 
+	nd=4;
+	dimension.resize(4);
+	boundary.resize(4);
+      }
     };
 
+
+
   namespace QCD {
 
     using namespace Grid;
@@ -89,13 +111,6 @@ namespace Grid {
 	header.boundary[d] = std::string("PERIODIC");
       }
     }
-    template<class GaugeField>
-    inline void GaugeStatistics(GaugeField & data,FieldMetaData &header)
-    {
-      // How to convert data precision etc...
-      header.link_trace=Grid::QCD::WilsonLoops<PeriodicGimplR>::linkTrace(data);
-      header.plaquette =Grid::QCD::WilsonLoops<PeriodicGimplR>::avgPlaquette(data);
-    }
 
     inline void MachineCharacteristics(FieldMetaData &header)
     {
@@ -133,7 +148,7 @@ namespace Grid {
 	s << "BOUNDARY_"<<i+1<<" = " << field.boundary[i] << std::endl;	\
       }									\
 									\
-      s << "CHECKSUM (NERSC) = "<< std::hex << std::setw(10) << field.checksum << std::dec<<std::endl; \
+      s << "CHECKSUM = "<< std::hex << std::setw(10) << field.checksum << std::dec<<std::endl; \
       s << "SCIDAC_CHECKSUMA = "<< std::hex << std::setw(10) << field.scidac_checksuma << std::dec<<std::endl; \
       s << "SCIDAC_CHECKSUMB = "<< std::hex << std::setw(10) << field.scidac_checksumb << std::dec<<std::endl; \
       s << "ENSEMBLE_ID = "     << field.ensemble_id      << std::endl;	\
@@ -146,6 +161,48 @@ namespace Grid {
       s << "FLOATING_POINT = "  << field.floating_point   << std::endl;	\
       s << "END_HEADER"         << std::endl;
 
+template<class vobj> inline void PrepareMetaData(Lattice<vobj> & field, FieldMetaData &header)
+{
+  GridBase *grid = field._grid;
+  std::string format = getFormatString<vobj>();
+   header.floating_point = format;
+   header.checksum = 0x0; // Nersc checksum unused in ILDG, Scidac
+   GridMetaData(grid,header); 
+   MachineCharacteristics(header);
+ }
+ inline void GaugeStatistics(Lattice<vLorentzColourMatrixF> & data,FieldMetaData &header)
+ {
+   // How to convert data precision etc...
+   header.link_trace=Grid::QCD::WilsonLoops<PeriodicGimplF>::linkTrace(data);
+   header.plaquette =Grid::QCD::WilsonLoops<PeriodicGimplF>::avgPlaquette(data);
+ }
+ inline void GaugeStatistics(Lattice<vLorentzColourMatrixD> & data,FieldMetaData &header)
+ {
+   // How to convert data precision etc...
+   header.link_trace=Grid::QCD::WilsonLoops<PeriodicGimplD>::linkTrace(data);
+   header.plaquette =Grid::QCD::WilsonLoops<PeriodicGimplD>::avgPlaquette(data);
+ }
+ template<> inline void PrepareMetaData<vLorentzColourMatrixF>(Lattice<vLorentzColourMatrixF> & field, FieldMetaData &header)
+ {
+   
+   GridBase *grid = field._grid;
+   std::string format = getFormatString<vLorentzColourMatrixF>();
+   header.floating_point = format;
+   header.checksum = 0x0; // Nersc checksum unused in ILDG, Scidac
+   GridMetaData(grid,header); 
+   GaugeStatistics(field,header);
+   MachineCharacteristics(header);
+ }
+ template<> inline void PrepareMetaData<vLorentzColourMatrixD>(Lattice<vLorentzColourMatrixD> & field, FieldMetaData &header)
+ {
+   GridBase *grid = field._grid;
+   std::string format = getFormatString<vLorentzColourMatrixD>();
+   header.floating_point = format;
+   header.checksum = 0x0; // Nersc checksum unused in ILDG, Scidac
+   GridMetaData(grid,header); 
+   GaugeStatistics(field,header);
+   MachineCharacteristics(header);
+ }
 
     //////////////////////////////////////////////////////////////////////
     // Utilities ; these are QCD aware
@@ -171,6 +228,48 @@ namespace Grid {
     typedef iLorentzColour2x3<ComplexF> LorentzColour2x3F;
     typedef iLorentzColour2x3<ComplexD> LorentzColour2x3D;
 
+/////////////////////////////////////////////////////////////////////////////////
+// Simple classes for precision conversion
+/////////////////////////////////////////////////////////////////////////////////
+template <class fobj, class sobj>
+struct BinarySimpleUnmunger {
+  typedef typename getPrecision<fobj>::real_scalar_type fobj_stype;
+  typedef typename getPrecision<sobj>::real_scalar_type sobj_stype;
+  
+  void operator()(sobj &in, fobj &out) {
+    // take word by word and transform accoding to the status
+    fobj_stype *out_buffer = (fobj_stype *)&out;
+    sobj_stype *in_buffer = (sobj_stype *)&in;
+    size_t fobj_words = sizeof(out) / sizeof(fobj_stype);
+    size_t sobj_words = sizeof(in) / sizeof(sobj_stype);
+    assert(fobj_words == sobj_words);
+    
+    for (unsigned int word = 0; word < sobj_words; word++)
+      out_buffer[word] = in_buffer[word];  // type conversion on the fly
+    
+  }
+};
+
+template <class fobj, class sobj>
+struct BinarySimpleMunger {
+  typedef typename getPrecision<fobj>::real_scalar_type fobj_stype;
+  typedef typename getPrecision<sobj>::real_scalar_type sobj_stype;
+
+  void operator()(fobj &in, sobj &out) {
+    // take word by word and transform accoding to the status
+    fobj_stype *in_buffer = (fobj_stype *)&in;
+    sobj_stype *out_buffer = (sobj_stype *)&out;
+    size_t fobj_words = sizeof(in) / sizeof(fobj_stype);
+    size_t sobj_words = sizeof(out) / sizeof(sobj_stype);
+    assert(fobj_words == sobj_words);
+    
+    for (unsigned int word = 0; word < sobj_words; word++)
+      out_buffer[word] = in_buffer[word];  // type conversion on the fly
+    
+  }
+};
+
+
     template<class fobj,class sobj>
     struct GaugeSimpleMunger{
       void operator()(fobj &in, sobj &out) {
@@ -220,6 +319,7 @@ namespace Grid {
 	}
       }
     };
-
   }
+
+
 }
diff --git a/lib/parallelIO/NerscIO.h b/lib/parallelIO/NerscIO.h
index cc37b537..786839f2 100644
--- a/lib/parallelIO/NerscIO.h
+++ b/lib/parallelIO/NerscIO.h
@@ -179,7 +179,7 @@ namespace Grid {
 	assert(0);
       }
 
-      GaugeStatistics<GaugeField>(Umu,clone);
+      GaugeStatistics(Umu,clone);
 
       std::cout<<GridLogMessage <<"NERSC Configuration "<<file<<" checksum "<<std::hex<<nersc_csum<< std::dec
 	       <<" header   "<<std::hex<<header.checksum<<std::dec <<std::endl;
@@ -233,7 +233,7 @@ namespace Grid {
 
 	GridMetaData(grid,header);
 	assert(header.nd==4);
-	GaugeStatistics<GaugeField>(Umu,header);
+	GaugeStatistics(Umu,header);
 	MachineCharacteristics(header);
 
 	int offset;