diff --git a/Grid/serialisation/BaseIO.h b/Grid/serialisation/BaseIO.h index 043762ea..2a780fb4 100644 --- a/Grid/serialisation/BaseIO.h +++ b/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(malloc(TotalNumElements * sizeof(Scalar))); + pCopyBuffer = new Scalar[TotalNumElements]; pWriteBuffer = pCopyBuffer; Scalar * pCopy = pCopyBuffer; std::array MyIndex; @@ -330,7 +330,7 @@ namespace Grid { } } upcast->template writeMultiDim(s, TotalDims, pWriteBuffer, TotalNumElements); - if( pCopyBuffer ) free( pCopyBuffer ); + if( pCopyBuffer ) delete [] pCopyBuffer; } template diff --git a/Grid/serialisation/XmlIO.h b/Grid/serialisation/XmlIO.h index af3d5d76..40d5f2bb 100644 --- a/Grid/serialisation/XmlIO.h +++ b/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( Dimensions.size() ); write("rank", Rank ); + std::vector 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 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(); } diff --git a/Grid/util/Eigen.h b/Grid/util/EigenUtil.h similarity index 50% rename from Grid/util/Eigen.h rename to Grid/util/EigenUtil.h index 9af510fe..c6f10276 100644 --- a/Grid/util/Eigen.h +++ b/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 -#define GRID_UTIL_EIGEN_H +#ifndef GRID_UTIL_EIGENUTIL_H +#define GRID_UTIL_EIGENUTIL_H #include #include @@ -35,25 +35,28 @@ namespace Grid { template typename std::enable_if::value, void>::type for_all_do_lambda( Lambda lambda, typename ETensor::Scalar &scalar, typename ETensor::Index &Seq, - std::array::Rank> &MyIndex) + const std::array &MyIndex, + const std::array::Rank> &DimGridTensor, + std::array::Rank> &GridTensorIndex) { - lambda( scalar, Seq++, MyIndex ); + lambda( scalar, Seq++, MyIndex, GridTensorIndex ); } // for_all helper function to call the lambda for container template typename std::enable_if::value, void>::type for_all_do_lambda( Lambda lambda, typename ETensor::Scalar &container, typename ETensor::Index &Seq, - std::array::Rank> &MyIndex) + const std::array &MyIndex, + const std::array::Rank> &DimGridTensor, + std::array::Rank> &GridTensorIndex) { using Traits = EigenIO::Traits; - 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-- ) - MyIndex[rank + i] = 0; + for( auto i = InnerRank - 1; i != -1 && ++GridTensorIndex[i] == DimGridTensor[i]; i-- ) + 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; - const std::size_t NumScalars = ET.size(); - assert( NumScalars > 0 ); using Index = typename ETensor::Index; + using Traits = EigenIO::Traits; + // 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 Dims; - for(auto i = 0; i < rank; i++ ) { - auto dim = ET.dimension(i); - assert( dim > 0 ); - Dims[i] = static_cast(dim); - assert( Dims[i] == dim ); // check we didn't lose anything in the conversion - ScalarElementCount *= Dims[i]; + std::array DimTensor, MyIndex; + for(int i = 0; i < rank; i++ ) { + DimTensor[i] = ET.dimension(i); + ScalarElementCount *= DimTensor[i]; + MyIndex[i] = 0; } - // Check that the number of containers is correct ... and we didn't lose anything in conversions - assert( NumScalars == ScalarElementCount ); - // If the Scalar is actually a container, add the inner Scalar's dimensions - size_t InnerScalarCount{1}; - for(auto i = 0; i < InnerRank; i++ ) { - auto dim = Traits::Dimension(i); - assert( dim > 0 ); - Dims[rank + i] = static_cast(dim); - assert( Dims[rank + i] == dim ); // check we didn't lose anything in the conversion - InnerScalarCount *= dim; + assert( NumScalars == ScalarElementCount && "EigenUtil: tensor size not product of dimensions" ); + // Save the GridTensor dimensions + const auto InnerRank{Traits::Rank}; + std::array DimGridTensor, GridTensorIndex; + for(int i = 0; i < InnerRank; i++ ) { + DimGridTensor[i] = Traits::Dimension(i); + GridTensorIndex[i] = 0; } - assert(Traits::count == InnerScalarCount); - std::array MyIndex; - for( auto &idx : MyIndex ) idx = 0; + // Now walk the tensor in memory order Index Seq = 0; Scalar * pScalar = ET.data(); - for( std::size_t j = 0; j < NumScalars; j++ ) { - for_all_do_lambda( lambda, * pScalar, Seq, MyIndex ); + for( Index j = 0; j < NumScalars; j++ ) { + for_all_do_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; - for( auto i = 0; i < rank + InnerRank ; i++ ) { - NewSeq *= Dims[i]; - NewSeq += MyIndex[i]; + Seq = 0; + for( auto i = 0; i < rank; i++ ) { + Seq *= DimTensor[i]; + Seq += MyIndex[i]; } - Seq = static_cast( 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 std::enable_if::value && !EigenIO::Traits::is_complex>::type @@ -124,11 +120,13 @@ namespace Grid { { using Traits = EigenIO::Traits; using scalar_type = typename Traits::scalar_type; - for_all( ET, [&](scalar_type &c, typename ETensor::Index n, const std::array &Dims ) { + using Index = typename ETensor::Index; + for_all( ET, [&](scalar_type &c, Index n, const std::array &TensorIndex, + const std::array &ScalarIndex ) { scalar_type x = Inc * static_cast(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; using scalar_type = typename Traits::scalar_type; - for_all( ET, [&](scalar_type &c, typename ETensor::Index n, const std::array &Dims ) { + using Index = typename ETensor::Index; + for_all( ET, [&](scalar_type &c, Index n, const std::array &TensorIndex, + const std::array &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 std::enable_if::value, void>::type - dump_tensor_func(T &t, const char * pName = nullptr) + dump_tensor(T &t, const char * pName = nullptr) { using Traits = EigenIO::Traits; + 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 &Dims ){ + for_all( t, [&](scalar_type &c, Index n, const std::array &TensorIndex, + const std::array &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 std::enable_if::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 ) - std::cout << pName; + if( 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 std::enable_if::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 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 diff --git a/tests/IO/Test_serialisation.cc b/tests/IO/Test_serialisation.cc index 56cdcb4b..229e7f12 100644 --- a/tests/IO/Test_serialisation.cc +++ b/tests/IO/Test_serialisation.cc @@ -29,7 +29,7 @@ Author: Michael Marshall *************************************************************************************/ /* END LEGAL */ #include -#include +#include using namespace Grid; using namespace Grid::QCD; @@ -158,8 +158,8 @@ public: template void EigenTensorTestSingle(const char * MyTypeName, typename EigenIO::Traits::scalar_type Flag, - unsigned short Precision, std::string &filename, const char * pszExtension, unsigned int &TestNum, - IndexTypes... otherDims) + unsigned short Precision, std::string &filename, const char * pszExtension, + unsigned int &TestNum, IndexTypes... otherDims) { using Traits = EigenIO::Traits; using scalar_type = typename Traits::scalar_type; @@ -187,15 +187,17 @@ void EigenTensorTest(const char * pszExtension, unsigned short Precision = 0) EigenTensorTestSingle(TEST_PARAMS( Tensor_9_4_2 )); { unsigned short Flag = 1; - TensorRank5UShort t; EigenTensorTestSingle(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, - const std::array &Dims ) { + using Index = typename TensorRank5UShortAlt::Index; + // NB: I can't call + for_all( t2, [&](unsigned short c, Index n, + const std::array &TensorIndex, + const std::array::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::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