diff --git a/.gitignore b/.gitignore index 5338acb9..40156f9d 100644 --- a/.gitignore +++ b/.gitignore @@ -88,6 +88,7 @@ Thumbs.db # build directory # ################### build*/* +Documentation/_build # IDE related files # ##################### diff --git a/Grid/serialisation/BaseIO.cc b/Grid/serialisation/BaseIO.cc new file mode 100644 index 00000000..9afc20b3 --- /dev/null +++ b/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 + +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 + +NAMESPACE_BEGIN(Grid) + +std::uint64_t EigenIO::EigenResizeCounter(0); + +NAMESPACE_END(Grid) diff --git a/Grid/serialisation/BaseIO.h b/Grid/serialisation/BaseIO.h index 49406201..25481301 100644 --- a/Grid/serialisation/BaseIO.h +++ b/Grid/serialisation/BaseIO.h @@ -9,6 +9,7 @@ Author: Antonin Portelli Author: Peter Boyle Author: Guido Cossu +Author: Michael Marshall 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 #ifndef GRID_SERIALISATION_ABSTRACT_READER_H #define GRID_SERIALISATION_ABSTRACT_READER_H +#include #include #include #include @@ -110,6 +112,10 @@ namespace Grid { template inline typename std::enable_if::value, typename Traits::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::value, void>::type Reader::Reshape(ETensor &t, const std::array &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(t.size()) * sizeof(typename ETensor::Scalar); //t.reshape( dims ); t.resize( dims ); + EigenIO::EigenResizeCounter += static_cast(t.size()) * sizeof(typename ETensor::Scalar); } template diff --git a/Grid/serialisation/Hdf5IO.cc b/Grid/serialisation/Hdf5IO.cc index 77396809..d3c7061e 100644 --- a/Grid/serialisation/Hdf5IO.cc +++ b/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 + Author: Peter Boyle + Author: Guido Cossu + Author: Michael Marshall + + 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 using namespace Grid; diff --git a/Grid/serialisation/Hdf5IO.h b/Grid/serialisation/Hdf5IO.h index 19537599..46cb07e1 100644 --- a/Grid/serialisation/Hdf5IO.h +++ b/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 + Author: Antonin Portelli + Author: Guido Cossu + Author: Michael Marshall + + 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 void writeDefault(const std::string &s, const U &x); template - typename std::enable_if>::is_number, void>::type + void writeRagged(const std::string &s, const std::vector &x); + template + typename std::enable_if>::value>::type writeDefault(const std::string &s, const std::vector &x); template - typename std::enable_if>::is_number, void>::type - writeDefault(const std::string &s, const std::vector &x); + typename std::enable_if>::value>::type + writeDefault(const std::string &s, const std::vector &x) { writeRagged(s, x); } template void writeMultiDim(const std::string &s, const std::vector & Dimensions, const U * pDataRowMajor, size_t NumElements); H5NS::Group & getGroup(void); @@ -64,11 +97,13 @@ namespace Grid template void readDefault(const std::string &s, U &output); template - typename std::enable_if>::is_number, void>::type + void readRagged(const std::string &s, std::vector &x); + template + typename std::enable_if>::value>::type readDefault(const std::string &s, std::vector &x); template - typename std::enable_if>::is_number, void>::type - readDefault(const std::string &s, std::vector &x); + typename std::enable_if>::value>::type + readDefault(const std::string &s, std::vector &x) { readRagged(s, x); } template void readMultiDim(const std::string &s, std::vector &buf, std::vector &dim); H5NS::Group & getGroup(void); @@ -176,24 +211,30 @@ namespace Grid } template - typename std::enable_if>::is_number, void>::type + typename std::enable_if>::value>::type Hdf5Writer::writeDefault(const std::string &s, const std::vector &x) { - // alias to element type - typedef typename element>::type Element; - - // flatten the vector and getting dimensions - Flatten> flat(x); - std::vector dim; - const auto &flatx = flat.getFlatVector(); - for (auto &d: flat.getDim()) - dim.push_back(d); - writeMultiDim(s, dim, &flatx[0], flatx.size()); + if (isRegularShape(x)) + { + // alias to element type + using Scalar = typename is_flattenable>::type; + + // flatten the vector and getting dimensions + Flatten> flat(x); + std::vector dim; + const auto &flatx = flat.getFlatVector(); + for (auto &d: flat.getDim()) + dim.push_back(d); + writeMultiDim(s, dim, &flatx[0], flatx.size()); + } + else + { + writeRagged(s, x); + } } template - typename std::enable_if>::is_number, void>::type - Hdf5Writer::writeDefault(const std::string &s, const std::vector &x) + void Hdf5Writer::writeRagged(const std::string &s, const std::vector &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 &buf, std::vector &dim) { // alias to element type - typedef typename element>::type Element; + using Scalar = typename is_flattenable>::type; // read the dimensions H5NS::DataSpace dataSpace; @@ -260,37 +301,44 @@ namespace Grid H5NS::DataSet dataSet; dataSet = group_.openDataSet(s); - dataSet.read(buf.data(), Hdf5Type::type()); + dataSet.read(buf.data(), Hdf5Type::type()); } else { H5NS::Attribute attribute; attribute = group_.openAttribute(s); - attribute.read(Hdf5Type::type(), buf.data()); + attribute.read(Hdf5Type::type(), buf.data()); } } template - typename std::enable_if>::is_number, void>::type + typename std::enable_if>::value>::type Hdf5Reader::readDefault(const std::string &s, std::vector &x) { - // alias to element type - typedef typename element>::type Element; + 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 + using Scalar = typename is_flattenable>::type; - std::vector dim; - std::vector buf; - readMultiDim( s, buf, dim ); + std::vector dim; + std::vector buf; + readMultiDim( s, buf, dim ); - // reconstruct the multidimensional vector - Reconstruct> r(buf, dim); - - x = r.getVector(); + // reconstruct the multidimensional vector + Reconstruct> r(buf, dim); + + x = r.getVector(); + } } template - typename std::enable_if>::is_number, void>::type - Hdf5Reader::readDefault(const std::string &s, std::vector &x) + void Hdf5Reader::readRagged(const std::string &s, std::vector &x) { uint64_t size; diff --git a/Grid/serialisation/MacroMagic.h b/Grid/serialisation/MacroMagic.h index 0495b91e..de456305 100644 --- a/Grid/serialisation/MacroMagic.h +++ b/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 \ -static inline void write(Writer &WR,const std::string &s, const cname &obj){ \ +static inline void write(::Grid::Writer &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 \ -static inline void read(Reader &RD,const std::string &s, cname &obj){ \ +static inline void read(::Grid::Reader &RD,const std::string &s, cname &obj){ \ if (!push(RD,s))\ {\ std::cout << ::Grid::GridLogWarning << "IO: Cannot open node '" << s << "'" << std::endl; \ diff --git a/Grid/serialisation/VectorUtils.h b/Grid/serialisation/VectorUtils.h index dd5ff0b8..8f490c64 100644 --- a/Grid/serialisation/VectorUtils.h +++ b/Grid/serialisation/VectorUtils.h @@ -9,7 +9,8 @@ Author: Antonin Portelli Author: Peter Boyle Author: paboyle - + Author: Michael Marshall + 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 @@ -236,21 +237,36 @@ namespace Grid { } } - // Vector element trait ////////////////////////////////////////////////////// - template - struct element + // is_flattenable::value is true if T is a std::vector<> which can be flattened ////////////////////// + template + struct is_flattenable : std::false_type { - typedef T type; - static constexpr bool is_number = false; + using type = T; + using grid_type = T; + static constexpr int vecRank = 0; + static constexpr bool isGridTensor = false; + static constexpr bool children_flattenable = std::is_arithmetic::value or is_complex::value; }; - + template - struct element> + struct is_flattenable::value>::type> : std::false_type { - typedef typename element::type type; - static constexpr bool is_number = std::is_arithmetic::value - or is_complex::value - or element::is_number; + using type = typename GridTypeMapper::scalar_type; + using grid_type = T; + static constexpr int vecRank = 0; + static constexpr bool isGridTensor = true; + static constexpr bool children_flattenable = true; + }; + + template + struct is_flattenable, typename std::enable_if::children_flattenable>::type> + : std::true_type + { + using type = typename is_flattenable::type; + using grid_type = typename is_flattenable::grid_type; + static constexpr bool isGridTensor = is_flattenable::isGridTensor; + static constexpr int vecRank = is_flattenable::vecRank + 1; + static constexpr bool children_flattenable = true; }; // Vector flattening utility class //////////////////////////////////////////// @@ -259,23 +275,30 @@ namespace Grid { class Flatten { public: - typedef typename element::type Element; + using Scalar = typename is_flattenable::type; + static constexpr bool isGridTensor = is_flattenable::isGridTensor; public: - explicit Flatten(const V &vector); - const V & getVector(void); - const std::vector & getFlatVector(void); - const std::vector & getDim(void); + explicit Flatten(const V &vector); + const V & getVector(void) const { return vector_; } + const std::vector & getFlatVector(void) const { return flatVector_; } + const std::vector & getDim(void) const { return dim_; } private: - void accumulate(const Element &e); - template - void accumulate(const W &v); - void accumulateDim(const Element &e); - template - void accumulateDim(const W &v); + template typename std::enable_if::value && !is_flattenable::isGridTensor>::type + accumulate(const W &e); + template typename std::enable_if::value && is_flattenable::isGridTensor>::type + accumulate(const W &e); + template typename std::enable_if< is_flattenable::value>::type + accumulate(const W &v); + template typename std::enable_if::value && !is_flattenable::isGridTensor>::type + accumulateDim(const W &e) {} // Innermost is a scalar - do nothing + template typename std::enable_if::value && is_flattenable::isGridTensor>::type + accumulateDim(const W &e); + template typename std::enable_if< is_flattenable::value>::type + accumulateDim(const W &v); private: - const V &vector_; - std::vector flatVector_; - std::vector dim_; + const V &vector_; + std::vector flatVector_; + std::vector dim_; }; // Class to reconstruct a multidimensional std::vector @@ -283,38 +306,57 @@ namespace Grid { class Reconstruct { public: - typedef typename element::type Element; + using Scalar = typename is_flattenable::type; + static constexpr bool isGridTensor = is_flattenable::isGridTensor; public: - Reconstruct(const std::vector &flatVector, + Reconstruct(const std::vector &flatVector, const std::vector &dim); - const V & getVector(void); - const std::vector & getFlatVector(void); - const std::vector & getDim(void); + const V & getVector(void) const { return vector_; } + const std::vector & getFlatVector(void) const { return flatVector_; } + const std::vector & getDim(void) const { return dim_; } private: - void fill(std::vector &v); - template - void fill(W &v); - void resize(std::vector &v, const unsigned int dim); - template - void resize(W &v, const unsigned int dim); + template typename std::enable_if::value && !is_flattenable::isGridTensor>::type + fill(W &v); + template typename std::enable_if::value && is_flattenable::isGridTensor>::type + fill(W &v); + template typename std::enable_if< is_flattenable::value>::type + fill(W &v); + template typename std::enable_if< is_flattenable::value && is_flattenable::vecRank==1>::type + resize(W &v, const unsigned int dim); + template typename std::enable_if< is_flattenable::value && (is_flattenable::vecRank>1)>::type + resize(W &v, const unsigned int dim); + template typename std::enable_if::isGridTensor>::type + checkInnermost(const W &e) {} // Innermost is a scalar - do nothing + template typename std::enable_if< is_flattenable::isGridTensor>::type + checkInnermost(const W &e); private: - V vector_; - const std::vector &flatVector_; - std::vector dim_; - size_t ind_{0}; - unsigned int dimInd_{0}; + V vector_; + const std::vector &flatVector_; + std::vector dim_; + size_t ind_{0}; + unsigned int dimInd_{0}; }; // Flatten class template implementation template - void Flatten::accumulate(const Element &e) + template typename std::enable_if::value && !is_flattenable::isGridTensor>::type + Flatten::accumulate(const W &e) { flatVector_.push_back(e); } template - template - void Flatten::accumulate(const W &v) + template typename std::enable_if::value && is_flattenable::isGridTensor>::type + Flatten::accumulate(const W &e) + { + for (const Scalar &x: e) { + flatVector_.push_back(x); + } + } + + template + template typename std::enable_if::value>::type + Flatten::accumulate(const W &v) { for (auto &e: v) { @@ -323,11 +365,17 @@ namespace Grid { } template - void Flatten::accumulateDim(const Element &e) {}; + template typename std::enable_if::value && is_flattenable::isGridTensor>::type + Flatten::accumulateDim(const W &e) + { + using Traits = GridTypeMapper::grid_type>; + for (int rank=0; rank < Traits::Rank; ++rank) + dim_.push_back(Traits::Dimension(rank)); + } template - template - void Flatten::accumulateDim(const W &v) + template typename std::enable_if::value>::type + Flatten::accumulateDim(const W &v) { dim_.push_back(v.size()); accumulateDim(v[0]); @@ -337,42 +385,36 @@ namespace Grid { Flatten::Flatten(const V &vector) : vector_(vector) { - accumulate(vector_); accumulateDim(vector_); - } - - template - const V & Flatten::getVector(void) - { - return vector_; - } - - template - const std::vector::Element> & - Flatten::getFlatVector(void) - { - return flatVector_; - } - - template - const std::vector & Flatten::getDim(void) - { - return dim_; + std::size_t TotalSize{ dim_[0] }; + for (int i = 1; i < dim_.size(); ++i) { + TotalSize *= dim_[i]; + } + flatVector_.reserve(TotalSize); + accumulate(vector_); } // Reconstruct class template implementation template - void Reconstruct::fill(std::vector &v) + template typename std::enable_if::value && !is_flattenable::isGridTensor>::type + Reconstruct::fill(W &v) + { + v = flatVector_[ind_++]; + } + + template + template typename std::enable_if::value && is_flattenable::isGridTensor>::type + Reconstruct::fill(W &v) { for (auto &e: v) { e = flatVector_[ind_++]; } } - + template - template - void Reconstruct::fill(W &v) + template typename std::enable_if::value>::type + Reconstruct::fill(W &v) { for (auto &e: v) { @@ -381,14 +423,15 @@ namespace Grid { } template - void Reconstruct::resize(std::vector &v, const unsigned int dim) + template typename std::enable_if::value && is_flattenable::vecRank==1>::type + Reconstruct::resize(W &v, const unsigned int dim) { v.resize(dim_[dim]); } template - template - void Reconstruct::resize(W &v, const unsigned int dim) + template typename std::enable_if::value && (is_flattenable::vecRank>1)>::type + Reconstruct::resize(W &v, const unsigned int dim) { v.resize(dim_[dim]); for (auto &e: v) @@ -398,34 +441,31 @@ namespace Grid { } template - Reconstruct::Reconstruct(const std::vector &flatVector, + template typename std::enable_if::isGridTensor>::type + Reconstruct::checkInnermost(const W &) + { + using Traits = GridTypeMapper::grid_type>; + const int gridRank{Traits::Rank}; + const int dimRank{static_cast(dim_.size())}; + assert(dimRank >= gridRank && "Tensor rank too low for Grid tensor"); + for (int i=0; i + Reconstruct::Reconstruct(const std::vector &flatVector, const std::vector &dim) : flatVector_(flatVector) , dim_(dim) { + checkInnermost(vector_); + assert(dim_.size() == is_flattenable::vecRank && "Tensor rank doesn't match nested std::vector rank"); resize(vector_, 0); fill(vector_); } - template - const V & Reconstruct::getVector(void) - { - return vector_; - } - - template - const std::vector::Element> & - Reconstruct::getFlatVector(void) - { - return flatVector_; - } - - template - const std::vector & Reconstruct::getDim(void) - { - return dim_; - } - // Vector IO utilities /////////////////////////////////////////////////////// // helper function to read space-separated values template @@ -459,6 +499,64 @@ namespace Grid { return os; } + + // In general, scalar types are considered "flattenable" (regularly shaped) + template + bool isRegularShapeHelper(const std::vector &, std::vector &, int, bool) + { + return true; + } + + template + bool isRegularShapeHelper(const std::vector> &v, std::vector &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 + bool isRegularShape(const T &t) { return true; } + + template + bool isRegularShape(const std::vector &v) { return !v.empty(); } + + // Return non-zero if all dimensions of this std::vector> are regularly shaped + template + bool isRegularShape(const std::vector> &v) + { + if (v.empty() || v[0].empty()) + return false; + // Make sure all of my rows are the same size + std::vector Dims; + Dims.reserve(is_flattenable::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 diff --git a/Grid/tensors/Tensor_class.h b/Grid/tensors/Tensor_class.h index 36becc49..be045ede 100644 --- a/Grid/tensors/Tensor_class.h +++ b/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"; diff --git a/Grid/threads/Accelerator.cc b/Grid/threads/Accelerator.cc index 32a6693e..52fadc0b 100644 --- a/Grid/threads/Accelerator.cc +++ b/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"); diff --git a/configure.ac b/configure.ac index 721d890e..406b0b74 100644 --- a/configure.ac +++ b/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) diff --git a/documentation/Grid.pdf b/documentation/Grid.pdf index 868c6db4..df3304eb 100644 Binary files a/documentation/Grid.pdf and b/documentation/Grid.pdf differ diff --git a/documentation/manual.rst b/documentation/manual.rst index d51f07c1..e545bdaf 100644 --- a/documentation/manual.rst +++ b/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 struct is_tensor : std::integral_constant, T>::value> {}; + // Is this an Eigen tensor of a supported scalar + template struct is_tensor_of_scalar : public std::false_type {}; + template struct is_tensor_of_scalar::value && is_scalar::value>::type> : public std::true_type {}; + // Is this an Eigen tensor of a supported container + template struct is_tensor_of_container : public std::false_type {}; + template struct is_tensor_of_container::value && isGridTensor::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 + void readMultiDim(const std::string &s, std::vector &buf, std::vector &dim); + template + void writeMultiDim(const std::string &s, const std::vector & 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 +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``, where ``T`` is an arithmetic or complex type, +the Hdf5Writer writes the data as an Hdf5 DataSet object. + +However, nested ``std::vector>`` 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 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`` 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 ----------------------- diff --git a/tests/IO/Test_serialisation.cc b/tests/IO/Test_serialisation.cc index 27fe589e..e1596ea6 100644 --- a/tests/IO/Test_serialisation.cc +++ b/tests/IO/Test_serialisation.cc @@ -48,7 +48,9 @@ public: std::vector, array, std::vector >, twodimarray, std::vector> > >, cmplx3darray, - SpinColourMatrix, scm + SpinColourMatrix, scm, + std::vector > >, ragged, + std::vector >, vscm ); myclass() {} myclass(int i) @@ -56,6 +58,10 @@ public: , twodimarray(3,std::vector(5, 1.23456)) , cmplx3darray(3,std::vector>>(5, std::vector>(7, std::complex(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(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(Counter, -Counter); + Counter++; + } + } } }; diff --git a/tests/Test_meson_field.cc b/tests/Test_meson_field.cc new file mode 100644 index 00000000..25d908d7 --- /dev/null +++ b/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 + +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 +#include + +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> t(&grid); + LatticeCoordinate(t, Tp); + std::vector seeds({1,2,3,4}); + GridParallelRNG pRNG(&grid); + pRNG.SeedFixedIntegers(seeds); + + // MesonField lhs and rhs vectors + std::vector phi(VDIM,&grid); + std::vector 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 Gmu = { + Gamma::Algebra::GammaX, + Gamma::Algebra::GammaY, + Gamma::Algebra::GammaZ, + Gamma::Algebra::GammaT + }; + + // momentum phases e^{ipx} + std::vector> 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 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 Mpp(momenta.size(),Gmu.size(),Nt,VDIM,VDIM); + Eigen::Tensor Mpr(momenta.size(),Gmu.size(),Nt,VDIM,VDIM); + Eigen::Tensor Mrr(momenta.size(),Gmu.size(),Nt,VDIM,VDIM); + + // timer + double start,stop; + + //execute meson field routine + start = usecond(); + A2Autils::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::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::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; +}