Merge branch 'release/v0.7.0'

Double precision compile fix

TODO

@@ -2,9 +2,9 @@ TODO:
 ---------------
 
 Peter's work list:
-1)- Precision conversion and sort out localConvert      <-- 
+2)- Precision conversion and sort out localConvert      <-- 
-2)- Remove DenseVector, DenseMatrix; Use Eigen instead. <-- 
+3)- Remove DenseVector, DenseMatrix; Use Eigen instead. <-- started 
-
+4)- Binary I/O speed up & x-strips
 -- Profile CG, BlockCG, etc... Flop count/rate -- PARTIAL, time but no flop/s yet
 -- Physical propagator interface
 -- Conserved currents
@@ -13,7 +13,6 @@ Peter's work list:
 -- HDCR resume
 
 Recent DONE 
--- Binary I/O speed up & x-strips                      <-- DONE
 -- Cut down the exterior overhead                      <-- DONE
 -- Interior legs from SHM comms                        <-- DONE
 -- Half-precision comms                                <-- DONE

benchmarks/Benchmark_comms.cc

@@ -31,32 +31,6 @@ using namespace std;
 using namespace Grid;
 using namespace Grid::QCD;
 
-struct time_statistics{
-  double mean;
-  double err;
-  double min;
-  double max;
-
-  void statistics(std::vector<double> v){
-      double sum = std::accumulate(v.begin(), v.end(), 0.0);
-      mean = sum / v.size();
-
-      std::vector<double> diff(v.size());
-      std::transform(v.begin(), v.end(), diff.begin(), [=](double x) { return x - mean; });
-      double sq_sum = std::inner_product(diff.begin(), diff.end(), diff.begin(), 0.0);
-      err = std::sqrt(sq_sum / (v.size()*(v.size() - 1)));
-
-      auto result = std::minmax_element(v.begin(), v.end());
-      min = *result.first;
-      max = *result.second;
-}
-};
-
-void header(){
-  std::cout <<GridLogMessage << " L  "<<"\t"<<" Ls  "<<"\t"
-            <<std::setw(11)<<"bytes"<<"MB/s uni (err/min/max)"<<"\t\t"<<"MB/s bidi (err/min/max)"<<std::endl;
-};
-
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
@@ -66,19 +40,15 @@ int main (int argc, char ** argv)
   int threads = GridThread::GetThreads();
   std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
 
-  int Nloop=100;
+  int Nloop=10;
   int nmu=0;
-  int maxlat=24;
   for(int mu=0;mu<Nd;mu++) if (mpi_layout[mu]>1) nmu++;
 
-  std::cout << GridLogMessage << "Number of iterations to average: "<< Nloop << std::endl;
-  std::vector<double> t_time(Nloop);
-  time_statistics timestat;
-
   std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
   std::cout<<GridLogMessage << "= Benchmarking concurrent halo exchange in "<<nmu<<" dimensions"<<std::endl;
   std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
-  header();
+  std::cout<<GridLogMessage << "  L  "<<"\t\t"<<" Ls  "<<"\t\t"<<"bytes"<<"\t\t"<<"MB/s uni"<<"\t\t"<<"MB/s bidi"<<std::endl;
+  int maxlat=24;
   for(int lat=4;lat<=maxlat;lat+=4){
     for(int Ls=8;Ls<=32;Ls*=2){
 
@@ -88,9 +58,6 @@ int main (int argc, char ** argv)
       				    lat*mpi_layout[3]});
 
       GridCartesian     Grid(latt_size,simd_layout,mpi_layout);
-      RealD Nrank = Grid._Nprocessors;
-      RealD Nnode = Grid.NodeCount();
-      RealD ppn = Nrank/Nnode;
 
       std::vector<std::vector<HalfSpinColourVectorD> > xbuf(8,std::vector<HalfSpinColourVectorD>(lat*lat*lat*Ls));
       std::vector<std::vector<HalfSpinColourVectorD> > rbuf(8,std::vector<HalfSpinColourVectorD>(lat*lat*lat*Ls));
@@ -98,8 +65,8 @@ int main (int argc, char ** argv)
       int ncomm;
       int bytes=lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD);
 
-      for(int i=0;i<Nloop;i++){
       double start=usecond();
+      for(int i=0;i<Nloop;i++){
 
 	std::vector<CartesianCommunicator::CommsRequest_t> requests;
 
@@ -135,24 +102,18 @@ int main (int argc, char ** argv)
 	}
 	Grid.SendToRecvFromComplete(requests);
 	Grid.Barrier();
-	double stop=usecond();
+
-	t_time[i] = stop-start; // microseconds
       }
+      double stop=usecond();
 
-      timestat.statistics(t_time);
+      double dbytes    = bytes;
-
+      double xbytes    = Nloop*dbytes*2.0*ncomm;
-      double dbytes    = bytes*ppn;
-      double xbytes    = dbytes*2.0*ncomm;
       double rbytes    = xbytes;
       double bidibytes = xbytes+rbytes;
 
-      std::cout<<GridLogMessage << std::setw(4) << lat<<"\t"<<Ls<<"\t"
+      double time = stop-start; // microseconds
-               <<std::setw(11) << bytes<< std::fixed << std::setprecision(1) << std::setw(7)
-               <<std::right<< xbytes/timestat.mean<<"  "<< xbytes*timestat.err/(timestat.mean*timestat.mean)<< " "
-               <<xbytes/timestat.max <<" "<< xbytes/timestat.min  
-               << "\t\t"<<std::setw(7)<< bidibytes/timestat.mean<< "  " << bidibytes*timestat.err/(timestat.mean*timestat.mean) << " "
-               << bidibytes/timestat.max << " " << bidibytes/timestat.min << std::endl;
 
+      std::cout<<GridLogMessage << lat<<"\t\t"<<Ls<<"\t\t"<<bytes<<"\t\t"<<xbytes/time<<"\t\t"<<bidibytes/time<<std::endl;
     }
   }    
 
@@ -160,7 +121,8 @@ int main (int argc, char ** argv)
   std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
   std::cout<<GridLogMessage << "= Benchmarking sequential halo exchange in "<<nmu<<" dimensions"<<std::endl;
   std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
-  header();
+  std::cout<<GridLogMessage << "  L  "<<"\t\t"<<" Ls  "<<"\t\t"<<"bytes"<<"\t\t"<<"MB/s uni"<<"\t\t"<<"MB/s bidi"<<std::endl;
+
 
   for(int lat=4;lat<=maxlat;lat+=4){
     for(int Ls=8;Ls<=32;Ls*=2){
@@ -168,9 +130,6 @@ int main (int argc, char ** argv)
       std::vector<int> latt_size  ({lat,lat,lat,lat});
 
       GridCartesian     Grid(latt_size,simd_layout,mpi_layout);
-      RealD Nrank = Grid._Nprocessors;
-      RealD Nnode = Grid.NodeCount();
-      RealD ppn = Nrank/Nnode;
 
       std::vector<std::vector<HalfSpinColourVectorD> > xbuf(8,std::vector<HalfSpinColourVectorD>(lat*lat*lat*Ls));
       std::vector<std::vector<HalfSpinColourVectorD> > rbuf(8,std::vector<HalfSpinColourVectorD>(lat*lat*lat*Ls));
@@ -179,8 +138,8 @@ int main (int argc, char ** argv)
       int ncomm;
       int bytes=lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD);
 
-      for(int i=0;i<Nloop;i++){
       double start=usecond();
+      for(int i=0;i<Nloop;i++){
     
 	ncomm=0;
 	for(int mu=0;mu<4;mu++){
@@ -219,34 +178,27 @@ int main (int argc, char ** argv)
 	  }
 	}
 	Grid.Barrier();
-	double stop=usecond();
-	t_time[i] = stop-start; // microseconds
-
       }
 
-      timestat.statistics(t_time);
+      double stop=usecond();
       
-      double dbytes    = bytes*ppn;
+      double dbytes    = bytes;
-      double xbytes    = dbytes*2.0*ncomm;
+      double xbytes    = Nloop*dbytes*2.0*ncomm;
       double rbytes    = xbytes;
       double bidibytes = xbytes+rbytes;
 
-    std::cout<<GridLogMessage << std::setw(4) << lat<<"\t"<<Ls<<"\t"
+      double time = stop-start;
-               <<std::setw(11) << bytes<< std::fixed << std::setprecision(1) << std::setw(7)
-               <<std::right<< xbytes/timestat.mean<<"  "<< xbytes*timestat.err/(timestat.mean*timestat.mean)<< " "
-               <<xbytes/timestat.max <<" "<< xbytes/timestat.min  
-               << "\t\t"<<std::setw(7)<< bidibytes/timestat.mean<< "  " << bidibytes*timestat.err/(timestat.mean*timestat.mean) << " "
-               << bidibytes/timestat.max << " " << bidibytes/timestat.min << std::endl;
 
-      
+      std::cout<<GridLogMessage << lat<<"\t\t"<<Ls<<"\t\t"<<bytes<<"\t\t"<<xbytes/time<<"\t\t"<<bidibytes/time<<std::endl;
     }
   }  
 
 
+  Nloop=10;
   std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
   std::cout<<GridLogMessage << "= Benchmarking concurrent STENCIL halo exchange in "<<nmu<<" dimensions"<<std::endl;
   std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
-  header();
+  std::cout<<GridLogMessage << "  L  "<<"\t\t"<<" Ls  "<<"\t\t"<<"bytes"<<"\t\t"<<"MB/s uni"<<"\t\t"<<"MB/s bidi"<<std::endl;
 
   for(int lat=4;lat<=maxlat;lat+=4){
     for(int Ls=8;Ls<=32;Ls*=2){
@@ -257,9 +209,6 @@ int main (int argc, char ** argv)
       				    lat*mpi_layout[3]});
 
       GridCartesian     Grid(latt_size,simd_layout,mpi_layout);
-      RealD Nrank = Grid._Nprocessors;
-      RealD Nnode = Grid.NodeCount();
-      RealD ppn = Nrank/Nnode;
 
       std::vector<HalfSpinColourVectorD *> xbuf(8);
       std::vector<HalfSpinColourVectorD *> rbuf(8);
@@ -267,115 +216,16 @@ int main (int argc, char ** argv)
       for(int d=0;d<8;d++){
 	xbuf[d] = (HalfSpinColourVectorD *)Grid.ShmBufferMalloc(lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD));
 	rbuf[d] = (HalfSpinColourVectorD *)Grid.ShmBufferMalloc(lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD));
-	bzero((void *)xbuf[d],lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD));
-	bzero((void *)rbuf[d],lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD));
       }
 
       int ncomm;
       int bytes=lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD);
 
-      double dbytes;
+      double start=usecond();
       for(int i=0;i<Nloop;i++){
-	double start=usecond();
-
-	dbytes=0;
-	ncomm=0;
 
 	std::vector<CartesianCommunicator::CommsRequest_t> requests;
 
-	for(int mu=0;mu<4;mu++){
-	
-
-	  if (mpi_layout[mu]>1 ) {
-	  
-	    ncomm++;
-	    int comm_proc=1;
-	    int xmit_to_rank;
-	    int recv_from_rank;
-	    Grid.ShiftedRanks(mu,comm_proc,xmit_to_rank,recv_from_rank);
-	    dbytes+=
-	      Grid.StencilSendToRecvFromBegin(requests,
-					      (void *)&xbuf[mu][0],
-					      xmit_to_rank,
-					      (void *)&rbuf[mu][0],
-					      recv_from_rank,
-					      bytes);
-	
-	    comm_proc = mpi_layout[mu]-1;
-	  
-	    Grid.ShiftedRanks(mu,comm_proc,xmit_to_rank,recv_from_rank);
-	    dbytes+=
-	      Grid.StencilSendToRecvFromBegin(requests,
-					      (void *)&xbuf[mu+4][0],
-					      xmit_to_rank,
-					      (void *)&rbuf[mu+4][0],
-					      recv_from_rank,
-					      bytes);
-	  
-	  }
-	}
-	Grid.StencilSendToRecvFromComplete(requests);
-	Grid.Barrier();
-	double stop=usecond();
-	t_time[i] = stop-start; // microseconds
-	
-      }
-
-      timestat.statistics(t_time);
-
-      dbytes=dbytes*ppn;
-      double xbytes    = dbytes*0.5;
-      double rbytes    = dbytes*0.5;
-      double bidibytes = dbytes;
-
-      std::cout<<GridLogMessage << std::setw(4) << lat<<"\t"<<Ls<<"\t"
-               <<std::setw(11) << bytes<< std::fixed << std::setprecision(1) << std::setw(7)
-               <<std::right<< xbytes/timestat.mean<<"  "<< xbytes*timestat.err/(timestat.mean*timestat.mean)<< " "
-               <<xbytes/timestat.max <<" "<< xbytes/timestat.min  
-               << "\t\t"<<std::setw(7)<< bidibytes/timestat.mean<< "  " << bidibytes*timestat.err/(timestat.mean*timestat.mean) << " "
-               << bidibytes/timestat.max << " " << bidibytes/timestat.min << std::endl;
-
-
-    }
-  }    
-
-
-  std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
-  std::cout<<GridLogMessage << "= Benchmarking sequential STENCIL halo exchange in "<<nmu<<" dimensions"<<std::endl;
-  std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
-  header();
-
-  for(int lat=4;lat<=maxlat;lat+=4){
-    for(int Ls=8;Ls<=32;Ls*=2){
-
-      std::vector<int> latt_size  ({lat*mpi_layout[0],
-      				    lat*mpi_layout[1],
-      				    lat*mpi_layout[2],
-      				    lat*mpi_layout[3]});
-
-      GridCartesian     Grid(latt_size,simd_layout,mpi_layout);
-      RealD Nrank = Grid._Nprocessors;
-      RealD Nnode = Grid.NodeCount();
-      RealD ppn = Nrank/Nnode;
-
-      std::vector<HalfSpinColourVectorD *> xbuf(8);
-      std::vector<HalfSpinColourVectorD *> rbuf(8);
-      Grid.ShmBufferFreeAll();
-      for(int d=0;d<8;d++){
-	xbuf[d] = (HalfSpinColourVectorD *)Grid.ShmBufferMalloc(lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD));
-	rbuf[d] = (HalfSpinColourVectorD *)Grid.ShmBufferMalloc(lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD));
-	bzero((void *)xbuf[d],lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD));
-	bzero((void *)rbuf[d],lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD));
-      }
-
-      int ncomm;
-      int bytes=lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD);
-      double dbytes;
-      for(int i=0;i<Nloop;i++){
-	double start=usecond();
-
-	std::vector<CartesianCommunicator::CommsRequest_t> requests;
-	dbytes=0;
 	ncomm=0;
 	for(int mu=0;mu<4;mu++){
 	
@@ -387,52 +237,123 @@ int main (int argc, char ** argv)
 	    int recv_from_rank;
 	    
 	    Grid.ShiftedRanks(mu,comm_proc,xmit_to_rank,recv_from_rank);
-	    dbytes+=
+	    Grid.StencilSendToRecvFromBegin(requests,
-	      Grid.StencilSendToRecvFromBegin(requests,
+					    (void *)&xbuf[mu][0],
-					      (void *)&xbuf[mu][0],
+					    xmit_to_rank,
-					      xmit_to_rank,
+					    (void *)&rbuf[mu][0],
-					      (void *)&rbuf[mu][0],
+					    recv_from_rank,
-					      recv_from_rank,
+					    bytes);
-					      bytes);
+	
+	    comm_proc = mpi_layout[mu]-1;
+	  
+	    Grid.ShiftedRanks(mu,comm_proc,xmit_to_rank,recv_from_rank);
+	    Grid.StencilSendToRecvFromBegin(requests,
+					    (void *)&xbuf[mu+4][0],
+					    xmit_to_rank,
+					    (void *)&rbuf[mu+4][0],
+					    recv_from_rank,
+					    bytes);
+	  
+	  }
+	}
+	Grid.StencilSendToRecvFromComplete(requests);
+	Grid.Barrier();
+
+      }
+      double stop=usecond();
+
+      double dbytes    = bytes;
+      double xbytes    = Nloop*dbytes*2.0*ncomm;
+      double rbytes    = xbytes;
+      double bidibytes = xbytes+rbytes;
+
+      double time = stop-start; // microseconds
+
+      std::cout<<GridLogMessage << lat<<"\t\t"<<Ls<<"\t\t"<<bytes<<"\t\t"<<xbytes/time<<"\t\t"<<bidibytes/time<<std::endl;
+    }
+  }    
+
+
+
+  Nloop=100;
+  std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
+  std::cout<<GridLogMessage << "= Benchmarking sequential STENCIL halo exchange in "<<nmu<<" dimensions"<<std::endl;
+  std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
+  std::cout<<GridLogMessage << "  L  "<<"\t\t"<<" Ls  "<<"\t\t"<<"bytes"<<"\t\t"<<"MB/s uni"<<"\t\t"<<"MB/s bidi"<<std::endl;
+
+  for(int lat=4;lat<=maxlat;lat+=4){
+    for(int Ls=8;Ls<=32;Ls*=2){
+
+      std::vector<int> latt_size  ({lat*mpi_layout[0],
+      				    lat*mpi_layout[1],
+      				    lat*mpi_layout[2],
+      				    lat*mpi_layout[3]});
+
+      GridCartesian     Grid(latt_size,simd_layout,mpi_layout);
+
+      std::vector<HalfSpinColourVectorD *> xbuf(8);
+      std::vector<HalfSpinColourVectorD *> rbuf(8);
+      Grid.ShmBufferFreeAll();
+      for(int d=0;d<8;d++){
+	xbuf[d] = (HalfSpinColourVectorD *)Grid.ShmBufferMalloc(lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD));
+	rbuf[d] = (HalfSpinColourVectorD *)Grid.ShmBufferMalloc(lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD));
+      }
+
+      int ncomm;
+      int bytes=lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD);
+
+      double start=usecond();
+      for(int i=0;i<Nloop;i++){
+
+	std::vector<CartesianCommunicator::CommsRequest_t> requests;
+
+	ncomm=0;
+	for(int mu=0;mu<4;mu++){
+	
+	  if (mpi_layout[mu]>1 ) {
+	  
+	    ncomm++;
+	    int comm_proc=1;
+	    int xmit_to_rank;
+	    int recv_from_rank;
+	    
+	    Grid.ShiftedRanks(mu,comm_proc,xmit_to_rank,recv_from_rank);
+	    Grid.StencilSendToRecvFromBegin(requests,
+					    (void *)&xbuf[mu][0],
+					    xmit_to_rank,
+					    (void *)&rbuf[mu][0],
+					    recv_from_rank,
+					    bytes);
 	    Grid.StencilSendToRecvFromComplete(requests);
 	    requests.resize(0);
 
 	    comm_proc = mpi_layout[mu]-1;
 	  
 	    Grid.ShiftedRanks(mu,comm_proc,xmit_to_rank,recv_from_rank);
-	    dbytes+=
+	    Grid.StencilSendToRecvFromBegin(requests,
-	      Grid.StencilSendToRecvFromBegin(requests,
+					    (void *)&xbuf[mu+4][0],
-					      (void *)&xbuf[mu+4][0],
+					    xmit_to_rank,
-					      xmit_to_rank,
+					    (void *)&rbuf[mu+4][0],
-					      (void *)&rbuf[mu+4][0],
+					    recv_from_rank,
-					      recv_from_rank,
+					    bytes);
-					      bytes);
 	    Grid.StencilSendToRecvFromComplete(requests);
 	    requests.resize(0);
 	  
 	  }
 	}
 	Grid.Barrier();
-	double stop=usecond();
+
-	t_time[i] = stop-start; // microseconds
-	
       }
+      double stop=usecond();
 
-      timestat.statistics(t_time);
+      double dbytes    = bytes;
+      double xbytes    = Nloop*dbytes*2.0*ncomm;
+      double rbytes    = xbytes;
+      double bidibytes = xbytes+rbytes;
 
-      dbytes=dbytes*ppn;
+      double time = stop-start; // microseconds
-      double xbytes    = dbytes*0.5;
-      double rbytes    = dbytes*0.5;
-      double bidibytes = dbytes;
 
-
+      std::cout<<GridLogMessage << lat<<"\t\t"<<Ls<<"\t\t"<<bytes<<"\t\t"<<xbytes/time<<"\t\t"<<bidibytes/time<<std::endl;
-      std::cout<<GridLogMessage << std::setw(4) << lat<<"\t"<<Ls<<"\t"
-               <<std::setw(11) << bytes<< std::fixed << std::setprecision(1) << std::setw(7)
-               <<std::right<< xbytes/timestat.mean<<"  "<< xbytes*timestat.err/(timestat.mean*timestat.mean)<< " "
-               <<xbytes/timestat.max <<" "<< xbytes/timestat.min  
-               << "\t\t"<<std::setw(7)<< bidibytes/timestat.mean<< "  " << bidibytes*timestat.err/(timestat.mean*timestat.mean) << " "
-               << bidibytes/timestat.max << " " << bidibytes/timestat.min << std::endl;
- 
     }
   }    
 

benchmarks/Benchmark_memory_bandwidth.cc

@@ -55,8 +55,8 @@ int main (int argc, char ** argv)
   std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
   std::cout<<GridLogMessage << "  L  "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s"<<"\t\t"<<"Gflop/s"<<"\t\t seconds"<<std::endl;
   std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
-  uint64_t lmax=64;
+  uint64_t lmax=44;
-#define NLOOP (100*lmax*lmax*lmax*lmax/vol)
+#define NLOOP (1*lmax*lmax*lmax*lmax/vol)
   for(int lat=4;lat<=lmax;lat+=4){
 
       std::vector<int> latt_size  ({lat*mpi_layout[0],lat*mpi_layout[1],lat*mpi_layout[2],lat*mpi_layout[3]});

benchmarks/Benchmark_su3.cc

@@ -35,9 +35,9 @@ using namespace Grid::QCD;
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
-#define LMAX (64)
+#define LMAX (32)
 
-  int Nloop=20;
+  int Nloop=200;
 
   std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
   std::vector<int> mpi_layout  = GridDefaultMpi();

configure.ac

@@ -27,7 +27,7 @@ AX_GXX_VERSION
 AC_DEFINE_UNQUOTED([GXX_VERSION],["$GXX_VERSION"],
       [version of g++ that will compile the code])
 
-CXXFLAGS="-g $CXXFLAGS"
+CXXFLAGS="-O3 $CXXFLAGS"
 
 
 ############### Checks for typedefs, structures, and compiler characteristics
@@ -184,10 +184,6 @@ AC_SEARCH_LIBS([limeCreateReader], [lime],
 In order to use ILGG file format please install or provide the correct path to your installation
 Info at: http://usqcd.jlab.org/usqcd-docs/c-lime/)])
 
-AC_SEARCH_LIBS([crc32], [z],
-               [AC_DEFINE([HAVE_ZLIB], [1], [Define to 1 if you have the `LIBZ' library])]
-               [have_zlib=true],
-	       [AC_MSG_ERROR(zlib library was not found in your system.)])
 
 AC_SEARCH_LIBS([H5Fopen], [hdf5_cpp],
                [AC_DEFINE([HAVE_HDF5], [1], [Define to 1 if you have the `HDF5' library])]

extras/Hadrons/Modules/MGauge/Load.cc

@@ -65,7 +65,7 @@ void TLoad::setup(void)
 // execution ///////////////////////////////////////////////////////////////////
 void TLoad::execute(void)
 {
-    FieldMetaData  header;
+    NerscField  header;
     std::string fileName = par().file + "."
                            + std::to_string(env().getTrajectory());
     
@@ -74,5 +74,5 @@ void TLoad::execute(void)
     LatticeGaugeField &U = *env().createLattice<LatticeGaugeField>(getName());
     NerscIO::readConfiguration(U, header, fileName);
     LOG(Message) << "NERSC header:" << std::endl;
-    dump_meta_data(header, LOG(Message));
+    dump_nersc_header(header, LOG(Message));
 }

lib/Grid.h

@@ -42,7 +42,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #include <Grid/GridQCDcore.h>
 #include <Grid/qcd/action/Action.h>
 #include <Grid/qcd/smearing/Smearing.h>
-#include <Grid/parallelIO/MetaData.h>
 #include <Grid/qcd/hmc/HMC_aggregate.h>
 
 #endif

lib/GridStd.h

@@ -7,7 +7,6 @@
 #include <cassert>
 #include <complex>
 #include <vector>
-#include <string>
 #include <iostream>
 #include <iomanip>
 #include <random>
@@ -19,7 +18,6 @@
 #include <ctime>
 #include <sys/time.h>
 #include <chrono>
-#include <zlib.h>
 
 ///////////////////
 // Grid config

lib/cartesian/Cartesian_base.h

@@ -50,6 +50,7 @@ public:
 
     GridBase(const std::vector<int> & processor_grid) : CartesianCommunicator(processor_grid) {};
 
+
     // Physics Grid information.
     std::vector<int> _simd_layout;// Which dimensions get relayed out over simd lanes.
     std::vector<int> _fdimensions;// (full) Global dimensions of array prior to cb removal
@@ -62,12 +63,13 @@ public:
     int _isites;
     int _fsites;                  // _isites*_osites = product(dimensions).
     int _gsites;
-    std::vector<int> _slice_block;// subslice information
+    std::vector<int> _slice_block;   // subslice information
     std::vector<int> _slice_stride;
     std::vector<int> _slice_nblock;
 
-    std::vector<int> _lstart;     // local start of array in gcoors _processor_coor[d]*_ldimensions[d]
+    // Might need these at some point
-    std::vector<int> _lend  ;     // local end of array in gcoors   _processor_coor[d]*_ldimensions[d]+_ldimensions_[d]-1
+    //    std::vector<int> _lstart;     // local start of array in gcoors. _processor_coor[d]*_ldimensions[d]
+    //    std::vector<int> _lend;       // local end of array in gcoors    _processor_coor[d]*_ldimensions[d]+_ldimensions_[d]-1
 
 public:
 
@@ -174,7 +176,6 @@ public:
     inline int gSites(void) const { return _isites*_osites*_Nprocessors; }; 
     inline int Nd    (void) const { return _ndimension;};
 
-    inline const std::vector<int> LocalStarts(void)             { return _lstart;    };
     inline const std::vector<int> &FullDimensions(void)         { return _fdimensions;};
     inline const std::vector<int> &GlobalDimensions(void)       { return _gdimensions;};
     inline const std::vector<int> &LocalDimensions(void)        { return _ldimensions;};

lib/cartesian/Cartesian_full.h

@@ -76,8 +76,6 @@ public:
         _ldimensions.resize(_ndimension);
         _rdimensions.resize(_ndimension);
         _simd_layout.resize(_ndimension);
-	_lstart.resize(_ndimension);
-	_lend.resize(_ndimension);
             
         _ostride.resize(_ndimension);
         _istride.resize(_ndimension);
@@ -96,10 +94,8 @@ public:
 	  // Use a reduced simd grid
 	  _ldimensions[d]= _gdimensions[d]/_processors[d];  //local dimensions
 	  _rdimensions[d]= _ldimensions[d]/_simd_layout[d]; //overdecomposition
-	  _lstart[d]     = _processor_coor[d]*_ldimensions[d];
+	  _osites *= _rdimensions[d];
-	  _lend[d]       = _processor_coor[d]*_ldimensions[d]+_ldimensions[d]-1;
+	  _isites *= _simd_layout[d];
-	  _osites  *= _rdimensions[d];
-	  _isites  *= _simd_layout[d];
                 
 	  // Addressing support
 	  if ( d==0 ) {

lib/cartesian/Cartesian_red_black.h

@@ -151,8 +151,6 @@ public:
       _ldimensions.resize(_ndimension);
       _rdimensions.resize(_ndimension);
       _simd_layout.resize(_ndimension);
-      _lstart.resize(_ndimension);
-      _lend.resize(_ndimension);
       
       _ostride.resize(_ndimension);
       _istride.resize(_ndimension);
@@ -171,8 +169,6 @@ public:
 	  _gdimensions[d] = _gdimensions[d]/2; // Remove a checkerboard
 	}
 	_ldimensions[d] = _gdimensions[d]/_processors[d];
-	_lstart[d]     = _processor_coor[d]*_ldimensions[d];
-	_lend[d]       = _processor_coor[d]*_ldimensions[d]+_ldimensions[d]-1;
 
 	// Use a reduced simd grid
 	_simd_layout[d] = simd_layout[d];

lib/communicator/Communicator_base.cc

@@ -60,7 +60,6 @@ void CartesianCommunicator::ShmBufferFreeAll(void) {
 /////////////////////////////////
 // Grid information queries
 /////////////////////////////////
-int                      CartesianCommunicator::Dimensions(void)         { return _ndimension; };
 int                      CartesianCommunicator::IsBoss(void)            { return _processor==0; };
 int                      CartesianCommunicator::BossRank(void)          { return 0; };
 int                      CartesianCommunicator::ThisRank(void)          { return _processor; };
@@ -92,7 +91,6 @@ void CartesianCommunicator::GlobalSumVector(ComplexD *c,int N)
 #if !defined( GRID_COMMS_MPI3) && !defined (GRID_COMMS_MPI3L)
 
 int                      CartesianCommunicator::NodeCount(void)    { return ProcessorCount();};
-int                      CartesianCommunicator::RankCount(void)    { return ProcessorCount();};
 
 double CartesianCommunicator::StencilSendToRecvFromBegin(std::vector<CommsRequest_t> &list,
 						       void *xmit,

lib/communicator/Communicator_base.h

@@ -148,7 +148,6 @@ class CartesianCommunicator {
   int  RankFromProcessorCoor(std::vector<int> &coor);
   void ProcessorCoorFromRank(int rank,std::vector<int> &coor);
   
-  int                      Dimensions(void)        ;
   int                      IsBoss(void)            ;
   int                      BossRank(void)          ;
   int                      ThisRank(void)          ;
@@ -156,7 +155,6 @@ class CartesianCommunicator {
   const std::vector<int> & ProcessorGrid(void)     ;
   int                      ProcessorCount(void)    ;
   int                      NodeCount(void)    ;
-  int                      RankCount(void)    ;
 
   ////////////////////////////////////////////////////////////////////////////////
   // very VERY rarely (Log, serial RNG) we need world without a grid
@@ -177,8 +175,6 @@ class CartesianCommunicator {
   void GlobalSumVector(ComplexF *c,int N);
   void GlobalSum(ComplexD &c);
   void GlobalSumVector(ComplexD *c,int N);
-  void GlobalXOR(uint32_t &);
-  void GlobalXOR(uint64_t &);
   
   template<class obj> void GlobalSum(obj &o){
     typedef typename obj::scalar_type scalar_type;

lib/communicator/Communicator_mpi.cc

@@ -83,14 +83,6 @@ void CartesianCommunicator::GlobalSum(uint64_t &u){
   int ierr=MPI_Allreduce(MPI_IN_PLACE,&u,1,MPI_UINT64_T,MPI_SUM,communicator);
   assert(ierr==0);
 }
-void CartesianCommunicator::GlobalXOR(uint32_t &u){
-  int ierr=MPI_Allreduce(MPI_IN_PLACE,&u,1,MPI_UINT32_T,MPI_BXOR,communicator);
-  assert(ierr==0);
-}
-void CartesianCommunicator::GlobalXOR(uint64_t &u){
-  int ierr=MPI_Allreduce(MPI_IN_PLACE,&u,1,MPI_UINT64_T,MPI_BXOR,communicator);
-  assert(ierr==0);
-}
 void CartesianCommunicator::GlobalSum(float &f){
   int ierr=MPI_Allreduce(MPI_IN_PLACE,&f,1,MPI_FLOAT,MPI_SUM,communicator);
   assert(ierr==0);

lib/communicator/Communicator_mpi3.cc

@@ -65,7 +65,6 @@ std::vector<int> CartesianCommunicator::MyGroup;
 std::vector<void *> CartesianCommunicator::ShmCommBufs;
 
 int CartesianCommunicator::NodeCount(void)    { return GroupSize;};
-int CartesianCommunicator::RankCount(void)    { return WorldSize;};
 
 
 #undef FORCE_COMMS
@@ -510,14 +509,6 @@ void CartesianCommunicator::GlobalSum(uint64_t &u){
   int ierr=MPI_Allreduce(MPI_IN_PLACE,&u,1,MPI_UINT64_T,MPI_SUM,communicator);
   assert(ierr==0);
 }
-void CartesianCommunicator::GlobalXOR(uint32_t &u){
-  int ierr=MPI_Allreduce(MPI_IN_PLACE,&u,1,MPI_UINT32_T,MPI_BXOR,communicator);
-  assert(ierr==0);
-}
-void CartesianCommunicator::GlobalXOR(uint64_t &u){
-  int ierr=MPI_Allreduce(MPI_IN_PLACE,&u,1,MPI_UINT64_T,MPI_BXOR,communicator);
-  assert(ierr==0);
-}
 void CartesianCommunicator::GlobalSum(float &f){
   int ierr=MPI_Allreduce(MPI_IN_PLACE,&f,1,MPI_FLOAT,MPI_SUM,communicator);
   assert(ierr==0);

lib/communicator/Communicator_none.cc

@@ -59,8 +59,6 @@ void CartesianCommunicator::GlobalSum(double &){}
 void CartesianCommunicator::GlobalSum(uint32_t &){}
 void CartesianCommunicator::GlobalSum(uint64_t &){}
 void CartesianCommunicator::GlobalSumVector(double *,int N){}
-void CartesianCommunicator::GlobalXOR(uint32_t &){}
-void CartesianCommunicator::GlobalXOR(uint64_t &){}
 
 void CartesianCommunicator::SendRecvPacket(void *xmit,
 					   void *recv,

lib/lattice/Lattice_reduction.h

@@ -328,8 +328,6 @@ static void sliceMaddVector(Lattice<vobj> &R,std::vector<RealD> &a,const Lattice
   typedef typename vobj::vector_type vector_type;
   typedef typename vobj::tensor_reduced tensor_reduced;
   
-  scalar_type zscale(scale);
-
   GridBase *grid  = X._grid;
 
   int Nsimd  =grid->Nsimd();
@@ -355,7 +353,7 @@ static void sliceMaddVector(Lattice<vobj> &R,std::vector<RealD> &a,const Lattice
       grid->iCoorFromIindex(icoor,l);
       int ldx =r+icoor[orthogdim]*rd;
       scalar_type *as =(scalar_type *)&av;
-      as[l] = scalar_type(a[ldx])*zscale;
+      as[l] = scalar_type(a[ldx])*scale;
     }
 
     tensor_reduced at; at=av;

lib/lattice/Lattice_transfer.h

@@ -551,10 +551,7 @@ void Replicate(Lattice<vobj> &coarse,Lattice<vobj> & fine)
 
 //Copy SIMD-vectorized lattice to array of scalar objects in lexicographic order
 template<typename vobj, typename sobj>
-typename std::enable_if<isSIMDvectorized<vobj>::value && !isSIMDvectorized<sobj>::value, void>::type 
+typename std::enable_if<isSIMDvectorized<vobj>::value && !isSIMDvectorized<sobj>::value, void>::type unvectorizeToLexOrdArray(std::vector<sobj> &out, const Lattice<vobj> &in){
-unvectorizeToLexOrdArray(std::vector<sobj> &out, const Lattice<vobj> &in)
-{
-
   typedef typename vobj::vector_type vtype;
   
   GridBase* in_grid = in._grid;
@@ -593,54 +590,6 @@ unvectorizeToLexOrdArray(std::vector<sobj> &out, const Lattice<vobj> &in)
     extract1(in_vobj, out_ptrs, 0);
   }
 }
-//Copy SIMD-vectorized lattice to array of scalar objects in lexicographic order
-template<typename vobj, typename sobj>
-typename std::enable_if<isSIMDvectorized<vobj>::value 
-                    && !isSIMDvectorized<sobj>::value, void>::type 
-vectorizeFromLexOrdArray( std::vector<sobj> &in, Lattice<vobj> &out)
-{
-
-  typedef typename vobj::vector_type vtype;
-  
-  GridBase* grid = out._grid;
-  assert(in.size()==grid->lSites());
-  
-  int ndim     = grid->Nd();
-  int nsimd    = vtype::Nsimd();
-
-  std::vector<std::vector<int> > icoor(nsimd);
-      
-  for(int lane=0; lane < nsimd; lane++){
-    icoor[lane].resize(ndim);
-    grid->iCoorFromIindex(icoor[lane],lane);
-  }
-  
-  parallel_for(uint64_t oidx = 0; oidx < grid->oSites(); oidx++){ //loop over outer index
-    //Assemble vector of pointers to output elements
-    std::vector<sobj*> ptrs(nsimd);
-
-    std::vector<int> ocoor(ndim);
-    grid->oCoorFromOindex(ocoor, oidx);
-
-    std::vector<int> lcoor(grid->Nd());
-      
-    for(int lane=0; lane < nsimd; lane++){
-
-      for(int mu=0;mu<ndim;mu++){
-	lcoor[mu] = ocoor[mu] + grid->_rdimensions[mu]*icoor[lane][mu];
-      }
-
-      int lex;
-      Lexicographic::IndexFromCoor(lcoor, lex, grid->_ldimensions);
-      ptrs[lane] = &in[lex];
-    }
-    
-    //pack from those ptrs
-    vobj vecobj;
-    merge1(vecobj, ptrs, 0);
-    out._odata[oidx] = vecobj; 
-  }
-}
 
 //Convert a Lattice from one precision to another
 template<class VobjOut, class VobjIn>
@@ -666,7 +615,7 @@ void precisionChange(Lattice<VobjOut> &out, const Lattice<VobjIn> &in){
   std::vector<SobjOut> in_slex_conv(in_grid->lSites());
   unvectorizeToLexOrdArray(in_slex_conv, in);
     
-  parallel_for(uint64_t out_oidx=0;out_oidx<out_grid->oSites();out_oidx++){
+  parallel_for(int out_oidx=0;out_oidx<out_grid->oSites();out_oidx++){
     std::vector<int> out_ocoor(ndim);
     out_grid->oCoorFromOindex(out_ocoor, out_oidx);
 

lib/lattice/Lattice_unary.h

@@ -66,21 +66,10 @@ namespace Grid {
     Lattice<obj> ret(rhs._grid);
     ret.checkerboard = rhs.checkerboard;
     conformable(ret,rhs);
-    obj unit(1.0);
     parallel_for(int ss=0;ss<rhs._grid->oSites();ss++){
-      //ret._odata[ss]=Exponentiate(rhs._odata[ss],alpha, Nexp);
+      ret._odata[ss]=Exponentiate(rhs._odata[ss],alpha, Nexp);
-      ret._odata[ss] = unit;
-      for(int i=Nexp; i>=1;--i)
-	ret._odata[ss] = unit + ret._odata[ss]*rhs._odata[ss]*(alpha/RealD(i));
-
     }
-
     return ret;
-
-    
-    
-
-    
   }
 
 

lib/parallelIO/IldgIO.h

@@ -27,7 +27,6 @@ directory
 #ifndef GRID_ILDG_IO_H
 #define GRID_ILDG_IO_H
 
-#ifdef HAVE_LIME
 #include <algorithm>
 #include <fstream>
 #include <iomanip>
@@ -38,672 +37,213 @@ directory
 #include <sys/utsname.h>
 #include <unistd.h>
 
-//C-Lime is a must have for this functionality
+#ifdef HAVE_LIME
-extern "C" {  
+
+extern "C" {  // for linkage
 #include "lime.h"
 }
 
 namespace Grid {
 namespace QCD {
 
-  /////////////////////////////////
+inline void ILDGGrid(GridBase *grid, ILDGField &header) {
-  // Encode word types as strings
+  assert(grid->_ndimension == 4);  // emit error if not
-  /////////////////////////////////
+  header.dimension.resize(4);
- template<class word> inline std::string ScidacWordMnemonic(void){ return std::string("unknown"); }
+  header.boundary.resize(4);
- template<> inline std::string ScidacWordMnemonic<double>  (void){ return std::string("D"); }
+  for (int d = 0; d < 4; d++) {
- template<> inline std::string ScidacWordMnemonic<float>   (void){ return std::string("F"); }
+    header.dimension[d] = grid->_fdimensions[d];
- template<> inline std::string ScidacWordMnemonic< int32_t>(void){ return std::string("I32_t"); }
+    // Read boundary conditions from ... ?
- template<> inline std::string ScidacWordMnemonic<uint32_t>(void){ return std::string("U32_t"); }
+    header.boundary[d] = std::string("periodic");
- template<> inline std::string ScidacWordMnemonic< int64_t>(void){ return std::string("I64_t"); }
+  }
- template<> inline std::string ScidacWordMnemonic<uint64_t>(void){ return std::string("U64_t"); }
+}
 
-  /////////////////////////////////////////
+inline void ILDGChecksum(uint32_t *buf, uint32_t buf_size_bytes,
-  // Encode a generic tensor as a string
+                         uint32_t &csum) {
-  /////////////////////////////////////////
+  BinaryIO::Uint32Checksum(buf, buf_size_bytes, csum);
- template<class vobj> std::string ScidacRecordTypeString(int &colors, int &spins, int & typesize,int &datacount) { 
+}
 
-   typedef typename getPrecision<vobj>::real_scalar_type stype;
+//////////////////////////////////////////////////////////////////////
+// Utilities ; these are QCD aware
+//////////////////////////////////////////////////////////////////////
+template <class GaugeField>
+inline void ILDGStatistics(GaugeField &data, ILDGField &header) {
+  // How to convert data precision etc...
+  header.link_trace = Grid::QCD::WilsonLoops<PeriodicGimplR>::linkTrace(data);
+  header.plaquette = Grid::QCD::WilsonLoops<PeriodicGimplR>::avgPlaquette(data);
+  // header.polyakov =
+}
 
-   int _ColourN       = indexRank<ColourIndex,vobj>();
+// Forcing QCD here
-   int _ColourScalar  =  isScalar<ColourIndex,vobj>();
+template <class fobj, class sobj>
-   int _ColourVector  =  isVector<ColourIndex,vobj>();
+struct ILDGMunger {
-   int _ColourMatrix  =  isMatrix<ColourIndex,vobj>();
+  void operator()(fobj &in, sobj &out, uint32_t &csum) {
-
+    for (int mu = 0; mu < 4; mu++) {
-   int _SpinN       = indexRank<SpinIndex,vobj>();
+      for (int i = 0; i < 3; i++) {
-   int _SpinScalar  =  isScalar<SpinIndex,vobj>();
+        for (int j = 0; j < 3; j++) {
-   int _SpinVector  =  isVector<SpinIndex,vobj>();
+          out(mu)()(i, j) = in(mu)()(i, j);
-   int _SpinMatrix  =  isMatrix<SpinIndex,vobj>();
+        }
-
-   int _LorentzN       = indexRank<LorentzIndex,vobj>();
-   int _LorentzScalar  =  isScalar<LorentzIndex,vobj>();
-   int _LorentzVector  =  isVector<LorentzIndex,vobj>();
-   int _LorentzMatrix  =  isMatrix<LorentzIndex,vobj>();
-
-   std::stringstream stream;
-
-   stream << "GRID_";
-   stream << ScidacWordMnemonic<stype>();
-
-   //   std::cout << " Lorentz N/S/V/M : " << _LorentzN<<" "<<_LorentzScalar<<"/"<<_LorentzVector<<"/"<<_LorentzMatrix<<std::endl;
-   //   std::cout << " Spin    N/S/V/M : " << _SpinN   <<" "<<_SpinScalar   <<"/"<<_SpinVector   <<"/"<<_SpinMatrix<<std::endl;
-   //   std::cout << " Colour  N/S/V/M : " << _ColourN <<" "<<_ColourScalar <<"/"<<_ColourVector <<"/"<<_ColourMatrix<<std::endl;
-
-   if ( _LorentzVector )   stream << "_LorentzVector"<<_LorentzN;
-   if ( _LorentzMatrix )   stream << "_LorentzMatrix"<<_LorentzN;
-
-   if ( _SpinVector )   stream << "_SpinVector"<<_SpinN;
-   if ( _SpinMatrix )   stream << "_SpinMatrix"<<_SpinN;
-
-   if ( _ColourVector )   stream << "_ColourVector"<<_ColourN;
-   if ( _ColourMatrix )   stream << "_ColourMatrix"<<_ColourN;
-
-   if ( _ColourScalar && _LorentzScalar && _SpinScalar )   stream << "_Complex";
-
-
-   typesize = sizeof(typename vobj::scalar_type);
-
-   if ( _ColourMatrix ) typesize*= _ColourN*_ColourN;
-   else                 typesize*= _ColourN;
-
-   if ( _SpinMatrix )   typesize*= _SpinN*_SpinN;
-   else                 typesize*= _SpinN;
-
-   colors    = _ColourN;
-   spins     = _SpinN;
-   datacount = _LorentzN;
-
-   return stream.str();
- }
- 
- template<class vobj> std::string ScidacRecordTypeString(Lattice<vobj> & lat,int &colors, int &spins, int & typesize,int &datacount) { 
-   return ScidacRecordTypeString<vobj>(colors,spins,typesize,datacount);
- };
-
-
- ////////////////////////////////////////////////////////////
- // Helper to fill out metadata
- ////////////////////////////////////////////////////////////
- template<class vobj> void ScidacMetaData(Lattice<vobj> & field,
-					  FieldMetaData &header,
-					  scidacRecord & _scidacRecord,
-					  scidacFile   & _scidacFile) 
- {
-   typedef typename getPrecision<vobj>::real_scalar_type stype;
-
-   /////////////////////////////////////
-   // Pull Grid's metadata
-   /////////////////////////////////////
-   PrepareMetaData(field,header);
-
-   /////////////////////////////////////
-   // Scidac Private File structure
-   /////////////////////////////////////
-   _scidacFile              = scidacFile(field._grid);
-
-   /////////////////////////////////////
-   // Scidac Private Record structure
-   /////////////////////////////////////
-   scidacRecord sr;
-   sr.datatype   = ScidacRecordTypeString(field,sr.colors,sr.spins,sr.typesize,sr.datacount);
-   sr.date       = header.creation_date;
-   sr.precision  = ScidacWordMnemonic<stype>();
-   sr.recordtype = GRID_IO_FIELD;
-
-   _scidacRecord = sr;
-
-   std::cout << GridLogMessage << "Build SciDAC datatype " <<sr.datatype<<std::endl;
- }
- 
- ///////////////////////////////////////////////////////
- // Scidac checksum
- ///////////////////////////////////////////////////////
- static int scidacChecksumVerify(scidacChecksum &scidacChecksum_,uint32_t scidac_csuma,uint32_t scidac_csumb)
- {
-   uint32_t scidac_checksuma = stoull(scidacChecksum_.suma,0,16);
-   uint32_t scidac_checksumb = stoull(scidacChecksum_.sumb,0,16);
-   if ( scidac_csuma !=scidac_checksuma) return 0;
-   if ( scidac_csumb !=scidac_checksumb) return 0;
-    return 1;
- }
-
-////////////////////////////////////////////////////////////////////////////////////
-// Lime, ILDG and Scidac I/O classes
-////////////////////////////////////////////////////////////////////////////////////
-class GridLimeReader : public BinaryIO {
- public:
-   ///////////////////////////////////////////////////
-   // FIXME: format for RNG? Now just binary out instead
-   ///////////////////////////////////////////////////
-
-   FILE       *File;
-   LimeReader *LimeR;
-   std::string filename;
-
-   /////////////////////////////////////////////
-   // Open the file
-   /////////////////////////////////////////////
-   void open(std::string &_filename) 
-   {
-     filename= _filename;
-     File = fopen(filename.c_str(), "r");
-     LimeR = limeCreateReader(File);
-   }
-   /////////////////////////////////////////////
-   // Close the file
-   /////////////////////////////////////////////
-   void close(void){
-     fclose(File);
-     //     limeDestroyReader(LimeR);
-   }
-
-  ////////////////////////////////////////////
-  // Read a generic lattice field and verify checksum
-  ////////////////////////////////////////////
-  template<class vobj>
-  void readLimeLatticeBinaryObject(Lattice<vobj> &field,std::string record_name)
-  {
-    typedef typename vobj::scalar_object sobj;
-    scidacChecksum scidacChecksum_;
-    uint32_t nersc_csum,scidac_csuma,scidac_csumb;
-
-    std::string format = getFormatString<vobj>();
-
-    while ( limeReaderNextRecord(LimeR) == LIME_SUCCESS ) { 
-
-      std::cout << GridLogMessage << limeReaderType(LimeR) <<std::endl;
-	
-      if ( strncmp(limeReaderType(LimeR), record_name.c_str(),strlen(record_name.c_str()) )  ) {
-
-
-	off_t offset= ftell(File);
-	BinarySimpleMunger<sobj,sobj> munge;
-	BinaryIO::readLatticeObject< sobj, sobj >(field, filename, munge, offset, format,nersc_csum,scidac_csuma,scidac_csumb);
-
-	/////////////////////////////////////////////
-	// Insist checksum is next record
-	/////////////////////////////////////////////
-	readLimeObject(scidacChecksum_,std::string("scidacChecksum"),record_name);
-
-	/////////////////////////////////////////////
-	// Verify checksums
-	/////////////////////////////////////////////
-	scidacChecksumVerify(scidacChecksum_,scidac_csuma,scidac_csumb);
-	return;
       }
     }
-  }
+    ILDGChecksum((uint32_t *)&in, sizeof(in), csum);
-  ////////////////////////////////////////////
+  };
-  // Read a generic serialisable object
-  ////////////////////////////////////////////
-  template<class serialisable_object>
-  void readLimeObject(serialisable_object &object,std::string object_name,std::string record_name)
-  {
-    std::string xmlstring;
-    // should this be a do while; can we miss a first record??
-    while ( limeReaderNextRecord(LimeR) == LIME_SUCCESS ) { 
-
-      uint64_t nbytes = limeReaderBytes(LimeR);//size of this record (configuration)
-
-      if ( strncmp(limeReaderType(LimeR), record_name.c_str(),strlen(record_name.c_str()) )  ) {
-	std::vector<char> xmlc(nbytes+1,'\0');
-	limeReaderReadData((void *)&xmlc[0], &nbytes, LimeR);    
-	XmlReader RD(&xmlc[0],"");
-	read(RD,object_name,object);
-	return;
-      }
-
-    }  
-    assert(0);
-  }
 };
 
-class GridLimeWriter : public BinaryIO {
+template <class fobj, class sobj>
+struct ILDGUnmunger {
+  void operator()(sobj &in, fobj &out, uint32_t &csum) {
+    for (int mu = 0; mu < 4; mu++) {
+      for (int i = 0; i < 3; i++) {
+        for (int j = 0; j < 3; j++) {
+          out(mu)()(i, j) = in(mu)()(i, j);
+        }
+      }
+    }
+    ILDGChecksum((uint32_t *)&out, sizeof(out), csum);
+  };
+};
+
+////////////////////////////////////////////////////////////////////////////////
+// Write and read from fstream; compute header offset for payload
+////////////////////////////////////////////////////////////////////////////////
+enum ILDGstate {ILDGread, ILDGwrite};
+
+class ILDGIO : public BinaryIO {
+  FILE *File;
+  LimeWriter *LimeW;
+  LimeRecordHeader *LimeHeader;
+  LimeReader *LimeR;
+  std::string filename;
+
+
  public:
-   ///////////////////////////////////////////////////
+  ILDGIO(std::string file, ILDGstate RW) {
-   // FIXME: format for RNG? Now just binary out instead
+      filename = file;
-   ///////////////////////////////////////////////////
+    if (RW == ILDGwrite){
+      File = fopen(file.c_str(), "w");
+      // check if opened correctly
 
-   FILE       *File;
+      LimeW = limeCreateWriter(File);
-   LimeWriter *LimeW;
+    } else {
-   std::string filename;
+      File = fopen(file.c_str(), "r");
+      // check if opened correctly
 
-   void open(std::string &_filename) { 
+      LimeR = limeCreateReader(File);
-     filename= _filename;
+    }
-     File = fopen(filename.c_str(), "w");
+  }
-     LimeW = limeCreateWriter(File); assert(LimeW != NULL );
+
-   }
+  ~ILDGIO() { fclose(File); }
-   /////////////////////////////////////////////
+
-   // Close the file
+  int createHeader(std::string message, int MB, int ME, size_t PayloadSize, LimeWriter* L){
-   /////////////////////////////////////////////
-   void close(void) {
-     fclose(File);
-     //  limeDestroyWriter(LimeW);
-   }
-  ///////////////////////////////////////////////////////
-  // Lime utility functions
-  ///////////////////////////////////////////////////////
-  int createLimeRecordHeader(std::string message, int MB, int ME, size_t PayloadSize)
-  {
     LimeRecordHeader *h;
     h = limeCreateHeader(MB, ME, const_cast<char *>(message.c_str()), PayloadSize);
-    assert(limeWriteRecordHeader(h, LimeW) >= 0);
+    int status = limeWriteRecordHeader(h, L);
+    if (status < 0) {
+      std::cerr << "ILDG Header error\n";
+      return status;
+    }
     limeDestroyHeader(h);
     return LIME_SUCCESS;
   }
-  ////////////////////////////////////////////
-  // Write a generic serialisable object
-  ////////////////////////////////////////////
-  template<class serialisable_object>
-  void writeLimeObject(int MB,int ME,serialisable_object &object,std::string object_name,std::string record_name)
-  {
-    std::string xmlstring;
-    {
-      XmlWriter WR("","");
-      write(WR,object_name,object);
-      xmlstring = WR.XmlString();
-    }
-    uint64_t nbytes = xmlstring.size();
-    int err;
-    LimeRecordHeader *h = limeCreateHeader(MB, ME,(char *)record_name.c_str(), nbytes); assert(h!= NULL);
 
-    err=limeWriteRecordHeader(h, LimeW);                    assert(err>=0);
+  unsigned int writeHeader(ILDGField &header) {
-    err=limeWriteRecordData(&xmlstring[0], &nbytes, LimeW); assert(err>=0);
+    // write header in LIME
-    err=limeWriterCloseRecord(LimeW);                       assert(err>=0);
+    n_uint64_t nbytes;
-    limeDestroyHeader(h);
+    int MB_flag = 1, ME_flag = 0;
-  }
-  ////////////////////////////////////////////
-  // Write a generic lattice field and csum
-  ////////////////////////////////////////////
-  template<class vobj>
-  void writeLimeLatticeBinaryObject(Lattice<vobj> &field,std::string record_name)
-  {
-    ////////////////////////////////////////////
-    // Create record header
-    ////////////////////////////////////////////
-    typedef typename vobj::scalar_object sobj;
-    int err;
-    uint32_t nersc_csum,scidac_csuma,scidac_csumb;
-    uint64_t PayloadSize = sizeof(sobj) * field._grid->_gsites;
-    createLimeRecordHeader(record_name, 0, 0, PayloadSize);
 
-    ////////////////////////////////////////////////////////////////////
+    char message[] = "ildg-format";
-    // NB: FILE and iostream are jointly writing disjoint sequences in the
+    nbytes = strlen(message);
-    // the same file through different file handles (integer units).
+    LimeHeader = limeCreateHeader(MB_flag, ME_flag, message, nbytes);
-    // 
+    limeWriteRecordHeader(LimeHeader, LimeW);
-    // These are both buffered, so why I think this code is right is as follows.
+    limeDestroyHeader(LimeHeader);
-    //
+    // save the xml header here
-    // i)  write record header to FILE *File, telegraphing the size. 
+    // use the xml_writer to c++ streams in pugixml
-    // ii) ftell reads the offset from FILE *File .
+    // and convert to char message
-    // iii) iostream / MPI Open independently seek this offset. Write sequence direct to disk.
+    limeWriteRecordData(message, &nbytes, LimeW);
-    //      Closes iostream and flushes.
+    limeWriterCloseRecord(LimeW);
-    // iv) fseek on FILE * to end of this disjoint section.
-    //  v) Continue writing scidac record.
-    ////////////////////////////////////////////////////////////////////
-    off_t offset = ftell(File);
-    std::string format = getFormatString<vobj>();
-    BinarySimpleMunger<sobj,sobj> munge;
-    BinaryIO::writeLatticeObject<vobj,sobj>(field, filename, munge, offset, format,nersc_csum,scidac_csuma,scidac_csumb);
-    err=limeWriterCloseRecord(LimeW);  assert(err>=0);
-    ////////////////////////////////////////
-    // Write checksum element, propagaing forward from the BinaryIO
-    // Always pair a checksum with a binary object, and close message
-    ////////////////////////////////////////
-    scidacChecksum checksum;
-    std::stringstream streama; streama << std::hex << scidac_csuma;
-    std::stringstream streamb; streamb << std::hex << scidac_csumb;
-    checksum.suma= streama.str();
-    checksum.sumb= streamb.str();
-    std::cout << GridLogMessage<<" writing scidac checksums "<<std::hex<<scidac_csuma<<"/"<<scidac_csumb<<std::dec<<std::endl;
-    writeLimeObject(0,1,checksum,std::string("scidacChecksum"    ),std::string(SCIDAC_CHECKSUM));
-  }
-};
 
-class ScidacWriter : public GridLimeWriter {
+    return 0;
- public:
+  }
-
+
-   template<class SerialisableUserFile>
+  unsigned int readHeader(ILDGField &header) {
-   void writeScidacFileRecord(GridBase *grid,SerialisableUserFile &_userFile)
+    return 0;
-   {
-     scidacFile    _scidacFile(grid);
-     writeLimeObject(1,0,_scidacFile,_scidacFile.SerialisableClassName(),std::string(SCIDAC_PRIVATE_FILE_XML));
-     writeLimeObject(0,1,_userFile,_userFile.SerialisableClassName(),std::string(SCIDAC_FILE_XML));
-   }
-  ////////////////////////////////////////////////
-  // Write generic lattice field in scidac format
-  ////////////////////////////////////////////////
-   template <class vobj, class userRecord>
-  void writeScidacFieldRecord(Lattice<vobj> &field,userRecord _userRecord) 
-  {
-    typedef typename vobj::scalar_object sobj;
-    uint64_t nbytes;
-    GridBase * grid = field._grid;
-
-    ////////////////////////////////////////
-    // fill the Grid header
-    ////////////////////////////////////////
-    FieldMetaData header;
-    scidacRecord  _scidacRecord;
-    scidacFile    _scidacFile;
-
-    ScidacMetaData(field,header,_scidacRecord,_scidacFile);
-
-    //////////////////////////////////////////////
-    // Fill the Lime file record by record
-    //////////////////////////////////////////////
-    writeLimeObject(1,0,header ,std::string("FieldMetaData"),std::string(GRID_FORMAT)); // Open message 
-    writeLimeObject(0,0,_userRecord,_userRecord.SerialisableClassName(),std::string(SCIDAC_RECORD_XML));
-    writeLimeObject(0,0,_scidacRecord,_scidacRecord.SerialisableClassName(),std::string(SCIDAC_PRIVATE_RECORD_XML));
-    writeLimeLatticeBinaryObject(field,std::string(ILDG_BINARY_DATA));      // Closes message with checksum
-  }
-};
-
-class IldgWriter : public ScidacWriter {
- public:
-
-  ///////////////////////////////////
-  // A little helper
-  ///////////////////////////////////
-  void writeLimeIldgLFN(std::string &LFN)
-  {
-    uint64_t PayloadSize = LFN.size();
-    int err;
-    createLimeRecordHeader(ILDG_DATA_LFN, 0 , 0, PayloadSize);
-    err=limeWriteRecordData(const_cast<char*>(LFN.c_str()), &PayloadSize,LimeW); assert(err>=0);
-    err=limeWriterCloseRecord(LimeW); assert(err>=0);
   }
 
-  ////////////////////////////////////////////////////////////////
-  // Special ILDG operations ; gauge configs only.
-  // Don't require scidac records EXCEPT checksum
-  // Use Grid MetaData object if present.
-  ////////////////////////////////////////////////////////////////
   template <class vsimd>
-  void writeConfiguration(Lattice<iLorentzColourMatrix<vsimd> > &Umu,int sequence,std::string LFN,std::string description) 
+  uint32_t readConfiguration(Lattice<iLorentzColourMatrix<vsimd> > &Umu) {
-  {
+    typedef Lattice<iLorentzColourMatrix<vsimd> > GaugeField;
-    GridBase * grid = Umu._grid;
+    typedef LorentzColourMatrixD sobjd;
+    typedef LorentzColourMatrixF sobjf;
+    typedef iLorentzColourMatrix<vsimd> itype;
+    typedef LorentzColourMatrix sobj;
+    GridBase *grid = Umu._grid;
+
+    ILDGField header;
+    readHeader(header);
+
+    // now just the conf, ignore the header
+    std::string format = std::string("IEEE64BIG");
+    do {limeReaderNextRecord(LimeR);}
+    while (strncmp(limeReaderType(LimeR), "ildg-binary-data",16));
+
+    n_uint64_t nbytes = limeReaderBytes(LimeR);//size of this record (configuration)
+
+
+    ILDGtype ILDGt(true, LimeR);
+    // this is special for double prec data, just for the moment
+    uint32_t csum = BinaryIO::readObjectParallel< itype, sobjd >(
+       Umu, filename, ILDGMunger<sobjd, sobj>(), 0, format, ILDGt);
+
+    // Check configuration 
+    // todo
+
+    return csum;
+  }
+
+  template <class vsimd>
+  uint32_t writeConfiguration(Lattice<iLorentzColourMatrix<vsimd> > &Umu, std::string format) {
     typedef Lattice<iLorentzColourMatrix<vsimd> > GaugeField;
     typedef iLorentzColourMatrix<vsimd> vobj;
     typedef typename vobj::scalar_object sobj;
+    typedef LorentzColourMatrixD fobj;
 
-    uint64_t nbytes;
+    ILDGField header;
+    // fill the header
+    header.floating_point = format;
 
-    ////////////////////////////////////////
+    ILDGUnmunger<fobj, sobj> munge;
-    // fill the Grid header
+    unsigned int offset = writeHeader(header);
-    ////////////////////////////////////////
-    FieldMetaData header;
-    scidacRecord  _scidacRecord;
-    scidacFile    _scidacFile;
 
-    ScidacMetaData(Umu,header,_scidacRecord,_scidacFile);
+    BinaryIO::Uint32Checksum<vobj, fobj>(Umu, munge, header.checksum);
 
-    std::string format = header.floating_point;
+    // Write data record header
-    header.ensemble_id    = description;
+    n_uint64_t PayloadSize = sizeof(fobj) * Umu._grid->_gsites;
-    header.ensemble_label = description;
+    createHeader("ildg-binary-data", 0, 1, PayloadSize, LimeW);
-    header.sequence_number = sequence;
-    header.ildg_lfn = LFN;
 
-    assert ( (format == std::string("IEEE32BIG"))  
+    ILDGtype ILDGt(true, LimeW);
-           ||(format == std::string("IEEE64BIG")) );
+    uint32_t csum = BinaryIO::writeObjectParallel<vobj, fobj>(
+       Umu, filename, munge, 0, header.floating_point, ILDGt);
 
-    //////////////////////////////////////////////////////
+    limeWriterCloseRecord(LimeW);
-    // Fill ILDG header data struct
-    //////////////////////////////////////////////////////
-    ildgFormat ildgfmt ;
-    ildgfmt.field     = std::string("su3gauge");
 
-    if ( format == std::string("IEEE32BIG") ) { 
+    // Last record
-      ildgfmt.precision = 32;
+    // the logical file name LNF
-    } else { 
+    // look into documentation on how to generate this string
-      ildgfmt.precision = 64;
+    std::string LNF = "empty"; 
-    }
-    ildgfmt.version = 1.0;
-    ildgfmt.lx = header.dimension[0];
-    ildgfmt.ly = header.dimension[1];
-    ildgfmt.lz = header.dimension[2];
-    ildgfmt.lt = header.dimension[3];
-    assert(header.nd==4);
-    assert(header.nd==header.dimension.size());
 
-    //////////////////////////////////////////////////////////////////////////////
-    // Fill the USQCD info field
-    //////////////////////////////////////////////////////////////////////////////
-    usqcdInfo info;
-    info.version=1.0;
-    info.plaq   = header.plaquette;
-    info.linktr = header.link_trace;
 
-    std::cout << GridLogMessage << " Writing config; IldgIO "<<std::endl;
+    PayloadSize = sizeof(LNF);
-    //////////////////////////////////////////////
+    createHeader("ildg-binary-lfn", 1 , 1, PayloadSize, LimeW);
-    // Fill the Lime file record by record
+    limeWriteRecordData(const_cast<char*>(LNF.c_str()), &PayloadSize, LimeW);
-    //////////////////////////////////////////////
+
-    writeLimeObject(1,0,header ,std::string("FieldMetaData"),std::string(GRID_FORMAT)); // Open message 
+    limeWriterCloseRecord(LimeW);
-    writeLimeObject(0,0,_scidacFile,_scidacFile.SerialisableClassName(),std::string(SCIDAC_PRIVATE_FILE_XML));
+
-    writeLimeObject(0,1,info,info.SerialisableClassName(),std::string(SCIDAC_FILE_XML));
+    return csum;
-    writeLimeObject(1,0,_scidacRecord,_scidacRecord.SerialisableClassName(),std::string(SCIDAC_PRIVATE_RECORD_XML));
-    writeLimeObject(0,0,info,info.SerialisableClassName(),std::string(SCIDAC_RECORD_XML));
-    writeLimeObject(0,0,ildgfmt,std::string("ildgFormat")   ,std::string(ILDG_FORMAT)); // rec
-    writeLimeIldgLFN(header.ildg_lfn);                                                 // rec
-    writeLimeLatticeBinaryObject(Umu,std::string(ILDG_BINARY_DATA));      // Closes message with checksum
-    //    limeDestroyWriter(LimeW);
-    fclose(File);
   }
+
+  // format for RNG? Now just binary out
 };
-
+}
-class IldgReader : public GridLimeReader {
+}
- public:
-
-  ////////////////////////////////////////////////////////////////
-  // Read either Grid/SciDAC/ILDG configuration
-  // Don't require scidac records EXCEPT checksum
-  // Use Grid MetaData object if present.
-  // Else use ILDG MetaData object if present.
-  // Else use SciDAC MetaData object if present.
-  ////////////////////////////////////////////////////////////////
-  template <class vsimd>
-  void readConfiguration(Lattice<iLorentzColourMatrix<vsimd> > &Umu, FieldMetaData &FieldMetaData_) {
-
-    typedef Lattice<iLorentzColourMatrix<vsimd> > GaugeField;
-    typedef typename GaugeField::vector_object  vobj;
-    typedef typename vobj::scalar_object sobj;
-
-    typedef LorentzColourMatrixF fobj;
-    typedef LorentzColourMatrixD dobj;
-
-    GridBase *grid = Umu._grid;
-
-    std::vector<int> dims = Umu._grid->FullDimensions();
-
-    assert(dims.size()==4);
-
-    // Metadata holders
-    ildgFormat     ildgFormat_    ;
-    std::string    ildgLFN_       ;
-    scidacChecksum scidacChecksum_; 
-    usqcdInfo      usqcdInfo_     ;
-
-    // track what we read from file
-    int found_ildgFormat    =0;
-    int found_ildgLFN       =0;
-    int found_scidacChecksum=0;
-    int found_usqcdInfo     =0;
-    int found_ildgBinary =0;
-    int found_FieldMetaData =0;
-
-    uint32_t nersc_csum;
-    uint32_t scidac_csuma;
-    uint32_t scidac_csumb;
-
-    // Binary format
-    std::string format;
-
-    //////////////////////////////////////////////////////////////////////////
-    // Loop over all records
-    // -- Order is poorly guaranteed except ILDG header preceeds binary section.
-    // -- Run like an event loop.
-    // -- Impose trust hierarchy. Grid takes precedence & look for ILDG, and failing
-    //    that Scidac. 
-    // -- Insist on Scidac checksum record.
-    //////////////////////////////////////////////////////////////////////////
-
-    while ( limeReaderNextRecord(LimeR) == LIME_SUCCESS ) { 
-
-      uint64_t nbytes = limeReaderBytes(LimeR);//size of this record (configuration)
-      
-      //////////////////////////////////////////////////////////////////
-      // If not BINARY_DATA read a string and parse
-      //////////////////////////////////////////////////////////////////
-      if ( strncmp(limeReaderType(LimeR), ILDG_BINARY_DATA,strlen(ILDG_BINARY_DATA) )  ) {
-	
-	// Copy out the string
-	std::vector<char> xmlc(nbytes+1,'\0');
-	limeReaderReadData((void *)&xmlc[0], &nbytes, LimeR);    
-	std::cout << GridLogMessage<< "Non binary record :" <<limeReaderType(LimeR) <<std::endl; //<<"\n"<<(&xmlc[0])<<std::endl;
-
-	//////////////////////////////////
-	// ILDG format record
-	if ( !strncmp(limeReaderType(LimeR), ILDG_FORMAT,strlen(ILDG_FORMAT)) ) { 
-
-	  XmlReader RD(&xmlc[0],"");
-	  read(RD,"ildgFormat",ildgFormat_);
-
-	  if ( ildgFormat_.precision == 64 ) format = std::string("IEEE64BIG");
-	  if ( ildgFormat_.precision == 32 ) format = std::string("IEEE32BIG");
-
-	  assert( ildgFormat_.lx == dims[0]);
-	  assert( ildgFormat_.ly == dims[1]);
-	  assert( ildgFormat_.lz == dims[2]);
-	  assert( ildgFormat_.lt == dims[3]);
-
-	  found_ildgFormat = 1;
-	}
-
-	if ( !strncmp(limeReaderType(LimeR), ILDG_DATA_LFN,strlen(ILDG_DATA_LFN)) ) {
-	  FieldMetaData_.ildg_lfn = std::string(&xmlc[0]);
-	  found_ildgLFN = 1;
-	}
-
-	if ( !strncmp(limeReaderType(LimeR), GRID_FORMAT,strlen(ILDG_FORMAT)) ) { 
-
-	  XmlReader RD(&xmlc[0],"");
-	  read(RD,"FieldMetaData",FieldMetaData_);
-
-	  format = FieldMetaData_.floating_point;
-
-	  assert(FieldMetaData_.dimension[0] == dims[0]);
-	  assert(FieldMetaData_.dimension[1] == dims[1]);
-	  assert(FieldMetaData_.dimension[2] == dims[2]);
-	  assert(FieldMetaData_.dimension[3] == dims[3]);
-
-	  found_FieldMetaData = 1;
-	}
-
-	if ( !strncmp(limeReaderType(LimeR), SCIDAC_RECORD_XML,strlen(SCIDAC_RECORD_XML)) ) { 
-	  XmlReader RD(&xmlc[0],"");
-	  read(RD,"usqcdInfo",usqcdInfo_);
-	  found_usqcdInfo = 1;
-	}
-
-	if ( !strncmp(limeReaderType(LimeR), SCIDAC_CHECKSUM,strlen(SCIDAC_CHECKSUM)) ) { 
-	  XmlReader RD(&xmlc[0],"");
-	  read(RD,"scidacChecksum",scidacChecksum_);
-	  found_scidacChecksum = 1;
-	}
-
-      } else {  
-	/////////////////////////////////
-	// Binary data
-	/////////////////////////////////
-	std::cout << GridLogMessage << "ILDG Binary record found : "  ILDG_BINARY_DATA << std::endl;
-	off_t offset= ftell(File);
-
-	if ( format == std::string("IEEE64BIG") ) {
-	  GaugeSimpleMunger<dobj, sobj> munge;
-	  BinaryIO::readLatticeObject< vobj, dobj >(Umu, filename, munge, offset, format,nersc_csum,scidac_csuma,scidac_csumb);
-	} else { 
-	  GaugeSimpleMunger<fobj, sobj> munge;
-	  BinaryIO::readLatticeObject< vobj, fobj >(Umu, filename, munge, offset, format,nersc_csum,scidac_csuma,scidac_csumb);
-	}
-
-	found_ildgBinary = 1;
-      }
-
-    }
-
-    //////////////////////////////////////////////////////
-    // Minimally must find binary segment and checksum
-    // Since this is an ILDG reader require ILDG format
-    //////////////////////////////////////////////////////
-    assert(found_ildgBinary);
-    assert(found_ildgFormat);
-    assert(found_scidacChecksum);
-
-    // Must find something with the lattice dimensions
-    assert(found_FieldMetaData||found_ildgFormat);
-
-    if ( found_FieldMetaData ) {
-
-      std::cout << GridLogMessage<<"Grid MetaData was record found: configuration was probably written by Grid ! Yay ! "<<std::endl;
-
-    } else { 
-
-      assert(found_ildgFormat);
-      assert ( ildgFormat_.field == std::string("su3gauge") );
-
-      ///////////////////////////////////////////////////////////////////////////////////////
-      // Populate our Grid metadata as best we can
-      ///////////////////////////////////////////////////////////////////////////////////////
-
-      std::ostringstream vers; vers << ildgFormat_.version;
-      FieldMetaData_.hdr_version = vers.str();
-      FieldMetaData_.data_type = std::string("4D_SU3_GAUGE_3X3");
-
-      FieldMetaData_.nd=4;
-      FieldMetaData_.dimension.resize(4);
-
-      FieldMetaData_.dimension[0] = ildgFormat_.lx ;
-      FieldMetaData_.dimension[1] = ildgFormat_.ly ;
-      FieldMetaData_.dimension[2] = ildgFormat_.lz ;
-      FieldMetaData_.dimension[3] = ildgFormat_.lt ;
-
-      if ( found_usqcdInfo ) { 
-	FieldMetaData_.plaquette = usqcdInfo_.plaq;
-	FieldMetaData_.link_trace= usqcdInfo_.linktr;
-	std::cout << GridLogMessage <<"This configuration was probably written by USQCD "<<std::endl;
-	std::cout << GridLogMessage <<"USQCD xml record Plaquette : "<<FieldMetaData_.plaquette<<std::endl;
-	std::cout << GridLogMessage <<"USQCD xml record LinkTrace : "<<FieldMetaData_.link_trace<<std::endl;
-      } else { 
-	FieldMetaData_.plaquette = 0.0;
-	FieldMetaData_.link_trace= 0.0;
-	std::cout << GridLogWarning << "This configuration is unsafe with no plaquette records that can verify it !!! "<<std::endl;
-      }
-    }
-
-    ////////////////////////////////////////////////////////////
-    // Really really want to mandate a scidac checksum
-    ////////////////////////////////////////////////////////////
-    if ( found_scidacChecksum ) {
-      FieldMetaData_.scidac_checksuma = stoull(scidacChecksum_.suma,0,16);
-      FieldMetaData_.scidac_checksumb = stoull(scidacChecksum_.sumb,0,16);
-      scidacChecksumVerify(scidacChecksum_,scidac_csuma,scidac_csumb);
-      assert( scidac_csuma ==FieldMetaData_.scidac_checksuma);
-      assert( scidac_csumb ==FieldMetaData_.scidac_checksumb);
-      std::cout << GridLogMessage<<"SciDAC checksums match " << std::endl;
-    } else { 
-      std::cout << GridLogWarning<<"SciDAC checksums not found. This is unsafe. " << std::endl;
-      assert(0); // Can I insist always checksum ?
-    }
-
-    if ( found_FieldMetaData || found_usqcdInfo ) {
-      FieldMetaData checker;
-      GaugeStatistics(Umu,checker);
-      assert(fabs(checker.plaquette  - FieldMetaData_.plaquette )<1.0e-5);
-      assert(fabs(checker.link_trace - FieldMetaData_.link_trace)<1.0e-5);
-      std::cout << GridLogMessage<<"Plaquette and link trace match " << std::endl;
-    }
-  }
- };
-
-}}
 
 //HAVE_LIME
 #endif

lib/parallelIO/IldgIOtypes.h

@@ -34,198 +34,47 @@ extern "C" { // for linkage
 
 namespace Grid {
 
-/////////////////////////////////////////////////////////////////////////////////
+struct ILDGtype {
-// Data representation of records that enter ILDG and SciDac formats
+  bool is_ILDG;
-/////////////////////////////////////////////////////////////////////////////////
+  LimeWriter* LW;
+  LimeReader* LR;
 
-#define GRID_FORMAT      "grid-format"
+  ILDGtype(bool is, LimeWriter* L) : is_ILDG(is), LW(L), LR(NULL) {}
-#define ILDG_FORMAT      "ildg-format"
+  ILDGtype(bool is, LimeReader* L) : is_ILDG(is), LW(NULL), LR(L) {}
-#define ILDG_BINARY_DATA "ildg-binary-data"
+  ILDGtype() : is_ILDG(false), LW(NULL), LR(NULL) {}
-#define ILDG_DATA_LFN    "ildg-data-lfn"
+};
-#define SCIDAC_CHECKSUM           "scidac-checksum"
-#define SCIDAC_PRIVATE_FILE_XML   "scidac-private-file-xml"
-#define SCIDAC_FILE_XML           "scidac-file-xml"
-#define SCIDAC_PRIVATE_RECORD_XML "scidac-private-record-xml"
-#define SCIDAC_RECORD_XML         "scidac-record-xml"
-#define SCIDAC_BINARY_DATA        "scidac-binary-data"
-// Unused SCIDAC records names; could move to support this functionality
-#define SCIDAC_SITELIST           "scidac-sitelist"
 
-  ////////////////////////////////////////////////////////////
+class ILDGField {
-  const int GRID_IO_SINGLEFILE = 0; // hardcode lift from QIO compat
-  const int GRID_IO_MULTIFILE  = 1; // hardcode lift from QIO compat
-  const int GRID_IO_FIELD      = 0; // hardcode lift from QIO compat
-  const int GRID_IO_GLOBAL     = 1; // hardcode lift from QIO compat
-  ////////////////////////////////////////////////////////////
-
-/////////////////////////////////////////////////////////////////////////////////
-// QIO uses mandatory "private" records fixed format
-// Private is in principle "opaque" however it can't be changed now because that would break existing 
-// file compatability, so should be correct to assume the undocumented but defacto file structure.
-/////////////////////////////////////////////////////////////////////////////////
-
-////////////////////////
-// Scidac private file xml
-// <?xml version="1.0" encoding="UTF-8"?><scidacFile><version>1.1</version><spacetime>4</spacetime><dims>16 16 16 32 </dims><volfmt>0</volfmt></scidacFile>
-////////////////////////
-struct scidacFile : Serializable {
  public:
-  GRID_SERIALIZABLE_CLASS_MEMBERS(scidacFile,
+  // header strings (not in order)
-                                  double, version,
+  std::vector<int> dimension;
-                                  int, spacetime,
+  std::vector<std::string> boundary;
-				  std::string, dims, // must convert to int
+  int data_start;
-                                  int, volfmt);
+  std::string hdr_version;
-
+  std::string storage_format;
-  std::vector<int> getDimensions(void) { 
+  // Checks on data
-    std::stringstream stream(dims);
+  double link_trace;
-    std::vector<int> dimensions;
+  double plaquette;
-    int n;
+  uint32_t checksum;
-    while(stream >> n){
+  unsigned int sequence_number;
-      dimensions.push_back(n);
+  std::string data_type;
-    }
+  std::string ensemble_id;
-    return dimensions;
+  std::string ensemble_label;
-  }
+  std::string creator;
-
+  std::string creator_hardware;
-  void setDimensions(std::vector<int> dimensions) { 
+  std::string creation_date;
-    char delimiter = ' ';
+  std::string archive_date;
-    std::stringstream stream;
+  std::string floating_point;
-    for(int i=0;i<dimensions.size();i++){ 
-      stream << dimensions[i];
-      if ( i != dimensions.size()-1) { 
-	stream << delimiter <<std::endl;
-      }
-    }
-    dims = stream.str();
-  }
-
-  // Constructor provides Grid
-  scidacFile() =default; // default constructor
-  scidacFile(GridBase * grid){
-    version      = 1.0;
-    spacetime    = grid->_ndimension;
-    setDimensions(grid->FullDimensions()); 
-    volfmt       = GRID_IO_SINGLEFILE;
-  }
-
 };
-
-///////////////////////////////////////////////////////////////////////
-// scidac-private-record-xml : example
-// <scidacRecord>
-// <version>1.1</version><date>Tue Jul 26 21:14:44 2011 UTC</date><recordtype>0</recordtype>
-// <datatype>QDP_D3_ColorMatrix</datatype><precision>D</precision><colors>3</colors><spins>4</spins>
-// <typesize>144</typesize><datacount>4</datacount>
-// </scidacRecord>
-///////////////////////////////////////////////////////////////////////
-
-struct scidacRecord : Serializable {
- public:
-  GRID_SERIALIZABLE_CLASS_MEMBERS(scidacRecord,
-                                  double, version,
-                                  std::string, date,
-				  int, recordtype,
-				  std::string, datatype,
-				  std::string, precision,
-				  int, colors,
-				  int, spins,
-				  int, typesize,
-				  int, datacount);
-
-  scidacRecord() { version =1.0; }
-
-};
-
-////////////////////////
-// ILDG format
-////////////////////////
-struct ildgFormat : Serializable {
-public:
-  GRID_SERIALIZABLE_CLASS_MEMBERS(ildgFormat,
-				  double, version,
-				  std::string, field,
-				  int, precision,
-				  int, lx,
-				  int, ly,
-				  int, lz,
-				  int, lt);
-  ildgFormat() { version=1.0; };
-};
-////////////////////////
-// USQCD info
-////////////////////////
-struct usqcdInfo : Serializable { 
- public:
-  GRID_SERIALIZABLE_CLASS_MEMBERS(usqcdInfo,
-				  double, version,
-				  double, plaq,
-				  double, linktr,
-				  std::string, info);
-  usqcdInfo() { 
-    version=1.0; 
-  };
-};
-////////////////////////
-// Scidac Checksum
-////////////////////////
-struct scidacChecksum : Serializable { 
- public:
-  GRID_SERIALIZABLE_CLASS_MEMBERS(scidacChecksum,
-				  double, version,
-				  std::string, suma,
-				  std::string, sumb);
-  scidacChecksum() { 
-    version=1.0; 
-  };
-};
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-// Type:           scidac-file-xml         <title>MILC ILDG archival gauge configuration</title>
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-// Type:           
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-////////////////////////
-// Scidac private file xml 
-// <?xml version="1.0" encoding="UTF-8"?><scidacFile><version>1.1</version><spacetime>4</spacetime><dims>16 16 16 32 </dims><volfmt>0</volfmt></scidacFile> 
-////////////////////////                                                                                                                                                                              
-
-#if 0
-////////////////////////////////////////////////////////////////////////////////////////
-// From http://www.physics.utah.edu/~detar/scidac/qio_2p3.pdf
-////////////////////////////////////////////////////////////////////////////////////////
-struct usqcdPropFile : Serializable { 
- public:
-  GRID_SERIALIZABLE_CLASS_MEMBERS(usqcdPropFile,
-				  double, version,
-				  std::string, type,
-				  std::string, info);
-  usqcdPropFile() { 
-    version=1.0; 
-  };
-};
-struct usqcdSourceInfo : Serializable { 
- public:
-  GRID_SERIALIZABLE_CLASS_MEMBERS(usqcdSourceInfo,
-				  double, version,
-				  std::string, info);
-  usqcdSourceInfo() { 
-    version=1.0; 
-  };
-};
-struct usqcdPropInfo : Serializable { 
- public:
-  GRID_SERIALIZABLE_CLASS_MEMBERS(usqcdPropInfo,
-				  double, version,
-				  int, spin,
-				  int, color,
-				  std::string, info);
-  usqcdPropInfo() { 
-    version=1.0; 
-  };
-};
-#endif
-
 }
+#else
+namespace Grid {
+
+struct ILDGtype {
+  bool is_ILDG;
+  ILDGtype() : is_ILDG(false) {}
+};
+}
+
 #endif
 #endif

lib/parallelIO/MetaData.h

@@ -1,325 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./lib/parallelIO/NerscIO.h
-
-    Copyright (C) 2015
-
-
-    Author: Peter Boyle <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 */
-
-#include <algorithm>
-#include <iostream>
-#include <iomanip>
-#include <fstream>
-#include <map>
-#include <unistd.h>
-#include <sys/utsname.h>
-#include <pwd.h>
-
-namespace Grid {
-
-  ///////////////////////////////////////////////////////
-  // 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:
-
-      GRID_SERIALIZABLE_CLASS_MEMBERS(FieldMetaData,
-				      int, nd,
-				      std::vector<int>, dimension,
-				      std::vector<std::string>, boundary,
-				      int, data_start,
-				      std::string, hdr_version,
-				      std::string, storage_format,
-				      double, link_trace,
-				      double, plaquette,
-				      uint32_t, checksum,
-				      uint32_t, scidac_checksuma,
-				      uint32_t, scidac_checksumb,
-				      unsigned int, sequence_number,
-				      std::string, data_type,
-				      std::string, ensemble_id,
-				      std::string, ensemble_label,
-				      std::string, ildg_lfn,
-				      std::string, creator,
-				      std::string, creator_hardware,
-				      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;
-
-
-    //////////////////////////////////////////////////////////////////////
-    // Bit and Physical Checksumming and QA of data
-    //////////////////////////////////////////////////////////////////////
-    inline void GridMetaData(GridBase *grid,FieldMetaData &header)
-    {
-      int nd = grid->_ndimension;
-      header.nd = nd;
-      header.dimension.resize(nd);
-      header.boundary.resize(nd);
-      for(int d=0;d<nd;d++) {
-	header.dimension[d] = grid->_fdimensions[d];
-      }
-      for(int d=0;d<nd;d++) {
-	header.boundary[d] = std::string("PERIODIC");
-      }
-    }
-
-    inline void MachineCharacteristics(FieldMetaData &header)
-    {
-      // Who
-      struct passwd *pw = getpwuid (getuid());
-      if (pw) header.creator = std::string(pw->pw_name); 
-
-      // When
-      std::time_t t = std::time(nullptr);
-      std::tm tm_ = *std::localtime(&t);
-      std::ostringstream oss; 
-      //      oss << std::put_time(&tm_, "%c %Z");
-      header.creation_date = oss.str();
-      header.archive_date  = header.creation_date;
-
-      // What
-      struct utsname name;  uname(&name);
-      header.creator_hardware = std::string(name.nodename)+"-";
-      header.creator_hardware+= std::string(name.machine)+"-";
-      header.creator_hardware+= std::string(name.sysname)+"-";
-      header.creator_hardware+= std::string(name.release);
-    }
-
-#define dump_meta_data(field, s)					\
-      s << "BEGIN_HEADER"      << std::endl;				\
-      s << "HDR_VERSION = "    << field.hdr_version    << std::endl;	\
-      s << "DATATYPE = "       << field.data_type      << std::endl;	\
-      s << "STORAGE_FORMAT = " << field.storage_format << std::endl;	\
-      for(int i=0;i<4;i++){						\
-	s << "DIMENSION_" << i+1 << " = " << field.dimension[i] << std::endl ; \
-      }									\
-      s << "LINK_TRACE = " << std::setprecision(10) << field.link_trace << std::endl; \
-      s << "PLAQUETTE  = " << std::setprecision(10) << field.plaquette  << std::endl; \
-      for(int i=0;i<4;i++){						\
-	s << "BOUNDARY_"<<i+1<<" = " << field.boundary[i] << 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;	\
-      s << "ENSEMBLE_LABEL = "  << field.ensemble_label   << std::endl;	\
-      s << "SEQUENCE_NUMBER = " << field.sequence_number  << std::endl;	\
-      s << "CREATOR = "         << field.creator          << std::endl;	\
-      s << "CREATOR_HARDWARE = "<< field.creator_hardware << std::endl;	\
-      s << "CREATION_DATE = "   << field.creation_date    << std::endl;	\
-      s << "ARCHIVE_DATE = "    << field.archive_date     << std::endl;	\
-      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
-    //////////////////////////////////////////////////////////////////////
-    inline void reconstruct3(LorentzColourMatrix & cm)
-    {
-      const int x=0;
-      const int y=1;
-      const int z=2;
-      for(int mu=0;mu<Nd;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,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
-      }
-    }
-
-    ////////////////////////////////////////////////////////////////////////////////
-    // Some data types for intermediate storage
-    ////////////////////////////////////////////////////////////////////////////////
-    template<typename vtype> using iLorentzColour2x3 = iVector<iVector<iVector<vtype, Nc>, 2>, Nd >;
-
-    typedef iLorentzColour2x3<Complex>  LorentzColour2x3;
-    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) {
-        for (int mu = 0; mu < Nd; mu++) {
-          for (int i = 0; i < Nc; i++) {
-          for (int j = 0; j < Nc; j++) {
-	    out(mu)()(i, j) = in(mu)()(i, j);
-	  }}
-        }
-      };
-    };
-
-    template <class fobj, class sobj>
-    struct GaugeSimpleUnmunger {
-
-      void operator()(sobj &in, fobj &out) {
-        for (int mu = 0; mu < Nd; mu++) {
-          for (int i = 0; i < Nc; i++) {
-          for (int j = 0; j < Nc; j++) {
-	    out(mu)()(i, j) = in(mu)()(i, j);
-	  }}
-        }
-      };
-    };
-
-    template<class fobj,class sobj>
-    struct Gauge3x2munger{
-      void operator() (fobj &in,sobj &out){
-	for(int mu=0;mu<Nd;mu++){
-	  for(int i=0;i<2;i++){
-	  for(int j=0;j<3;j++){
-	    out(mu)()(i,j) = in(mu)(i)(j);
-	  }}
-	}
-	reconstruct3(out);
-      }
-    };
-
-    template<class fobj,class sobj>
-    struct Gauge3x2unmunger{
-      void operator() (sobj &in,fobj &out){
-	for(int mu=0;mu<Nd;mu++){
-	  for(int i=0;i<2;i++){
-	  for(int j=0;j<3;j++){
-	    out(mu)(i)(j) = in(mu)()(i,j);
-	  }}
-	}
-      }
-    };
-  }
-
-
-}

lib/parallelIO/NerscIO.h

@@ -30,11 +30,182 @@
 #ifndef GRID_NERSC_IO_H
 #define GRID_NERSC_IO_H
 
+#include <algorithm>
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+#include <map>
+
+#include <unistd.h>
+#include <sys/utsname.h>
+#include <pwd.h>
+
 namespace Grid {
   namespace QCD {
 
     using namespace Grid;
 
+    ////////////////////////////////////////////////////////////////////////////////
+    // Some data types for intermediate storage
+    ////////////////////////////////////////////////////////////////////////////////
+    template<typename vtype> using iLorentzColour2x3 = iVector<iVector<iVector<vtype, Nc>, 2>, 4 >;
+
+    typedef iLorentzColour2x3<Complex>  LorentzColour2x3;
+    typedef iLorentzColour2x3<ComplexF> LorentzColour2x3F;
+    typedef iLorentzColour2x3<ComplexD> LorentzColour2x3D;
+
+    ////////////////////////////////////////////////////////////////////////////////
+    // header specification/interpretation
+    ////////////////////////////////////////////////////////////////////////////////
+    class NerscField {
+    public:
+      // header strings (not in order)
+      int dimension[4];
+      std::string boundary[4]; 
+      int data_start;
+      std::string hdr_version;
+      std::string storage_format;
+      // Checks on data
+      double link_trace;
+      double plaquette;
+      uint32_t checksum;
+      unsigned int sequence_number;
+      std::string data_type;
+      std::string ensemble_id ;
+      std::string ensemble_label ;
+      std::string creator ;
+      std::string creator_hardware ;
+      std::string creation_date ;
+      std::string archive_date ;
+      std::string floating_point;
+    };
+
+    //////////////////////////////////////////////////////////////////////
+    // Bit and Physical Checksumming and QA of data
+    //////////////////////////////////////////////////////////////////////
+
+    inline void NerscGrid(GridBase *grid,NerscField &header)
+    {
+      assert(grid->_ndimension==4);
+      for(int d=0;d<4;d++) {
+	header.dimension[d] = grid->_fdimensions[d];
+      }
+      for(int d=0;d<4;d++) {
+	header.boundary[d] = std::string("PERIODIC");
+      }
+    }
+    template<class GaugeField>
+    inline void NerscStatistics(GaugeField & data,NerscField &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 NerscMachineCharacteristics(NerscField &header)
+    {
+      // Who
+      struct passwd *pw = getpwuid (getuid());
+      if (pw) header.creator = std::string(pw->pw_name); 
+
+      // When
+      std::time_t t = std::time(nullptr);
+      std::tm tm = *std::localtime(&t);
+      std::ostringstream oss; 
+      //  oss << std::put_time(&tm, "%c %Z");
+      header.creation_date = oss.str();
+      header.archive_date  = header.creation_date;
+
+      // What
+      struct utsname name;  uname(&name);
+      header.creator_hardware = std::string(name.nodename)+"-";
+      header.creator_hardware+= std::string(name.machine)+"-";
+      header.creator_hardware+= std::string(name.sysname)+"-";
+      header.creator_hardware+= std::string(name.release);
+
+    }
+    //////////////////////////////////////////////////////////////////////
+    // Utilities ; these are QCD aware
+    //////////////////////////////////////////////////////////////////////
+    inline void NerscChecksum(uint32_t *buf,uint32_t buf_size_bytes,uint32_t &csum)
+    {
+      BinaryIO::Uint32Checksum(buf,buf_size_bytes,csum);
+    }
+    inline void reconstruct3(LorentzColourMatrix & cm)
+    {
+      const int x=0;
+      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,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
+      }
+    }
+
+    template<class fobj,class sobj>
+    struct NerscSimpleMunger{
+      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 j = 0; j < Nc; j++) {
+              out(mu)()(i, j) = in(mu)()(i, j);
+            }
+          }
+        }
+        NerscChecksum((uint32_t *)&in, sizeof(in), csum);
+      };
+    };
+
+    template <class fobj, class sobj>
+    struct NerscSimpleUnmunger {
+      void operator()(sobj &in, fobj &out, uint32_t &csum) {
+        for (int mu = 0; mu < Nd; mu++) {
+          for (int i = 0; i < Nc; i++) {
+            for (int j = 0; j < Nc; j++) {
+              out(mu)()(i, j) = in(mu)()(i, j);
+            }
+          }
+        }
+        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); 
+
+	for(int mu=0;mu<4;mu++){
+	  for(int i=0;i<2;i++){
+	    for(int j=0;j<3;j++){
+	      out(mu)()(i,j) = in(mu)(i)(j);
+	    }}
+	}
+	reconstruct3(out);
+      }
+    };
+
+    template<class fobj,class sobj>
+    struct Nersc3x2unmunger{
+
+      void operator() (sobj &in,fobj &out,uint32_t &csum){
+
+
+	for(int mu=0;mu<4;mu++){
+	  for(int i=0;i<2;i++){
+	    for(int j=0;j<3;j++){
+	      out(mu)(i)(j) = in(mu)()(i,j);
+	    }}
+	}
+
+	NerscChecksum((uint32_t *)&out,sizeof(out),csum); 
+
+      }
+    };
+
+
     ////////////////////////////////////////////////////////////////////////////////
     // Write and read from fstream; comput header offset for payload
     ////////////////////////////////////////////////////////////////////////////////
@@ -45,17 +216,42 @@ namespace Grid {
 	std::ofstream fout(file,std::ios::out);
       }
   
-      static inline unsigned int writeHeader(FieldMetaData &field,std::string file)
+#define dump_nersc_header(field, s)					\
+      s << "BEGIN_HEADER"      << std::endl;				\
+      s << "HDR_VERSION = "    << field.hdr_version    << std::endl;	\
+      s << "DATATYPE = "       << field.data_type      << std::endl;	\
+      s << "STORAGE_FORMAT = " << field.storage_format << std::endl;	\
+      for(int i=0;i<4;i++){						\
+	s << "DIMENSION_" << i+1 << " = " << field.dimension[i] << std::endl ; \
+      }									\
+      s << "LINK_TRACE = " << std::setprecision(10) << field.link_trace << std::endl; \
+      s << "PLAQUETTE  = " << std::setprecision(10) << field.plaquette  << std::endl; \
+      for(int i=0;i<4;i++){						\
+	s << "BOUNDARY_"<<i+1<<" = " << field.boundary[i] << std::endl;	\
+      }									\
+									\
+      s << "CHECKSUM = "<< std::hex << std::setw(10) << field.checksum << std::dec<<std::endl; \
+      s << "ENSEMBLE_ID = "     << field.ensemble_id      << std::endl;	\
+      s << "ENSEMBLE_LABEL = "  << field.ensemble_label   << std::endl;	\
+      s << "SEQUENCE_NUMBER = " << field.sequence_number  << std::endl;	\
+      s << "CREATOR = "         << field.creator          << std::endl;	\
+      s << "CREATOR_HARDWARE = "<< field.creator_hardware << std::endl;	\
+      s << "CREATION_DATE = "   << field.creation_date    << std::endl;	\
+      s << "ARCHIVE_DATE = "    << field.archive_date     << std::endl;	\
+      s << "FLOATING_POINT = "  << field.floating_point   << std::endl;	\
+      s << "END_HEADER"         << std::endl;
+  
+      static inline unsigned int writeHeader(NerscField &field,std::string file)
       {
       std::ofstream fout(file,std::ios::out|std::ios::in);
       fout.seekp(0,std::ios::beg);
-      dump_meta_data(field, fout);
+      dump_nersc_header(field, fout);
       field.data_start = fout.tellp();
       return field.data_start;
     }
 
       // for the header-reader
-      static inline int readHeader(std::string file,GridBase *grid,  FieldMetaData &field)
+      static inline int readHeader(std::string file,GridBase *grid,  NerscField &field)
       {
       int offset=0;
       std::map<std::string,std::string> header;
@@ -127,21 +323,21 @@ namespace Grid {
       return field.data_start;
     }
 
-    /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+      /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-    // Now the meat: the object readers
+      // Now the meat: the object readers
-    /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+      /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+#define PARALLEL_READ
+#define PARALLEL_WRITE
 
-    template<class vsimd>
+      template<class vsimd>
-    static inline void readConfiguration(Lattice<iLorentzColourMatrix<vsimd> > &Umu,
+      static inline void readConfiguration(Lattice<iLorentzColourMatrix<vsimd> > &Umu,NerscField& header,std::string file)
-					 FieldMetaData& header,
+      {
-					 std::string file)
-    {
       typedef Lattice<iLorentzColourMatrix<vsimd> > GaugeField;
 
       GridBase *grid = Umu._grid;
       int offset = readHeader(file,Umu._grid,header);
 
-      FieldMetaData clone(header);
+      NerscField clone(header);
 
       std::string format(header.floating_point);
 
@@ -150,78 +346,76 @@ namespace Grid {
       int ieee64big = (format == std::string("IEEE64BIG"));
       int ieee64    = (format == std::string("IEEE64"));
 
-      uint32_t nersc_csum,scidac_csuma,scidac_csumb;
+      uint32_t csum;
       // depending on datatype, set up munger;
       // munger is a function of <floating point, Real, data_type>
       if ( header.data_type == std::string("4D_SU3_GAUGE") ) {
-	if ( ieee32 || ieee32big ) {
+      if ( ieee32 || ieee32big ) {
-	  BinaryIO::readLatticeObject<iLorentzColourMatrix<vsimd>, LorentzColour2x3F> 
+#ifdef PARALLEL_READ
-	    (Umu,file,Gauge3x2munger<LorentzColour2x3F,LorentzColourMatrix>(), offset,format,
+	csum=BinaryIO::readObjectParallel<iLorentzColourMatrix<vsimd>, LorentzColour2x3F> 
-	     nersc_csum,scidac_csuma,scidac_csumb);
+	  (Umu,file,Nersc3x2munger<LorentzColour2x3F,LorentzColourMatrix>(), offset,format);
-	}
+#else
-	if ( ieee64 || ieee64big ) {
+	csum=BinaryIO::readObjectSerial<iLorentzColourMatrix<vsimd>, LorentzColour2x3F> 
-	  BinaryIO::readLatticeObject<iLorentzColourMatrix<vsimd>, LorentzColour2x3D> 
+	  (Umu,file,Nersc3x2munger<LorentzColour2x3F,LorentzColourMatrix>(), offset,format);
-	    (Umu,file,Gauge3x2munger<LorentzColour2x3D,LorentzColourMatrix>(),offset,format,
+#endif
-	     nersc_csum,scidac_csuma,scidac_csumb);
+      }
-	}
+      if ( ieee64 || ieee64big ) {
+#ifdef PARALLEL_READ
+	csum=BinaryIO::readObjectParallel<iLorentzColourMatrix<vsimd>, LorentzColour2x3D> 
+	  (Umu,file,Nersc3x2munger<LorentzColour2x3D,LorentzColourMatrix>(),offset,format);
+#else 
+	csum=BinaryIO::readObjectSerial<iLorentzColourMatrix<vsimd>, LorentzColour2x3D> 
+	  (Umu,file,Nersc3x2munger<LorentzColour2x3D,LorentzColourMatrix>(),offset,format);
+#endif
+      }
       } else if ( header.data_type == std::string("4D_SU3_GAUGE_3x3") ) {
 	if ( ieee32 || ieee32big ) {
-	  BinaryIO::readLatticeObject<iLorentzColourMatrix<vsimd>,LorentzColourMatrixF>
+#ifdef PARALLEL_READ
-	    (Umu,file,GaugeSimpleMunger<LorentzColourMatrixF,LorentzColourMatrix>(),offset,format,
+	  csum=BinaryIO::readObjectParallel<iLorentzColourMatrix<vsimd>,LorentzColourMatrixF>
-	     nersc_csum,scidac_csuma,scidac_csumb);
+	    (Umu,file,NerscSimpleMunger<LorentzColourMatrixF,LorentzColourMatrix>(),offset,format);
+#else
+	  csum=BinaryIO::readObjectSerial<iLorentzColourMatrix<vsimd>,LorentzColourMatrixF>
+	    (Umu,file,NerscSimpleMunger<LorentzColourMatrixF,LorentzColourMatrix>(),offset,format);
+#endif
 	}
 	if ( ieee64 || ieee64big ) {
-	  BinaryIO::readLatticeObject<iLorentzColourMatrix<vsimd>,LorentzColourMatrixD>
+#ifdef PARALLEL_READ
-	    (Umu,file,GaugeSimpleMunger<LorentzColourMatrixD,LorentzColourMatrix>(),offset,format,
+	  csum=BinaryIO::readObjectParallel<iLorentzColourMatrix<vsimd>,LorentzColourMatrixD>
-	     nersc_csum,scidac_csuma,scidac_csumb);
+	    (Umu,file,NerscSimpleMunger<LorentzColourMatrixD,LorentzColourMatrix>(),offset,format);
+#else
+	  csum=BinaryIO::readObjectSerial<iLorentzColourMatrix<vsimd>,LorentzColourMatrixD>
+	    (Umu,file,NerscSimpleMunger<LorentzColourMatrixD,LorentzColourMatrix>(),offset,format);
+#endif
 	}
       } else {
 	assert(0);
       }
 
-      GaugeStatistics(Umu,clone);
+      NerscStatistics<GaugeField>(Umu,clone);
 
-      std::cout<<GridLogMessage <<"NERSC Configuration "<<file<<" checksum "<<std::hex<<nersc_csum<< std::dec
+      std::cout<<GridLogMessage <<"NERSC Configuration "<<file<<" checksum "<<std::hex<<            csum<< std::dec
 	       <<" header   "<<std::hex<<header.checksum<<std::dec <<std::endl;
       std::cout<<GridLogMessage <<"NERSC Configuration "<<file<<" plaquette "<<clone.plaquette
 	       <<" header    "<<header.plaquette<<std::endl;
       std::cout<<GridLogMessage <<"NERSC Configuration "<<file<<" link_trace "<<clone.link_trace
 	       <<" header    "<<header.link_trace<<std::endl;
-
-      if ( fabs(clone.plaquette -header.plaquette ) >=  1.0e-5 ) { 
-	std::cout << " Plaquette mismatch "<<std::endl;
-	std::cout << Umu[0]<<std::endl;
-	std::cout << Umu[1]<<std::endl;
-      }
-      if ( nersc_csum != header.checksum ) { 
-	std::cerr << " checksum mismatch " << std::endl;
-	std::cerr << " plaqs " << clone.plaquette << " " << header.plaquette << std::endl;
-	std::cerr << " trace " << clone.link_trace<< " " << header.link_trace<< std::endl;
-	std::cerr << " nersc_csum  " <<std::hex<< nersc_csum << " " << header.checksum<< std::dec<< std::endl;
-	exit(0);
-      }
       assert(fabs(clone.plaquette -header.plaquette ) < 1.0e-5 );
       assert(fabs(clone.link_trace-header.link_trace) < 1.0e-6 );
-      assert(nersc_csum == header.checksum );
+      assert(csum == header.checksum );
-      
+
       std::cout<<GridLogMessage <<"NERSC Configuration "<<file<< " and plaquette, link trace, and checksum agree"<<std::endl;
-    }
+      }
 
       template<class vsimd>
-      static inline void writeConfiguration(Lattice<iLorentzColourMatrix<vsimd> > &Umu,
+      static inline void writeConfiguration(Lattice<iLorentzColourMatrix<vsimd> > &Umu,std::string file, int two_row,int bits32)
-					    std::string file, 
-					    int two_row,
-					    int bits32)
       {
 	typedef Lattice<iLorentzColourMatrix<vsimd> > GaugeField;
 
 	typedef iLorentzColourMatrix<vsimd> vobj;
 	typedef typename vobj::scalar_object sobj;
 
-	FieldMetaData header;
-	///////////////////////////////////////////
 	// Following should become arguments
-	///////////////////////////////////////////
+	NerscField header;
 	header.sequence_number = 1;
 	header.ensemble_id     = "UKQCD";
 	header.ensemble_label  = "DWF";
@@ -231,32 +425,45 @@ namespace Grid {
   
 	GridBase *grid = Umu._grid;
 
-	GridMetaData(grid,header);
+	NerscGrid(grid,header);
-	assert(header.nd==4);
+	NerscStatistics<GaugeField>(Umu,header);
-	GaugeStatistics(Umu,header);
+	NerscMachineCharacteristics(header);
-	MachineCharacteristics(header);
 
+	uint32_t csum;
 	int offset;
   
 	truncate(file);
 
-	// Sod it -- always write 3x3 double
+	if ( two_row ) { 
-	header.floating_point = std::string("IEEE64BIG");
-	header.data_type      = std::string("4D_SU3_GAUGE_3x3");
-	GaugeSimpleUnmunger<fobj3D,sobj> munge;
-	offset = writeHeader(header,file);
 
-	uint32_t nersc_csum,scidac_csuma,scidac_csumb;
+	  header.floating_point = std::string("IEEE64BIG");
-	BinaryIO::writeLatticeObject<vobj,fobj3D>(Umu,file,munge,offset,header.floating_point,
+	  header.data_type      = std::string("4D_SU3_GAUGE");
-								  nersc_csum,scidac_csuma,scidac_csumb);
+	  Nersc3x2unmunger<fobj2D,sobj> munge;
-	header.checksum = nersc_csum;
+	  BinaryIO::Uint32Checksum<vobj,fobj2D>(Umu, munge,header.checksum);
-	writeHeader(header,file);
+	  offset = writeHeader(header,file);
+#ifdef PARALLEL_WRITE
+	  csum=BinaryIO::writeObjectParallel<vobj,fobj2D>(Umu,file,munge,offset,header.floating_point);
+#else
+	  csum=BinaryIO::writeObjectSerial<vobj,fobj2D>(Umu,file,munge,offset,header.floating_point);
+#endif
+	} else { 
+	  header.floating_point = std::string("IEEE64BIG");
+	  header.data_type      = std::string("4D_SU3_GAUGE_3x3");
+	  NerscSimpleUnmunger<fobj3D,sobj> munge;
+	  BinaryIO::Uint32Checksum<vobj,fobj3D>(Umu, munge,header.checksum);
+	  offset = writeHeader(header,file);
+#ifdef PARALLEL_WRITE
+	  csum=BinaryIO::writeObjectParallel<vobj,fobj3D>(Umu,file,munge,offset,header.floating_point);
+#else
+	  csum=BinaryIO::writeObjectSerial<vobj,fobj3D>(Umu,file,munge,offset,header.floating_point);
+#endif
+	}
 
-	std::cout<<GridLogMessage <<"Written NERSC Configuration on "<< file << " checksum "
+	std::cout<<GridLogMessage <<"Written NERSC Configuration on "<< file << " checksum "<<std::hex<<csum<< std::dec<<" plaq "<< header.plaquette <<std::endl;
-		 <<std::hex<<header.checksum
-		 <<std::dec<<" plaq "<< header.plaquette <<std::endl;
 
       }
+
+
       ///////////////////////////////
       // RNG state
       ///////////////////////////////
@@ -265,19 +472,19 @@ namespace Grid {
 	typedef typename GridParallelRNG::RngStateType RngStateType;
 
 	// Following should become arguments
-	FieldMetaData header;
+	NerscField header;
 	header.sequence_number = 1;
 	header.ensemble_id     = "UKQCD";
 	header.ensemble_label  = "DWF";
 
 	GridBase *grid = parallel._grid;
 
-	GridMetaData(grid,header);
+	NerscGrid(grid,header);
-	assert(header.nd==4);
 	header.link_trace=0.0;
 	header.plaquette=0.0;
-	MachineCharacteristics(header);
+	NerscMachineCharacteristics(header);
 
+	uint32_t csum;
 	int offset;
   
 #ifdef RNG_RANLUX
@@ -295,19 +502,15 @@ namespace Grid {
 
 	truncate(file);
 	offset = writeHeader(header,file);
-	uint32_t nersc_csum,scidac_csuma,scidac_csumb;
+	csum=BinaryIO::writeRNGSerial(serial,parallel,file,offset);
-	BinaryIO::writeRNG(serial,parallel,file,offset,nersc_csum,scidac_csuma,scidac_csumb);
+	header.checksum = csum;
-	header.checksum = nersc_csum;
 	offset = writeHeader(header,file);
 
-	std::cout<<GridLogMessage 
+	std::cout<<GridLogMessage <<"Written NERSC RNG STATE "<<file<< " checksum "<<std::hex<<csum<<std::dec<<std::endl;
-		 <<"Written NERSC RNG STATE "<<file<< " checksum "
-		 <<std::hex<<header.checksum
-		 <<std::dec<<std::endl;
 
       }
     
-      static inline void readRNGState(GridSerialRNG &serial,GridParallelRNG & parallel,FieldMetaData& header,std::string file)
+      static inline void readRNGState(GridSerialRNG &serial,GridParallelRNG & parallel,NerscField& header,std::string file)
       {
 	typedef typename GridParallelRNG::RngStateType RngStateType;
 
@@ -315,7 +518,7 @@ namespace Grid {
 
 	int offset = readHeader(file,grid,header);
 
-	FieldMetaData clone(header);
+	NerscField clone(header);
 
 	std::string format(header.floating_point);
 	std::string data_type(header.data_type);
@@ -335,19 +538,15 @@ namespace Grid {
 
 	// depending on datatype, set up munger;
 	// munger is a function of <floating point, Real, data_type>
-	uint32_t nersc_csum,scidac_csuma,scidac_csumb;
+	uint32_t csum=BinaryIO::readRNGSerial(serial,parallel,file,offset);
-	BinaryIO::readRNG(serial,parallel,file,offset,nersc_csum,scidac_csuma,scidac_csumb);
 
-	if ( nersc_csum != header.checksum ) { 
+	assert(csum == header.checksum );
-	  std::cerr << "checksum mismatch "<<std::hex<< nersc_csum <<" "<<header.checksum<<std::dec<<std::endl;
-	  exit(0);
-	}
-	assert(nersc_csum == header.checksum );
 
 	std::cout<<GridLogMessage <<"Read NERSC RNG file "<<file<< " format "<< data_type <<std::endl;
       }
 
     };
 
+
   }}
 #endif

lib/qcd/action/fermion/FermionOperatorImpl.h

@@ -644,16 +644,19 @@ class StaggeredImpl : public PeriodicGaugeImpl<GaugeImplTypes<S, Representation:
 
     INHERIT_GIMPL_TYPES(Gimpl);
       
+    template <typename vtype> using iImplScalar            = iScalar<iScalar<iScalar<vtype> > >;
     template <typename vtype> using iImplSpinor            = iScalar<iScalar<iVector<vtype, Dimension> > >;
     template <typename vtype> using iImplHalfSpinor        = iScalar<iScalar<iVector<vtype, Dimension> > >;
     template <typename vtype> using iImplDoubledGaugeField = iVector<iScalar<iMatrix<vtype, Dimension> >, Nds>;
     template <typename vtype> using iImplPropagator        = iScalar<iScalar<iMatrix<vtype, Dimension> > >;
     
+    typedef iImplScalar<Simd>            SiteComplex;
     typedef iImplSpinor<Simd>            SiteSpinor;
     typedef iImplHalfSpinor<Simd>        SiteHalfSpinor;
     typedef iImplDoubledGaugeField<Simd> SiteDoubledGaugeField;
     typedef iImplPropagator<Simd>        SitePropagator;
     
+    typedef Lattice<SiteComplex>           ComplexField;
     typedef Lattice<SiteSpinor>            FermionField;
     typedef Lattice<SiteDoubledGaugeField> DoubledGaugeField;
     typedef Lattice<SitePropagator> PropagatorField;
@@ -772,6 +775,7 @@ class StaggeredImpl : public PeriodicGaugeImpl<GaugeImplTypes<S, Representation:
 
     INHERIT_GIMPL_TYPES(Gimpl);
 
+    template <typename vtype> using iImplScalar            = iScalar<iScalar<iScalar<vtype> > >;
     template <typename vtype> using iImplSpinor            = iScalar<iScalar<iVector<vtype, Dimension> > >;
     template <typename vtype> using iImplHalfSpinor        = iScalar<iScalar<iVector<vtype, Dimension> > >;
     template <typename vtype> using iImplDoubledGaugeField = iVector<iScalar<iMatrix<vtype, Dimension> >, Nds>;
@@ -788,10 +792,12 @@ class StaggeredImpl : public PeriodicGaugeImpl<GaugeImplTypes<S, Representation:
     typedef Lattice<SiteDoubledGaugeField> DoubledGaugeField;
     typedef Lattice<SitePropagator> PropagatorField;
     
+    typedef iImplScalar<Simd>            SiteComplex;
     typedef iImplSpinor<Simd>            SiteSpinor;
     typedef iImplHalfSpinor<Simd>        SiteHalfSpinor;
 
     
+    typedef Lattice<SiteComplex>           ComplexField;
     typedef Lattice<SiteSpinor>            FermionField;
     
     typedef SimpleCompressor<SiteSpinor> Compressor;

lib/qcd/action/gauge/GaugeImplTypes.h

@@ -40,15 +40,12 @@ namespace QCD {
   typedef typename GImpl::Simd Simd;                \
   typedef typename GImpl::LinkField GaugeLinkField; \
   typedef typename GImpl::Field GaugeField;         \
-  typedef typename GImpl::ComplexField ComplexField;\
   typedef typename GImpl::SiteField SiteGaugeField; \
-  typedef typename GImpl::SiteComplex SiteComplex;  \
   typedef typename GImpl::SiteLink SiteGaugeLink;
 
-#define INHERIT_FIELD_TYPES(Impl)		    \
+#define INHERIT_FIELD_TYPES(Impl)             \
-  typedef typename Impl::Simd Simd;		    \
+  typedef typename Impl::Simd Simd;           \
-  typedef typename Impl::ComplexField ComplexField; \
+  typedef typename Impl::SiteField SiteField; \
-  typedef typename Impl::SiteField SiteField;	    \
   typedef typename Impl::Field Field;
 
 // hardcodes the exponential approximation in the template
@@ -56,15 +53,12 @@ template <class S, int Nrepresentation = Nc, int Nexp = 12 > class GaugeImplType
 public:
   typedef S Simd;
 
-  template <typename vtype> using iImplScalar     = iScalar<iScalar<iScalar<vtype> > >;
+  template <typename vtype> using iImplGaugeLink  = iScalar<iScalar<iMatrix<vtype, Nrepresentation>>>;
-  template <typename vtype> using iImplGaugeLink  = iScalar<iScalar<iMatrix<vtype, Nrepresentation> > >;
+  template <typename vtype> using iImplGaugeField = iVector<iScalar<iMatrix<vtype, Nrepresentation>>, Nd>;
-  template <typename vtype> using iImplGaugeField = iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nd>;
 
-  typedef iImplScalar<Simd>     SiteComplex;
   typedef iImplGaugeLink<Simd>  SiteLink;
   typedef iImplGaugeField<Simd> SiteField;
 
-  typedef Lattice<SiteComplex> ComplexField;
   typedef Lattice<SiteLink>  LinkField; 
   typedef Lattice<SiteField> Field;
 

lib/qcd/hmc/checkpointers/BinaryCheckpointer.h

@@ -62,50 +62,36 @@ class BinaryHmcCheckpointer : public BaseHmcCheckpointer<Impl> {
     fout.close();
   }
 
-  void TrajectoryComplete(int traj, Field &U, GridSerialRNG &sRNG, GridParallelRNG &pRNG) {
+  void TrajectoryComplete(int traj, Field &U, GridSerialRNG &sRNG,
-
+                          GridParallelRNG &pRNG) {
     if ((traj % Params.saveInterval) == 0) {
       std::string config, rng;
       this->build_filenames(traj, Params, config, rng);
 
-      uint32_t nersc_csum;
+      BinaryIO::BinarySimpleUnmunger<sobj_double, sobj> munge;
-      uint32_t scidac_csuma;
-      uint32_t scidac_csumb;
-      
-      BinarySimpleUnmunger<sobj_double, sobj> munge;
       truncate(rng);
-      BinaryIO::writeRNG(sRNG, pRNG, rng, 0,nersc_csum,scidac_csuma,scidac_csumb);
+      BinaryIO::writeRNGSerial(sRNG, pRNG, rng, 0);
       truncate(config);
-
+      uint32_t csum = BinaryIO::writeObjectParallel<vobj, sobj_double>(
-      BinaryIO::writeLatticeObject<vobj, sobj_double>(U, config, munge, 0, Params.format,
+          U, config, munge, 0, Params.format);
-						      nersc_csum,scidac_csuma,scidac_csumb);
 
       std::cout << GridLogMessage << "Written Binary Configuration " << config
-                << " checksum " << std::hex 
+                << " checksum " << std::hex << csum << std::dec << std::endl;
-		<< nersc_csum   <<"/"
-		<< scidac_csuma   <<"/"
-		<< scidac_csumb 
-		<< std::dec << std::endl;
     }
-
   };
 
-  void CheckpointRestore(int traj, Field &U, GridSerialRNG &sRNG, GridParallelRNG &pRNG) {
+  void CheckpointRestore(int traj, Field &U, GridSerialRNG &sRNG,
+                         GridParallelRNG &pRNG) {
     std::string config, rng;
     this->build_filenames(traj, Params, config, rng);
 
-    BinarySimpleMunger<sobj_double, sobj> munge;
+    BinaryIO::BinarySimpleMunger<sobj_double, sobj> munge;
+    BinaryIO::readRNGSerial(sRNG, pRNG, rng, 0);
+    uint32_t csum = BinaryIO::readObjectParallel<vobj, sobj_double>(
+        U, config, munge, 0, Params.format);
 
-    uint32_t nersc_csum;
-    uint32_t scidac_csuma;
-    uint32_t scidac_csumb;
-    BinaryIO::readRNG(sRNG, pRNG, rng, 0,nersc_csum,scidac_csuma,scidac_csumb);
-    BinaryIO::readLatticeObject<vobj, sobj_double>(U, config, munge, 0, Params.format,
-						   nersc_csum,scidac_csuma,scidac_csumb);
-    
     std::cout << GridLogMessage << "Read Binary Configuration " << config
-              << " checksums " << std::hex << nersc_csum<<"/"<<scidac_csuma<<"/"<<scidac_csumb 
+              << " checksum " << std::hex << csum << std::dec << std::endl;
-	      << std::dec << std::endl;
   };
 };
 }

lib/qcd/hmc/checkpointers/ILDGCheckpointer.h

@@ -54,9 +54,9 @@ class ILDGHmcCheckpointer : public BaseHmcCheckpointer<Implementation> {
 
     // check here that the format is valid
     int ieee32big = (Params.format == std::string("IEEE32BIG"));
-    int ieee32    = (Params.format == std::string("IEEE32"));
+    int ieee32 = (Params.format == std::string("IEEE32"));
     int ieee64big = (Params.format == std::string("IEEE64BIG"));
-    int ieee64    = (Params.format == std::string("IEEE64"));
+    int ieee64 = (Params.format == std::string("IEEE64"));
 
     if (!(ieee64big || ieee32 || ieee32big || ieee64)) {
       std::cout << GridLogError << "Unrecognized file format " << Params.format
@@ -74,20 +74,13 @@ class ILDGHmcCheckpointer : public BaseHmcCheckpointer<Implementation> {
     if ((traj % Params.saveInterval) == 0) {
       std::string config, rng;
       this->build_filenames(traj, Params, config, rng);
-      
+
-      uint32_t nersc_csum,scidac_csuma,scidac_csumb;
+      ILDGIO IO(config, ILDGwrite);
-      BinaryIO::writeRNG(sRNG, pRNG, rng, 0,nersc_csum,scidac_csuma,scidac_csumb);
+      BinaryIO::writeRNGSerial(sRNG, pRNG, rng, 0);
-      IldgWriter _IldgWriter;
+      uint32_t csum = IO.writeConfiguration(U, Params.format);
-      _IldgWriter.open(config);
-      _IldgWriter.writeConfiguration(U, traj, config, config);
-      _IldgWriter.close();
 
       std::cout << GridLogMessage << "Written ILDG Configuration on " << config
-                << " checksum " << std::hex 
+                << " checksum " << std::hex << csum << std::dec << std::endl;
-		<< nersc_csum<<"/"
-		<< scidac_csuma<<"/"
-		<< scidac_csumb
-		<< std::dec << std::endl;
     }
   };
 
@@ -96,21 +89,12 @@ class ILDGHmcCheckpointer : public BaseHmcCheckpointer<Implementation> {
     std::string config, rng;
     this->build_filenames(traj, Params, config, rng);
 
-    uint32_t nersc_csum,scidac_csuma,scidac_csumb;
+    ILDGIO IO(config, ILDGread);
-    BinaryIO::readRNG(sRNG, pRNG, rng, 0,nersc_csum,scidac_csuma,scidac_csumb);
+    BinaryIO::readRNGSerial(sRNG, pRNG, rng, 0);
-
+    uint32_t csum = IO.readConfiguration(U);  // format from the header
-    FieldMetaData header;
-    IldgReader _IldgReader;
-    _IldgReader.open(config);
-    _IldgReader.readConfiguration(config,U,header);  // format from the header
-    _IldgReader.close();
 
     std::cout << GridLogMessage << "Read ILDG Configuration from " << config
-              << " checksum " << std::hex 
+              << " checksum " << std::hex << csum << std::dec << std::endl;
-	      << nersc_csum<<"/"
-	      << scidac_csuma<<"/"
-	      << scidac_csumb
-	      << std::dec << std::endl;
   };
 };
 }

lib/qcd/hmc/checkpointers/NerscCheckpointer.h

@@ -70,7 +70,7 @@ class NerscHmcCheckpointer : public BaseHmcCheckpointer<Gimpl> {
     std::string config, rng;
     this->build_filenames(traj, Params, config, rng);
 
-    FieldMetaData header;
+    NerscField header;
     NerscIO::readRNGState(sRNG, pRNG, header, rng);
     NerscIO::readConfiguration(U, header, config);
   };

lib/qcd/utils/Utils.h

@@ -12,4 +12,7 @@
 #include <Grid/qcd/utils/SUnAdjoint.h>
 #include <Grid/qcd/utils/SUnTwoIndex.h>
 
+
+
+
 #endif

lib/qcd/utils/WilsonLoops.h

@@ -73,7 +73,7 @@ public:
   //////////////////////////////////////////////////
   // trace of directed plaquette oriented in mu,nu plane
   //////////////////////////////////////////////////
-  static void traceDirPlaquette(ComplexField &plaq,
+  static void traceDirPlaquette(LatticeComplex &plaq,
                                 const std::vector<GaugeMat> &U, const int mu,
                                 const int nu) {
     GaugeMat sp(U[0]._grid);
@@ -83,9 +83,9 @@ public:
   //////////////////////////////////////////////////
   // sum over all planes of plaquette
   //////////////////////////////////////////////////
-  static void sitePlaquette(ComplexField &Plaq,
+  static void sitePlaquette(LatticeComplex &Plaq,
                             const std::vector<GaugeMat> &U) {
-    ComplexField sitePlaq(U[0]._grid);
+    LatticeComplex sitePlaq(U[0]._grid);
     Plaq = zero;
     for (int mu = 1; mu < Nd; mu++) {
       for (int nu = 0; nu < mu; nu++) {
@@ -104,11 +104,11 @@ public:
       U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
     }
 
-    ComplexField Plaq(Umu._grid);
+    LatticeComplex Plaq(Umu._grid);
 
     sitePlaquette(Plaq, U);
-    auto Tp = sum(Plaq);
+    TComplex Tp = sum(Plaq);
-    auto p = TensorRemove(Tp);
+    Complex p = TensorRemove(Tp);
     return p.real();
   }
 
@@ -129,15 +129,15 @@ public:
   static RealD linkTrace(const GaugeLorentz &Umu) {
     std::vector<GaugeMat> U(Nd, Umu._grid);
 
-    ComplexField Tr(Umu._grid);
+    LatticeComplex Tr(Umu._grid);
     Tr = zero;
     for (int mu = 0; mu < Nd; mu++) {
       U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
       Tr = Tr + trace(U[mu]);
     }
 
-    auto Tp = sum(Tr);
+    TComplex Tp = sum(Tr);
-    auto p = TensorRemove(Tp);
+    Complex p = TensorRemove(Tp);
 
     double vol = Umu._grid->gSites();
 
@@ -330,8 +330,8 @@ public:
 
     double coeff = 8.0/(32.0*M_PI*M_PI);
 
-    ComplexField qfield = coeff*trace(Bx*Ex + By*Ey + Bz*Ez);
+    LatticeComplex qfield = coeff*trace(Bx*Ex + By*Ey + Bz*Ez);
-    auto Tq = sum(qfield);
+    TComplex Tq = sum(qfield);
     return TensorRemove(Tq).real();
   }
 
@@ -350,16 +350,16 @@ public:
                adj(Gimpl::CovShiftForward(
                    U[nu], nu, Gimpl::CovShiftForward(U[nu], nu, U[mu])));
   }
-  static void traceDirRectangle(ComplexField &rect,
+  static void traceDirRectangle(LatticeComplex &rect,
                                 const std::vector<GaugeMat> &U, const int mu,
                                 const int nu) {
     GaugeMat sp(U[0]._grid);
     dirRectangle(sp, U, mu, nu);
     rect = trace(sp);
   }
-  static void siteRectangle(ComplexField &Rect,
+  static void siteRectangle(LatticeComplex &Rect,
                             const std::vector<GaugeMat> &U) {
-    ComplexField siteRect(U[0]._grid);
+    LatticeComplex siteRect(U[0]._grid);
     Rect = zero;
     for (int mu = 1; mu < Nd; mu++) {
       for (int nu = 0; nu < mu; nu++) {
@@ -379,12 +379,12 @@ public:
       U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
     }
 
-    ComplexField Rect(Umu._grid);
+    LatticeComplex Rect(Umu._grid);
 
     siteRectangle(Rect, U);
 
-    auto Tp = sum(Rect);
+    TComplex Tp = sum(Rect);
-    auto p = TensorRemove(Tp);
+    Complex p = TensorRemove(Tp);
     return p.real();
   }
   //////////////////////////////////////////////////

lib/serialisation/MacroMagic.h

@@ -110,12 +110,11 @@ THE SOFTWARE.
 
 #define GRID_MACRO_MEMBER(A,B)        A B;
 #define GRID_MACRO_COMP_MEMBER(A,B) result = (result and (lhs. B == rhs. B));
-#define GRID_MACRO_OS_WRITE_MEMBER(A,B) os<< #A <<" " #B << " = " << obj. B << " ; " <<std::endl;
+#define GRID_MACRO_OS_WRITE_MEMBER(A,B) os<< #A <<" "#B <<" = "<< obj. B <<" ; " <<std::endl;
 #define GRID_MACRO_READ_MEMBER(A,B) Grid::read(RD,#B,obj. B);
 #define GRID_MACRO_WRITE_MEMBER(A,B) Grid::write(WR,#B,obj. B);
 
 #define GRID_SERIALIZABLE_CLASS_MEMBERS(cname,...)\
-  std::string SerialisableClassName(void) {return std::string(#cname);}	\
 GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_MEMBER,__VA_ARGS__))\
 template <typename T>\
 static inline void write(Writer<T> &WR,const std::string &s, const cname &obj){ \

lib/serialisation/XmlIO.cc

@@ -32,21 +32,16 @@ using namespace Grid;
 using namespace std;
 
 // Writer implementation ///////////////////////////////////////////////////////
-XmlWriter::XmlWriter(const string &fileName, string toplev) : fileName_(fileName)
+XmlWriter::XmlWriter(const string &fileName)
+: fileName_(fileName)
 {
-  if ( toplev == std::string("") ) {
+  node_ = doc_.append_child();
-    node_=doc_;
+  node_.set_name("grid");
-  } else { 
-    node_=doc_.append_child();
-    node_.set_name(toplev.c_str());
-  }
 }
 
 XmlWriter::~XmlWriter(void)
 {
-  if ( fileName_ != std::string("") ) { 
+  doc_.save_file(fileName_.c_str(), "  ");
-    doc_.save_file(fileName_.c_str(), "  ");
-  }
 }
 
 void XmlWriter::push(const string &s)
@@ -58,44 +53,21 @@ void XmlWriter::pop(void)
 {
   node_ = node_.parent();
 }
-std::string XmlWriter::XmlString(void)
-{
-  std::ostringstream oss; 
-  doc_.save(oss);
-  return oss.str();
-}
-
-XmlReader::XmlReader(const char *xmlstring,string toplev) : fileName_("")
-{
-  pugi::xml_parse_result result;
-  result = doc_.load_string(xmlstring);
-  if ( !result ) {
-    cerr << "XML error description: " << result.description() << "\n";
-    cerr << "XML error offset     : " << result.offset        << "\n";
-    abort();
-  }
-  if ( toplev == std::string("") ) {
-    node_ = doc_;
-  } else { 
-    node_ = doc_.child(toplev.c_str());
-  }
-}
 
 // Reader implementation ///////////////////////////////////////////////////////
-XmlReader::XmlReader(const string &fileName,string toplev) : fileName_(fileName)
+XmlReader::XmlReader(const string &fileName)
+: fileName_(fileName)
 {
-  pugi::xml_parse_result result;
+  pugi::xml_parse_result result = doc_.load_file(fileName_.c_str());
-  result = doc_.load_file(fileName_.c_str());
+  
-  if ( !result ) {
+  if ( !result )
+  {
     cerr << "XML error description: " << result.description() << "\n";
     cerr << "XML error offset     : " << result.offset        << "\n";
     abort();
   }
-  if ( toplev == std::string("") ) {
+  
-    node_ = doc_;
+  node_ = doc_.child("grid");
-  } else { 
-    node_ = doc_.child(toplev.c_str());
-  }
 }
 
 bool XmlReader::push(const string &s)

lib/serialisation/XmlIO.h

@@ -44,9 +44,10 @@ namespace Grid
 {
   
   class XmlWriter: public Writer<XmlWriter>
-  {    
+  {
+    
   public:
-    XmlWriter(const std::string &fileName,std::string toplev = std::string("grid") );
+    XmlWriter(const std::string &fileName);
     virtual ~XmlWriter(void);
     void push(const std::string &s);
     void pop(void);
@@ -54,7 +55,6 @@ namespace Grid
     void writeDefault(const std::string &s, const U &x);
     template <typename U>
     void writeDefault(const std::string &s, const std::vector<U> &x);
-    std::string XmlString(void);
   private:
     pugi::xml_document doc_;
     pugi::xml_node     node_;
@@ -64,8 +64,7 @@ namespace Grid
   class XmlReader: public Reader<XmlReader>
   {
   public:
-    XmlReader(const char *xmlstring,std::string toplev = std::string("grid") );
+    XmlReader(const std::string &fileName);
-    XmlReader(const std::string &fileName,std::string toplev = std::string("grid") );
     virtual ~XmlReader(void) = default;
     bool push(const std::string &s);
     void pop(void);
@@ -119,7 +118,7 @@ namespace Grid
     std::string buf;
     
     readDefault(s, buf);
-    //    std::cout << s << "   " << buf << std::endl;
+    std::cout << s << "   " << buf << std::endl;
     fromString(output, buf);
   }
   

lib/simd/Grid_generic.h

@@ -281,8 +281,8 @@ namespace Optimization {
 
   struct PrecisionChange {
     static inline vech StoH (const vecf &a,const vecf &b) {
-      vech ret;
 #ifdef USE_FP16
+      vech ret;
       vech *ha = (vech *)&a;
       vech *hb = (vech *)&b;
       const int nf = W<float>::r;
@@ -493,8 +493,6 @@ namespace Optimization {
     
     return a;
   }
-
-  #undef acc  // EIGEN compatibility
 }
 
 //////////////////////////////////////////////////////////////////////////////////////

lib/simd/Grid_vector_types.h

@@ -327,16 +327,18 @@ class Grid_simd {
   // provides support
   ///////////////////////////////////////
 
+#if (__GNUC__ == 5 ) || ( ( __GNUC__ == 6 ) && __GNUC_MINOR__ < 3 )
+#pragma GCC push_options 
+#pragma GCC optimize ("O0") 
+#endif
   template <class functor>
   friend inline Grid_simd SimdApply(const functor &func, const Grid_simd &v) {
     Grid_simd ret;
     Grid_simd::conv_t conv;
-    Grid_simd::scalar_type s;
+
-    
     conv.v = v.v;
     for (int i = 0; i < Nsimd(); i++) {
-      s = conv.s[i];
+      conv.s[i] = func(conv.s[i]);
-      conv.s[i] = func(s);
     }
     ret.v = conv.v;
     return ret;
@@ -348,18 +350,18 @@ class Grid_simd {
     Grid_simd ret;
     Grid_simd::conv_t cx;
     Grid_simd::conv_t cy;
-    Grid_simd::scalar_type sx,sy;
 
     cx.v = x.v;
     cy.v = y.v;
     for (int i = 0; i < Nsimd(); i++) {
-      sx = cx.s[i];
+      cx.s[i] = func(cx.s[i], cy.s[i]);
-      sy = cy.s[i];
-      cx.s[i] = func(sx,sy);
     }
     ret.v = cx.v;
     return ret;
   }
+#if (__GNUC__ == 5 ) || ( ( __GNUC__ == 6 ) && __GNUC_MINOR__ < 3 )
+#pragma GCC pop_options
+#endif
   ///////////////////////
   // Exchange 
   // Al Ah , Bl Bh -> Al Bl Ah,Bh
@@ -421,6 +423,7 @@ class Grid_simd {
   
 };  // end of Grid_simd class definition
 
+
 inline void permute(ComplexD &y,ComplexD b, int perm) {  y=b; }
 inline void permute(ComplexF &y,ComplexF b, int perm) {  y=b; }
 inline void permute(RealD &y,RealD b, int perm) {  y=b; }
@@ -830,6 +833,8 @@ inline void precisionChange(vComplexD *out,vComplexF *in,int nvec){ precisionCha
 inline void precisionChange(vComplexD *out,vComplexH *in,int nvec){ precisionChange((vRealD *)out,(vRealH *)in,nvec);}
 inline void precisionChange(vComplexF *out,vComplexH *in,int nvec){ precisionChange((vRealF *)out,(vRealH *)in,nvec);}
 
+
+
 // Check our vector types are of an appropriate size.
 #if defined QPX
 static_assert(2*sizeof(SIMD_Ftype) == sizeof(SIMD_Dtype), "SIMD vector lengths incorrect");
@@ -844,14 +849,21 @@ static_assert(sizeof(SIMD_Ftype) == sizeof(SIMD_Itype), "SIMD vector lengths inc
 /////////////////////////////////////////
 template <typename T>
 struct is_simd : public std::false_type {};
-template <> struct is_simd<vRealF>     : public std::true_type {};
+template <>
-template <> struct is_simd<vRealD>     : public std::true_type {};
+struct is_simd<vRealF> : public std::true_type {};
-template <> struct is_simd<vComplexF>  : public std::true_type {};
+template <>
-template <> struct is_simd<vComplexD>  : public std::true_type {};
+struct is_simd<vRealD> : public std::true_type {};
-template <> struct is_simd<vInteger>   : public std::true_type {};
+template <>
+struct is_simd<vComplexF> : public std::true_type {};
+template <>
+struct is_simd<vComplexD> : public std::true_type {};
+template <>
+struct is_simd<vInteger> : public std::true_type {};
 
-template <typename T> using IfSimd    = Invoke<std::enable_if<is_simd<T>::value, int> >;
+template <typename T>
-template <typename T> using IfNotSimd = Invoke<std::enable_if<!is_simd<T>::value, unsigned> >;
+using IfSimd = Invoke<std::enable_if<is_simd<T>::value, int> >;
+template <typename T>
+using IfNotSimd = Invoke<std::enable_if<!is_simd<T>::value, unsigned> >;
 }
 
 #endif

lib/simd/Grid_vector_unops.h

@@ -179,6 +179,13 @@ inline Grid_simd<S, V> div(const Grid_simd<S, V> &r, Integer y) {
 ////////////////////////////////////////////////////////////////////////////
 // Allows us to assign into **conformable** real vectors from complex
 ////////////////////////////////////////////////////////////////////////////
+//  template < class S, class V >
+//  inline auto ComplexRemove(const Grid_simd<S,V> &c) ->
+//  Grid_simd<Grid_simd<S,V>::Real,V> {
+//    Grid_simd<Grid_simd<S,V>::Real,V> ret;
+//    ret.v = c.v;
+//    return ret;
+//  }
 template <class scalar>
 struct AndFunctor {
   scalar operator()(const scalar &x, const scalar &y) const { return x & y; }

lib/tensors/Tensor_index.h

@@ -47,28 +47,6 @@ template<int Level>
 class TensorIndexRecursion {
 
  public:
-
-  ////////////////////////////////////////////////////
-  // Type Queries
-  ////////////////////////////////////////////////////
-  template<class vtype>       static inline int indexRank(const iScalar<vtype> tmp)  { return TensorIndexRecursion<Level-1>::indexRank(tmp._internal);  }
-  template<class vtype,int N> static inline int indexRank(const iVector<vtype,N> tmp){ return TensorIndexRecursion<Level-1>::indexRank(tmp._internal[0]);  }
-  template<class vtype,int N> static inline int indexRank(const iMatrix<vtype,N> tmp){ return TensorIndexRecursion<Level-1>::indexRank(tmp._internal[0][0]);  }
-
-  template<class vtype>       static inline int isScalar(const iScalar<vtype> tmp)  { return TensorIndexRecursion<Level-1>::isScalar(tmp._internal);  }
-  template<class vtype,int N> static inline int isScalar(const iVector<vtype,N> tmp){ return TensorIndexRecursion<Level-1>::isScalar(tmp._internal[0]);  }
-  template<class vtype,int N> static inline int isScalar(const iMatrix<vtype,N> tmp){ return TensorIndexRecursion<Level-1>::isScalar(tmp._internal[0][0]);  }
-
-  template<class vtype>       static inline int isVector(const iScalar<vtype> tmp)  { return TensorIndexRecursion<Level-1>::isVector(tmp._internal);  }
-  template<class vtype,int N> static inline int isVector(const iVector<vtype,N> tmp){ return TensorIndexRecursion<Level-1>::isVector(tmp._internal[0]);  }
-  template<class vtype,int N> static inline int isVector(const iMatrix<vtype,N> tmp){ return TensorIndexRecursion<Level-1>::isVector(tmp._internal[0][0]);  }
-
-  template<class vtype>       static inline int isMatrix(const iScalar<vtype> tmp)  { return TensorIndexRecursion<Level-1>::isMatrix(tmp._internal);  }
-  template<class vtype,int N> static inline int isMatrix(const iVector<vtype,N> tmp){ return TensorIndexRecursion<Level-1>::isMatrix(tmp._internal[0]);  }
-  template<class vtype,int N> static inline int isMatrix(const iMatrix<vtype,N> tmp){ return TensorIndexRecursion<Level-1>::isMatrix(tmp._internal[0][0]);  }
-  ////////////////////////////////////////////////////
-  // Trace
-  ////////////////////////////////////////////////////
   template<class vtype>
   static auto traceIndex(const iScalar<vtype> arg) ->  iScalar<decltype(TensorIndexRecursion<Level-1>::traceIndex(arg._internal))> 
   {
@@ -237,24 +215,6 @@ class TensorIndexRecursion {
 template<>
 class TensorIndexRecursion<0> {
  public:
-  ////////////////////////////////////////////////////
-  // Type Queries
-  ////////////////////////////////////////////////////
-  template<class vtype>       static inline int indexRank(const iScalar<vtype> tmp)  { return 1; }
-  template<class vtype,int N> static inline int indexRank(const iVector<vtype,N> tmp){ return N; }
-  template<class vtype,int N> static inline int indexRank(const iMatrix<vtype,N> tmp){ return N; }
-
-  template<class vtype>       static inline int isScalar(const iScalar<vtype> tmp)  { return true;}
-  template<class vtype,int N> static inline int isScalar(const iVector<vtype,N> tmp){ return false;}
-  template<class vtype,int N> static inline int isScalar(const iMatrix<vtype,N> tmp){ return false;}
-
-  template<class vtype>       static inline int isVector(const iScalar<vtype> tmp)  { return false;}
-  template<class vtype,int N> static inline int isVector(const iVector<vtype,N> tmp){ return true;}
-  template<class vtype,int N> static inline int isVector(const iMatrix<vtype,N> tmp){ return false;}
-
-  template<class vtype>       static inline int isMatrix(const iScalar<vtype> tmp)  { return false;}
-  template<class vtype,int N> static inline int isMatrix(const iVector<vtype,N> tmp){ return false;}
-  template<class vtype,int N> static inline int isMatrix(const iMatrix<vtype,N> tmp){ return true;}
 
   /////////////////////////////////////////
   // Ends recursion for trace (scalar/vector/matrix)
@@ -342,26 +302,6 @@ class TensorIndexRecursion<0> {
 ////////////////////////////////////////////////////////////////////////////////////////////////////////
 // External wrappers
 ////////////////////////////////////////////////////////////////////////////////////////////////////////
-template<int Level,class vtype> inline int indexRank(void)
-{
-  vtype tmp;
-  return TensorIndexRecursion<Level>::indexRank(tmp);
-}
-template<int Level,class vtype> inline int isScalar(void)
-{
-  vtype tmp;
-  return TensorIndexRecursion<Level>::isScalar(tmp);
-}
-template<int Level,class vtype> inline int isVector(void)
-{
-  vtype tmp;
-  return TensorIndexRecursion<Level>::isVector(tmp);
-}
-template<int Level,class vtype> inline int isMatrix(void)
-{
-  vtype tmp;
-  return TensorIndexRecursion<Level>::isMatrix(tmp);
-}
 
 template<int Level,class vtype> inline auto traceIndex (const vtype &arg) -> RemoveCRV(TensorIndexRecursion<Level>::traceIndex(arg))
 {

lib/tensors/Tensor_traits.h

@@ -281,8 +281,8 @@ namespace Grid {
   template<typename T>
   class getPrecision{
   public:
-    //get the vector_obj (i.e. a grid Tensor) if its a Lattice<vobj>, do nothing otherwise (i.e. if fundamental or grid Tensor)
+    typedef typename getVectorType<T>::type vector_obj; //get the vector_obj (i.e. a grid Tensor) if its a Lattice<vobj>, do nothing otherwise (i.e. if fundamental or grid Tensor)
-    typedef typename getVectorType<T>::type vector_obj; 
+  
     typedef typename GridTypeMapper<vector_obj>::scalar_type scalar_type; //get the associated scalar type. Works on fundamental and tensor types
     typedef typename GridTypeMapper<scalar_type>::Realified real_scalar_type; //remove any std::complex wrapper, should get us to the fundamental type
 

tests/IO/Test_ildg_io.cc

@@ -1,99 +0,0 @@
-    /*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./tests/Test_nersc_io.cc
-
-    Copyright (C) 2015
-
-Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-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 */
-#include <Grid/Grid.h>
-
-using namespace std;
-using namespace Grid;
-using namespace Grid::QCD;
-
-
-int main (int argc, char ** argv)
-{
-  Grid_init(&argc,&argv);
-
-  std::cout <<GridLogMessage<< " main "<<std::endl;
-
-  std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
-  std::vector<int> mpi_layout  = GridDefaultMpi();
-  //std::vector<int> latt_size  ({48,48,48,96});
-  //std::vector<int> latt_size  ({32,32,32,32});
-  std::vector<int> latt_size  ({16,16,16,32});
-  std::vector<int> clatt_size  ({4,4,4,8});
-  int orthodir=3;
-  int orthosz =latt_size[orthodir];
-    
-  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;
-
-  std::cout <<GridLogMessage<< " seeding... "<<std::endl;
-  pRNGa.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
-  sRNGa.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
-  std::cout <<GridLogMessage<< " ...done "<<std::endl;
-
-  LatticeGaugeField Umu(&Fine);
-  LatticeGaugeField Umu_diff(&Fine);
-  LatticeGaugeField Umu_saved(&Fine);
-
-  std::vector<LatticeColourMatrix> U(4,&Fine);
-  
-  SU3::HotConfiguration(pRNGa,Umu);
-
-
-  FieldMetaData header;
-
-  std::cout <<GridLogMessage<<"**************************************"<<std::endl;
-  std::cout <<GridLogMessage<<"** Writing out  ILDG conf    *********"<<std::endl;
-  std::cout <<GridLogMessage<<"**************************************"<<std::endl;
-  std::string file("./ckpoint_ildg.4000");
-  IldgWriter _IldgWriter;
-  _IldgWriter.open(file);
-  _IldgWriter.writeConfiguration(Umu,4000,std::string("dummy_ildg_LFN"),std::string("dummy_config"));
-  _IldgWriter.close();
-
-  Umu_saved = Umu;
-  std::cout <<GridLogMessage<<"**************************************"<<std::endl;
-  std::cout <<GridLogMessage<<"** Reading back ILDG conf    *********"<<std::endl;
-  std::cout <<GridLogMessage<<"**************************************"<<std::endl;
-  IldgReader _IldgReader;
-  _IldgReader.open(file);
-  _IldgReader.readConfiguration(Umu,header);
-  _IldgReader.close();
-  Umu_diff = Umu - Umu_saved;
-
-  std::cout <<GridLogMessage<< "norm2 Gauge Diff = "<<norm2(Umu_diff)<<std::endl;
-
-  Grid_finalize();
-}

tests/IO/Test_ildg_read.cc

@@ -1,115 +0,0 @@
-    /*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./tests/Test_nersc_io.cc
-
-    Copyright (C) 2015
-
-Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-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 */
-#include <Grid/Grid.h>
-
-using namespace std;
-using namespace Grid;
-using namespace Grid::QCD;
-
-
-int main (int argc, char ** argv)
-{
-  Grid_init(&argc,&argv);
-
-
-  std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
-  std::vector<int> mpi_layout  = GridDefaultMpi();
-  std::vector<int> latt_size = GridDefaultLatt();
-  int orthodir=3;
-  int orthosz =latt_size[orthodir];
-    
-  GridCartesian     Fine(latt_size,simd_layout,mpi_layout);
-
-  LatticeGaugeField Umu(&Fine);
-  std::vector<LatticeColourMatrix> U(4,&Fine);
-  
-  FieldMetaData header;
-  std::string file("./ildg.file");
-  IldgReader IR;
-  IR.open(file);
-  IR.readConfiguration(Umu,header);
-  IR.close();
-
-  for(int mu=0;mu<Nd;mu++){
-    U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
-  }
-
-  // Painful ; fix syntactical niceness
-  LatticeComplex LinkTrace(&Fine);
-  LinkTrace=zero;
-  for(int mu=0;mu<Nd;mu++){
-    LinkTrace = LinkTrace + trace(U[mu]);
-  }
-
-  // (1+2+3)=6 = N(N-1)/2 terms
-  LatticeComplex Plaq(&Fine);
-
-  Plaq = zero;
-
-  for(int mu=1;mu<Nd;mu++){
-    for(int nu=0;nu<mu;nu++){
-      Plaq = Plaq + trace(U[mu]*Cshift(U[nu],mu,1)*adj(Cshift(U[mu],nu,1))*adj(U[nu]));
-    }
-  }
-
-  double vol = Fine.gSites();
-  Complex PlaqScale(1.0/vol/6.0/3.0);
-  std::cout<<GridLogMessage <<"PlaqScale" << PlaqScale<<std::endl;
-
-  std::vector<TComplex> Plaq_T(orthosz);
-  sliceSum(Plaq,Plaq_T,Nd-1);
-  int Nt = Plaq_T.size();
-
-  TComplex Plaq_T_sum; 
-  Plaq_T_sum=zero;
-  for(int t=0;t<Nt;t++){
-    Plaq_T_sum = Plaq_T_sum+Plaq_T[t];
-    Complex Pt=TensorRemove(Plaq_T[t]);
-    std::cout<<GridLogMessage << "sliced ["<<t<<"]" <<Pt*PlaqScale*Real(Nt)<<std::endl;
-  }
-
-  {
-    Complex Pt = TensorRemove(Plaq_T_sum);
-    std::cout<<GridLogMessage << "total " <<Pt*PlaqScale<<std::endl;
-  }  
-
-
-  TComplex Tp = sum(Plaq);
-  Complex p  = TensorRemove(Tp);
-  std::cout<<GridLogMessage << "calculated plaquettes " <<p*PlaqScale<<std::endl;
-
-
-  Complex LinkTraceScale(1.0/vol/4.0/3.0);
-  TComplex Tl = sum(LinkTrace);
-  Complex l  = TensorRemove(Tl);
-  std::cout<<GridLogMessage << "calculated link trace " <<l*LinkTraceScale<<std::endl;
-
-  Grid_finalize();
-}

tests/IO/Test_nersc_io.cc

@@ -38,13 +38,10 @@ int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  std::cout <<GridLogMessage<< " main "<<std::endl;
 
   std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
   std::vector<int> mpi_layout  = GridDefaultMpi();
-  //std::vector<int> latt_size  ({48,48,48,96});
+  std::vector<int> latt_size  ({16,16,16,16});
-  //std::vector<int> latt_size  ({32,32,32,32});
-  std::vector<int> latt_size  ({16,16,16,32});
   std::vector<int> clatt_size  ({4,4,4,8});
   int orthodir=3;
   int orthosz =latt_size[orthodir];
@@ -52,32 +49,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;
 
-  std::cout <<GridLogMessage<< " seeding... "<<std::endl;
   pRNGa.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
   sRNGa.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
-  std::cout <<GridLogMessage<< " ...done "<<std::endl;
+  
-
   std::string rfile("./ckpoint_rng.4000");
-  FieldMetaData rngheader;
   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 <<GridLogMessage<< " difference between restored randoms and orig "<<norm2( tmpa ) <<" / "<< norm2(tmpb)<<std::endl;
+  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 <<GridLogMessage<< " serial RNG numbers "<<a<<" "<<b<<std::endl;
+  std::cout << " serial RNG numbers "<<a<<" "<<b<<std::endl;
+
 
   LatticeGaugeField Umu(&Fine);
   LatticeGaugeField Umu_diff(&Fine);
@@ -85,20 +80,15 @@ int main (int argc, char ** argv)
 
   std::vector<LatticeColourMatrix> U(4,&Fine);
   
-  SU3::HotConfiguration(pRNGa,Umu);
+  SU3::ColdConfiguration(pRNGa,Umu);
 
-  FieldMetaData header;
+  NerscField header;
   std::string file("./ckpoint_lat.4000");
 
   int precision32 = 0;
   int tworow      = 0;
   NerscIO::writeConfiguration(Umu,file,tworow,precision32);
-  Umu_saved = Umu;
   NerscIO::readConfiguration(Umu,header,file);
-  Umu_diff = Umu - Umu_saved;
-  //std::cout << "Umu_save "<<Umu_saved[0]<<std::endl;
-  //std::cout << "Umu_read "<<Umu[0]<<std::endl;
-  std::cout <<GridLogMessage<< "norm2 Gauge Diff = "<<norm2(Umu_diff)<<std::endl;
 
   for(int mu=0;mu<Nd;mu++){
     U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
@@ -125,6 +115,7 @@ int main (int argc, char ** argv)
 #endif
   double vol = Fine.gSites();
   Complex PlaqScale(1.0/vol/6.0/3.0);
+  std::cout<<GridLogMessage <<"PlaqScale" << PlaqScale<<std::endl;
 
   std::vector<TComplex> Plaq_T(orthosz);
   sliceSum(Plaq,Plaq_T,Nd-1);
@@ -148,6 +139,7 @@ int main (int argc, char ** argv)
   Complex p  = TensorRemove(Tp);
   std::cout<<GridLogMessage << "calculated plaquettes " <<p*PlaqScale<<std::endl;
 
+
   Complex LinkTraceScale(1.0/vol/4.0/3.0);
   TComplex Tl = sum(LinkTrace);
   Complex l  = TensorRemove(Tl);

tests/IO/Test_nersc_read.cc

@@ -50,7 +50,7 @@ int main (int argc, char ** argv)
   LatticeGaugeField Umu(&Fine);
   std::vector<LatticeColourMatrix> U(4,&Fine);
   
-  FieldMetaData header;
+  NerscField header;
   std::string file("./ckpoint_lat");
   NerscIO::readConfiguration(Umu,header,file);
 

tests/IO/Test_serialisation.cc

@@ -31,7 +31,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 
 
 using namespace Grid;
-using namespace Grid::QCD;
 
 GRID_SERIALIZABLE_ENUM(myenum, undef, red, 1, blue, 2, green, 3);
   
@@ -45,8 +44,8 @@ public:
                           double, y,
                           bool , b,
                           std::vector<double>, array,
-                          std::vector<std::vector<double> >, twodimarray,
+                          std::vector<std::vector<double>>, twodimarray,
-                          std::vector<std::vector<std::vector<Complex> > >, cmplx3darray
+                          std::vector<std::vector<std::vector<Complex>>>, cmplx3darray
                           );
   myclass() {}
   myclass(int i)
@@ -238,7 +237,7 @@ int main(int argc,char **argv)
     std::cout << "Loaded (JSON) -----------------" << std::endl;
     std::cout << jcopy1 << std::endl << jveccopy1 << std::endl;
   }
-  
+
 /* 
   // This is still work in progress
   {

tests/Test_dwf_mixedcg_prec_halfcomms.cc

@@ -1,110 +0,0 @@
-    /*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./tests/Test_dwf_cg_prec.cc
-
-    Copyright (C) 2015
-
-Author: Peter Boyle <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 */
-#include <Grid/Grid.h>
-
-using namespace std;
-using namespace Grid;
-using namespace Grid::QCD;
-
-template<class d>
-struct scal {
-  d internal;
-};
-
-  Gamma::Algebra Gmu [] = {
-    Gamma::Algebra::GammaX,
-    Gamma::Algebra::GammaY,
-    Gamma::Algebra::GammaZ,
-    Gamma::Algebra::GammaT
-  };
-
-int main (int argc, char ** argv)
-{
-  Grid_init(&argc,&argv);
-
-  const int Ls=24;
-
-  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexD::Nsimd()),GridDefaultMpi());
-  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
-  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
-  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
-
-  GridCartesian         * UGrid_f   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
-  GridRedBlackCartesian * UrbGrid_f = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid_f);
-  GridCartesian         * FGrid_f   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid_f);
-  GridRedBlackCartesian * FrbGrid_f = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid_f);
-  
-  std::vector<int> seeds4({1,2,3,4});
-  std::vector<int> seeds5({5,6,7,8});
-  GridParallelRNG          RNG5(FGrid);  RNG5.SeedFixedIntegers(seeds5);
-  GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds4);
-
-  LatticeFermionD    src(FGrid); random(RNG5,src);
-  LatticeFermionD result(FGrid); result=zero;
-  LatticeGaugeFieldD Umu(UGrid);
-  LatticeGaugeFieldF Umu_f(UGrid_f); 
-  
-  SU3::HotConfiguration(RNG4,Umu);
-
-  precisionChange(Umu_f,Umu);
-  
-  RealD mass=0.1;
-  RealD M5=1.8;
-  DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
-  DomainWallFermionFH Ddwf_f(Umu_f,*FGrid_f,*FrbGrid_f,*UGrid_f,*UrbGrid_f,mass,M5);
-
-
-  LatticeFermionD    src_o(FrbGrid);
-  LatticeFermionD result_o(FrbGrid);
-  LatticeFermionD result_o_2(FrbGrid);
-  pickCheckerboard(Odd,src_o,src);
-  result_o.checkerboard = Odd;
-  result_o = zero;
-  result_o_2.checkerboard = Odd;
-  result_o_2 = zero;
-
-  SchurDiagMooeeOperator<DomainWallFermionD,LatticeFermionD> HermOpEO(Ddwf);
-  SchurDiagMooeeOperator<DomainWallFermionFH,LatticeFermionF> HermOpEO_f(Ddwf_f);
-
-  std::cout << "Starting mixed CG" << std::endl;
-  MixedPrecisionConjugateGradient<LatticeFermionD,LatticeFermionF> mCG(1.0e-8, 10000, 50, FrbGrid_f, HermOpEO_f, HermOpEO);
-  mCG.InnerTolerance = 3.0e-5;
-  mCG(src_o,result_o);
-
-  std::cout << "Starting regular CG" << std::endl;
-  ConjugateGradient<LatticeFermionD> CG(1.0e-8,10000);
-  CG(HermOpEO,src_o,result_o_2);
-
-  LatticeFermionD diff_o(FrbGrid);
-  RealD diff = axpy_norm(diff_o, -1.0, result_o, result_o_2);
-
-  std::cout << "Diff between mixed and regular CG: " << diff << std::endl;
-
-  
-  Grid_finalize();
-}

tests/forces/Test_contfrac_force.cc

@@ -139,7 +139,7 @@ int main (int argc, char ** argv)
 
   }
 
-  ComplexD dSpred    = sum(dS);
+  Complex dSpred    = sum(dS);
 
   std::cout << GridLogMessage << " S      "<<S<<std::endl;
   std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;

tests/forces/Test_dwf_force.cc

@@ -150,7 +150,7 @@ int main (int argc, char ** argv)
 
   }
 
-  ComplexD dSpred    = sum(dS);
+  Complex dSpred    = sum(dS);
 
   std::cout << GridLogMessage << " S      "<<S<<std::endl;
   std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;

tests/forces/Test_dwf_gpforce.cc

@@ -194,9 +194,9 @@ int main (int argc, char ** argv)
 
   }
 
-  ComplexD dSpred    = sum(dS);
+  Complex dSpred    = sum(dS);
-  ComplexD dSm       = sum(dSmom);
+  Complex dSm       = sum(dSmom);
-  ComplexD dSm2      = sum(dSmom2);
+  Complex dSm2      = sum(dSmom2);
 
 
   std::cout << GridLogMessage <<"Initial mom hamiltonian is "<< Hmom <<std::endl;

tests/forces/Test_gp_rect_force.cc

@@ -113,7 +113,7 @@ int main (int argc, char ** argv)
     dS = dS - trace(mommu*UdSdUmu)*dt*2.0;
 
   }
-  ComplexD dSpred    = sum(dS);
+  Complex dSpred    = sum(dS);
 
   std::cout << GridLogMessage << " S      "<<S<<std::endl;
   std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;

tests/forces/Test_gpdwf_force.cc

@@ -143,7 +143,7 @@ int main (int argc, char ** argv)
     dS = dS+trace(mommu*forcemu)*dt;
   }
 
-  ComplexD dSpred    = sum(dS);
+  Complex dSpred    = sum(dS);
 
   // From TwoFlavourPseudoFermion:
   //////////////////////////////////////////////////////

tests/forces/Test_gpwilson_force.cc

@@ -143,7 +143,7 @@ int main (int argc, char ** argv)
     dS = dS+trace(mommu*forcemu)*dt;
   }
 
-  ComplexD dSpred    = sum(dS);
+  Complex dSpred    = sum(dS);
 
   std::cout << GridLogMessage << " S      "<<S<<std::endl;
   std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;

tests/forces/Test_laplacian_force.cc

@@ -128,7 +128,7 @@ int main (int argc, char ** argv)
     dS = dS + trace(mommu*UdSdUmu)*dt*2.0;
   }
 
-  ComplexD dSpred    = sum(dS);
+  Complex dSpred    = sum(dS);
 
   std::cout << GridLogMessage << " S      "<<S<<std::endl;
   std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;

tests/forces/Test_mobius_force.cc

@@ -141,7 +141,7 @@ int main (int argc, char ** argv)
 
   }
 
-  ComplexD dSpred    = sum(dS);
+  Complex dSpred    = sum(dS);
 
   std::cout << GridLogMessage << " -- S         "<<S<<std::endl;
   std::cout << GridLogMessage << " -- Sprime    "<<Sprime<<std::endl;

tests/forces/Test_partfrac_force.cc

@@ -141,7 +141,7 @@ int main (int argc, char ** argv)
 
   }
 
-  ComplexD dSpred    = sum(dS);
+  Complex dSpred    = sum(dS);
 
   std::cout << GridLogMessage << " S      "<<S<<std::endl;
   std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;

tests/forces/Test_rect_force.cc

@@ -112,7 +112,7 @@ int main (int argc, char ** argv)
     dS = dS - trace(mommu*UdSdUmu)*dt*2.0;
 
   }
-  ComplexD dSpred    = sum(dS);
+  Complex dSpred    = sum(dS);
 
   std::cout << GridLogMessage << " S      "<<S<<std::endl;
   std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;

tests/forces/Test_wilson_force.cc

@@ -178,9 +178,9 @@ int main (int argc, char ** argv)
 
   }
 
-  ComplexD dSpred    = sum(dS);
+  Complex dSpred    = sum(dS);
-  ComplexD dSm       = sum(dSmom);
+  Complex dSm       = sum(dSmom);
-  ComplexD dSm2      = sum(dSmom2);
+  Complex dSm2      = sum(dSmom2);
 
 
   std::cout << GridLogMessage <<"Initial mom hamiltonian is "<< Hmom <<std::endl;

tests/forces/Test_zmobius_force.cc

@@ -155,7 +155,7 @@ int main (int argc, char ** argv)
 
   }
 
-  ComplexD dSpred    = sum(dS);
+  Complex dSpred    = sum(dS);
 
   std::cout << GridLogMessage << " S      "<<S<<std::endl;
   std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;