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

.gitignore

@@ -88,6 +88,7 @@ Thumbs.db
 # build directory #
 ###################
 build*/*
+Documentation/_build
 
 # IDE related files #
 #####################

Grid/serialisation/BaseIO.cc

@@ -0,0 +1,35 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/serialisation/BaseIO.h
+
+Copyright (C) 2015
+
+Author: Michael Marshall <michael.marshall@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/GridCore.h>
+
+NAMESPACE_BEGIN(Grid)
+
+std::uint64_t EigenIO::EigenResizeCounter(0);
+
+NAMESPACE_END(Grid)

Grid/serialisation/BaseIO.h

@@ -9,6 +9,7 @@
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 Author: Guido Cossu <guido.cossu@ed.ac.uk>
+Author: Michael Marshall <michael.marshall@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
@@ -30,6 +31,7 @@ Author: Guido Cossu <guido.cossu@ed.ac.uk>
 #ifndef GRID_SERIALISATION_ABSTRACT_READER_H
 #define GRID_SERIALISATION_ABSTRACT_READER_H
 
+#include <atomic>
 #include <type_traits>
 #include <Grid/tensors/Tensors.h>
 #include <Grid/serialisation/VectorUtils.h>
@@ -110,6 +112,10 @@ namespace Grid {
     template <typename ET>
     inline typename std::enable_if<is_tensor_of_container<ET>::value, typename Traits<ET>::scalar_type *>::type
     getFirstScalar(ET &eigenTensor) { return eigenTensor.data()->begin(); }
+
+    // Counter for resized EigenTensors (poor man's substitute for allocator)
+    // Defined in BinaryIO.cc
+    extern std::uint64_t EigenResizeCounter;
   }
 
   // Abstract writer/reader classes ////////////////////////////////////////////
@@ -497,8 +503,14 @@ namespace Grid {
   typename std::enable_if<EigenIO::is_tensor_variable<ETensor>::value, void>::type
   Reader<T>::Reshape(ETensor &t, const std::array<typename ETensor::Index, ETensor::NumDimensions> &dims )
   {
+#ifdef GRID_OMP
+    // The memory counter is the reason this must be done from the primary thread
+    assert(omp_in_parallel()==0 && "Deserialisation which resizes Eigen tensor must happen from primary thread");
+#endif
+    EigenIO::EigenResizeCounter -= static_cast<uint64_t>(t.size()) * sizeof(typename ETensor::Scalar);
     //t.reshape( dims );
     t.resize( dims );
+    EigenIO::EigenResizeCounter += static_cast<uint64_t>(t.size()) * sizeof(typename ETensor::Scalar);
   }
 
   template <typename T>

Grid/serialisation/Hdf5IO.cc

@@ -1,3 +1,34 @@
+/*************************************************************************************
+ 
+ Grid physics library, www.github.com/paboyle/Grid
+ 
+ Source file: ./Grid/serialisation/VectorUtils.h
+ 
+ Copyright (C) 2015
+ 
+ Author: Antonin Portelli <antonin.portelli@me.com>
+ Author: Peter Boyle <paboyle@ed.ac.uk>
+ Author: Guido Cossu <guido.cossu@ed.ac.uk>
+ Author: Michael Marshall <michael.marshall@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 Grid;

Grid/serialisation/Hdf5IO.h

@@ -1,3 +1,34 @@
+/*************************************************************************************
+ 
+ Grid physics library, www.github.com/paboyle/Grid
+ 
+ Source file: ./Grid/serialisation/VectorUtils.h
+ 
+ Copyright (C) 2015
+ 
+ Author: Peter Boyle <paboyle@ed.ac.uk>
+ Author: Antonin Portelli <antonin.portelli@me.com>
+ Author: Guido Cossu <guido.cossu@ed.ac.uk>
+ Author: Michael Marshall <michael.marshall@ed.ac.uk>
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+ 
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ GNU General Public License for more details.
+ 
+ You should have received a copy of the GNU General Public License along
+ with this program; if not, write to the Free Software Foundation, Inc.,
+ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ 
+ See the full license in the file "LICENSE" in the top level distribution directory
+ *************************************************************************************/
+/*  END LEGAL */
+
 #ifndef GRID_SERIALISATION_HDF5_H
 #define GRID_SERIALISATION_HDF5_H
 
@@ -34,11 +65,13 @@ namespace Grid
     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
+    void writeRagged(const std::string &s, const std::vector<U> &x);
+    template <typename U>
+    typename std::enable_if<is_flattenable<std::vector<U>>::value>::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
+    typename std::enable_if<!is_flattenable<std::vector<U>>::value>::type
-    writeDefault(const std::string &s, const std::vector<U> &x);
+    writeDefault(const std::string &s, const std::vector<U> &x) { writeRagged(s, x); }
     template <typename U>
     void writeMultiDim(const std::string &s, const std::vector<size_t> & Dimensions, const U * pDataRowMajor, size_t NumElements);
     H5NS::Group & getGroup(void);
@@ -64,11 +97,13 @@ namespace Grid
     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
+    void readRagged(const std::string &s, std::vector<U> &x);
+    template <typename U>
+    typename std::enable_if<is_flattenable<std::vector<U>>::value>::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
+    typename std::enable_if<!is_flattenable<std::vector<U>>::value>::type
-    readDefault(const std::string &s, std::vector<U> &x);
+    readDefault(const std::string &s, std::vector<U> &x) { readRagged(s, x); }
     template <typename U>
     void readMultiDim(const std::string &s, std::vector<U> &buf, std::vector<size_t> &dim);
     H5NS::Group & getGroup(void);
@@ -176,11 +211,13 @@ namespace Grid
   }
 
   template <typename U>
-  typename std::enable_if<element<std::vector<U>>::is_number, void>::type
+  typename std::enable_if<is_flattenable<std::vector<U>>::value>::type
   Hdf5Writer::writeDefault(const std::string &s, const std::vector<U> &x)
+  {
+    if (isRegularShape(x))
     {
       // alias to element type
-    typedef typename element<std::vector<U>>::type Element;
+      using Scalar = typename is_flattenable<std::vector<U>>::type;
       
       // flatten the vector and getting dimensions
       Flatten<std::vector<U>> flat(x);
@@ -188,12 +225,16 @@ namespace Grid
       const auto           &flatx = flat.getFlatVector();
       for (auto &d: flat.getDim())
         dim.push_back(d);
-    writeMultiDim<Element>(s, dim, &flatx[0], flatx.size());
+      writeMultiDim<Scalar>(s, dim, &flatx[0], flatx.size());
+    }
+    else
+    {
+      writeRagged(s, x);
+    }
   }
   
   template <typename U>
-  typename std::enable_if<!element<std::vector<U>>::is_number, void>::type
+  void Hdf5Writer::writeRagged(const std::string &s, const std::vector<U> &x)
-  Hdf5Writer::writeDefault(const std::string &s, const std::vector<U> &x)
   {
     push(s);
     writeSingleAttribute(x.size(), HDF5_GRID_GUARD "vector_size",
@@ -229,7 +270,7 @@ namespace Grid
   void Hdf5Reader::readMultiDim(const std::string &s, std::vector<U> &buf, std::vector<size_t> &dim)
   {
     // alias to element type
-    typedef typename element<std::vector<U>>::type Element;
+    using Scalar = typename is_flattenable<std::vector<U>>::type;
     
     // read the dimensions
     H5NS::DataSpace       dataSpace;
@@ -260,26 +301,33 @@ namespace Grid
       H5NS::DataSet dataSet;
       
       dataSet = group_.openDataSet(s);
-      dataSet.read(buf.data(), Hdf5Type<Element>::type());
+      dataSet.read(buf.data(), Hdf5Type<Scalar>::type());
     }
     else
     {
       H5NS::Attribute attribute;
       
       attribute = group_.openAttribute(s);
-      attribute.read(Hdf5Type<Element>::type(), buf.data());
+      attribute.read(Hdf5Type<Scalar>::type(), buf.data());
     }
   }
 
   template <typename U>
-  typename std::enable_if<element<std::vector<U>>::is_number, void>::type
+  typename std::enable_if<is_flattenable<std::vector<U>>::value>::type
   Hdf5Reader::readDefault(const std::string &s, std::vector<U> &x)
+  {
+    if (H5Lexists        (group_.getId(), s.c_str(), H5P_DEFAULT) > 0
+     && H5Aexists_by_name(group_.getId(), s.c_str(), HDF5_GRID_GUARD "vector_size", H5P_DEFAULT ) > 0)
+    {
+      readRagged(s, x);
+    }
+    else
     {
       // alias to element type
-    typedef typename element<std::vector<U>>::type Element;
+      using Scalar = typename is_flattenable<std::vector<U>>::type;
 
       std::vector<size_t>   dim;
-    std::vector<Element>  buf;
+      std::vector<Scalar>   buf;
       readMultiDim( s, buf, dim );
 
       // reconstruct the multidimensional vector
@@ -287,10 +335,10 @@ namespace Grid
 
       x = r.getVector();
     }
+  }
   
   template <typename U>
-  typename std::enable_if<!element<std::vector<U>>::is_number, void>::type
+  void Hdf5Reader::readRagged(const std::string &s, std::vector<U> &x)
-  Hdf5Reader::readDefault(const std::string &s, std::vector<U> &x)
   {
     uint64_t size;
     

Grid/serialisation/MacroMagic.h

@@ -118,13 +118,13 @@ static inline std::string SerialisableClassName(void) {return std::string(#cname
 static constexpr bool isEnum = false; \
 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){ \
+static inline void write(::Grid::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){	\
+static inline void read(::Grid::Reader<T> &RD,const std::string &s, cname &obj){	\
   if (!push(RD,s))\
   {\
     std::cout << ::Grid::GridLogWarning << "IO: Cannot open node '" << s << "'" << std::endl; \

Grid/serialisation/VectorUtils.h

@@ -9,6 +9,7 @@
  Author: Antonin Portelli <antonin.portelli@me.com>
  Author: Peter Boyle <paboyle@ph.ed.ac.uk>
  Author: paboyle <paboyle@ph.ed.ac.uk>
+ Author: Michael Marshall <michael.marshall@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
@@ -236,21 +237,36 @@ namespace Grid {
     }
   }
 
-  // Vector element trait //////////////////////////////////////////////////////  
+  // is_flattenable<T>::value is true if T is a std::vector<> which can be flattened //////////////////////
-  template <typename T>
+  template <typename T, typename V = void>
-  struct element
+  struct is_flattenable : std::false_type
   {
-    typedef T type;
+    using type      = T;
-    static constexpr bool is_number = false;
+    using grid_type = T;
+    static constexpr int vecRank = 0;
+    static constexpr bool isGridTensor = false;
+    static constexpr bool children_flattenable = std::is_arithmetic<T>::value or is_complex<T>::value;
   };
 
   template <typename T>
-  struct element<std::vector<T>>
+  struct is_flattenable<T, typename std::enable_if<isGridTensor<T>::value>::type> : std::false_type
   {
-    typedef typename element<T>::type type;
+    using type      = typename GridTypeMapper<T>::scalar_type;
-    static constexpr bool is_number = std::is_arithmetic<T>::value
+    using grid_type = T;
-                                      or is_complex<T>::value
+    static constexpr int vecRank = 0;
-                                      or element<T>::is_number;
+    static constexpr bool isGridTensor = true;
+    static constexpr bool children_flattenable = true;
+  };
+
+  template <typename T>
+  struct is_flattenable<std::vector<T>, typename std::enable_if<is_flattenable<T>::children_flattenable>::type>
+  : std::true_type
+  {
+    using type      = typename is_flattenable<T>::type;
+    using grid_type = typename is_flattenable<T>::grid_type;
+    static constexpr bool isGridTensor = is_flattenable<T>::isGridTensor;
+    static constexpr int vecRank = is_flattenable<T>::vecRank + 1;
+    static constexpr bool children_flattenable = true;
   };
   
   // Vector flattening utility class ////////////////////////////////////////////
@@ -259,22 +275,29 @@ namespace Grid {
   class Flatten
   {
   public:
-    typedef typename element<V>::type Element;
+    using Scalar  = typename is_flattenable<V>::type;
+    static constexpr bool isGridTensor = is_flattenable<V>::isGridTensor;
   public:
     explicit                    Flatten(const V &vector);
-    const V &                    getVector(void);
+    const V &                   getVector(void)     const { return vector_; }
-    const std::vector<Element> & getFlatVector(void);
+    const std::vector<Scalar> & getFlatVector(void) const { return flatVector_; }
-    const std::vector<size_t>  & getDim(void);
+    const std::vector<size_t> & getDim(void)        const { return dim_; }
   private:
-    void accumulate(const Element &e);
+    template <typename W> typename std::enable_if<!is_flattenable<W>::value && !is_flattenable<W>::isGridTensor>::type
-    template <typename W>
+    accumulate(const W &e);
-    void accumulate(const W &v);
+    template <typename W> typename std::enable_if<!is_flattenable<W>::value &&  is_flattenable<W>::isGridTensor>::type
-    void accumulateDim(const Element &e);
+    accumulate(const W &e);
-    template <typename W>
+    template <typename W> typename std::enable_if< is_flattenable<W>::value>::type
-    void accumulateDim(const W &v);
+    accumulate(const W &v);
+    template <typename W> typename std::enable_if<!is_flattenable<W>::value && !is_flattenable<W>::isGridTensor>::type
+    accumulateDim(const W &e) {} // Innermost is a scalar - do nothing
+    template <typename W> typename std::enable_if<!is_flattenable<W>::value &&  is_flattenable<W>::isGridTensor>::type
+    accumulateDim(const W &e);
+    template <typename W> typename std::enable_if< is_flattenable<W>::value>::type
+    accumulateDim(const W &v);
   private:
     const V             &vector_;
-    std::vector<Element> flatVector_;
+    std::vector<Scalar> flatVector_;
     std::vector<size_t> dim_;
   };
   
@@ -283,23 +306,32 @@ namespace Grid {
   class Reconstruct
   {
   public:
-    typedef typename element<V>::type Element;
+    using Scalar  = typename is_flattenable<V>::type;
+    static constexpr bool isGridTensor = is_flattenable<V>::isGridTensor;
   public:
-    Reconstruct(const std::vector<Element> &flatVector,
+    Reconstruct(const std::vector<Scalar> &flatVector,
                 const std::vector<size_t> &dim);
-    const V &                    getVector(void);
+    const V &                   getVector(void)     const { return vector_; }
-    const std::vector<Element> & getFlatVector(void);
+    const std::vector<Scalar> & getFlatVector(void) const { return flatVector_; }
-    const std::vector<size_t>  & getDim(void);
+    const std::vector<size_t> & getDim(void)        const { return dim_; }
   private:
-    void fill(std::vector<Element> &v);
+    template <typename W> typename std::enable_if<!is_flattenable<W>::value && !is_flattenable<W>::isGridTensor>::type
-    template <typename W>
+    fill(W &v);
-    void fill(W &v);
+    template <typename W> typename std::enable_if<!is_flattenable<W>::value &&  is_flattenable<W>::isGridTensor>::type
-    void resize(std::vector<Element> &v, const unsigned int dim);
+    fill(W &v);
-    template <typename W>
+    template <typename W> typename std::enable_if< is_flattenable<W>::value>::type
-    void resize(W &v, const unsigned int dim);
+    fill(W &v);
+    template <typename W> typename std::enable_if< is_flattenable<W>::value &&  is_flattenable<W>::vecRank==1>::type
+    resize(W &v, const unsigned int dim);
+    template <typename W> typename std::enable_if< is_flattenable<W>::value && (is_flattenable<W>::vecRank>1)>::type
+    resize(W &v, const unsigned int dim);
+    template <typename W> typename std::enable_if<!is_flattenable<W>::isGridTensor>::type
+    checkInnermost(const W &e) {} // Innermost is a scalar - do nothing
+    template <typename W> typename std::enable_if< is_flattenable<W>::isGridTensor>::type
+    checkInnermost(const W &e);
   private:
     V                         vector_;
-    const std::vector<Element> &flatVector_;
+    const std::vector<Scalar> &flatVector_;
     std::vector<size_t>       dim_;
     size_t                    ind_{0};
     unsigned int              dimInd_{0};
@@ -307,14 +339,24 @@ namespace Grid {
 
   // Flatten class template implementation
   template <typename V>
-  void Flatten<V>::accumulate(const Element &e)
+  template <typename W> typename std::enable_if<!is_flattenable<W>::value && !is_flattenable<W>::isGridTensor>::type
+  Flatten<V>::accumulate(const W &e)
   {
     flatVector_.push_back(e);
   }
   
   template <typename V>
-  template <typename W>
+  template <typename W> typename std::enable_if<!is_flattenable<W>::value && is_flattenable<W>::isGridTensor>::type
-  void Flatten<V>::accumulate(const W &v)
+  Flatten<V>::accumulate(const W &e)
+  {
+    for (const Scalar &x: e) {
+      flatVector_.push_back(x);
+    }
+  }
+
+  template <typename V>
+  template <typename W> typename std::enable_if<is_flattenable<W>::value>::type
+  Flatten<V>::accumulate(const W &v)
   {
     for (auto &e: v)
     {
@@ -323,11 +365,17 @@ namespace Grid {
   }
   
   template <typename V>
-  void Flatten<V>::accumulateDim(const Element &e) {};
+  template <typename W> typename std::enable_if<!is_flattenable<W>::value && is_flattenable<W>::isGridTensor>::type
+  Flatten<V>::accumulateDim(const W &e)
+  {
+    using Traits = GridTypeMapper<typename is_flattenable<W>::grid_type>;
+    for (int rank=0; rank < Traits::Rank; ++rank)
+      dim_.push_back(Traits::Dimension(rank));
+  }
   
   template <typename V>
-  template <typename W>
+  template <typename W> typename std::enable_if<is_flattenable<W>::value>::type
-  void Flatten<V>::accumulateDim(const W &v)
+  Flatten<V>::accumulateDim(const W &v)
   {
     dim_.push_back(v.size());
     accumulateDim(v[0]);
@@ -337,32 +385,26 @@ namespace Grid {
   Flatten<V>::Flatten(const V &vector)
   : vector_(vector)
   {
-    accumulate(vector_);
     accumulateDim(vector_);
+    std::size_t TotalSize{ dim_[0] };
+    for (int i = 1; i < dim_.size(); ++i) {
+      TotalSize *= dim_[i];
     }
-  
+    flatVector_.reserve(TotalSize);
-  template <typename V>
+    accumulate(vector_);
-  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)
+  template <typename W> typename std::enable_if<!is_flattenable<W>::value && !is_flattenable<W>::isGridTensor>::type
+  Reconstruct<V>::fill(W &v)
+  {
+    v = flatVector_[ind_++];
+  }
+  
+  template <typename V>
+  template <typename W> typename std::enable_if<!is_flattenable<W>::value &&  is_flattenable<W>::isGridTensor>::type
+  Reconstruct<V>::fill(W &v)
   {
     for (auto &e: v)
     {
@@ -371,8 +413,8 @@ namespace Grid {
   }
 
   template <typename V>
-  template <typename W>
+  template <typename W> typename std::enable_if<is_flattenable<W>::value>::type
-  void Reconstruct<V>::fill(W &v)
+  Reconstruct<V>::fill(W &v)
   {
     for (auto &e: v)
     {
@@ -381,14 +423,15 @@ namespace Grid {
   }
   
   template <typename V>
-  void Reconstruct<V>::resize(std::vector<Element> &v, const unsigned int dim)
+  template <typename W> typename std::enable_if<is_flattenable<W>::value && is_flattenable<W>::vecRank==1>::type
+  Reconstruct<V>::resize(W &v, const unsigned int dim)
   {
     v.resize(dim_[dim]);
   }
   
   template <typename V>
-  template <typename W>
+  template <typename W> typename std::enable_if<is_flattenable<W>::value && (is_flattenable<W>::vecRank>1)>::type
-  void Reconstruct<V>::resize(W &v, const unsigned int dim)
+  Reconstruct<V>::resize(W &v, const unsigned int dim)
   {
     v.resize(dim_[dim]);
     for (auto &e: v)
@@ -398,34 +441,31 @@ namespace Grid {
   }
   
   template <typename V>
-  Reconstruct<V>::Reconstruct(const std::vector<Element> &flatVector,
+  template <typename W> typename std::enable_if<is_flattenable<W>::isGridTensor>::type
+  Reconstruct<V>::checkInnermost(const W &)
+  {
+    using Traits = GridTypeMapper<typename is_flattenable<W>::grid_type>;
+    const int gridRank{Traits::Rank};
+    const int dimRank{static_cast<int>(dim_.size())};
+    assert(dimRank >= gridRank && "Tensor rank too low for Grid tensor");
+    for (int i=0; i<gridRank; ++i) {
+      assert(dim_[dimRank - gridRank + i] == Traits::Dimension(i) && "Tensor dimension doesn't match Grid tensor");
+    }
+    dim_.resize(dimRank - gridRank);
+  }
+
+  template <typename V>
+  Reconstruct<V>::Reconstruct(const std::vector<Scalar> &flatVector,
                               const std::vector<size_t> &dim)
   : flatVector_(flatVector)
   , dim_(dim)
   {
+    checkInnermost(vector_);
+    assert(dim_.size() == is_flattenable<V>::vecRank && "Tensor rank doesn't match nested std::vector rank");
     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_;
-  }
-
   // Vector IO utilities ///////////////////////////////////////////////////////
   // helper function to read space-separated values
   template <typename T>
@@ -459,6 +499,64 @@ namespace Grid {
     
     return os;
   }
+
+  // In general, scalar types are considered "flattenable" (regularly shaped)
+  template <typename T>
+  bool isRegularShapeHelper(const std::vector<T> &, std::vector<std::size_t> &, int, bool)
+  {
+    return true;
+  }
+
+  template <typename T>
+  bool isRegularShapeHelper(const std::vector<std::vector<T>> &v, std::vector<std::size_t> &Dims, int Depth, bool bFirst)
+  {
+    if( bFirst)
+    {
+      assert( Dims.size() == Depth     && "Bug: Delete this message after testing" );
+      Dims.push_back(v[0].size());
+      if (!Dims[Depth])
+        return false;
+    }
+    else
+    {
+      assert( Dims.size() >= Depth + 1 && "Bug: Delete this message after testing" );
+    }
+    for (std::size_t i = 0; i < v.size(); ++i)
+    {
+      if (v[i].size() != Dims[Depth] || !isRegularShapeHelper(v[i], Dims, Depth + 1, bFirst && i==0))
+      {
+        return false;
+      }
+    }
+    return true;
+  }
+
+  template <typename T>
+  bool isRegularShape(const T &t) { return true; }
+
+  template <typename T>
+  bool isRegularShape(const std::vector<T> &v) { return !v.empty(); }
+
+  // Return non-zero if all dimensions of this std::vector<std::vector<T>> are regularly shaped
+  template <typename T>
+  bool isRegularShape(const std::vector<std::vector<T>> &v)
+  {
+    if (v.empty() || v[0].empty())
+      return false;
+    // Make sure all of my rows are the same size
+    std::vector<std::size_t> Dims;
+    Dims.reserve(is_flattenable<T>::vecRank);
+    Dims.push_back(v.size());
+    Dims.push_back(v[0].size());
+    for (std::size_t i = 0; i < Dims[0]; ++i)
+    {
+      if (v[i].size() != Dims[1] || !isRegularShapeHelper(v[i], Dims, 2, i==0))
+      {
+        return false;
+      }
+    }
+    return true;
+  }
 }
 
 // helper function to read space-separated values

Grid/tensors/Tensor_class.h

@@ -417,7 +417,7 @@ public:
       stream << "{";
       for (int j = 0; j < N; j++) {
 	stream << o._internal[i][j];
-	if (i < N - 1) stream << ",";
+	if (j < N - 1) stream << ",";
       }
       stream << "}";
       if (i != N - 1) stream << "\n\t\t";

Grid/threads/Accelerator.cc

@@ -84,11 +84,11 @@ void acceleratorInit(void)
     printf("AcceleratorCudaInit: using default device \n");
     printf("AcceleratorCudaInit: assume user either uses a) IBM jsrun, or \n");
     printf("AcceleratorCudaInit: b) invokes through a wrapping script to set CUDA_VISIBLE_DEVICES, UCX_NET_DEVICES, and numa binding \n");
-    printf("AcceleratorCudaInit: Configure options --enable-summit, --enable-select-gpu=no \n");
+    printf("AcceleratorCudaInit: Configure options --enable-setdevice=no \n");
   }
 #else
   printf("AcceleratorCudaInit: rank %d setting device to node rank %d\n",world_rank,rank);
-  printf("AcceleratorCudaInit: Configure options --enable-select-gpu=yes \n");
+  printf("AcceleratorCudaInit: Configure options --enable-setdevice=yes \n");
   cudaSetDevice(rank);
 #endif
   if ( world_rank == 0 )  printf("AcceleratorCudaInit: ================================================\n");

configure.ac

@@ -390,6 +390,7 @@ case ${CXXTEST} in
     CXXFLAGS="$CXXFLAGS -Xcompiler -fno-strict-aliasing --expt-extended-lambda --expt-relaxed-constexpr"
     if test $ac_openmp = yes; then
        CXXFLAGS="$CXXFLAGS -Xcompiler -fopenmp"
+       LDFLAGS="$LDFLAGS -Xcompiler -fopenmp"
     fi
     ;;
   hipcc)

documentation/manual.rst

@@ -1787,7 +1787,7 @@ Hdf5Writer     Hdf5Reader       HDF5
 
 Write interfaces, similar to the XML facilities in QDP++ are presented. However,
 the serialisation routines are automatically generated by the macro, and a virtual
-reader adn writer interface enables writing to any of a number of formats.
+reader and writer interface enables writing to any of a number of formats.
 
 **Example**::
 
@@ -1814,6 +1814,91 @@ reader adn writer interface enables writing to any of a number of formats.
    }
 
 
+Eigen tensor support -- added 2019H1
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The Serialisation library was expanded in 2019 to support de/serialisation of
+Eigen tensors. De/serialisation of existing types was not changed. Data files
+without Eigen tensors remain compatible with earlier versions of Grid and other readers.
+Conversely, data files containing serialised Eigen tensors is a breaking change.
+
+Eigen tensor serialisation support was added to BaseIO, which was modified to provide a Traits class
+to recognise Eigen tensors with elements that are either: primitive scalars (arithmetic and complex types);
+or Grid tensors.
+
+**Traits determining de/serialisable scalars**::
+
+  // Is this an Eigen tensor
+  template<typename T> struct is_tensor : std::integral_constant<bool,
+    std::is_base_of<Eigen::TensorBase<T, Eigen::ReadOnlyAccessors>, T>::value> {};
+  // Is this an Eigen tensor of a supported scalar
+  template<typename T, typename V = void> struct is_tensor_of_scalar : public std::false_type {};
+  template<typename T> struct is_tensor_of_scalar<T, typename std::enable_if<is_tensor<T>::value && is_scalar<typename T::Scalar>::value>::type> : public std::true_type {};
+  // Is this an Eigen tensor of a supported container
+  template<typename T, typename V = void> struct is_tensor_of_container : public std::false_type {};
+  template<typename T> struct is_tensor_of_container<T, typename std::enable_if<is_tensor<T>::value && isGridTensor<typename T::Scalar>::value>::type> : public std::true_type {};
+
+
+Eigen tensors are regular, multidimensional objects, and each Reader/Writer
+was extended to support this new datatype. Where the Eigen tensor contains
+a Grid tensor, the dimensions of the data written are the dimensions of the
+Eigen tensor plus the dimensions of the underlying Grid scalar. Dimensions
+of size 1 are preserved.
+
+**New Reader/Writer methods for multi-dimensional data**::
+
+  template <typename U>
+  void readMultiDim(const std::string &s, std::vector<U> &buf, std::vector<size_t> &dim);
+  template <typename U>
+  void writeMultiDim(const std::string &s, const std::vector<size_t> & Dimensions, const U * pDataRowMajor, size_t NumElements);
+
+
+On readback, the Eigen tensor rank must match the data being read, but the tensor
+dimensions will be resized if necessary. Resizing is not possible for Eigen::TensorMap<T>
+because these tensors use a buffer provided at construction, and this buffer cannot be changed.
+Deserialisation failures cause Grid to assert.
+
+
+HDF5 Optimisations -- added June 2021
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Grid serialisation is intended to be light, deterministic and provide a layer of abstraction over
+multiple file formats. HDF5 excels at handling multi-dimensional data, and the Grid HDF5Reader/HDF5Writer exploits this.
+When serialising nested ``std::vector<T>``, where ``T`` is an arithmetic or complex type,
+the Hdf5Writer writes the data as an Hdf5 DataSet object.
+
+However, nested ``std::vector<std::vector<...T>>`` might be "ragged", i.e. not necessarily regular. E.g. a 3d nested
+``std::vector`` might contain 2 rows, the first being a 2x2 block and the second row being a 1 x 2 block.
+A bug existed whereby this was not checked on write, so nested, ragged vectors
+were written as a regular dataset, with a buffer under/overrun and jumbled contents.
+
+Clearly this was not used in production, as the bug went undetected until now. Fixing this bug
+is an opportunity to further optimise the HDF5 file format.
+
+The goals of this change are to:
+
+* Make changes to the Hdf5 file format only -- i.e. do not impact other file formats
+
+* Implement file format changes in such a way that they are transparent to the Grid reader
+
+* Correct the bug for ragged vectors of numeric / complex types
+
+* Extend the support of nested std::vector<T> to arbitrarily nested Grid tensors
+
+
+The trait class ``element`` has been redefined to ``is_flattenable``, which is a trait class for
+potentially "flattenable" objects. These are (possibly nested) ``std::vector<T>`` where ``T`` is
+an arithmetic, complex or Grid tensor type. Flattenable objects are tested on write
+(with the function ``isRegularShape``) to see whether they actually are regular.
+
+Flattenable, regular objects are written to a multidimensional HDF5 DataSet.
+Otherwise, an Hdf5 sub group is created with the object "name", and each element of the outer dimension is
+recursively written to as object "name_n", where n is a 0-indexed number.
+
+On readback (by Grid)), the presence of a subgroup containing the attribute ``Grid_vector_size`` triggers a
+"ragged read", otherwise a read from a DataSet is attempted.
+
+
 Data parallel field IO
 -----------------------
 

tests/IO/Test_serialisation.cc

@@ -48,7 +48,9 @@ public:
                           std::vector<double>, array,
                           std::vector<std::vector<double> >, twodimarray,
 			  std::vector<std::vector<std::vector<std::complex<double>> > >, cmplx3darray,
-			  SpinColourMatrix, scm
+			  SpinColourMatrix, scm,
+                          std::vector<std::vector<std::vector<int> > >, ragged,
+                          std::vector<std::vector<SpinColourMatrix> >, vscm
                           );
   myclass() {}
   myclass(int i)
@@ -56,6 +58,10 @@ public:
   , twodimarray(3,std::vector<double>(5, 1.23456))
   , cmplx3darray(3,std::vector<std::vector<std::complex<double>>>(5, std::vector<std::complex<double>>(7, std::complex<double>(1.2, 3.4))))
   , ve(2, myenum::blue)
+  , ragged( {{{i+1},{i+2,i+3}}, // ragged
+            {{i+4,i+5,i+6,i+7},{i+8,i+9,i+10,i+11},{i+12,i+13,i+14,i+15}}, // block
+            {{i+16,i+17},{i+18,i+19,i+20}}} ) //ragged
+  , vscm(3, std::vector<SpinColourMatrix>(5))
   {
     e=myenum::red;
     x=i;
@@ -68,6 +74,13 @@ public:
     scm()(0, 2)(1, 1) = 6.336;
     scm()(2, 1)(2, 2) = 7.344;
     scm()(1, 1)(2, 0) = 8.3534;
+    int Counter = i;
+    for( auto & v : vscm ) {
+      for( auto & j : v ) {
+        j = std::complex<double>(Counter, -Counter);
+        Counter++;
+      }
+    }
   }
 };
 

tests/Test_meson_field.cc

@@ -0,0 +1,148 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: tests/core/Test_meson_field.cc
+
+Copyright (C) 2015-2018
+
+Author: Felix Erben <felix.erben@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
+*************************************************************************************/
+
+#include <Grid/Grid.h>
+#include <Grid/qcd/utils/A2Autils.h>
+
+using namespace Grid;
+
+const int TSRC = 0;  //timeslice where rho is nonzero
+const int VDIM = 5; //length of each vector
+
+typedef typename DomainWallFermionR::ComplexField ComplexField;
+typedef typename DomainWallFermionR::FermionField FermionField;
+
+int main(int argc, char *argv[])
+{
+  // initialization
+  Grid_init(&argc, &argv);
+  std::cout << GridLogMessage << "Grid initialized" << std::endl;
+
+  // Lattice and rng setup 
+  Coordinate latt_size   = GridDefaultLatt();
+  Coordinate simd_layout = GridDefaultSimd(4, vComplex::Nsimd());
+  Coordinate mpi_layout  = GridDefaultMpi();
+  GridCartesian    grid(latt_size,simd_layout,mpi_layout);
+  int Nt = GridDefaultLatt()[Tp];
+  Lattice<iScalar<vInteger>> t(&grid);
+  LatticeCoordinate(t, Tp);
+  std::vector<int> seeds({1,2,3,4});
+  GridParallelRNG          pRNG(&grid);
+  pRNG.SeedFixedIntegers(seeds);
+
+  // MesonField lhs and rhs vectors
+  std::vector<FermionField> phi(VDIM,&grid);
+  std::vector<FermionField> rho(VDIM,&grid);
+  FermionField rho_tmp(&grid);
+  std::cout << GridLogMessage << "Initialising random meson fields" << std::endl;
+  for (unsigned int i = 0; i < VDIM; ++i){
+    random(pRNG,phi[i]);
+    random(pRNG,rho_tmp); //ideally only nonzero on t=0
+    rho[i] = where((t==TSRC), rho_tmp, 0.*rho_tmp); //ideally only nonzero on t=0
+  }
+  std::cout << GridLogMessage << "Meson fields initialised, rho non-zero only for t = " << TSRC << std::endl;
+
+  // Gamma matrices used in the contraction
+  std::vector<Gamma::Algebra> Gmu = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
+  };
+
+  // momentum phases e^{ipx}
+  std::vector<std::vector<double>> momenta = {
+	  {0.,0.,0.},
+	  {1.,0.,0.},
+	  {1.,1.,0.},
+	  {1.,1.,1.},
+	  {2.,0.,0.}
+  };
+
+  std::cout << GridLogMessage << "Meson fields will be created for " << Gmu.size() << " Gamma matrices and " << momenta.size() << " momenta." << std::endl;
+
+  std::cout << GridLogMessage << "Computing complex phases" << std::endl;
+  std::vector<ComplexField> phases(momenta.size(),&grid);
+  ComplexField coor(&grid);
+  Complex Ci(0.0,1.0);
+  for (unsigned int j = 0; j < momenta.size(); ++j)
+  {
+    phases[j] = Zero();
+    for(unsigned int mu = 0; mu < momenta[j].size(); mu++)
+    {
+      LatticeCoordinate(coor, mu);
+      phases[j] = phases[j] + momenta[j][mu]/GridDefaultLatt()[mu]*coor;
+    }
+    phases[j] = exp((Real)(2*M_PI)*Ci*phases[j]);
+  }
+  std::cout << GridLogMessage << "Computing complex phases done." << std::endl;
+
+  Eigen::Tensor<ComplexD,5, Eigen::RowMajor> Mpp(momenta.size(),Gmu.size(),Nt,VDIM,VDIM);
+  Eigen::Tensor<ComplexD,5, Eigen::RowMajor> Mpr(momenta.size(),Gmu.size(),Nt,VDIM,VDIM);
+  Eigen::Tensor<ComplexD,5, Eigen::RowMajor> Mrr(momenta.size(),Gmu.size(),Nt,VDIM,VDIM);
+
+  // timer
+  double start,stop;
+
+  //execute meson field routine
+  start = usecond();
+  A2Autils<WilsonImplR>::MesonField(Mpp,&phi[0],&phi[0],Gmu,phases,Tp);
+  stop = usecond();
+  std::cout << GridLogMessage << "M(phi,phi) created, execution time " << stop-start << " us" << std::endl;
+  start = usecond();
+  /* Ideally, for this meson field we could pass TSRC (even better a list of timeslices)
+   * to the routine so that all the compnents which are predictably equal to zero are not computed. */
+  A2Autils<WilsonImplR>::MesonField(Mpr,&phi[0],&rho[0],Gmu,phases,Tp);
+  stop = usecond();
+  std::cout << GridLogMessage << "M(phi,rho) created, execution time " << stop-start << " us" << std::endl;
+  start = usecond();
+  A2Autils<WilsonImplR>::MesonField(Mrr,&rho[0],&rho[0],Gmu,phases,Tp);
+  stop = usecond();
+  std::cout << GridLogMessage << "M(rho,rho) created, execution time " << stop-start << " us" << std::endl;
+
+  std::string FileName = "Meson_Fields";
+#ifdef HAVE_HDF5
+  using Default_Reader = Grid::Hdf5Reader;
+  using Default_Writer = Grid::Hdf5Writer;
+  FileName.append(".h5");
+#else
+  using Default_Reader = Grid::BinaryReader;
+  using Default_Writer = Grid::BinaryWriter;
+  FileName.append(".bin");
+#endif
+
+  Default_Writer w(FileName);
+  write(w,"phi_phi",Mpp);
+  write(w,"phi_rho",Mpr);
+  write(w,"rho_rho",Mrr);
+
+  // epilogue
+  std::cout << GridLogMessage << "Grid is finalizing now" << std::endl;
+  Grid_finalize();
+  
+  return EXIT_SUCCESS;
+}