Hdf5 writing of scalar (i.e. no Grid subtypes) Eigen::Tensor works. But issues when adding Eigen::Tensor to serialisable object.

Grid/serialisation/BaseIO.h

@@ -55,6 +55,12 @@ namespace Grid {
     template <typename U>
     typename std::enable_if<!std::is_base_of<Serializable, U>::value, void>::type
     write(const std::string& s, const U &output);
+    template <typename U>
+    void write(const std::string &s, const iScalar<U> &output);
+    template <typename U, int N>
+    void write(const std::string &s, const iVector<U, N> &output);
+    template <typename U, int N>
+    void write(const std::string &s, const iMatrix<U, N> &output);
     template <typename Scalar_, int NumIndices_, int Options_, typename IndexType_>
     typename std::enable_if<std::is_arithmetic<Scalar_>::value || Grid::is_complex<Scalar_>::value, void>::type
     write(const std::string &s, const Eigen::Tensor<Scalar_, NumIndices_, Options_, IndexType_> &output);
@@ -64,13 +70,7 @@ namespace Grid {
     void write(const std::string &s, const Eigen::Tensor<iVector<U, N>, NumIndices_, Options_, IndexType_> &output);
     template<typename U, int N, int NumIndices_, int Options_, typename IndexType_>
     void write(const std::string &s, const Eigen::Tensor<iMatrix<U, N>, NumIndices_, Options_, IndexType_> &output);
-    template <typename U>
+
-    void write(const std::string &s, const iScalar<U> &output);
-    template <typename U, int N>
-    void write(const std::string &s, const iVector<U, N> &output);
-    template <typename U, int N>
-    void write(const std::string &s, const iMatrix<U, N> &output);
-    
     void         scientificFormat(const bool set);
     bool         isScientific(void);
     void         setPrecision(const unsigned int prec);
@@ -152,43 +152,6 @@ namespace Grid {
     upcast->writeDefault(s, output);
   }
 
-  // Eigen::Tensors of arithmetic/complex base type
-  template <typename T>
-  template <typename Scalar_, int NumIndices_, int Options_, typename IndexType_>
-  typename std::enable_if<std::is_arithmetic<Scalar_>::value || Grid::is_complex<Scalar_>::value, void>::type
-  Writer<T>::write(const std::string &s, const Eigen::Tensor<Scalar_, NumIndices_, Options_, IndexType_> &output)
-  {
-    //upcast->writeDefault(s, tensorToVec(output));
-    std::cout << "I really should add code to write Eigen::Tensor (arithmetic/complex) ..." << std::endl;
-  }
-
-  // Eigen::Tensors of iScalar<U>
-  template <typename T>
-  template<typename U, int NumIndices_, int Options_, typename IndexType_>
-  void Writer<T>::write(const std::string &s, const Eigen::Tensor<iScalar<U>, NumIndices_, Options_, IndexType_> &output)
-  {
-    //upcast->writeDefault(s, tensorToVec(output));
-    std::cout << "I really should add code to write Eigen::Tensor (iScalar) ..." << std::endl;
-  }
-
-  // Eigen::Tensors of iVector<U, N>
-  template <typename T>
-  template<typename U, int N, int NumIndices_, int Options_, typename IndexType_>
-  void Writer<T>::write(const std::string &s, const Eigen::Tensor<iVector<U, N>, NumIndices_, Options_, IndexType_> &output)
-  {
-    //upcast->writeDefault(s, tensorToVec(output));
-    std::cout << "I really should add code to write Eigen::Tensor (iVector) ..." << std::endl;
-  }
-
-  // Eigen::Tensors of iMatrix<U, N>
-  template <typename T>
-  template<typename U, int N, int NumIndices_, int Options_, typename IndexType_>
-  void Writer<T>::write(const std::string &s, const Eigen::Tensor<iMatrix<U, N>, NumIndices_, Options_, IndexType_> &output)
-  {
-    //upcast->writeDefault(s, tensorToVec(output));
-    std::cout << "I really should add code to write Eigen::Tensor (iMatrix) ..." << std::endl;
-  }
-
 
   template <typename T>
   template <typename U>
@@ -211,6 +174,83 @@ namespace Grid {
     upcast->writeDefault(s, tensorToVec(output));
   }
 
+  // Eigen::Tensors of arithmetic/complex base type
+  template <typename T>
+  template <typename Scalar_, int NumIndices_, int Options_, typename IndexType_>
+  typename std::enable_if<std::is_arithmetic<Scalar_>::value || Grid::is_complex<Scalar_>::value, void>::type
+  Writer<T>::write(const std::string &s, const Eigen::Tensor<Scalar_, NumIndices_, Options_, IndexType_> &output)
+  {
+    typedef Eigen::Tensor<Scalar_, NumIndices_, Options_, IndexType_> Tensor;
+    const typename Tensor::Index NumElements{output.size()};
+    assert( NumElements > 0 );
+    if( NumElements == 1 )
+      upcast->writeDefault(s, * output.data());
+    else {
+      // Create a single, flat vector to hold all the data
+      std::vector<Scalar_> flat(NumElements);
+      // We're not interested in trivial dimensions, i.e. dimensions = 1
+      const typename Tensor::Dimensions & DimsOriginal{output.dimensions()};
+      assert(DimsOriginal.size() == NumIndices_);
+      unsigned int TrivialDimCount{0};
+      for(auto i : DimsOriginal ) {
+        if( i <= 1 ) {
+          TrivialDimCount++;
+          assert( i == 1 ); // Not expecting dimension to be <= 0
+        }
+      }
+      const unsigned int ReducedDimCount{NumIndices_ - TrivialDimCount};
+      assert( ReducedDimCount > 0 ); // NB: We've already checked this is not a scalar
+      // Save a flat vector of the non-trivial dimensions
+      std::vector<size_t> ReducedDims(ReducedDimCount);
+      unsigned int ui = 0;
+      for(auto i : DimsOriginal ) {
+        if( i > 1 ) {
+          ReducedDims[ui] = static_cast<size_t>(i);
+          assert( ReducedDims[ui] == i ); // check we didn't lose anything in the conversion
+          ui++;
+        }
+      }
+      // Now copy all the data to my flat vector
+      // Regardless of the Eigen::Tensor storage order, the copy will be Row Major
+      std::array<typename Tensor::Index, NumIndices_> MyIndex;
+      for( int i = 0 ; i < NumIndices_ ; i++ ) MyIndex[i] = 0;
+      for( typename Tensor::Index n = 0; n < NumElements; n++ ) {
+        flat[n] = output( MyIndex );
+        // Now increment the index
+        for( int i = NumIndices_ - 1; i >= 0 && ++MyIndex[i] == DimsOriginal[i]; i-- )
+          MyIndex[i] = 0;
+      }
+      upcast->template writeMultiDim<Scalar_>(s, ReducedDims, flat);
+    }
+  }
+  
+  // Eigen::Tensors of iScalar<U>
+  template <typename T>
+  template<typename U, int NumIndices_, int Options_, typename IndexType_>
+  void Writer<T>::write(const std::string &s, const Eigen::Tensor<iScalar<U>, NumIndices_, Options_, IndexType_> &output)
+  {
+    //upcast->writeDefault(s, tensorToVec(output));
+    std::cout << "I really should add code to write Eigen::Tensor (iScalar) ..." << std::endl;
+  }
+  
+  // Eigen::Tensors of iVector<U, N>
+  template <typename T>
+  template<typename U, int N, int NumIndices_, int Options_, typename IndexType_>
+  void Writer<T>::write(const std::string &s, const Eigen::Tensor<iVector<U, N>, NumIndices_, Options_, IndexType_> &output)
+  {
+    //upcast->writeDefault(s, tensorToVec(output));
+    std::cout << "I really should add code to write Eigen::Tensor (iVector) ..." << std::endl;
+  }
+  
+  // Eigen::Tensors of iMatrix<U, N>
+  template <typename T>
+  template<typename U, int N, int NumIndices_, int Options_, typename IndexType_>
+  void Writer<T>::write(const std::string &s, const Eigen::Tensor<iMatrix<U, N>, NumIndices_, Options_, IndexType_> &output)
+  {
+    //upcast->writeDefault(s, tensorToVec(output));
+    std::cout << "I really should add code to write Eigen::Tensor (iMatrix) ..." << std::endl;
+  }
+
   template <typename T>
   void Writer<T>::scientificFormat(const bool set)
   {

Grid/serialisation/Hdf5IO.h

@@ -38,6 +38,8 @@ namespace Grid
     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>
+    void writeMultiDim(const std::string &s, const std::vector<size_t> & Dimensions, const std::vector<U> & DataRowMajor);
     H5NS::Group & getGroup(void);
   private:
     template <typename U>
@@ -101,6 +103,35 @@ namespace Grid
   template <>
   void Hdf5Writer::writeDefault(const std::string &s, const std::string &x);
   
+  template <typename U>
+  void Hdf5Writer::writeMultiDim(const std::string &s, const std::vector<size_t> & Dimensions, const std::vector<U> & DataRowMajor)
+  {
+    // Hdf5 needs the dimensions as hsize_t
+    std::vector<hsize_t> dim;
+    for (auto &d: Dimensions)
+      dim.push_back(d);
+    // write to file
+    H5NS::DataSpace dataSpace(dim.size(), dim.data());
+    
+    if (DataRowMajor.size() > dataSetThres_)
+    {
+      H5NS::DataSet           dataSet;
+      H5NS::DSetCreatPropList plist;
+      
+      plist.setChunk(dim.size(), dim.data());
+      plist.setFletcher32();
+      dataSet = group_.createDataSet(s, Hdf5Type<U>::type(), dataSpace, plist);
+      dataSet.write(&DataRowMajor[0], Hdf5Type<U>::type());
+    }
+    else
+    {
+      H5NS::Attribute attribute;
+      
+      attribute = group_.createAttribute(s, Hdf5Type<U>::type(), dataSpace);
+      attribute.write(Hdf5Type<U>::type(), &DataRowMajor[0]);
+    }
+  }
+
   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)
@@ -110,34 +141,11 @@ namespace Grid
     
     // flatten the vector and getting dimensions
     Flatten<std::vector<U>> flat(x);
-    std::vector<hsize_t> dim;
+    std::vector<size_t> dim;
     const auto           &flatx = flat.getFlatVector();
-    
     for (auto &d: flat.getDim())
-    {
       dim.push_back(d);
-    }
+    writeMultiDim<Element>(s, dim, flatx);
-    
-    // write to file
-    H5NS::DataSpace dataSpace(dim.size(), dim.data());
-    
-    if (flatx.size() > dataSetThres_)
-    {
-      H5NS::DataSet           dataSet;
-      H5NS::DSetCreatPropList plist;
-      
-      plist.setChunk(dim.size(), dim.data());
-      plist.setFletcher32();
-      dataSet = group_.createDataSet(s, Hdf5Type<Element>::type(), dataSpace, plist);
-      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>

Grid/serialisation/VectorUtils.h

@@ -53,6 +53,17 @@ namespace Grid {
     return os;
   }
   
+  // std::vector<std:vector<...>> nested to specified Rank //////////////////////////////////
+  template<typename T, unsigned int Rank>
+  struct NestedStdVector {
+    typedef typename std::vector<typename NestedStdVector<T, Rank - 1>::type> type;
+  };
+  
+  template<typename T>
+  struct NestedStdVector<T,0> {
+    typedef T type;
+  };
+  
   // Grid scalar tensors to nested std::vectors //////////////////////////////////
   template <typename T>
   struct TensorToVec

tests/hadrons/Test_hadrons_distil.cc

@@ -441,18 +441,21 @@ bool DebugEigenTest()
 }
 
 typedef iMatrix<Complex,7> OddBall;
-typedef Eigen::Tensor<int, 3, Eigen::RowMajor> TensorInt;
-typedef Eigen::Tensor<std::complex<double>, 3, Eigen::RowMajor> TensorComplex;
 typedef Eigen::Tensor<OddBall, 3, Eigen::RowMajor> TensorOddBall;
 
+typedef int TestScalar;
+//typedef std::complex<double> TestScalar;
+typedef Eigen::Tensor<TestScalar, 3, Eigen::RowMajor> TestTensor;
+
 // From Test_serialisation.cc
 class myclass: Serializable {
 public:
   GRID_SERIALIZABLE_CLASS_MEMBERS(myclass
-                                  //, OddBall, critter
+                                  //, TestTensor, critter
                                   , SpinColourVector, scv
                                   , SpinColourMatrix, scm
                                   );
+  //myclass() : critter(7,3,2) {}
 };
 
 template <typename W, typename R, typename O>
@@ -483,11 +486,31 @@ bool DebugIOTest(void) {
   ioTest<Hdf5Writer, Hdf5Reader, SpinColourMatrix>("iotest_matrix.h5", scm, "SpinColourMatrix");
   SpinColourVector scv;
   ioTest<Hdf5Writer, Hdf5Reader, SpinColourVector>("iotest_vector.h5", scv, "SpinColourVector");
+
+  TestTensor t(3,6,2);
+  TestScalar Val{1};
+  const TestScalar Inc{1};
+  for( int i = 0 ; i < 3 ; i++)
+    for( int j = 0 ; j < 6 ; j++)
+      for( int k = 0 ; k < 2 ; k++) {
+        t(i,j,k) = Val;
+        Val += Inc;
+      }
+  ioTest<Hdf5Writer, Hdf5Reader, TestTensor>("iotest_tensor.h5", t, "eigen_tensor_instance_name");
+
+  TestTensor t2(t);
+  auto bEqual = (t == t2).all();
+  std::cout << "(t2 == t) = " << bEqual << std::endl;
+  //if( * bEqual.data() != 0 )
+    //std::cout << "t2 == t" << std::endl;
+  //else
+    //std::cout << "t2 != t" << std::endl;
+
   myclass o;
   ioTest<Hdf5Writer, Hdf5Reader, myclass>("iotest_object.h5", o, "myclass_object_instance_name");
-  TensorInt t(3,6,2);
-  ioTest<Hdf5Writer, Hdf5Reader, TensorInt>("iotest_tensor.h5", t, "eigen_tensor_instance_name");
 
+  std::cout << "Wow!" << std::endl;
+  
   return true;
 }
 #endif