Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-10-23 09:13:38

0001 /*
0002  *  Copyright (c), 2017, Adrien Devresse <adrien.devresse@epfl.ch>
0003  *
0004  *  Distributed under the Boost Software License, Version 1.0.
0005  *    (See accompanying file LICENSE_1_0.txt or copy at
0006  *          http://www.boost.org/LICENSE_1_0.txt)
0007  *
0008  */
0009 #pragma once
0010 
0011 #include <cstdlib>
0012 #include <vector>
0013 
0014 #include "H5_definitions.hpp"
0015 #include "H5Utils.hpp"
0016 
0017 #include "../H5PropertyList.hpp"
0018 #include "h5s_wrapper.hpp"
0019 
0020 namespace HighFive {
0021 
0022 class ElementSet {
0023   public:
0024     ///
0025     /// \brief Create a list of points of N-dimension for selection.
0026     ///
0027     /// \param list List of continuous coordinates (e.g.: in 2 dimensions space
0028     /// `ElementSet{1, 2, 3 ,4}` creates points `(1, 2)` and `(3, 4)`).
0029     explicit ElementSet(std::initializer_list<std::size_t> list);
0030     ///
0031     /// \brief Create a list of points of N-dimension for selection.
0032     ///
0033     /// \param list List of N-dim points.
0034     explicit ElementSet(std::initializer_list<std::vector<std::size_t>> list);
0035     ///
0036     /// \brief Create a list of points of N-dimension for selection.
0037     ///
0038     /// \param element_ids List of continuous coordinates (e.g.: in 2 dimensions space
0039     /// `ElementSet{1, 2, 3 ,4}` creates points `(1, 2)` and `(3, 4)`).
0040     explicit ElementSet(const std::vector<std::size_t>& element_ids);
0041     ///
0042     /// \brief Create a list of points of N-dimension for selection.
0043     ///
0044     /// \param element_ids List of N-dim points.
0045     explicit ElementSet(const std::vector<std::vector<std::size_t>>& element_ids);
0046 
0047   private:
0048     std::vector<std::size_t> _ids;
0049 
0050     template <typename Derivate>
0051     friend class SliceTraits;
0052 };
0053 
0054 namespace detail {
0055 
0056 template <class To, class From>
0057 inline std::vector<To> convertSizeVector(const std::vector<From>& from) {
0058     std::vector<To> to(from.size());
0059     std::copy(from.cbegin(), from.cend(), to.begin());
0060 
0061     return to;
0062 }
0063 }  // namespace detail
0064 
0065 inline std::vector<hsize_t> toHDF5SizeVector(const std::vector<size_t>& from) {
0066     return detail::convertSizeVector<hsize_t>(from);
0067 }
0068 
0069 inline std::vector<size_t> toSTLSizeVector(const std::vector<hsize_t>& from) {
0070     return detail::convertSizeVector<size_t>(from);
0071 }
0072 
0073 struct RegularHyperSlab {
0074     RegularHyperSlab() = default;
0075 
0076     RegularHyperSlab(std::vector<size_t> offset_,
0077                      std::vector<size_t> count_ = {},
0078                      std::vector<size_t> stride_ = {},
0079                      std::vector<size_t> block_ = {})
0080         : offset(toHDF5SizeVector(offset_))
0081         , count(toHDF5SizeVector(count_))
0082         , stride(toHDF5SizeVector(stride_))
0083         , block(toHDF5SizeVector(block_)) {}
0084 
0085     static RegularHyperSlab fromHDF5Sizes(std::vector<hsize_t> offset_,
0086                                           std::vector<hsize_t> count_ = {},
0087                                           std::vector<hsize_t> stride_ = {},
0088                                           std::vector<hsize_t> block_ = {}) {
0089         RegularHyperSlab slab;
0090         slab.offset = offset_;
0091         slab.count = count_;
0092         slab.stride = stride_;
0093         slab.block = block_;
0094 
0095         return slab;
0096     }
0097 
0098     size_t rank() const {
0099         return std::max(std::max(offset.size(), count.size()),
0100                         std::max(stride.size(), block.size()));
0101     }
0102 
0103     /// Dimensions when all gaps are removed.
0104     std::vector<size_t> packedDims() const {
0105         auto n_dims = rank();
0106         auto dims = std::vector<size_t>(n_dims, 0);
0107 
0108         for (size_t i = 0; i < n_dims; ++i) {
0109             dims[i] = count[i] * (block.empty() ? 1 : block[i]);
0110         }
0111 
0112         return dims;
0113     }
0114 
0115     std::vector<hsize_t> offset;
0116     std::vector<hsize_t> count;
0117     std::vector<hsize_t> stride;
0118     std::vector<hsize_t> block;
0119 };
0120 
0121 class HyperSlab {
0122   public:
0123     HyperSlab() {
0124         selects.emplace_back(RegularHyperSlab{}, Op::None);
0125     };
0126 
0127     explicit HyperSlab(const RegularHyperSlab& sel) {
0128         selects.emplace_back(sel, Op::Set);
0129     }
0130 
0131     HyperSlab operator|(const RegularHyperSlab& sel) const {
0132         auto ret = *this;
0133         ret |= sel;
0134         return ret;
0135     }
0136 
0137     HyperSlab& operator|=(const RegularHyperSlab& sel) {
0138         selects.emplace_back(sel, Op::Or);
0139         return *this;
0140     }
0141 
0142     HyperSlab operator&(const RegularHyperSlab& sel) const {
0143         auto ret = *this;
0144         ret &= sel;
0145         return ret;
0146     }
0147 
0148     HyperSlab& operator&=(const RegularHyperSlab& sel) {
0149         selects.emplace_back(sel, Op::And);
0150         return *this;
0151     }
0152 
0153     HyperSlab operator^(const RegularHyperSlab& sel) const {
0154         auto ret = *this;
0155         ret ^= sel;
0156         return ret;
0157     }
0158 
0159     HyperSlab& operator^=(const RegularHyperSlab& sel) {
0160         selects.emplace_back(sel, Op::Xor);
0161         return *this;
0162     }
0163 
0164     HyperSlab& notA(const RegularHyperSlab& sel) {
0165         selects.emplace_back(sel, Op::NotA);
0166         return *this;
0167     }
0168 
0169     HyperSlab& notB(const RegularHyperSlab& sel) {
0170         selects.emplace_back(sel, Op::NotB);
0171         return *this;
0172     }
0173 
0174     DataSpace apply(const DataSpace& space_) const {
0175         return apply_impl(space_);
0176     }
0177 
0178   private:
0179     enum class Op {
0180         Noop,
0181         Set,
0182         Or,
0183         And,
0184         Xor,
0185         NotB,
0186         NotA,
0187         Append,
0188         Prepend,
0189         Invalid,
0190         None,
0191     };
0192 
0193     H5S_seloper_t convert(Op op) const {
0194         switch (op) {
0195         case Op::Noop:
0196             return H5S_SELECT_NOOP;
0197         case Op::Set:
0198             return H5S_SELECT_SET;
0199         case Op::Or:
0200             return H5S_SELECT_OR;
0201         case Op::And:
0202             return H5S_SELECT_AND;
0203         case Op::Xor:
0204             return H5S_SELECT_XOR;
0205         case Op::NotB:
0206             return H5S_SELECT_NOTB;
0207         case Op::NotA:
0208             return H5S_SELECT_NOTA;
0209         case Op::Append:
0210             return H5S_SELECT_APPEND;
0211         case Op::Prepend:
0212             return H5S_SELECT_PREPEND;
0213         case Op::Invalid:
0214             return H5S_SELECT_INVALID;
0215         default:
0216             throw DataSpaceException("Invalid HyperSlab operation.");
0217         }
0218     }
0219 
0220     struct Select_: public RegularHyperSlab {
0221         Select_(const RegularHyperSlab& sel, Op op_)
0222             : RegularHyperSlab(sel)
0223             , op(op_) {}
0224 
0225         Op op;
0226     };
0227 
0228     std::vector<Select_> selects;
0229 
0230   protected:
0231     DataSpace select_none(const DataSpace& outer_space) const {
0232         auto space = outer_space.clone();
0233         detail::h5s_select_none(space.getId());
0234         return space;
0235     }
0236 
0237     void select_hyperslab(DataSpace& space, const Select_& sel) const {
0238         detail::h5s_select_hyperslab(space.getId(),
0239                                      convert(sel.op),
0240                                      sel.offset.empty() ? nullptr : sel.offset.data(),
0241                                      sel.stride.empty() ? nullptr : sel.stride.data(),
0242                                      sel.count.empty() ? nullptr : sel.count.data(),
0243                                      sel.block.empty() ? nullptr : sel.block.data());
0244     }
0245 
0246 #if H5_VERSION_GE(1, 10, 6)
0247     /// The length of a stream of `Op::Or` starting at `begin`.
0248     size_t detect_streak(Select_ const* begin, Select_ const* end, Op op) const {
0249         assert(op == Op::Or);
0250         auto it = std::find_if(begin, end, [op](const Select_& sel) { return sel.op != op; });
0251         return static_cast<size_t>(it - begin);
0252     }
0253 
0254     DataSpace combine_selections(const DataSpace& left_space,
0255                                  Op op,
0256                                  const DataSpace& right_space) const {
0257         assert(op == Op::Or);
0258 
0259         auto left_type = detail::h5s_get_select_type(left_space.getId());
0260         auto right_type = detail::h5s_get_select_type(right_space.getId());
0261 
0262         // Since HDF5 doesn't allow `combine_selections` with a None
0263         // selection, we need to avoid the issue:
0264         if (left_type == H5S_SEL_NONE) {
0265             return right_space;
0266         } else if (right_type == H5S_SEL_NONE) {
0267             return left_space;
0268         } else if (left_type == H5S_SEL_ALL) {
0269             return left_space;
0270         } else if (right_type == H5S_SEL_ALL) {
0271             return right_space;
0272         } else {
0273             return detail::make_data_space(
0274                 detail::h5s_combine_select(left_space.getId(), convert(op), right_space.getId()));
0275         }
0276     }
0277 
0278     /// Reduce a sequence of `Op::Or` efficiently.
0279     ///
0280     /// The issue is that `H5Sselect_hyperslab` runs in time that linear of the
0281     /// number of block in the existing selection. Therefore, a loop that adds
0282     /// slab-by-slab has quadratic runtime in the number of slabs.
0283     ///
0284     /// Fortunately, `H5Scombine_select` doesn't suffer from the same problem.
0285     /// However, it's only available in 1.10.6 and newer.
0286     ///
0287     /// The solution is to use divide-and-conquer to reduce (long) streaks of
0288     /// `Op::Or` in what seems to be log-linear time.
0289     DataSpace reduce_streak(const DataSpace& outer_space,
0290                             Select_ const* begin,
0291                             Select_ const* end,
0292                             Op op) const {
0293         assert(op == Op::Or);
0294 
0295         if (begin == end) {
0296             throw std::runtime_error("Broken logic in 'DataSpace::reduce_streak'.");
0297         }
0298 
0299         std::ptrdiff_t distance = end - begin;
0300         if (distance == 1) {
0301             auto space = select_none(outer_space);
0302             select_hyperslab(space, *begin);
0303             return space;
0304         }
0305 
0306         Select_ const* mid = begin + distance / 2;
0307         auto right_space = reduce_streak(outer_space, begin, mid, op);
0308         auto left_space = reduce_streak(outer_space, mid, end, op);
0309 
0310         return combine_selections(left_space, op, right_space);
0311     }
0312 
0313     DataSpace apply_impl(const DataSpace& space_) const {
0314         auto space = space_.clone();
0315         auto n_selects = selects.size();
0316         for (size_t i = 0; i < n_selects; ++i) {
0317             auto begin = selects.data() + i;
0318             auto end = selects.data() + n_selects;
0319 
0320             auto n_ors = detect_streak(begin, end, Op::Or);
0321 
0322             if (n_ors > 1) {
0323                 auto right_space = reduce_streak(space_, begin, begin + n_ors, Op::Or);
0324                 space = combine_selections(space, Op::Or, right_space);
0325                 i += n_ors - 1;
0326             } else if (selects[i].op == Op::None) {
0327                 detail::h5s_select_none(space.getId());
0328             } else {
0329                 select_hyperslab(space, selects[i]);
0330             }
0331         }
0332         return space;
0333     }
0334 #else
0335     DataSpace apply_impl(const DataSpace& space_) const {
0336         auto space = space_.clone();
0337         for (const auto& sel: selects) {
0338             if (sel.op == Op::None) {
0339                 detail::h5s_select_none(space.getId());
0340             } else {
0341                 select_hyperslab(space, sel);
0342             }
0343         }
0344         return space;
0345     }
0346 #endif
0347 };
0348 
0349 template <typename Derivate>
0350 class SliceTraits {
0351   public:
0352     ///
0353     /// \brief Select an \p hyperslab in the current Slice/Dataset.
0354     ///
0355     /// HyperSlabs can be either regular or irregular. Irregular hyperslabs are typically generated
0356     /// by taking the union of regular hyperslabs. An irregular hyperslab, in general, does not fit
0357     /// nicely into a multi-dimensional array, but only a subset of such an array.
0358     ///
0359     /// Therefore, the only memspaces supported for general hyperslabs are one-dimensional arrays.
0360     Selection select(const HyperSlab& hyperslab) const;
0361 
0362     ///
0363     /// \brief Select an \p hyperslab in the current Slice/Dataset.
0364     ///
0365     /// If the selection can be read into a simple, multi-dimensional dataspace,
0366     /// then this overload enable specifying the shape of the memory dataspace
0367     /// with `memspace`. Note, that simple implies no offsets, strides or
0368     /// number of blocks, just the size of the block in each dimension.
0369     Selection select(const HyperSlab& hyperslab, const DataSpace& memspace) const;
0370 
0371     ///
0372     /// \brief Select a region in the current Slice/Dataset of \p count points at
0373     /// \p offset separated by \p stride. If strides are not provided they will
0374     /// default to 1 in all dimensions.
0375     ///
0376     /// vector offset and count have to be from the same dimension
0377     ///
0378     Selection select(const std::vector<size_t>& offset,
0379                      const std::vector<size_t>& count,
0380                      const std::vector<size_t>& stride = {},
0381                      const std::vector<size_t>& block = {}) const;
0382 
0383     ///
0384     /// \brief Select a set of columns in the last dimension of this dataset.
0385     ///
0386     /// The column indices must be smaller than the dimension size.
0387     ///
0388     Selection select(const std::vector<size_t>& columns) const;
0389 
0390     ///
0391     /// \brief Select a region in the current Slice/Dataset out of a list of elements.
0392     ///
0393     Selection select(const ElementSet& elements) const;
0394 
0395     template <typename T>
0396     T read(const DataTransferProps& xfer_props = DataTransferProps()) const;
0397 
0398     ///
0399     /// Read the entire dataset into a buffer
0400     ///
0401     /// An exception is raised is if the numbers of dimension of the buffer and
0402     /// of the dataset are different.
0403     ///
0404     /// The array type can be a N-pointer or a N-vector. For plain pointers
0405     /// not dimensionality checking will be performed, it is the user's
0406     /// responsibility to ensure that the right amount of space has been
0407     /// allocated.
0408     template <typename T>
0409     void read(T& array, const DataTransferProps& xfer_props = DataTransferProps()) const;
0410 
0411     ///
0412     /// Read the entire dataset into a raw buffer
0413     ///
0414     /// \deprecated Use `read_raw` instead.
0415     ///
0416     /// No dimensionality checks will be performed, it is the user's
0417     /// responsibility to ensure that the right amount of space has been
0418     /// allocated.
0419     /// \param array: A buffer containing enough space for the data
0420     /// \param dtype: The datatype of elements of the in memory buffer.
0421     /// \param xfer_props: Data Transfer properties
0422     template <typename T>
0423     void read(T* array,
0424               const DataType& dtype,
0425               const DataTransferProps& xfer_props = DataTransferProps()) const;
0426 
0427     ///
0428     /// Read the entire dataset into a raw buffer
0429     ///
0430     /// \deprecated Use `read_raw` instead.
0431     ///
0432     /// Same as `read(T*, const DataType&, const DataTransferProps&)`. However,
0433     /// this overload deduces the HDF5 datatype of the element of `array` from
0434     /// `T`. Note, that the file datatype is already fixed.
0435     ///
0436     /// \param array: A buffer containing enough space for the data
0437     /// \param xfer_props: Data Transfer properties
0438     template <typename T>
0439     void read(T* array, const DataTransferProps& xfer_props = DataTransferProps()) const;
0440 
0441     ///
0442     /// Read the entire dataset into a raw buffer
0443     ///
0444     /// No dimensionality checks will be performed, it is the user's
0445     /// responsibility to ensure that the right amount of space has been
0446     /// allocated.
0447     /// \param array: A buffer containing enough space for the data
0448     /// \param dtype: The type of the data, in case it cannot be automatically guessed
0449     /// \param xfer_props: Data Transfer properties
0450     template <typename T>
0451     void read_raw(T* array,
0452                   const DataType& dtype,
0453                   const DataTransferProps& xfer_props = DataTransferProps()) const;
0454 
0455     ///
0456     /// Read the entire dataset into a raw buffer
0457     ///
0458     /// Same as `read(T*, const DataType&, const DataTransferProps&)`. However,
0459     /// this overload deduces the HDF5 datatype of the element of `array` from
0460     /// `T`. Note, that the file datatype is already fixed.
0461     ///
0462     /// \param array: A buffer containing enough space for the data
0463     /// \param xfer_props: Data Transfer properties
0464     template <typename T>
0465     void read_raw(T* array, const DataTransferProps& xfer_props = DataTransferProps()) const;
0466 
0467 
0468     ///
0469     /// Write the integrality N-dimension buffer to this dataset
0470     /// An exception is raised is if the numbers of dimension of the buffer and
0471     /// of the dataset are different
0472     ///
0473     /// The array type can be a N-pointer or a N-vector ( e.g int** integer two
0474     /// dimensional array )
0475     template <typename T>
0476     void write(const T& buffer, const DataTransferProps& xfer_props = DataTransferProps());
0477 
0478     ///
0479     /// Write from a raw pointer into this dataset.
0480     ///
0481     /// No dimensionality checks will be performed, it is the user's
0482     /// responsibility to ensure that the buffer holds the right amount of
0483     /// elements. For n-dimensional matrices the buffer layout follows H5
0484     /// default conventions.
0485     ///
0486     /// Note, this is the shallowest wrapper around `H5Dwrite` and should
0487     /// be used if full control is needed. Generally prefer `write`.
0488     ///
0489     /// \param buffer: A buffer containing the data to be written
0490     /// \param dtype: The datatype of `buffer`, i.e. the memory data type.
0491     /// \param xfer_props: The HDF5 data transfer properties, e.g. collective MPI-IO.
0492     template <typename T>
0493     void write_raw(const T* buffer,
0494                    const DataType& mem_datatype,
0495                    const DataTransferProps& xfer_props = DataTransferProps());
0496 
0497     ///
0498     /// Write from a raw pointer into this dataset.
0499     ///
0500     /// Same as `write_raw(const T*, const DataTransferProps&)`. However, this
0501     /// overload attempts to guess the data type of `buffer`, i.e. the memory
0502     /// datatype. Note that the file datatype is already fixed.
0503     ///
0504     template <typename T>
0505     void write_raw(const T* buffer, const DataTransferProps& xfer_props = DataTransferProps());
0506 };
0507 
0508 }  // namespace HighFive