Merge branch 'develop' of github.com:edbennett/Grid into develop

.gitignore

@@ -9,6 +9,7 @@
 ################
 *~
 *#
+*.sublime-*
 
 # Precompiled Headers #
 #######################
@@ -103,4 +104,11 @@ lib/fftw/*
 # libtool macros #
 ##################
 m4/lt*
-m4/libtool.m4
+m4/libtool.m4
+
+# Buck files #
+##############
+.buck*
+buck-out
+BUCK
+make-bin-BUCK.sh

configure.ac

@@ -99,6 +99,13 @@ case ${ac_MKL} in
         AC_DEFINE([USE_MKL], [1], [Define to 1 if you use the Intel MKL]);;
 esac
 
+############### HDF5
+AC_ARG_WITH([hdf5],
+    [AS_HELP_STRING([--with-hdf5=prefix],
+    [try this for a non-standard install prefix of the HDF5 library])],
+    [AM_CXXFLAGS="-I$with_hdf5/include $AM_CXXFLAGS"]
+    [AM_LDFLAGS="-L$with_hdf5/lib $AM_LDFLAGS"])
+
 ############### first-touch
 AC_ARG_ENABLE([numa],
     [AC_HELP_STRING([--enable-numa=yes|no|prefix], [enable first touch numa opt])], 
@@ -145,6 +152,12 @@ AC_SEARCH_LIBS([fftw_execute], [fftw3],
                [AC_DEFINE([HAVE_FFTW], [1], [Define to 1 if you have the `FFTW' library])]
                [have_fftw=true])
 
+AC_SEARCH_LIBS([H5Fopen], [hdf5_cpp],
+               [AC_DEFINE([HAVE_HDF5], [1], [Define to 1 if you have the `HDF5' library])]
+               [have_hdf5=true]
+               [LIBS="${LIBS} -lhdf5"], [], [-lhdf5])
+AM_CONDITIONAL(BUILD_HDF5, [ test "${have_hdf5}X" == "trueX" ])
+
 CXXFLAGS=$CXXFLAGS_CPY
 LDFLAGS=$LDFLAGS_CPY
 
@@ -306,7 +319,7 @@ AM_CONDITIONAL(BUILD_COMMS_MPI3L, [ test "${comms_type}X" == "mpi3lX" ] )
 AM_CONDITIONAL(BUILD_COMMS_NONE,  [ test "${comms_type}X" == "noneX" ])
 
 ############### RNG selection
-AC_ARG_ENABLE([rng],[AC_HELP_STRING([--enable-rng=ranlux48|mt19937],\
+AC_ARG_ENABLE([rng],[AC_HELP_STRING([--enable-rng=ranlux48|mt19937|sitmo],\
 	            [Select Random Number Generator to be used])],\
 	            [ac_RNG=${enable_rng}],[ac_RNG=ranlux48])
 
@@ -317,6 +330,9 @@ case ${ac_RNG} in
      mt19937)
       AC_DEFINE([RNG_MT19937],[1],[RNG_MT19937] )
      ;;
+     sitmo)
+      AC_DEFINE([RNG_SITMO],[1],[RNG_SITMO] )
+     ;;
      *)
       AC_MSG_ERROR([${ac_RNG} unsupported --enable-rng option]); 
      ;;
@@ -410,6 +426,7 @@ RNG choice                  : ${ac_RNG}
 GMP                         : `if test "x$have_gmp" = xtrue; then echo yes; else echo no; fi`
 LAPACK                      : ${ac_LAPACK}
 FFTW                        : `if test "x$have_fftw" = xtrue; then echo yes; else echo no; fi`
+HDF5                        : `if test "x$have_hdf5" = xtrue; then echo yes; else echo no; fi`
 build DOXYGEN documentation : `if test "$DX_FLAG_doc" = '1'; then echo yes; else echo no; fi`
 ----- BUILD FLAGS -------------------------------------
 CXXFLAGS:

extras/Hadrons/Application.cc

@@ -42,7 +42,6 @@ using namespace Hadrons;
  ******************************************************************************/
 // constructors ////////////////////////////////////////////////////////////////
 Application::Application(void)
-: env_(Environment::getInstance())
 {
     LOG(Message) << "Modules available:" << std::endl;
     auto list = ModuleFactory::getInstance().getBuilderList();
@@ -74,11 +73,17 @@ Application::Application(const std::string parameterFileName)
     parameterFileName_ = parameterFileName;
 }
 
+// environment shortcut ////////////////////////////////////////////////////////
+Environment & Application::env(void) const
+{
+    return Environment::getInstance();
+}
+
 // access //////////////////////////////////////////////////////////////////////
 void Application::setPar(const Application::GlobalPar &par)
 {
     par_ = par;
-    env_.setSeed(strToVec<int>(par_.seed));
+    env().setSeed(strToVec<int>(par_.seed));
 }
 
 const Application::GlobalPar & Application::getPar(void)
@@ -89,7 +94,7 @@ const Application::GlobalPar & Application::getPar(void)
 // execute /////////////////////////////////////////////////////////////////////
 void Application::run(void)
 {
-    if (!parameterFileName_.empty() and (env_.getNModule() == 0))
+    if (!parameterFileName_.empty() and (env().getNModule() == 0))
     {
         parseParameterFile(parameterFileName_);
     }
@@ -124,7 +129,7 @@ void Application::parseParameterFile(const std::string parameterFileName)
     do
     {
         read(reader, "id", id);
-        env_.createModule(id.name, id.type, reader);
+        env().createModule(id.name, id.type, reader);
     } while (reader.nextElement("module"));
     pop(reader);
     pop(reader);
@@ -134,7 +139,7 @@ void Application::saveParameterFile(const std::string parameterFileName)
 {
     XmlWriter          writer(parameterFileName);
     ObjectId           id;
-    const unsigned int nMod = env_.getNModule();
+    const unsigned int nMod = env().getNModule();
     
     LOG(Message) << "Saving application to '" << parameterFileName << "'..." << std::endl;
     write(writer, "parameters", getPar());
@@ -142,10 +147,10 @@ void Application::saveParameterFile(const std::string parameterFileName)
     for (unsigned int i = 0; i < nMod; ++i)
     {
         push(writer, "module");
-        id.name = env_.getModuleName(i);
+        id.name = env().getModuleName(i);
-        id.type = env_.getModule(i)->getRegisteredName();
+        id.type = env().getModule(i)->getRegisteredName();
         write(writer, "id", id);
-        env_.getModule(i)->saveParameters(writer, "options");
+        env().getModule(i)->saveParameters(writer, "options");
         pop(writer);
     }
     pop(writer);
@@ -164,10 +169,10 @@ auto memPeak = [this](const std::vector<unsigned int> &program)\
     \
     msg = HadronsLogMessage.isActive();\
     HadronsLogMessage.Active(false);\
-    env_.dryRun(true);\
+    env().dryRun(true);\
-    memPeak = env_.executeProgram(program);\
+    memPeak = env().executeProgram(program);\
-    env_.dryRun(false);\
+    env().dryRun(false);\
-    env_.freeAll();\
+    env().freeAll();\
     HadronsLogMessage.Active(true);\
     \
     return memPeak;\
@@ -179,7 +184,7 @@ void Application::schedule(void)
     
     // build module dependency graph
     LOG(Message) << "Building module graph..." << std::endl;
-    auto graph = env_.makeModuleGraph();
+    auto graph = env().makeModuleGraph();
     auto con = graph.getConnectedComponents();
     
     // constrained topological sort using a genetic algorithm
@@ -256,7 +261,7 @@ void Application::saveSchedule(const std::string filename)
                  << std::endl;
     for (auto address: program_)
     {
-        program.push_back(env_.getModuleName(address));
+        program.push_back(env().getModuleName(address));
     }
     write(writer, "schedule", program);
 }
@@ -274,7 +279,7 @@ void Application::loadSchedule(const std::string filename)
     program_.clear();
     for (auto &name: program)
     {
-        program_.push_back(env_.getModuleAddress(name));
+        program_.push_back(env().getModuleAddress(name));
     }
     scheduled_ = true;
     memPeak_   = memPeak(program_);
@@ -291,7 +296,7 @@ void Application::printSchedule(void)
     for (unsigned int i = 0; i < program_.size(); ++i)
     {
         LOG(Message) << std::setw(4) << i + 1 << ": "
-                     << env_.getModuleName(program_[i]) << std::endl;
+                     << env().getModuleName(program_[i]) << std::endl;
     }
 }
 
@@ -304,9 +309,9 @@ void Application::configLoop(void)
     {
         LOG(Message) << BIG_SEP << " Starting measurement for trajectory " << t
                      << " " << BIG_SEP << std::endl;
-        env_.setTrajectory(t);
+        env().setTrajectory(t);
-        env_.executeProgram(program_);
+        env().executeProgram(program_);
-        env_.freeAll();
     }
     LOG(Message) << BIG_SEP << " End of measurement " << BIG_SEP << std::endl;
+    env().freeAll();
 }

extras/Hadrons/Application.hpp

@@ -98,11 +98,13 @@ public:
     void printSchedule(void);
     // loop on configurations
     void configLoop(void);
+private:
+    // environment shortcut
+    Environment & env(void) const;
 private:
     long unsigned int         locVol_;
     std::string               parameterFileName_{""};
     GlobalPar                 par_;
-    Environment               &env_;
     std::vector<unsigned int> program_;
     Environment::Size         memPeak_;
     bool                      scheduled_{false};
@@ -115,14 +117,14 @@ private:
 template <typename M>
 void Application::createModule(const std::string name)
 {
-    env_.createModule<M>(name);
+    env().createModule<M>(name);
 }
 
 template <typename M>
 void Application::createModule(const std::string name,
                                const typename M::Par &par)
 {
-    env_.createModule<M>(name, par);
+    env().createModule<M>(name, par);
 }
 
 END_HADRONS_NAMESPACE

extras/Hadrons/Environment.cc

@@ -41,8 +41,9 @@ using namespace Hadrons;
 // constructor /////////////////////////////////////////////////////////////////
 Environment::Environment(void)
 {
+    nd_ = GridDefaultLatt().size();
     grid4d_.reset(SpaceTimeGrid::makeFourDimGrid(
-        GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()),
+        GridDefaultLatt(), GridDefaultSimd(nd_, vComplex::Nsimd()),
         GridDefaultMpi()));
     gridRb4d_.reset(SpaceTimeGrid::makeFourDimRedBlackGrid(grid4d_.get()));
     auto loc = getGrid()->LocalDimensions();
@@ -126,6 +127,11 @@ GridRedBlackCartesian * Environment::getRbGrid(const unsigned int Ls) const
     }
 }
 
+unsigned int Environment::getNd(void) const
+{
+    return nd_;
+}
+
 // random number generator /////////////////////////////////////////////////////
 void Environment::setSeed(const std::vector<int> &seed)
 {
@@ -410,12 +416,19 @@ Environment::Size Environment::executeProgram(const std::vector<std::string> &p)
 // general memory management ///////////////////////////////////////////////////
 void Environment::addObject(const std::string name, const int moduleAddress)
 {
-    ObjInfo info;
+    if (!hasObject(name))
-    
+    {
-    info.name   = name;
+        ObjInfo info;
-    info.module = moduleAddress;
+        
-    object_.push_back(std::move(info));
+        info.name   = name;
-    objectAddress_[name] = static_cast<unsigned int>(object_.size() - 1);
+        info.module = moduleAddress;
+        object_.push_back(std::move(info));
+        objectAddress_[name] = static_cast<unsigned int>(object_.size() - 1);
+    }
+    else
+    {
+        HADRON_ERROR("object '" + name + "' already exists");
+    }
 }
 
 void Environment::registerObject(const unsigned int address,
@@ -444,6 +457,10 @@ void Environment::registerObject(const unsigned int address,
 void Environment::registerObject(const std::string name,
                                  const unsigned int size, const unsigned int Ls)
 {
+    if (!hasObject(name))
+    {
+        addObject(name);
+    }
     registerObject(getObjectAddress(name), size, Ls);
 }
 
@@ -573,6 +590,30 @@ bool Environment::hasRegisteredObject(const std::string name) const
     }
 }
 
+bool Environment::hasCreatedObject(const unsigned int address) const
+{
+    if (hasObject(address))
+    {
+        return (object_[address].data != nullptr);
+    }
+    else
+    {
+        return false;
+    }
+}
+
+bool Environment::hasCreatedObject(const std::string name) const
+{
+    if (hasObject(name))
+    {
+        return hasCreatedObject(getObjectAddress(name));
+    }
+    else
+    {
+        return false;
+    }
+}
+
 bool Environment::isObject5d(const unsigned int address) const
 {
     return (getObjectLs(address) > 1);

extras/Hadrons/Environment.hpp

@@ -106,6 +106,7 @@ public:
     void                    createGrid(const unsigned int Ls);
     GridCartesian *         getGrid(const unsigned int Ls = 1) const;
     GridRedBlackCartesian * getRbGrid(const unsigned int Ls = 1) const;
+    unsigned int            getNd(void) const;
     // random number generator
     void                    setSeed(const std::vector<int> &seed);
     GridParallelRNG *       get4dRng(void) const;
@@ -137,7 +138,7 @@ public:
     Size                    executeProgram(const std::vector<std::string> &p);
     // general memory management
     void                    addObject(const std::string name,
-                                      const int moduleAddress);
+                                      const int moduleAddress = -1);
     void                    registerObject(const unsigned int address,
                                            const unsigned int size,
                                            const unsigned int Ls = 1);
@@ -176,6 +177,8 @@ public:
     bool                    hasObject(const std::string name) const;
     bool                    hasRegisteredObject(const unsigned int address) const;
     bool                    hasRegisteredObject(const std::string name) const;
+    bool                    hasCreatedObject(const unsigned int address) const;
+    bool                    hasCreatedObject(const std::string name) const;
     bool                    isObject5d(const unsigned int address) const;
     bool                    isObject5d(const std::string name) const;
     Environment::Size       getTotalSize(void) const;
@@ -198,6 +201,7 @@ private:
     std::map<unsigned int, GridPt>         grid5d_;
     GridRbPt                               gridRb4d_;
     std::map<unsigned int, GridRbPt>       gridRb5d_;
+    unsigned int                           nd_;
     // random number generator
     RngPt                                  rng4d_;
     // module and related maps

extras/Hadrons/GeneticScheduler.hpp

@@ -134,7 +134,7 @@ void GeneticScheduler<T>::nextGeneration(void)
     LOG(Debug) << "Starting population:\n" << *this << std::endl;
     
     // random mutations
-    PARALLEL_FOR_LOOP
+    //PARALLEL_FOR_LOOP
     for (unsigned int i = 0; i < par_.popSize; ++i)
     {
         doMutation();
@@ -142,7 +142,7 @@ void GeneticScheduler<T>::nextGeneration(void)
     LOG(Debug) << "After mutations:\n" << *this << std::endl;
     
     // mating
-    PARALLEL_FOR_LOOP
+    //PARALLEL_FOR_LOOP
     for (unsigned int i = 0; i < par_.popSize/2; ++i)
     {
         doCrossover();
@@ -166,7 +166,7 @@ void GeneticScheduler<T>::initPopulation(void)
     {
         auto p = graph_.topoSort(gen_);
         
-        population_.emplace(func_(p), p);
+        population_.insert(std::make_pair(func_(p), p));
     }
 }
 
@@ -180,8 +180,8 @@ void GeneticScheduler<T>::doCrossover(void)
     crossover(c1, c2, p1, p2);
     PARALLEL_CRITICAL
     {
-        population_.emplace(func_(c1), c1);
+        population_.insert(std::make_pair(func_(c1), c1));
-        population_.emplace(func_(c2), c2);
+        population_.insert(std::make_pair(func_(c2), c2));
     }
 }
 
@@ -200,7 +200,7 @@ void GeneticScheduler<T>::doMutation(void)
         mutation(m, it->second);
         PARALLEL_CRITICAL
         {
-            population_.emplace(func_(m), m);
+            population_.insert(std::make_pair(func_(m), m));
         }
     }
 }

extras/Hadrons/Modules/MSource/SeqGamma.hpp

@@ -147,7 +147,7 @@ void TSeqGamma<FImpl>::execute(void)
     g  = makeGammaProd(par().gamma);
     p  = strToVec<Real>(par().mom);
     ph = zero;
-    for(unsigned int mu = 0; mu < Nd; mu++)
+    for(unsigned int mu = 0; mu < env().getNd(); mu++)
     {
         LatticeCoordinate(coor, mu);
         ph = ph + p[mu]*coor;

lib/Grid.h

@@ -59,13 +59,13 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 ///////////////////
 // Grid headers
 ///////////////////
-#include <Grid/serialisation/Serialisation.h>
 #include "Config.h"
 #include <Grid/Timer.h>
 #include <Grid/PerfCount.h>
 #include <Grid/Log.h>
 #include <Grid/AlignedAllocator.h>
 #include <Grid/Simd.h>
+#include <Grid/serialisation/Serialisation.h>
 #include <Grid/Threads.h>
 #include <Grid/Lexicographic.h>
 #include <Grid/Init.h>

lib/Init.cc

@@ -41,7 +41,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #include <signal.h>
 #include <iostream>
 #include <iterator>
-#include <Grid.h>
+#include <Grid/Grid.h>
 #include <algorithm>
 #include <iterator>
 #include <cstdlib>

lib/Log.cc

@@ -29,7 +29,7 @@ See the full license in the file "LICENSE" in the top level distribution
 directory
 *************************************************************************************/
 /*  END LEGAL */
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 #include <cxxabi.h>
 

lib/Makefile.am

@@ -1,4 +1,5 @@
 extra_sources=
+extra_headers=
 if BUILD_COMMS_MPI
   extra_sources+=communicator/Communicator_mpi.cc
   extra_sources+=communicator/Communicator_base.cc
@@ -24,6 +25,12 @@ if BUILD_COMMS_NONE
   extra_sources+=communicator/Communicator_base.cc
 endif
 
+if BUILD_HDF5
+  extra_sources+=serialisation/Hdf5IO.cc 
+  extra_headers+=serialisation/Hdf5IO.h
+  extra_headers+=serialisation/Hdf5Type.h
+endif
+
 #
 # Libraries
 #
@@ -32,6 +39,9 @@ include Eigen.inc
 
 lib_LIBRARIES = libGrid.a
 
-libGrid_a_SOURCES              = $(CCFILES) $(extra_sources)
+CCFILES += $(extra_sources)
+HFILES  += $(extra_headers)
+
+libGrid_a_SOURCES              = $(CCFILES)
 libGrid_adir                   = $(pkgincludedir)
 nobase_dist_pkginclude_HEADERS = $(HFILES) $(eigen_files) Config.h

lib/PerfCount.cc

@@ -26,8 +26,8 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     *************************************************************************************/
     /*  END LEGAL */
 
-#include <Grid.h>
+#include <Grid/Grid.h>
-#include <PerfCount.h>
+#include <Grid/PerfCount.h>
 
 namespace Grid {
 

lib/Stat.cc

@@ -1,6 +1,6 @@
-#include <Grid.h>
+#include <Grid/Grid.h>
-#include <PerfCount.h>
+#include <Grid/PerfCount.h>
-#include <Stat.h>
+#include <Grid/Stat.h>
 
 
 namespace Grid { 

lib/algorithms/approx/MultiShiftFunction.cc

@@ -25,7 +25,7 @@ Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 namespace Grid {
 double MultiShiftFunction::approx(double x)

lib/algorithms/approx/Remez.cc

@@ -20,7 +20,7 @@
 #include<iomanip>
 #include<cassert>
 
-#include<algorithms/approx/Remez.h>
+#include<Grid/algorithms/approx/Remez.h>
 
 // Constructor
 AlgRemez::AlgRemez(double lower, double upper, long precision) 

lib/algorithms/iterative/Matrix.h

@@ -36,7 +36,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #include <iomanip>
 #include <complex>
 #include <typeinfo>
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 
 /** Sign function **/

lib/communicator/Communicator_base.cc

@@ -25,7 +25,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include "Grid.h"
+#include <Grid/Grid.h>
 namespace Grid {
 
 ///////////////////////////////////////////////////////////////

lib/communicator/Communicator_mpi.cc

@@ -25,7 +25,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include "Grid.h"
+#include <Grid/Grid.h>
 #include <mpi.h>
 
 namespace Grid {

lib/communicator/Communicator_mpi3.cc

@@ -25,7 +25,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include "Grid.h"
+#include <Grid/Grid.h>
 #include <mpi.h>
 
 namespace Grid {

lib/communicator/Communicator_none.cc

@@ -25,7 +25,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include "Grid.h"
+#include <Grid/Grid.h>
 namespace Grid {
 
 ///////////////////////////////////////////////////////////////////////////////////////////////////

lib/communicator/Communicator_shmem.cc

@@ -25,7 +25,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include "Grid.h"
+#include <Grid/Grid.h>
 #include <mpp/shmem.h>
 
 namespace Grid {

lib/lattice/Lattice_rng.h

@@ -30,6 +30,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #define GRID_LATTICE_RNG_H
 
 #include <random>
+#include <Grid/sitmo_rng/sitmo_prng_engine.hpp>
 
 namespace Grid {
 
@@ -114,10 +115,14 @@ namespace Grid {
     typedef uint64_t      RngStateType;
     typedef std::ranlux48 RngEngine;
     static const int RngStateCount = 15;
-#else
+#elif RNG_MT19937 
     typedef std::mt19937 RngEngine;
     typedef uint32_t     RngStateType;
     static const int     RngStateCount = std::mt19937::state_size;
+#elif RNG_SITMO
+    typedef sitmo::prng_engine 	RngEngine;
+    typedef uint64_t    	RngStateType;
+    static const int    	RngStateCount = 4;
 #endif
     std::vector<RngEngine>                             _generators;
     std::vector<std::uniform_real_distribution<RealD>> _uniform;

lib/lattice/Lattice_transfer.h

@@ -386,6 +386,7 @@ void InsertSlice(Lattice<vobj> &lowDim,Lattice<vobj> & higherDim,int slice, int
   }
 
   // the above should guarantee that the operations are local
+  // Guido: check the threading here
   //PARALLEL_FOR_LOOP
   for(int idx=0;idx<lg->lSites();idx++){
     std::vector<int> lcoor(nl);

lib/pugixml/pugixml.cc

@@ -14,7 +14,7 @@
 #ifndef SOURCE_PUGIXML_CPP
 #define SOURCE_PUGIXML_CPP
 
-#include <pugixml/pugixml.h>
+#include <Grid/pugixml/pugixml.h>
 
 #include <stdlib.h>
 #include <stdio.h>

lib/qcd/action/fermion/CayleyFermion5D.cc

@@ -30,7 +30,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     /*  END LEGAL */
 
 #include <Grid/Eigen/Dense>
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 
 namespace Grid {

lib/qcd/action/fermion/CayleyFermion5Dcache.cc

@@ -29,7 +29,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     *************************************************************************************/
     /*  END LEGAL */
 
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 
 namespace Grid {

lib/qcd/action/fermion/CayleyFermion5Ddense.cc

@@ -30,7 +30,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     /*  END LEGAL */
 
 #include <Grid/Eigen/Dense>
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 
 namespace Grid {

lib/qcd/action/fermion/CayleyFermion5Dssp.cc

@@ -29,7 +29,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     *************************************************************************************/
     /*  END LEGAL */
 
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 
 namespace Grid {

lib/qcd/action/fermion/CayleyFermion5Dvec.cc

@@ -30,7 +30,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     /*  END LEGAL */
 
 
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 
 namespace Grid {

lib/qcd/action/fermion/ContinuedFractionFermion5D.cc

@@ -26,7 +26,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 namespace Grid {
   namespace QCD {

lib/qcd/action/fermion/PartialFractionFermion5D.cc

@@ -26,7 +26,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include <Grid.h>
+#include <Grid/Grid.h>
 namespace Grid {
   namespace QCD {
 

lib/qcd/action/fermion/WilsonFermion.cc

@@ -29,7 +29,7 @@ See the full license in the file "LICENSE" in the top level distribution
 directory
 *************************************************************************************/
 /*  END LEGAL */
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 namespace Grid {
 namespace QCD {

lib/qcd/action/fermion/WilsonFermion5D.cc

@@ -29,8 +29,8 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include <Grid.h>
+#include <Grid/Grid.h>
-#include <PerfCount.h>
+#include <Grid/PerfCount.h>
 
 namespace Grid {
 namespace QCD {

lib/qcd/action/fermion/WilsonKernels.cc

@@ -28,7 +28,7 @@ See the full license in the file "LICENSE" in the top level distribution
 directory
 *************************************************************************************/
 /*  END LEGAL */
-#include <Grid.h>
+#include <Grid/Grid.h>
 namespace Grid {
 namespace QCD {
 

lib/qcd/action/fermion/WilsonKernelsAsm.cc

@@ -30,7 +30,7 @@ Author: Guido Cossu <guido.cossu@ed.ac.uk>
 *************************************************************************************/
 /*  END LEGAL */
 
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 
 namespace Grid {

lib/qcd/action/fermion/WilsonKernelsHand.cc

@@ -26,7 +26,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 #define REGISTER
 

lib/qcd/action/fermion/WilsonTMFermion.cc

@@ -25,7 +25,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 namespace Grid {
 namespace QCD {

lib/qcd/hmc/HMC.cc

@@ -27,7 +27,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 namespace Grid{
   namespace QCD{

lib/qcd/spin/Dirac.cc

@@ -25,7 +25,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 namespace Grid {
 

lib/qcd/utils/SpaceTimeGrid.cc

@@ -25,7 +25,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 namespace Grid { 
   namespace QCD {

lib/serialisation/BaseIO.h

@@ -32,6 +32,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #include <type_traits>
 
 namespace Grid {
+  // Vector IO utilities ///////////////////////////////////////////////////////
   // helper function to read space-separated values
   template <typename T>
   std::vector<T> strToVec(const std::string s)
@@ -67,6 +68,77 @@ namespace Grid {
     return os;
   }
   
+  // Vector element trait //////////////////////////////////////////////////////  
+  template <typename T>
+  struct element
+  {
+    typedef T type;
+    static constexpr bool is_number = false;
+  };
+  
+  template <typename T>
+  struct element<std::vector<T>>
+  {
+    typedef typename element<T>::type type;
+    static constexpr bool is_number = std::is_arithmetic<T>::value
+                                      or is_complex<T>::value
+                                      or element<T>::is_number;
+  };
+  
+  // Vector flatening utility class ////////////////////////////////////////////
+  // Class to flatten a multidimensional std::vector
+  template <typename V>
+  class Flatten
+  {
+  public:
+    typedef typename element<V>::type Element;
+  public:
+    explicit                     Flatten(const V &vector);
+    const V &                    getVector(void);
+    const std::vector<Element> & getFlatVector(void);
+    const std::vector<size_t>  & getDim(void);
+  private:
+    void accumulate(const Element &e);
+    template <typename W>
+    void accumulate(const W &v);
+    void accumulateDim(const Element &e);
+    template <typename W>
+    void accumulateDim(const W &v);
+  private:
+    const V              &vector_;
+    std::vector<Element> flatVector_;
+    std::vector<size_t>  dim_;
+  };
+  
+  
+  // Class to reconstruct a multidimensional std::vector
+  template <typename V>
+  class Reconstruct
+  {
+  public:
+    typedef typename element<V>::type Element;
+  public:
+    Reconstruct(const std::vector<Element> &flatVector,
+                const std::vector<size_t> &dim);
+    const V &                    getVector(void);
+    const std::vector<Element> & getFlatVector(void);
+    const std::vector<size_t>  & getDim(void);
+  private:
+    void fill(std::vector<Element> &v);
+    template <typename W>
+    void fill(W &v);
+    void resize(std::vector<Element> &v, const unsigned int dim);
+    template <typename W>
+    void resize(W &v, const unsigned int dim);
+  private:
+    V                          vector_;
+    const std::vector<Element> &flatVector_;
+    std::vector<size_t>        dim_;
+    size_t                     ind_{0};
+    unsigned int               dimInd_{0};
+  };
+  
+  // Abstract writer/reader classes ////////////////////////////////////////////
   // static polymorphism implemented using CRTP idiom
   class Serializable;
   
@@ -132,7 +204,128 @@ namespace Grid {
     }
   };
   
-  // Generic writer interface
+  // Flatten class template implementation /////////////////////////////////////
+  template <typename V>
+  void Flatten<V>::accumulate(const Element &e)
+  {
+    flatVector_.push_back(e);
+  }
+  
+  template <typename V>
+  template <typename W>
+  void Flatten<V>::accumulate(const W &v)
+  {
+    for (auto &e: v)
+    {
+      accumulate(e);
+    }
+  }
+  
+  template <typename V>
+  void Flatten<V>::accumulateDim(const Element &e) {};
+  
+  template <typename V>
+  template <typename W>
+  void Flatten<V>::accumulateDim(const W &v)
+  {
+    dim_.push_back(v.size());
+    accumulateDim(v[0]);
+  }
+  
+  template <typename V>
+  Flatten<V>::Flatten(const V &vector)
+  : vector_(vector)
+  {
+    accumulate(vector_);
+    accumulateDim(vector_);
+  }
+  
+  template <typename V>
+  const V & Flatten<V>::getVector(void)
+  {
+    return vector_;
+  }
+  
+  template <typename V>
+  const std::vector<typename Flatten<V>::Element> &
+  Flatten<V>::getFlatVector(void)
+  {
+    return flatVector_;
+  }
+  
+  template <typename V>
+  const std::vector<size_t> & Flatten<V>::getDim(void)
+  {
+    return dim_;
+  }
+  
+  // Reconstruct class template implementation /////////////////////////////////
+  template <typename V>
+  void Reconstruct<V>::fill(std::vector<Element> &v)
+  {
+    for (auto &e: v)
+    {
+      e = flatVector_[ind_++];
+    }
+  }
+  
+  template <typename V>
+  template <typename W>
+  void Reconstruct<V>::fill(W &v)
+  {
+    for (auto &e: v)
+    {
+      fill(e);
+    }
+  }
+  
+  template <typename V>
+  void Reconstruct<V>::resize(std::vector<Element> &v, const unsigned int dim)
+  {
+    v.resize(dim_[dim]);
+  }
+  
+  template <typename V>
+  template <typename W>
+  void Reconstruct<V>::resize(W &v, const unsigned int dim)
+  {
+    v.resize(dim_[dim]);
+    for (auto &e: v)
+    {
+      resize(e, dim + 1);
+    }
+  }
+  
+  template <typename V>
+  Reconstruct<V>::Reconstruct(const std::vector<Element> &flatVector,
+                              const std::vector<size_t> &dim)
+  : flatVector_(flatVector)
+  , dim_(dim)
+  {
+    resize(vector_, 0);
+    fill(vector_);
+  }
+  
+  template <typename V>
+  const V & Reconstruct<V>::getVector(void)
+  {
+    return vector_;
+  }
+  
+  template <typename V>
+  const std::vector<typename Reconstruct<V>::Element> &
+  Reconstruct<V>::getFlatVector(void)
+  {
+    return flatVector_;
+  }
+  
+  template <typename V>
+  const std::vector<size_t> & Reconstruct<V>::getDim(void)
+  {
+    return dim_;
+  }
+  
+  // Generic writer interface //////////////////////////////////////////////////
   template <typename T>
   inline void push(Writer<T> &w, const std::string &s)
   {
@@ -217,7 +410,7 @@ namespace Grid {
     upcast->writeDefault(s, output);
   }
   
-  // Reader template implementation ////////////////////////////////////////////
+  // Reader template implementation
   template <typename T>
   Reader<T>::Reader(void)
   {

lib/serialisation/BinaryIO.cc

@@ -26,7 +26,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 using namespace Grid;
 using namespace std;

lib/serialisation/Hdf5IO.cc

@@ -0,0 +1,103 @@
+#include <Grid/Grid.h>
+
+using namespace Grid;
+#ifndef H5_NO_NAMESPACE
+using namespace H5NS;
+#endif
+
+// Writer implementation ///////////////////////////////////////////////////////
+Hdf5Writer::Hdf5Writer(const std::string &fileName)
+: fileName_(fileName)
+, file_(fileName.c_str(), H5F_ACC_TRUNC)
+{
+  group_ = file_.openGroup("/");
+  writeSingleAttribute(dataSetThres_, HDF5_GRID_GUARD "dataset_threshold",
+                       Hdf5Type<unsigned int>::type());
+}
+
+void Hdf5Writer::push(const std::string &s)
+{
+  group_ = group_.createGroup(s);
+  path_.push_back(s);
+}
+
+void Hdf5Writer::pop(void)
+{
+  path_.pop_back();
+  if (path_.empty())
+  {
+    group_ = file_.openGroup("/");
+  }
+  else
+  {
+    auto binOp = [](const std::string &a, const std::string &b)->std::string
+    {
+      return a + "/" + b;
+    };
+    
+    group_ = group_.openGroup(std::accumulate(path_.begin(), path_.end(),
+                                              std::string(""), binOp));
+  }
+}
+
+template <>
+void Hdf5Writer::writeDefault(const std::string &s, const std::string &x)
+{
+  StrType     strType(PredType::C_S1, x.size());
+  
+  writeSingleAttribute(*(x.data()), s, strType);
+}
+
+void Hdf5Writer::writeDefault(const std::string &s, const char *x)
+{
+  std::string sx(x);
+  
+  writeDefault(s, sx);
+}
+
+// Reader implementation ///////////////////////////////////////////////////////
+Hdf5Reader::Hdf5Reader(const std::string &fileName)
+: fileName_(fileName)
+, file_(fileName.c_str(), H5F_ACC_RDONLY)
+{
+  group_ = file_.openGroup("/");
+  readSingleAttribute(dataSetThres_, HDF5_GRID_GUARD "dataset_threshold",
+                      Hdf5Type<unsigned int>::type());
+}
+
+void Hdf5Reader::push(const std::string &s)
+{
+  group_ = group_.openGroup(s);
+  path_.push_back(s);
+}
+
+void Hdf5Reader::pop(void)
+{
+  path_.pop_back();
+  if (path_.empty())
+  {
+    group_ = file_.openGroup("/");
+  }
+  else
+  {
+    auto binOp = [](const std::string &a, const std::string &b)->std::string
+    {
+      return a + "/" + b;
+    };
+    
+    group_ = group_.openGroup(std::accumulate(path_.begin(), path_.end(),
+                                              std::string(""), binOp));
+  }
+}
+
+template <>
+void Hdf5Reader::readDefault(const std::string &s, std::string &x)
+{
+  Attribute attribute;
+  
+  attribute       = group_.openAttribute(s);
+  StrType strType = attribute.getStrType();
+  
+  x.resize(strType.getSize());
+  attribute.read(strType, &(x[0]));
+}

lib/serialisation/Hdf5IO.h

@@ -0,0 +1,242 @@
+#ifndef GRID_SERIALISATION_HDF5_H
+#define GRID_SERIALISATION_HDF5_H
+
+#include <stack>
+#include <string>
+#include <vector>
+#include <H5Cpp.h>
+#include "Hdf5Type.h"
+
+#ifndef H5_NO_NAMESPACE
+#define H5NS H5
+#endif
+
+// default thresold above which datasets are used instead of attributes
+#ifndef HDF5_DEF_DATASET_THRES
+#define HDF5_DEF_DATASET_THRES 6u
+#endif
+
+// name guard for Grid metadata
+#define HDF5_GRID_GUARD "_Grid_"
+
+namespace Grid
+{
+  class Hdf5Writer: public Writer<Hdf5Writer>
+  {
+  public:
+    Hdf5Writer(const std::string &fileName);
+    virtual ~Hdf5Writer(void) = default;
+    void push(const std::string &s);
+    void pop(void);
+    void writeDefault(const std::string &s, const char *x);
+    template <typename U>
+    void writeDefault(const std::string &s, const U &x);
+    template <typename U>
+    typename std::enable_if<element<std::vector<U>>::is_number, void>::type
+    writeDefault(const std::string &s, const std::vector<U> &x);
+    template <typename U>
+    typename std::enable_if<!element<std::vector<U>>::is_number, void>::type
+    writeDefault(const std::string &s, const std::vector<U> &x);
+  private:
+    template <typename U>
+    void writeSingleAttribute(const U &x, const std::string &name,
+                              const H5NS::DataType &type);
+  private:
+    std::string              fileName_;
+    std::vector<std::string> path_;
+    H5NS::H5File             file_;
+    H5NS::Group              group_;
+    unsigned int             dataSetThres_{HDF5_DEF_DATASET_THRES};
+  };
+  
+  class Hdf5Reader: public Reader<Hdf5Reader>
+  {
+  public:
+    Hdf5Reader(const std::string &fileName);
+    virtual ~Hdf5Reader(void) = default;
+    void push(const std::string &s);
+    void pop(void);
+    template <typename U>
+    void readDefault(const std::string &s, U &output);
+    template <typename U>
+    typename std::enable_if<element<std::vector<U>>::is_number, void>::type
+    readDefault(const std::string &s, std::vector<U> &x);
+    template <typename U>
+    typename std::enable_if<!element<std::vector<U>>::is_number, void>::type
+    readDefault(const std::string &s, std::vector<U> &x);
+  private:
+    template <typename U>
+    void readSingleAttribute(U &x, const std::string &name,
+                             const H5NS::DataType &type);
+  private:
+    std::string              fileName_;
+    std::vector<std::string> path_;
+    H5NS::H5File             file_;
+    H5NS::Group              group_;
+    unsigned int             dataSetThres_;
+  };
+  
+  // Writer template implementation ////////////////////////////////////////////
+  template <typename U>
+  void Hdf5Writer::writeSingleAttribute(const U &x, const std::string &name,
+                                        const H5NS::DataType &type)
+  {
+    H5NS::Attribute attribute;
+    hsize_t         attrDim = 1;
+    H5NS::DataSpace attrSpace(1, &attrDim);
+    
+    attribute = group_.createAttribute(name, type, attrSpace);
+    attribute.write(type, &x);
+  }
+  
+  template <typename U>
+  void Hdf5Writer::writeDefault(const std::string &s, const U &x)
+  {
+    writeSingleAttribute(x, s, Hdf5Type<U>::type());
+  }
+  
+  template <>
+  void Hdf5Writer::writeDefault(const std::string &s, const std::string &x);
+  
+  template <typename U>
+  typename std::enable_if<element<std::vector<U>>::is_number, void>::type
+  Hdf5Writer::writeDefault(const std::string &s, const std::vector<U> &x)
+  {
+    // alias to element type
+    typedef typename element<std::vector<U>>::type Element;
+    
+    // flatten the vector and getting dimensions
+    Flatten<std::vector<U>> flat(x);
+    std::vector<hsize_t> dim;
+    const auto           &flatx = flat.getFlatVector();
+    
+    for (auto &d: flat.getDim())
+    {
+      dim.push_back(d);
+    }
+    
+    // write to file
+    H5NS::DataSpace dataSpace(dim.size(), dim.data());
+    
+    if (flatx.size() > dataSetThres_)
+    {
+      H5NS::DataSet dataSet;
+      
+      dataSet = group_.createDataSet(s, Hdf5Type<Element>::type(), dataSpace);
+      dataSet.write(flatx.data(), Hdf5Type<Element>::type());
+    }
+    else
+    {
+      H5NS::Attribute attribute;
+      
+      attribute = group_.createAttribute(s, Hdf5Type<Element>::type(), dataSpace);
+      attribute.write(Hdf5Type<Element>::type(), flatx.data());
+    }
+  }
+  
+  template <typename U>
+  typename std::enable_if<!element<std::vector<U>>::is_number, void>::type
+  Hdf5Writer::writeDefault(const std::string &s, const std::vector<U> &x)
+  {
+    push(s);
+    writeSingleAttribute(x.size(), HDF5_GRID_GUARD "vector_size",
+                         Hdf5Type<uint64_t>::type());
+    for (hsize_t i = 0; i < x.size(); ++i)
+    {
+      write(s + "_" + std::to_string(i), x[i]);
+    }
+    pop();
+  }
+  
+  // Reader template implementation ////////////////////////////////////////////
+  template <typename U>
+  void Hdf5Reader::readSingleAttribute(U &x, const std::string &name,
+                                       const H5NS::DataType &type)
+  {
+    H5NS::Attribute attribute;
+    
+    attribute = group_.openAttribute(name);
+    attribute.read(type, &x);
+  }
+  
+  template <typename U>
+  void Hdf5Reader::readDefault(const std::string &s, U &output)
+  {
+    readSingleAttribute(output, s, Hdf5Type<U>::type());
+  }
+  
+  template <>
+  void Hdf5Reader::readDefault(const std::string &s, std::string &x);
+  
+  template <typename U>
+  typename std::enable_if<element<std::vector<U>>::is_number, void>::type
+  Hdf5Reader::readDefault(const std::string &s, std::vector<U> &x)
+  {
+    // alias to element type
+    typedef typename element<std::vector<U>>::type Element;
+    
+    // read the dimensions
+    H5NS::DataSpace       dataSpace;
+    std::vector<hsize_t>  hdim;
+    std::vector<size_t>   dim;
+    hsize_t               size = 1;
+    
+    if (group_.attrExists(s))
+    {
+      dataSpace = group_.openAttribute(s).getSpace();
+    }
+    else
+    {
+      dataSpace = group_.openDataSet(s).getSpace();
+    }
+    hdim.resize(dataSpace.getSimpleExtentNdims());
+    dataSpace.getSimpleExtentDims(hdim.data());
+    for (auto &d: hdim)
+    {
+      dim.push_back(d);
+      size *= d;
+    }
+    
+    // read the flat vector
+    std::vector<Element> buf(size);
+
+    if (size > dataSetThres_)
+    {
+      H5NS::DataSet dataSet;
+      
+      dataSet = group_.openDataSet(s);
+      dataSet.read(buf.data(), Hdf5Type<Element>::type());
+    }
+    else
+    {
+      H5NS::Attribute attribute;
+      
+      attribute = group_.openAttribute(s);
+      attribute.read(Hdf5Type<Element>::type(), buf.data());
+    }
+    
+    // reconstruct the multidimensional vector
+    Reconstruct<std::vector<U>> r(buf, dim);
+    
+    x = r.getVector();
+  }
+  
+  template <typename U>
+  typename std::enable_if<!element<std::vector<U>>::is_number, void>::type
+  Hdf5Reader::readDefault(const std::string &s, std::vector<U> &x)
+  {
+    uint64_t size;
+    
+    push(s);
+    readSingleAttribute(size, HDF5_GRID_GUARD "vector_size",
+                        Hdf5Type<uint64_t>::type());
+    x.resize(size);
+    for (hsize_t i = 0; i < x.size(); ++i)
+    {
+      read(s + "_" + std::to_string(i), x[i]);
+    }
+    pop();
+  }
+}
+
+#endif

lib/serialisation/Hdf5Type.h

@@ -0,0 +1,77 @@
+#ifndef GRID_SERIALISATION_HDF5_TYPE_H
+#define GRID_SERIALISATION_HDF5_TYPE_H
+
+#include <H5Cpp.h>
+#include <complex>
+#include <memory>
+
+#ifndef H5_NO_NAMESPACE
+#define H5NS H5
+#endif
+
+#define HDF5_NATIVE_TYPE(predType, cType)\
+template <>\
+class Hdf5Type<cType>\
+{\
+public:\
+  static inline const H5NS::DataType & type(void)\
+  {\
+    return H5NS::PredType::predType;\
+  }\
+  static constexpr bool isNative = true;\
+};
+
+#define DEFINE_HDF5_NATIVE_TYPES \
+HDF5_NATIVE_TYPE(NATIVE_B8,      bool);\
+HDF5_NATIVE_TYPE(NATIVE_CHAR,    char);\
+HDF5_NATIVE_TYPE(NATIVE_SCHAR,   signed char);\
+HDF5_NATIVE_TYPE(NATIVE_UCHAR,   unsigned char);\
+HDF5_NATIVE_TYPE(NATIVE_SHORT,   short);\
+HDF5_NATIVE_TYPE(NATIVE_USHORT,  unsigned short);\
+HDF5_NATIVE_TYPE(NATIVE_INT,     int);\
+HDF5_NATIVE_TYPE(NATIVE_UINT,    unsigned int);\
+HDF5_NATIVE_TYPE(NATIVE_LONG,    long);\
+HDF5_NATIVE_TYPE(NATIVE_ULONG,   unsigned long);\
+HDF5_NATIVE_TYPE(NATIVE_LLONG,   long long);\
+HDF5_NATIVE_TYPE(NATIVE_ULLONG,  unsigned long long);\
+HDF5_NATIVE_TYPE(NATIVE_FLOAT,   float);\
+HDF5_NATIVE_TYPE(NATIVE_DOUBLE,  double);\
+HDF5_NATIVE_TYPE(NATIVE_LDOUBLE, long double);
+
+namespace Grid
+{
+  template <typename T> class Hdf5Type
+  {
+  public:
+    static constexpr bool isNative = false;
+  };
+  
+  DEFINE_HDF5_NATIVE_TYPES;
+  
+  template <typename R>
+  class Hdf5Type<std::complex<R>>
+  {
+  public:
+    static inline const H5NS::DataType & type(void)
+    {
+      if (typePtr_ == nullptr)
+      {
+        typePtr_.reset(new H5NS::CompType(sizeof(std::complex<R>)));
+        typePtr_->insertMember("re", 0,         Hdf5Type<R>::type());
+        typePtr_->insertMember("im", sizeof(R), Hdf5Type<R>::type());
+      }
+
+      return *typePtr_;
+    }
+    static constexpr bool isNative = false;
+  private:
+    static std::unique_ptr<H5NS::CompType> typePtr_;
+  };
+  
+  template <typename R>
+  std::unique_ptr<H5NS::CompType> Hdf5Type<std::complex<R>>::typePtr_ = nullptr;
+}
+
+#undef HDF5_NATIVE_TYPE
+
+#endif /* GRID_SERIALISATION_HDF5_TYPE_H */

lib/serialisation/MacroMagic.h

@@ -109,38 +109,36 @@ 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_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,...)		\
+#define GRID_SERIALIZABLE_CLASS_MEMBERS(cname,...)\
-\
+GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_MEMBER,__VA_ARGS__))\
-\
-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){ \
   push(WR,s);\
   GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_WRITE_MEMBER,__VA_ARGS__))	\
   pop(WR);\
-} \
+}\
-\
-\
 template <typename T>\
 static inline void read(Reader<T> &RD,const std::string &s, cname &obj){	\
   push(RD,s);\
   GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_READ_MEMBER,__VA_ARGS__))	\
   pop(RD);\
-} \
+}\
-\
-\
 friend inline std::ostream & operator << (std::ostream &os, const cname &obj ) { \
   os<<"class "<<#cname<<" {"<<std::endl;\
   GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_OS_WRITE_MEMBER,__VA_ARGS__))	\
     os<<"}";								\
   return os;\
-};  
+}\
+friend inline bool operator==(const cname &lhs, const cname &rhs) {\
+  bool result = true;\
+  GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_COMP_MEMBER,__VA_ARGS__))\
+  return result;\
+}
 
 #define GRID_ENUM_TYPE(obj) std::remove_reference<decltype(obj)>::type
 #define GRID_MACRO_ENUMVAL(A,B) A = B,
@@ -149,7 +147,7 @@ friend inline std::ostream & operator << (std::ostream &os, const cname &obj ) {
 #define GRID_MACRO_ENUMCASEIO(A,B) case GRID_ENUM_TYPE(obj)::A: os << #A; break;
 
 #define GRID_SERIALIZABLE_ENUM(name,undefname,...)\
-class name: public Serializable\
+class name: public Grid::Serializable\
 {\
 public:\
   enum EnumType\
@@ -161,7 +159,7 @@ public:\
   name(void): value_(undefname) {};\
   name(EnumType value): value_(value) {};\
   template <typename T>\
-  static inline void write(Writer<T> &WR,const std::string &s, const name &obj)\
+  static inline void write(Grid::Writer<T> &WR,const std::string &s, const name &obj)\
   {\
     switch (obj.value_)\
     {\
@@ -171,7 +169,7 @@ public:\
   }\
   \
   template <typename T>\
-  static inline void read(Reader<T> &RD,const std::string &s, name &obj)\
+  static inline void read(Grid::Reader<T> &RD,const std::string &s, name &obj)\
   {\
     std::string buf;\
     Grid::read(RD, s, buf);\

lib/serialisation/Serialisation.h

@@ -36,6 +36,9 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #include "BinaryIO.h"
 #include "TextIO.h"
 #include "XmlIO.h"
+#ifdef HAVE_HDF5
+#include "Hdf5IO.h"
+#endif
 //////////////////////////////////////////
 // Todo:
 //////////////////////////////////////////

lib/serialisation/TextIO.cc

@@ -26,7 +26,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 using namespace Grid;
 using namespace std;

lib/serialisation/XmlIO.cc

@@ -26,7 +26,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 using namespace Grid;
 using namespace std;

lib/sitmo_rng/sitmo_prng_engine.hpp

@@ -0,0 +1,390 @@
+//  Copyright (c) 2012-2016 M.A. (Thijs) van den Berg, http://sitmo.com/
+//
+//  Use, modification and distribution are subject to the MIT Software License. 
+//  
+//  The MIT License (MIT)
+//  Permission is hereby granted, free of charge, to any person obtaining a copy
+//  of this software and associated documentation files (the "Software"), to deal
+//  in the Software without restriction, including without limitation the rights
+//  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+//  copies of the Software, and to permit persons to whom the Software is
+//  furnished to do so, subject to the following conditions:
+//  
+//  The above copyright notice and this permission notice shall be included in
+//  all copies or substantial portions of the Software.
+//  
+//  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+//  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+//  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+//  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+//  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+//  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+//  THE SOFTWARE.
+
+// version history:
+// version 1,  6 Sep 2012
+// version 2, 10 Dec 2013
+//      bug fix in the discard() routine, it was discarding to many elements
+//      added the version() method
+// version 3...5, 13 Dec 2013
+//      fixed type-conversion earning
+//      fixed potential issues with constructor template matching
+// version 6, 4 March 2016
+//      made min() max() constexpr for C+11 compiler (thanks to James Joseph Balamuta)
+
+#ifndef SITMO_PRNG_ENGINE_HPP
+#define SITMO_PRNG_ENGINE_HPP
+#include <iostream>
+
+#ifdef __GNUC__
+    #include <stdint.h>                 // respecting the C99 standard.
+#endif
+#ifdef _MSC_VER
+    typedef unsigned __int64 uint64_t;  // Visual Studio 6.0(VC6) and newer..
+    typedef unsigned __int32 uint32_t;
+#endif
+
+// Double mixing function
+#define MIX2(x0,x1,rx,z0,z1,rz) \
+    x0 += x1; \
+    z0 += z1; \
+    x1 = (x1 << rx) | (x1 >> (64-rx)); \
+    z1 = (z1 << rz) | (z1 >> (64-rz)); \
+    x1 ^= x0; \
+    z1 ^= z0; 
+
+
+// Double mixing function with key adition
+#define MIXK(x0,x1,rx,z0,z1,rz,k0,k1,l0,l1) \
+    x1 += k1; \
+    z1 += l1; \
+    x0 += x1+k0; \
+    z0 += z1+l0; \
+    x1 = (x1 << rx) | (x1 >> (64-rx)); \
+    z1 = (z1 << rz) | (z1 >> (64-rz)); \
+    x1 ^= x0; \
+    z1 ^= z0; \
+
+
+namespace sitmo {
+
+    // enable_if for C__98 compilers
+    template<bool C, typename T = void>
+    struct sitmo_enable_if { typedef T type; };
+
+    template<typename T>
+    struct sitmo_enable_if<false, T> { };
+
+    // SFINAE check for the existence of a "void generate(int*,int*)"member function
+    template<typename T>
+    struct has_generate_template
+    {
+        typedef char (&Two)[2];;
+        template<typename F, void (F::*)(int *, int *)> struct helper {};
+        template<typename C> static char test(helper<C, &C::template generate<int*> >*);
+        template<typename C> static Two test(...);
+        static bool const value = sizeof(test<T>(0)) == sizeof(char);
+    };
+
+
+class prng_engine
+{
+public:
+    // "req" are requirements as stated in the C++ 11 draft n3242=11-0012
+    //
+    // req: 26.5.1.3 Uniform random number generator requirements, p.906, table 116, row 1
+    typedef uint32_t result_type;
+
+    // req: 26.5.1.3 Uniform random number generator requirements, p.906, table 116, row 3 & 4
+#if __cplusplus <= 199711L
+    static result_type (min)() { return 0; }
+    static result_type (max)() { return 0xFFFFFFFF; }
+#else
+    static constexpr result_type (min)() { return 0; }
+    static constexpr result_type (max)() { return 0xFFFFFFFF; }
+#endif    
+
+    // -------------------------------------------------
+    // Constructors
+    // -------------------------------------------------
+
+    // req: 26.5.1.4 Random number engine requirements, p.907 table 117, row 1
+    // Creates an engine with the same initial state as all other
+    // default-constructed engines of type E.
+    prng_engine() 
+    {
+        seed();
+    }
+
+    // req: 26.5.1.4 Random number engine requirements, p.907 table 117, row 2
+    // Creates an engine that compares equal to x.
+    prng_engine(const prng_engine& x)
+    {
+        for (unsigned short i=0; i<4; ++i) {
+            _s[i] = x._s[i];
+            _k[i] = x._k[i];
+            _o[i] = x._o[i];
+        }
+        _o_counter = x._o_counter;
+    }
+    
+    
+    // req: 26.5.1.4 Random number engine requirements, p.907 table 117, row 3
+    // Creates an engine with initial O(size of state) state determined by s.
+    prng_engine(uint32_t s) 
+    {
+        seed(s);
+    }
+    
+        // req: 26.5.1.4 Random number engine requirements, p.908 table 117, row 4
+    // Creates an engine with an initial state that depends on a sequence
+    // produced by one call to q.generate.
+    template<class Seq> 
+    prng_engine(Seq& q, typename sitmo_enable_if< has_generate_template<Seq>::value >::type* = 0 )
+    {
+        seed(q);
+    }
+
+    // -------------------------------------------------
+    // Seeding
+    // -------------------------------------------------
+    
+    // req: 26.5.1.4 Random number engine requirements, p.908 table 117, row 5
+    void seed()
+    {
+        for (unsigned short i=0; i<4; ++i) {
+            _k[i] = 0;
+            _s[i] = 0;
+        }
+        _o_counter = 0;
+
+        _o[0] = 0x09218ebde6c85537;
+        _o[1] = 0x55941f5266d86105;
+        _o[2] = 0x4bd25e16282434dc;
+        _o[3] = 0xee29ec846bd2e40b;
+    }
+
+    // req: 26.5.1.4 Random number engine requirements, p.908 table 117, row 6
+    // s needs to be of return_type, which is uint32_t
+    void seed(uint32_t s)
+    { 
+        for (unsigned short i=0; i<4; ++i) {
+            _k[i] = 0;
+            _s[i] = 0;
+        }
+        _k[0] = s; 
+        _o_counter = 0;
+        
+        encrypt_counter();
+    }
+    
+    // req: 26.5.1.4 Random number engine requirements, p.908 table 117, row 7
+    template<class Seq> 
+    void seed(Seq& q, typename sitmo_enable_if< has_generate_template<Seq>::value >::type* = 0 )
+    {
+        typename Seq::result_type w[8];
+        q.generate(&w[0], &w[8]);
+
+        for (unsigned short i=0; i<4; ++i) {
+            _k[i] = ( static_cast<uint64_t>(w[2*i]) << 32) | w[2*i+1];
+            _s[i] = 0;
+        }
+        _o_counter = 0;
+        
+        encrypt_counter();
+    }
+    
+    // req: 26.5.1.4 Random number engine requirements, p.908 table 117, row 8
+    // Advances e’s state ei to ei+1 = TA(ei) and returns GA(ei).
+    uint32_t operator()()
+    {
+        // can we return a value from the current block?
+        if (_o_counter < 8) {
+            unsigned short _o_index = _o_counter >> 1;
+            _o_counter++;
+            if (_o_counter&1) 
+                return _o[_o_index] & 0xFFFFFFFF; 
+            else
+                return _o[_o_index] >> 32;
+        }
+        
+        // generate a new block and return the first 32 bits
+        inc_counter();
+        encrypt_counter();
+        _o_counter = 1; // the next call
+        return _o[0] & 0xFFFFFFFF;   // this call
+    }
+    
+    // -------------------------------------------------
+    // misc
+    // -------------------------------------------------
+    
+    // req: 26.5.1.4 Random number engine requirements, p.908 table 117, row 9
+    // Advances e’s state ei to ei+z by any means equivalent to z
+    // consecutive calls e().
+    void discard(uint64_t z)
+    {
+        // check if we stay in the current block
+        if (z < 8 - _o_counter) {
+            _o_counter += static_cast<unsigned short>(z);
+            return;
+        }
+
+        // we will have to generate a new block...
+        z -= (8 - _o_counter);  // discard the remainder of the current blok
+        _o_counter = z % 8;     // set the pointer in the correct element in the new block
+        z -= _o_counter;        // update z
+        z >>= 3;                // the number of buffers is elements/8
+        ++z;                    // and one more because we crossed the buffer line
+        inc_counter(z);
+        encrypt_counter();
+    }
+
+   // -------------------------------------------------
+   // IO
+   // -------------------------------------------------
+   template<class CharT, class Traits>
+   friend std::basic_ostream<CharT,Traits>&
+   operator<<(std::basic_ostream<CharT,Traits>& os, const prng_engine& s) {
+       for (unsigned short i=0; i<4; ++i)
+           os << s._k[i] << ' ' << s._s[i] << ' ' << s._o[i] << ' ';
+       os << s._o_counter;
+       return os;
+   }
+
+   template<class CharT, class Traits>
+   friend std::basic_istream<CharT,Traits>&
+   operator>>(std::basic_istream<CharT,Traits>& is, prng_engine& s) {
+       for (unsigned short i=0; i<4; ++i)
+           is >> s._k[i] >> s._s[i] >> s._o[i];
+       is >> s._o_counter;
+       return is;
+   } 
+    // req: 26.5.1.4 Random number engine requirements, p.908 table 117, row 10
+    // This operator is an equivalence relation. With Sx and Sy as the infinite 
+    // sequences of values that would be generated by repeated future calls to 
+    // x() and y(), respectively, returns true if Sx = Sy; else returns false.
+    bool operator==(const prng_engine& y) 
+    {
+        if (_o_counter != y._o_counter) return false;
+        for (unsigned short i=0; i<4; ++i) {
+            if (_s[i] != y._s[i]) return false;
+            if (_k[i] != y._k[i]) return false;
+            if (_o[i] != y._o[i]) return false;
+        }
+        return true;
+    }
+    
+    // req: 26.5.1.4 Random number engine requirements, p.908 table 117, row 11
+    bool operator!=(const prng_engine& y) 
+    {
+        return !(*this == y);
+    }
+    
+    // Extra function to set the key
+    void set_key(uint64_t k0=0, uint64_t k1=0, uint64_t k2=0, uint64_t k3=0)
+    {
+        _k[0] = k0; _k[1] = k1; _k[2] = k2; _k[3] = k3;
+        encrypt_counter();
+    }
+    
+    // set the counter
+    void set_counter(uint64_t s0=0, uint64_t s1=0, uint64_t s2=0, uint64_t s3=0, unsigned short o_counter=0)
+    {
+        _s[0] = s0; 
+        _s[1] = s1; 
+        _s[2] = s2; 
+        _s[3] = s3;
+        _o_counter = o_counter % 8;
+        encrypt_counter();
+    }
+    
+    
+    // versioning
+    uint32_t version()
+    {
+        return 5;
+    }
+    
+private:
+    void encrypt_counter()
+    {
+        uint64_t b[4];
+        uint64_t k[5];
+
+        for (unsigned short i=0; i<4; ++i) b[i] = _s[i];
+        for (unsigned short i=0; i<4; ++i) k[i] = _k[i];
+
+        k[4] = 0x1BD11BDAA9FC1A22 ^ k[0] ^ k[1] ^ k[2] ^ k[3];
+
+        MIXK(b[0], b[1], 14,   b[2], b[3], 16,   k[0], k[1], k[2], k[3]);
+        MIX2(b[0], b[3], 52,   b[2], b[1], 57);
+        MIX2(b[0], b[1], 23,   b[2], b[3], 40);
+        MIX2(b[0], b[3],  5,   b[2], b[1], 37);
+        MIXK(b[0], b[1], 25,   b[2], b[3], 33,   k[1], k[2], k[3], k[4]+1);
+        MIX2(b[0], b[3], 46,   b[2], b[1], 12);
+        MIX2(b[0], b[1], 58,   b[2], b[3], 22);
+        MIX2(b[0], b[3], 32,   b[2], b[1], 32);
+
+        MIXK(b[0], b[1], 14,   b[2], b[3], 16,   k[2], k[3], k[4], k[0]+2);
+        MIX2(b[0], b[3], 52,   b[2], b[1], 57);
+        MIX2(b[0], b[1], 23,   b[2], b[3], 40);
+        MIX2(b[0], b[3],  5,   b[2], b[1], 37);
+        MIXK(b[0], b[1], 25,   b[2], b[3], 33,   k[3], k[4], k[0], k[1]+3);
+
+        MIX2(b[0], b[3], 46,   b[2], b[1], 12);
+        MIX2(b[0], b[1], 58,   b[2], b[3], 22);
+        MIX2(b[0], b[3], 32,   b[2], b[1], 32);
+
+        MIXK(b[0], b[1], 14,   b[2], b[3], 16,   k[4], k[0], k[1], k[2]+4);
+        MIX2(b[0], b[3], 52,   b[2], b[1], 57);
+        MIX2(b[0], b[1], 23,   b[2], b[3], 40);
+        MIX2(b[0], b[3],  5,   b[2], b[1], 37);
+
+        for (unsigned int i=0; i<4; ++i) _o[i] = b[i] + k[i];
+        _o[3] += 5;
+    }
+    
+    void inc_counter()
+    {
+        ++_s[0]; 
+        if (_s[0] != 0) return;
+        
+        ++_s[1]; 
+        if (_s[1] != 0) return;
+        
+        ++_s[2]; 
+        if (_s[2] != 0) return;
+        
+        ++_s[3];
+    }
+    
+    void inc_counter(uint64_t z)
+    {
+        if (z > 0xFFFFFFFFFFFFFFFF - _s[0]) {   // check if we will overflow the first 64 bit int
+            ++_s[1];
+            if (_s[1] == 0) {
+                ++_s[2];
+                if (_s[2] == 0) {
+                    ++_s[3];
+                }
+            }
+        }
+        _s[0] += z;
+    }
+    
+private:
+    uint64_t _k[4];             // key
+    uint64_t _s[4];             // state (counter)
+    uint64_t _o[4];             // cipher output    4 * 64 bit = 256 bit output
+    unsigned short _o_counter;  // output chunk counter, the 256 random bits in _o 
+                                // are returned in eight 32 bit chunks
+};
+
+
+} // namespace sitmo
+
+#undef MIXK
+#undef MIX2
+
+#endif

lib/stencil/Lebesgue.cc

@@ -26,7 +26,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include <Grid.h>
+#include <Grid/Grid.h>
 #include <algorithm>
 
 namespace Grid {

lib/stencil/Stencil_common.cc

@@ -26,7 +26,7 @@ Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include "Grid.h"
+#include <Grid/Grid.h>
 
 namespace Grid {
 }

scripts/filelist

@@ -4,9 +4,8 @@ home=`pwd`
 
 # library Make.inc
 cd $home/lib
-HFILES=`find . -type f -name '*.h' -not -path '*/Old/*' -not -path '*/Eigen/*'`
+HFILES=`find . -type f -name '*.h' -not -name '*Hdf5*' -not -path '*/Old/*' -not -path '*/Eigen/*'`
-HFILES="$HFILES"
+CCFILES=`find . -type f -name '*.cc' -not -name '*Communicator*.cc' -not -name '*Hdf5*'`
-CCFILES=`find . -type f -name '*.cc' -not  -name '*ommunicator*.cc'`
 echo HFILES=$HFILES > Make.inc
 echo >> Make.inc
 echo CCFILES=$CCFILES >> Make.inc

tests/IO/Test_serialisation.cc

@@ -28,130 +28,152 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     /*  END LEGAL */
 #include <Grid/Grid.h>
 
-namespace Grid {
-  
-  GRID_SERIALIZABLE_ENUM(myenum, undef, red, 1, blue, 2, green, 3);
-    
-  class myclass: Serializable {
-  public:
-    
-    GRID_SERIALIZABLE_CLASS_MEMBERS(myclass,
-                            myenum, e,
-                            std::vector<myenum>, ve,
-                            std::string, name,
-                            int, x,
-                            double, y,
-                            bool , b,
-                            std::vector<double>, array,
-                            std::vector<std::vector<double>>, twodimarray,
-                            );
-    
-    myclass() {}
-    myclass(int i)
-    : array(4,5.1), twodimarray(3,std::vector<double>(2,1.23456)), ve(2, myenum::blue)
-    {
-      e=myenum::red;
-      x=i;
-      y=2*i;
-      b=true;
-      name="bother said pooh";
-    }
-  };
-  
-}
-
 using namespace Grid;
 
-int16_t i16 = 1;
+GRID_SERIALIZABLE_ENUM(myenum, undef, red, 1, blue, 2, green, 3);
+  
+class myclass: Serializable {
+public:
+  GRID_SERIALIZABLE_CLASS_MEMBERS(myclass,
+                          myenum, e,
+                          std::vector<myenum>, ve,
+                          std::string, name,
+                          int, x,
+                          double, y,
+                          bool , b,
+                          std::vector<double>, array,
+                          std::vector<std::vector<double>>, twodimarray,
+                          std::vector<std::vector<std::vector<Complex>>>, cmplx3darray
+                          );
+  myclass() {}
+  myclass(int i)
+  : array(4,5.1)
+  , twodimarray(3,std::vector<double>(5, 1.23456))
+  , cmplx3darray(3,std::vector<std::vector<Complex>>(5, std::vector<Complex>(7, Complex(1.2, 3.4))))
+  , ve(2, myenum::blue)
+  {
+    e=myenum::red;
+    x=i;
+    y=2*i;
+    b=true;
+    name="bother said pooh";
+  }
+};
+
+int16_t  i16 = 1;
 uint16_t u16 = 2;
-int32_t i32 = 3;
+int32_t  i32 = 3;
 uint32_t u32 = 4;
-int64_t i64 = 5;
+int64_t  i64 = 5;
 uint64_t u64 = 6;
-float    f = M_PI;
+float    f   = M_PI;
-double   d = 2*M_PI;
+double   d   = 2*M_PI;
-bool     b = false;
+bool     b   = false;
+
+template <typename W, typename R, typename O>
+void ioTest(const std::string &filename, const O &object, const std::string &name)
+{
+  // writer needs to be destroyed so that writing physically happens
+  {
+    W writer(filename);
+    
+    write(writer, "testobject", object);
+  }
+  
+  R    reader(filename);
+  O    buf;
+  bool good;
+  
+  read(reader, "testobject", buf);
+  good = (object == buf);
+  std::cout << name << " IO test: " << (good ? "success" : "failure");
+  std::cout << std::endl;
+  if (!good) exit(EXIT_FAILURE);
+}
 
 int main(int argc,char **argv)
 {
-  {
+  std::cout << "==== basic IO" << std::endl;
-    XmlWriter WR("bother.xml");
+  XmlWriter WR("bother.xml");
-    
+  
-    // test basic type writing
+  // test basic type writing
-    push(WR,"BasicTypes");
+  std::cout << "-- basic writing to 'bother.xml'..." << std::endl;
-    write(WR,std::string("i16"),i16);
+  push(WR,"BasicTypes");
-    write(WR,"u16",u16);
+  write(WR,std::string("i16"),i16);
-    write(WR,"i32",i32);
+  write(WR,"u16",u16);
-    write(WR,"u32",u32);
+  write(WR,"i32",i32);
-    write(WR,"i64",i64);
+  write(WR,"u32",u32);
-    write(WR,"u64",u64);
+  write(WR,"i64",i64);
-    write(WR,"f",f);
+  write(WR,"u64",u64);
-    write(WR,"d",d);
+  write(WR,"f",f);
-    write(WR,"b",b);
+  write(WR,"d",d);
-    pop(WR);
+  write(WR,"b",b);
-    
+  pop(WR);
-    // test serializable class writing
+  
-    myclass obj(1234); // non-trivial constructor
+  // test serializable class writing
-    write(WR,"obj",obj);
+  myclass              obj(1234); // non-trivial constructor
-    WR.write("obj2", obj);
+  std::vector<myclass> vec;
-    std::cout << obj << std::endl;
+  
-    
+  std::cout << "-- serialisable class writing to 'bother.xml'..." << std::endl;
-    std::vector<myclass> vec;
+  write(WR,"obj",obj);
-    vec.push_back(myclass(1234));
+  WR.write("obj2", obj);
-    vec.push_back(myclass(5678));
+  vec.push_back(myclass(1234));
-    vec.push_back(myclass(3838));
+  vec.push_back(myclass(5678));
-    write(WR, "objvec", vec);
+  vec.push_back(myclass(3838));
-  };
+  write(WR, "objvec", vec);
+  std::cout << "-- serialisable class writing to std::cout:" << std::endl;
+  std::cout << obj << std::endl;
+  std::cout << "-- serialisable class comparison:" << std::endl;
+  std::cout << "vec[0] == obj: " << ((vec[0] == obj) ? "true" : "false") << std::endl;
+  std::cout << "vec[1] == obj: " << ((vec[1] == obj) ? "true" : "false") << std::endl;
   
   // read tests
-  myclass copy1, copy2, copy3;
+  std::cout << "\n==== IO self-consistency tests" << std::endl;
-  std::vector<myclass> veccopy1, veccopy2, veccopy3;
   //// XML
-  {
+  ioTest<XmlWriter, XmlReader>("iotest.xml", obj, "XML    (object)           ");
-    XmlReader RD("bother.xml");
+  ioTest<XmlWriter, XmlReader>("iotest.xml", vec, "XML    (vector of objects)");
-    read(RD,"obj",copy1);
-    read(RD,"objvec", veccopy1);
-    std::cout << "Loaded (XML) -----------------" << std::endl;
-    std::cout << copy1 << std::endl << veccopy1 << std::endl;
-  }
   //// binary
-  {
+  ioTest<BinaryWriter, BinaryReader>("iotest.bin", obj, "binary (object)           ");
-    BinaryWriter BWR("bother.bin");
+  ioTest<BinaryWriter, BinaryReader>("iotest.bin", vec, "binary (vector of objects)");
-    write(BWR,"discard",copy1 );
-    write(BWR,"discard",veccopy1 );
-  }
-  {
-    BinaryReader BRD("bother.bin");
-    read (BRD,"discard",copy2 );
-    read (BRD,"discard",veccopy2 );
-    std::cout << "Loaded (bin) -----------------" << std::endl;
-    std::cout << copy2 << std::endl << veccopy2 << std::endl;
-  }
   //// text
-  {
+  ioTest<TextWriter, TextReader>("iotest.dat", obj, "text   (object)           ");
-    TextWriter TWR("bother.txt");
+  ioTest<TextWriter, TextReader>("iotest.dat", vec, "text   (vector of objects)");
-    write(TWR,"discard",copy1 );
+  //// HDF5
-    write(TWR,"discard",veccopy1 );
+#ifdef HAVE_HDF5
-  }
+  ioTest<Hdf5Writer, Hdf5Reader>("iotest.h5", obj, "HDF5   (object)           ");
-  {
+  ioTest<Hdf5Writer, Hdf5Reader>("iotest.h5", vec, "HDF5   (vector of objects)");
-    TextReader TRD("bother.txt");
+#endif
-    read (TRD,"discard",copy3 );
-    read (TRD,"discard",veccopy3 );
-    std::cout << "Loaded (txt) -----------------" << std::endl;
-    std::cout << copy3 << std::endl << veccopy3 << std::endl;
-  }
   
-  std::vector<int> iv = strToVec<int>("1 2 2 4");
+  std::cout << "\n==== vector flattening/reconstruction" << std::endl;
-  std::vector<std::string> sv = strToVec<std::string>("bli bla blu");
+  typedef std::vector<std::vector<std::vector<double>>> vec3d;
   
-  for (auto &e: iv)
+  vec3d dv, buf;
+  double d = 0.;
+  
+  dv.resize(4);
+  for (auto &v1: dv)
   {
-    std::cout << e << " ";
+    v1.resize(3);
+    for (auto &v2: v1)
+    {
+      v2.resize(5);
+      for (auto &x: v2)
+      {
+        x = d++;
+      }
+    }
   }
-  std::cout << std::endl;
+  std::cout << "original 3D vector:" << std::endl;
-  for (auto &e: sv)
+  std::cout << dv << std::endl;
-  {
+  
-    std::cout << e << " ";
+  Flatten<vec3d> flatdv(dv);
-  }
+  
-  std::cout << std::endl;
+  std::cout << "\ndimensions:" << std::endl;
+  std::cout << flatdv.getDim() << std::endl;
+  std::cout << "\nflattened vector:" << std::endl;
+  std::cout << flatdv.getFlatVector() << std::endl;
+  
+  Reconstruct<vec3d> rec(flatdv.getFlatVector(), flatdv.getDim());
+  std::cout << "\nreconstructed vector:" << std::endl;
+  std::cout << flatdv.getVector() << std::endl;
 }

tests/hadrons/Makefile.am

@@ -1,3 +1,3 @@
-AM_LDFLAGS += -L../../extras/hadrons
+AM_LDFLAGS += -L../../extras/Hadrons
 
 include Make.inc

tests/qdpxx/Test_qdpxx_loops_staples.cc

@@ -25,7 +25,7 @@ Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 double calc_grid_p      (Grid::QCD::LatticeGaugeField & lat);
 double calc_chroma_p    (Grid::QCD::LatticeGaugeField & lat);

tests/qdpxx/Test_qdpxx_munprec.cc

@@ -26,7 +26,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 int    Ls=8;
 double M5=1.6;