Tensor serialisation is fully functional

Grid/serialisation/BaseIO.h

@@ -316,7 +316,7 @@ namespace Grid {
       pWriteBuffer = getFirstScalar( output );
     } else {
       // Regardless of the Eigen::Tensor storage order, the copy will be Row Major
-      pCopyBuffer = static_cast<Scalar *>(malloc(TotalNumElements * sizeof(Scalar)));
+      pCopyBuffer = new Scalar[TotalNumElements];
       pWriteBuffer = pCopyBuffer;
       Scalar * pCopy = pCopyBuffer;
       std::array<Index, TensorRank> MyIndex;
@@ -330,7 +330,7 @@ namespace Grid {
       }
     }
     upcast->template writeMultiDim<Scalar>(s, TotalDims, pWriteBuffer, TotalNumElements);
-    if( pCopyBuffer ) free( pCopyBuffer );
+    if( pCopyBuffer ) delete [] pCopyBuffer;
   }
 
   template <typename T>

Grid/serialisation/XmlIO.h

@@ -138,15 +138,30 @@ namespace Grid
   {
     push(s);
     size_t count = 1;
-    const size_t Rank = Dimensions.size();
+    const int Rank = static_cast<int>( Dimensions.size() );
     write("rank", Rank );
+    std::vector<size_t> MyIndex( Rank );
     for( auto d : Dimensions ) {
       write("dim", d);
       count *= d;
     }
     assert( count == NumElements && "XmlIO : element count doesn't match dimensions" );
-    while (NumElements--)
+    static const char sName[] = "tensor";
+    for( int i = 0 ; i < Rank ; i++ ) {
+      MyIndex[i] = 0;
+      push(sName);
+    }
+    while (NumElements--) {
       write("elem", *pDataRowMajor++);
+      int i;
+      for( i = Rank - 1 ; i != -1 && ++MyIndex[i] == Dimensions[i] ; i-- )
+        MyIndex[i] = 0;
+      int Rollover = Rank - 1 - i;
+      for( i = 0 ; i < Rollover ; i++ )
+        pop();
+      for( i = 0 ; NumElements && i < Rollover ; i++ )
+        push(sName);
+    }
     pop();
   }
 
@@ -189,7 +204,9 @@ namespace Grid
       std::cout << GridLogWarning << "XML: cannot open node '" << s << "'";
       std::cout << std::endl;
     } else {
-      size_t Rank;
+      static const char sName[] = "tensor";
+      static const char sNameDone[] = "tensor-done";
+      int Rank;
       read("rank", Rank);
       dim.resize( Rank );
       size_t NumElements = 1;
@@ -200,10 +217,27 @@ namespace Grid
         NumElements *= d;
       }
       buf.resize( NumElements );
+      std::vector<size_t> MyIndex( Rank );
+      for( int i = 0 ; i < Rank ; i++ ) {
+        MyIndex[i] = 0;
+        push(sName);
+      }
+
       for( auto &x : buf )
       {
+        NumElements--;
         read("elem", x);
         node_.child("elem").set_name("elem-done");
+        int i;
+        for( i = Rank - 1 ; i != -1 && ++MyIndex[i] == dim[i] ; i-- )
+          MyIndex[i] = 0;
+        int Rollover = Rank - 1 - i;
+        for( i = 0 ; i < Rollover ; i++ ) {
+          node_.set_name(sNameDone);
+          pop();
+        }
+        for( i = 0 ; NumElements && i < Rollover ; i++ )
+          push(sName);
       }
       pop();
     }

Grid/util/EigenUtil.h

@@ -2,7 +2,7 @@
  
  Grid physics library, www.github.com/paboyle/Grid
  
- Source file: Grid/util/Eigen.h
+ Source file: Grid/util/EigenUtil.h
  
  Copyright (C) 2019
  
@@ -25,8 +25,8 @@
  See the full license in the file "LICENSE" in the top level distribution directory
  *************************************************************************************/
 /*  END LEGAL */
-#ifndef GRID_UTIL_EIGEN_H
+#ifndef GRID_UTIL_EIGENUTIL_H
-#define GRID_UTIL_EIGEN_H
+#define GRID_UTIL_EIGENUTIL_H
 #include <Grid/tensors/Tensor_traits.h>
 #include <Grid/Eigen/unsupported/CXX11/Tensor>
 
@@ -35,25 +35,28 @@ namespace Grid {
   template <typename ETensor, typename Lambda>
   typename std::enable_if<EigenIO::is_tensor_of_scalar<ETensor>::value, void>::type
   for_all_do_lambda( Lambda lambda, typename ETensor::Scalar &scalar, typename ETensor::Index &Seq,
-                    std::array<std::size_t, ETensor::NumIndices + EigenIO::Traits<ETensor>::Rank> &MyIndex)
+      const std::array<typename ETensor::Index, ETensor::NumIndices> &MyIndex,
+      const std::array<int, EigenIO::Traits<ETensor>::Rank> &DimGridTensor,
+            std::array<int, EigenIO::Traits<ETensor>::Rank> &GridTensorIndex)
   {
-    lambda( scalar, Seq++, MyIndex );
+    lambda( scalar, Seq++, MyIndex, GridTensorIndex );
   }
   
   // for_all helper function to call the lambda for container
   template <typename ETensor, typename Lambda>
   typename std::enable_if<EigenIO::is_tensor_of_container<ETensor>::value, void>::type
   for_all_do_lambda( Lambda lambda, typename ETensor::Scalar &container, typename ETensor::Index &Seq,
-                    std::array<std::size_t, ETensor::NumIndices + EigenIO::Traits<ETensor>::Rank> &MyIndex)
+      const std::array<typename ETensor::Index, ETensor::NumIndices> &MyIndex,
+      const std::array<int, EigenIO::Traits<ETensor>::Rank> &DimGridTensor,
+            std::array<int, EigenIO::Traits<ETensor>::Rank> &GridTensorIndex)
   {
     using Traits = EigenIO::Traits<ETensor>;
-    const auto rank{ETensor::NumIndices};
     const auto InnerRank = Traits::Rank;
     for( typename Traits::scalar_type &Source : container ) {
-      lambda(Source, Seq++, MyIndex );
+      lambda(Source, Seq++, MyIndex, GridTensorIndex );
       // Now increment SubIndex
-      for( auto i = InnerRank - 1; i != -1 && ++MyIndex[rank + i] == Traits::Dimension(i); i-- )
+      for( auto i = InnerRank - 1; i != -1 && ++GridTensorIndex[i] == DimGridTensor[i]; i-- )
-        MyIndex[rank + i] = 0;
+        GridTensorIndex[i] = 0;
     }
   }
   
@@ -65,58 +68,51 @@ namespace Grid {
   for_all( ETensor &ET, Lambda lambda )
   {
     using Scalar = typename ETensor::Scalar; // This could be a Container - we'll check later
-    using Traits = EigenIO::Traits<ETensor>;
-    const std::size_t NumScalars = ET.size();
-    assert( NumScalars > 0 );
     using Index = typename ETensor::Index;
+    using Traits = EigenIO::Traits<ETensor>;
+    // Check that the number of elements in the container is the product of tensor dimensions
+    const Index NumScalars = ET.size();
+    assert( NumScalars > 0 && "EigenUtil: tensor has no elements" );
     Index ScalarElementCount{1};
-    const auto InnerRank = Traits::Rank;
     const auto rank{ETensor::NumIndices};
-    std::array<std::size_t, rank + InnerRank> Dims;
+    std::array<Index, rank> DimTensor, MyIndex;
-    for(auto i = 0; i < rank; i++ ) {
+    for(int i = 0; i < rank; i++ ) {
-      auto dim = ET.dimension(i);
+      DimTensor[i] = ET.dimension(i);
-      assert( dim > 0 );
+      ScalarElementCount *= DimTensor[i];
-      Dims[i] = static_cast<std::size_t>(dim);
+      MyIndex[i] = 0;
-      assert( Dims[i] == dim ); // check we didn't lose anything in the conversion
-      ScalarElementCount *= Dims[i];
     }
-    // Check that the number of containers is correct ... and we didn't lose anything in conversions
+    assert( NumScalars == ScalarElementCount && "EigenUtil: tensor size not product of dimensions" );
-    assert( NumScalars == ScalarElementCount );
+    // Save the GridTensor dimensions
-    // If the Scalar is actually a container, add the inner Scalar's dimensions
+    const auto InnerRank{Traits::Rank};
-    size_t InnerScalarCount{1};
+    std::array<int, InnerRank> DimGridTensor, GridTensorIndex;
-    for(auto i = 0; i < InnerRank; i++ ) {
+    for(int i = 0; i < InnerRank; i++ ) {
-      auto dim = Traits::Dimension(i);
+      DimGridTensor[i] = Traits::Dimension(i);
-      assert( dim > 0 );
+      GridTensorIndex[i] = 0;
-      Dims[rank + i] = static_cast<std::size_t>(dim);
-      assert( Dims[rank + i] == dim ); // check we didn't lose anything in the conversion
-      InnerScalarCount *= dim;
     }
-    assert(Traits::count == InnerScalarCount);
+    // Now walk the tensor in memory order
-    std::array<std::size_t, rank + InnerRank> MyIndex;
-    for( auto &idx : MyIndex ) idx = 0;
     Index Seq = 0;
     Scalar * pScalar = ET.data();
-    for( std::size_t j = 0; j < NumScalars; j++ ) {
+    for( Index j = 0; j < NumScalars; j++ ) {
-      for_all_do_lambda<ETensor, Lambda>( lambda, * pScalar, Seq, MyIndex );
+      for_all_do_lambda<ETensor, Lambda>( lambda, * pScalar, Seq, MyIndex, DimGridTensor, GridTensorIndex );
       // Now increment the index to pass to the lambda (bearing in mind we're walking in memory order)
       if( ETensor::Options & Eigen::RowMajor ) {
-        for( auto i = rank - 1; i != -1 && ++MyIndex[i] == Dims[i]; i-- )
+        for( auto i = rank - 1; i != -1 && ++MyIndex[i] == DimTensor[i]; i-- )
           MyIndex[i] = 0;
       } else {
-        for( auto i = 0; i < rank && ++MyIndex[i] == Dims[i]; i++ )
+        for( auto i = 0; i < rank && ++MyIndex[i] == DimTensor[i]; i++ )
           MyIndex[i] = 0;
-        size_t NewSeq = 0;
+        Seq = 0;
-        for( auto i = 0; i < rank + InnerRank ; i++ ) {
+        for( auto i = 0; i < rank; i++ ) {
-          NewSeq *= Dims[i];
+          Seq *= DimTensor[i];
-          NewSeq += MyIndex[i];
+          Seq += MyIndex[i];
         }
-        Seq = static_cast<Index>( NewSeq );
+        Seq *= Traits::count;
       }
       pScalar++;
     }
   }
   
-  // Sequential initialisation of tensors
+  // Sequential initialisation of tensors up to specified precision
   // Would have preferred to define template variables for this, but that's c++ 17
   template <typename ETensor>
   typename std::enable_if<EigenIO::is_tensor<ETensor>::value && !EigenIO::Traits<ETensor>::is_complex>::type
@@ -124,11 +120,13 @@ namespace Grid {
   {
     using Traits = EigenIO::Traits<ETensor>;
     using scalar_type = typename Traits::scalar_type;
-    for_all( ET, [&](scalar_type &c, typename ETensor::Index n, const std::array<size_t, ETensor::NumIndices + Traits::Rank> &Dims ) {
+    using Index = typename ETensor::Index;
+    for_all( ET, [&](scalar_type &c, Index n, const std::array<Index, ETensor::NumIndices> &TensorIndex,
+                     const std::array<int, Traits::Rank> &ScalarIndex ) {
       scalar_type x = Inc * static_cast<scalar_type>(n);
       if( Precision ) {
         std::stringstream s;
-        s << std::scientific << std::setprecision(Precision) << x;
+        s << std::setprecision(Precision) << x;
         s >> x;
       }
       c = x;
@@ -141,7 +139,9 @@ namespace Grid {
   {
     using Traits = EigenIO::Traits<ETensor>;
     using scalar_type = typename Traits::scalar_type;
-    for_all( ET, [&](scalar_type &c, typename ETensor::Index n, const std::array<size_t, ETensor::NumIndices + Traits::Rank> &Dims ) {
+    using Index = typename ETensor::Index;
+    for_all( ET, [&](scalar_type &c, Index n, const std::array<Index, ETensor::NumIndices> &TensorIndex,
+                     const std::array<int, Traits::Rank> &ScalarIndex ) {
       auto re = Inc.real();
       auto im = Inc.imag();
       re *= n;
@@ -159,24 +159,28 @@ namespace Grid {
   }
   
   // Helper to dump a tensor
-#ifdef DEBUG
-#define dump_tensor(args...) dump_tensor_func(args)
   template <typename T>
   typename std::enable_if<EigenIO::is_tensor<T>::value, void>::type
-  dump_tensor_func(T &t, const char * pName = nullptr)
+  dump_tensor(T &t, const char * pName = nullptr)
   {
     using Traits = EigenIO::Traits<T>;
+    using scalar_type = typename Traits::scalar_type;
+    using Index = typename T::Index;
     const auto rank{T::NumIndices};
     const auto &dims = t.dimensions();
-    std::cout << "Dumping rank " << rank << ((T::Options & Eigen::RowMajor) ? ", row" : ", column") << "-major tensor ";
+    std::cout << "Dumping rank " << rank + Traits::Rank << ((T::Options & Eigen::RowMajor) ? ", row" : ", column") << "-major tensor ";
     if( pName )
       std::cout << pName;
-    for( auto i = 0 ; i < rank; i++ ) std::cout << "[" << dims[i] << "]";
+    for( int i = 0 ; i < rank; i++ ) std::cout << "[" << dims[i] << "]";
+    for( int i = 0 ; i < Traits::Rank; i++ ) std::cout << "(" << Traits::Dimension(i) << ")";
     std::cout << " in memory order:" << std::endl;
-    for_all( t, [&](typename Traits::scalar_type &c, typename T::Index index, const std::array<size_t, T::NumIndices + Traits::Rank> &Dims ){
+    for_all( t, [&](scalar_type &c, Index n, const std::array<Index, rank> &TensorIndex,
+                    const std::array<int, Traits::Rank> &ScalarIndex ){
       std::cout << "  ";
-      for( auto dim : Dims )
+      for( auto dim : TensorIndex )
         std::cout << "[" << dim << "]";
+      for( auto dim : ScalarIndex )
+        std::cout << "(" << dim << ")";
       std::cout << " = " << c << std::endl;
     } );
     std::cout << "========================================" << std::endl;
@@ -184,64 +188,11 @@ namespace Grid {
   
   template <typename T>
   typename std::enable_if<!EigenIO::is_tensor<T>::value, void>::type
-  dump_tensor_func(T &t, const char * pName = nullptr)
+  dump_tensor(T &t, const char * pName = nullptr)
   {
     std::cout << "Dumping non-tensor object ";
-    if( pName )
+    if( pName ) std::cout << pName;
-      std::cout << pName;
     std::cout << "=" << t;
   }
-  
-  // Helper to dump a tensor in memory order
-  // Kind of superfluous given the above ... just keeping in case I need to fall back to this
-#define DumpMemoryOrder(args...) DumpMemoryOrder_func(args)
-  template <typename T>
-  typename std::enable_if<EigenIO::is_tensor_of_scalar<T>::value, void>::type
-  DumpMemoryOrder_func(T &t, const char * pName = nullptr)
-  {
-    const auto rank = t.rank();
-    const auto &dims = t.dimensions();
-    std::cout << "Dumping rank " << rank << ((T::Options & Eigen::RowMajor) ? ", row" : ", column") << "-major tensor ";
-    if( pName )
-      std::cout << pName;
-    for( auto d : dims ) std::cout << "[" << d << "]";
-    std::cout << " in memory order:" << std::endl;
-    const typename T::Scalar * p = t.data();
-    const auto size = t.size();
-    const typename T::Scalar * pEnd = p + size;
-    if( rank <= 2 ) {
-      for( unsigned int i = 0 ; i < t.size() ; i++ )
-        std::cout << "[" << i << "]=" << *p++ << " ";
-      std::cout << std::endl;
-    } else {
-      const auto innersize = dims[rank-2] * dims[rank-1];
-      using Index = typename T::Index;
-      std::vector<Index> idx(rank - 2);
-      for( auto &i : idx ) i = 0;
-      Index idxCounter = 0;
-      while( p < pEnd ) {
-        if( T::Options & Eigen::RowMajor ) {
-          if( pName )
-            std::cout << pName;
-          idxCounter = 0;
-          for(auto i = 0 ; i < rank - 2 ; i++)
-            std::cout << "[" << idx[i] << "]:";
-        }
-        for( unsigned int i = 0 ; i < innersize ; i++ )
-          std::cout << " [" << idxCounter++ << "]=" << *p++;
-        if( T::Options & Eigen::RowMajor )
-          std::cout << std::endl;
-        // Now increment MyIndex
-        for( auto i = rank - 3; i != -1 && ++idx[i] == dims[i]; i-- )
-          idx[i] = 0;
-      }
-      if( ! ( T::Options & Eigen::RowMajor ) )
-        std::cout << std::endl;
-    }
-  }
-#else
-#define dump_tensor(args...)
-#define DumpMemoryOrder(args...)
-#endif
 }
 #endif

tests/IO/Test_serialisation.cc

@@ -29,7 +29,7 @@ Author: Michael Marshall <michael.marshall@ed.ac.uk>
     *************************************************************************************/
     /*  END LEGAL */
 #include <Grid/Grid.h>
-#include <Grid/util/Eigen.h>
+#include <Grid/util/EigenUtil.h>
 
 using namespace Grid;
 using namespace Grid::QCD;
@@ -158,8 +158,8 @@ public:
 
 template <typename WTR_, typename RDR_, typename T, typename... IndexTypes>
 void EigenTensorTestSingle(const char * MyTypeName, typename EigenIO::Traits<T>::scalar_type Flag,
-                           unsigned short Precision, std::string &filename, const char * pszExtension, unsigned int &TestNum,
+                           unsigned short Precision, std::string &filename, const char * pszExtension,
-                           IndexTypes... otherDims)
+                           unsigned int &TestNum, IndexTypes... otherDims)
 {
   using Traits = EigenIO::Traits<T>;
   using scalar_type = typename Traits::scalar_type;
@@ -187,15 +187,17 @@ void EigenTensorTest(const char * pszExtension, unsigned short Precision = 0)
   EigenTensorTestSingle<WTR_, RDR_, Tensor_9_4_2>(TEST_PARAMS( Tensor_9_4_2 ));
   {
     unsigned short Flag = 1;
-    TensorRank5UShort t;
     EigenTensorTestSingle<WTR_, RDR_, TensorRank5UShort>(TEST_PARAMS( TensorRank5UShort ));
     std::cout << "    Testing alternate memory order read ... ";
     TensorRank5UShortAlt t2;
     RDR_ reader(filename);
     read(reader, "TensorRank5UShort", t2);
     bool good = true;
-    for_all( t2, [&](unsigned short c, unsigned short n,
+    using Index = typename TensorRank5UShortAlt::Index;
-                     const std::array<size_t, TensorRank5UShortAlt::NumIndices> &Dims ) {
+    // NB: I can't call
+    for_all( t2, [&](unsigned short c, Index n,
+                     const std::array<Index, TensorRank5UShortAlt::NumIndices> &TensorIndex,
+                     const std::array<int, EigenIO::Traits<TensorRank5UShortAlt>::Rank> &GridTensorIndex ){
       good = good && ( c == n );
     } );
     if (!good) {
@@ -240,20 +242,6 @@ void tensorConvTestFn(GridSerialRNG &rng, const std::string label)
 
 int main(int argc,char **argv)
 {
-  {
-    LSCTensor Bingo;
-    constexpr Complex Flag{1,-3.1415927};
-    Complex z{0};
-    SpinColourVector * pV = Bingo.data();
-    for( std::size_t i = Bingo.size(); i--; ) {
-      for( typename GridTypeMapper<SpinColourVector>::scalar_type &s : *pV++ ) {
-        s = z;
-        z += Flag;
-      }
-    }
-    dump_tensor( Bingo );
-  }
-
   Grid_init(&argc,&argv);
   std::cout << std::boolalpha << "==== basic IO" << std::endl; // display true / false for boolean