diff --git a/packages/HighFive/.gitrepo b/packages/HighFive/.gitrepo
index ba30619099f608abd928af1efa5324e3a37389c9..0512eb59369fa0287eb768cdb9215c29160d65a0 100644
--- a/packages/HighFive/.gitrepo
+++ b/packages/HighFive/.gitrepo
@@ -6,7 +6,7 @@
 [subrepo]
 	remote = git@github.com:BlueBrain/HighFive.git
 	branch = master
-	commit = b74fabbea7bc4422a6891b2f5a27305b8dd8cd1b
-	parent = d64d966d026b0bc07529a0328b8607036ecc24d6
+	commit = 20fbd4c4d328bc18a342eeef5ffd29f17de084e0
+	parent = 35908346b20b524d4d2d9f8f6220800c27595413
 	method = merge
 	cmdver = 0.4.6
diff --git a/packages/HighFive/deps/catch2 b/packages/HighFive/deps/catch2
index 3f0283de7a9c43200033da996ff9093be3ac84dc..4e8d92bf02f7d1c8006a0e7a5ecabd8e62d98502 160000
--- a/packages/HighFive/deps/catch2
+++ b/packages/HighFive/deps/catch2
@@ -1 +1 @@
-Subproject commit 3f0283de7a9c43200033da996ff9093be3ac84dc
+Subproject commit 4e8d92bf02f7d1c8006a0e7a5ecabd8e62d98502
diff --git a/packages/HighFive/doc/developer_guide.md b/packages/HighFive/doc/developer_guide.md
index 0f7af3e01695ea8d307f5f4207c1ff9bc63b1bfd..4537bef20bb87a6d78dbaa4d78d3ae2d387b4741 100644
--- a/packages/HighFive/doc/developer_guide.md
+++ b/packages/HighFive/doc/developer_guide.md
@@ -165,7 +165,7 @@ void check_write(...) {
 ### Test Organization
 #### Multi-Dimensional Arrays
 All tests for reading/writing whole multi-dimensional arrays to datasets or
-attributes belong in `tests/unit/tests_high_five_multi_dimensional.cpp`. This
+attributes belong in `tests/unit/test_all_types.cpp`. This
 includes write/read cycles; checking all the generic edges cases, e.g. empty
 arrays and mismatching sizes; and checking non-reallocation.
 
diff --git a/packages/HighFive/include/highfive/H5Easy.hpp b/packages/HighFive/include/highfive/H5Easy.hpp
index e793fd8753179263fb235129d717739abeb362fd..013c44135edf2c49680d962f1bcf6aa07d44b9dd 100644
--- a/packages/HighFive/include/highfive/H5Easy.hpp
+++ b/packages/HighFive/include/highfive/H5Easy.hpp
@@ -30,6 +30,7 @@
 #ifdef H5_USE_XTENSOR
 #include <xtensor/xarray.hpp>
 #include <xtensor/xtensor.hpp>
+#include "xtensor.hpp"
 #endif
 
 // optionally enable Eigen plug-in and load the library
@@ -41,6 +42,7 @@
 
 #ifdef H5_USE_EIGEN
 #include <Eigen/Eigen>
+#include "eigen.hpp"
 #endif
 
 // optionally enable OpenCV plug-in and load the library
@@ -52,6 +54,7 @@
 
 #ifdef H5_USE_OPENCV
 #include <opencv2/opencv.hpp>
+#include "experimental/opencv.hpp"
 #endif
 
 #include "H5File.hpp"
@@ -393,8 +396,5 @@ inline T loadAttribute(const File& file, const std::string& path, const std::str
 
 #include "h5easy_bits/H5Easy_Eigen.hpp"
 #include "h5easy_bits/H5Easy_misc.hpp"
-#include "h5easy_bits/H5Easy_opencv.hpp"
 #include "h5easy_bits/H5Easy_public.hpp"
 #include "h5easy_bits/H5Easy_scalar.hpp"
-#include "h5easy_bits/H5Easy_vector.hpp"
-#include "h5easy_bits/H5Easy_xtensor.hpp"
diff --git a/packages/HighFive/include/highfive/bits/H5Slice_traits.hpp b/packages/HighFive/include/highfive/bits/H5Slice_traits.hpp
index e2b481a36a6a242154c0b97eab87228a3bd27951..6812a091449ee5aa7c3343a439b123171fa56f8c 100644
--- a/packages/HighFive/include/highfive/bits/H5Slice_traits.hpp
+++ b/packages/HighFive/include/highfive/bits/H5Slice_traits.hpp
@@ -231,6 +231,85 @@ class HyperSlab {
     std::vector<Select_> selects;
 };
 
+///
+/// \brief Selects the Cartesian product of slices.
+///
+/// Given a one-dimensional dataset one might want to select the union of
+/// multiple, non-overlapping slices. For example,
+///
+///     using Slice = std::array<size_t, 2>;
+///     using Slices = std::vector<Slice>;
+///     auto slices = Slices{{0, 2}, {4, 10}};
+///     dset.select(ProductSet(slices);
+///
+/// to select elements `0`, `1` and `4`, ..., `9` (inclusive).
+///
+/// For a two-dimensional array one which to select the row specified above,
+/// but only columns `2`, `3` and `4`:
+///
+///    dset.select(ProductSet(slices, Slice{2, 5}));
+///    // Analogues with the roles of columns and rows reversed:
+///    dset.select(ProductSet(Slice{2, 5}, slices));
+///
+/// One can generalize once more and allow the unions of slices in both x- and
+/// y-dimension:
+///
+///     auto yslices = Slices{{1, 5}, {7, 8}};
+///     auto xslices = Slices{{0, 3}, {6, 8}};
+///     dset.select(ProductSet(yslices, xslices));
+///
+/// which selects the following from a 11x8 dataset:
+///
+///     . . . . . . . .
+///     x x x . . . x x
+///     x x x . . . x x
+///     x x x . . . x x
+///     x x x . . . x x
+///     . . . . . . . .
+///     . . . . . . . .
+///     x x x . . . x x
+///     . . . . . . . .
+///     . . . . . . . .
+///     . . . . . . . .
+///
+/// Final twist, the selection along and axis may be discrete indices, from
+/// which a vector of, possibly single-element, slices can be constructed. The
+/// corresponding types are `std::vector<size_t>` and `size_t` for multiple or
+/// just a single values. Note, looping over rows or columns one-by-one can be
+/// a very serious performance problem. In particular,
+///
+///     // Avoid:
+///     for(auto i : indices) {
+///         dset.select(i).read<double>();
+///     }
+///
+///     // Use:
+///     std::vector<size_t> tmp(indices.begin(), indices.end());
+///     dset.select(tmp).read<std::vector<double>>();
+///
+/// obvious omit the copy if `indices` already has the correct type.
+///
+/// The solution works analogous in higher dimensions. A selection `sk` along
+/// axis `k` can be interpreted as a subset `S_k` of the natural numbers. The
+/// index `i` is in `S_k` if it's selected by `sk`.  The `ProductSet` of `s0`,
+/// ..., `sN` selects the Cartesian product `S_0 x ... x S_N`.
+///
+/// Note that the selections along each axis must be sorted and non-overlapping.
+///
+class ProductSet {
+  public:
+    template <class... Slices>
+    explicit ProductSet(const Slices&... slices);
+
+  private:
+    HyperSlab slab;
+    std::vector<size_t> shape;
+
+    template <typename Derivate>
+    friend class SliceTraits;
+};
+
+
 template <typename Derivate>
 class SliceTraits {
   public:
@@ -277,6 +356,8 @@ class SliceTraits {
     ///
     Selection select(const ElementSet& elements) const;
 
+    Selection select(const ProductSet& product_set) const;
+
     template <typename T>
     T read(const DataTransferProps& xfer_props = DataTransferProps()) const;
 
diff --git a/packages/HighFive/include/highfive/bits/H5Slice_traits_misc.hpp b/packages/HighFive/include/highfive/bits/H5Slice_traits_misc.hpp
index 711d57ac2489c2ffb8d4448da3c6d7dcf1e035ac..c1bc9c9ea1751217336dbf591767061859a8b673 100644
--- a/packages/HighFive/include/highfive/bits/H5Slice_traits_misc.hpp
+++ b/packages/HighFive/include/highfive/bits/H5Slice_traits_misc.hpp
@@ -66,6 +66,153 @@ inline ElementSet::ElementSet(const std::vector<std::vector<std::size_t>>& eleme
     }
 }
 
+namespace detail {
+class HyperCube {
+  public:
+    HyperCube(size_t rank)
+        : offset(rank)
+        , count(rank) {}
+
+    void cross(const std::array<size_t, 2>& range, size_t axis) {
+        offset[axis] = range[0];
+        count[axis] = range[1] - range[0];
+    }
+
+    RegularHyperSlab asSlab() {
+        return RegularHyperSlab(offset, count);
+    }
+
+  private:
+    std::vector<size_t> offset;
+    std::vector<size_t> count;
+};
+
+inline void build_hyper_slab(HyperSlab& slab, size_t /* axis */, HyperCube& cube) {
+    slab |= cube.asSlab();
+}
+
+template <class... Slices>
+inline void build_hyper_slab(HyperSlab& slab,
+                             size_t axis,
+                             HyperCube& cube,
+                             const std::array<size_t, 2>& slice,
+                             const Slices&... higher_slices) {
+    cube.cross(slice, axis);
+    build_hyper_slab(slab, axis + 1, cube, higher_slices...);
+}
+
+template <class... Slices>
+inline void build_hyper_slab(HyperSlab& slab,
+                             size_t axis,
+                             HyperCube& cube,
+                             const std::vector<std::array<size_t, 2>>& slices,
+                             const Slices&... higher_slices) {
+    for (const auto& slice: slices) {
+        build_hyper_slab(slab, axis, cube, slice, higher_slices...);
+    }
+}
+
+template <class... Slices>
+inline void build_hyper_slab(HyperSlab& slab,
+                             size_t axis,
+                             HyperCube& cube,
+                             const std::vector<size_t>& ids,
+                             const Slices&... higher_slices) {
+    for (const auto& id: ids) {
+        auto slice = std::array<size_t, 2>{id, id + 1};
+        build_hyper_slab(slab, axis, cube, slice, higher_slices...);
+    }
+}
+
+template <class... Slices>
+inline void build_hyper_slab(HyperSlab& slab,
+                             size_t axis,
+                             HyperCube& cube,
+                             size_t id,
+                             const Slices&... higher_slices) {
+    auto slice = std::array<size_t, 2>{id, id + 1};
+    build_hyper_slab(slab, axis, cube, slice, higher_slices...);
+}
+
+inline void compute_squashed_shape(size_t /* axis */, std::vector<size_t>& /* shape */) {
+    // assert(axis == shape.size());
+}
+
+template <class... Slices>
+inline void compute_squashed_shape(size_t axis,
+                                   std::vector<size_t>& shape,
+                                   const std::array<size_t, 2>& slice,
+                                   const Slices&... higher_slices);
+
+template <class... Slices>
+inline void compute_squashed_shape(size_t axis,
+                                   std::vector<size_t>& shape,
+                                   const std::vector<size_t>& points,
+                                   const Slices&... higher_slices);
+
+template <class... Slices>
+inline void compute_squashed_shape(size_t axis,
+                                   std::vector<size_t>& shape,
+                                   size_t point,
+                                   const Slices&... higher_slices);
+
+template <class... Slices>
+inline void compute_squashed_shape(size_t axis,
+                                   std::vector<size_t>& shape,
+                                   const std::vector<std::array<size_t, 2>>& slices,
+                                   const Slices&... higher_slices);
+
+template <class... Slices>
+inline void compute_squashed_shape(size_t axis,
+                                   std::vector<size_t>& shape,
+                                   const std::array<size_t, 2>& slice,
+                                   const Slices&... higher_slices) {
+    shape[axis] = slice[1] - slice[0];
+    compute_squashed_shape(axis + 1, shape, higher_slices...);
+}
+
+template <class... Slices>
+inline void compute_squashed_shape(size_t axis,
+                                   std::vector<size_t>& shape,
+                                   const std::vector<size_t>& points,
+                                   const Slices&... higher_slices) {
+    shape[axis] = points.size();
+    compute_squashed_shape(axis + 1, shape, higher_slices...);
+}
+
+template <class... Slices>
+inline void compute_squashed_shape(size_t axis,
+                                   std::vector<size_t>& shape,
+                                   const std::vector<std::array<size_t, 2>>& slices,
+                                   const Slices&... higher_slices) {
+    shape[axis] = 0;
+    for (const auto& slice: slices) {
+        shape[axis] += slice[1] - slice[0];
+    }
+    compute_squashed_shape(axis + 1, shape, higher_slices...);
+}
+
+template <class... Slices>
+inline void compute_squashed_shape(size_t axis,
+                                   std::vector<size_t>& shape,
+                                   size_t /* point */,
+                                   const Slices&... higher_slices) {
+    shape[axis] = 1;
+    compute_squashed_shape(axis + 1, shape, higher_slices...);
+}
+}  // namespace detail
+
+template <class... Slices>
+inline ProductSet::ProductSet(const Slices&... slices) {
+    auto rank = sizeof...(slices);
+    detail::HyperCube cube(rank);
+    detail::build_hyper_slab(slab, 0, cube, slices...);
+
+    shape = std::vector<size_t>(rank, size_t(0));
+    detail::compute_squashed_shape(0, shape, slices...);
+}
+
+
 template <typename Derivate>
 inline Selection SliceTraits<Derivate>::select(const HyperSlab& hyperslab,
                                                const DataSpace& memspace) const {
@@ -157,6 +304,11 @@ inline Selection SliceTraits<Derivate>::select(const ElementSet& elements) const
     return detail::make_selection(DataSpace(num_elements), space, details::get_dataset(slice));
 }
 
+template <typename Derivate>
+inline Selection SliceTraits<Derivate>::select(const ProductSet& product_set) const {
+    return this->select(product_set.slab, DataSpace(product_set.shape));
+}
+
 
 template <typename Derivate>
 template <typename T>
diff --git a/packages/HighFive/include/highfive/h5easy_bits/H5Easy_Eigen.hpp b/packages/HighFive/include/highfive/h5easy_bits/H5Easy_Eigen.hpp
index a2e23bc0fa7c3ed5fdf3f09988127d0cc66ca281..794725436d5a63703e07a1680c3152d2d4e31419 100644
--- a/packages/HighFive/include/highfive/h5easy_bits/H5Easy_Eigen.hpp
+++ b/packages/HighFive/include/highfive/h5easy_bits/H5Easy_Eigen.hpp
@@ -14,66 +14,48 @@
 
 #ifdef H5_USE_EIGEN
 
+#include "../eigen.hpp"
+
 namespace H5Easy {
 
 namespace detail {
 
 template <typename T>
 struct io_impl<T, typename std::enable_if<std::is_base_of<Eigen::DenseBase<T>, T>::value>::type> {
-    // abbreviate row-major <-> col-major conversions
-    template <typename S>
-    struct types {
-        using row_major = Eigen::Ref<
-            const Eigen::Array<typename std::decay<T>::type::Scalar,
-                               std::decay<T>::type::RowsAtCompileTime,
-                               std::decay<T>::type::ColsAtCompileTime,
-                               std::decay<T>::type::ColsAtCompileTime == 1 ? Eigen::ColMajor
-                                                                           : Eigen::RowMajor,
-                               std::decay<T>::type::MaxRowsAtCompileTime,
-                               std::decay<T>::type::MaxColsAtCompileTime>,
-            0,
-            Eigen::InnerStride<1>>;
-
-        using col_major =
-            Eigen::Map<Eigen::Array<typename std::decay<T>::type::Scalar,
-                                    std::decay<T>::type::RowsAtCompileTime,
-                                    std::decay<T>::type::ColsAtCompileTime,
-                                    std::decay<T>::type::ColsAtCompileTime == 1 ? Eigen::ColMajor
-                                                                                : Eigen::RowMajor,
-                                    std::decay<T>::type::MaxRowsAtCompileTime,
-                                    std::decay<T>::type::MaxColsAtCompileTime>>;
-    };
-
-    // return the shape of Eigen::DenseBase<T> object as size 1 or 2 "std::vector<size_t>"
-    inline static std::vector<size_t> shape(const T& data) {
+    using EigenIndex = Eigen::DenseIndex;
+
+    // When creating a dataset for an Eigen object, the shape of the dataset is
+    // 1D for vectors. (legacy reasons)
+    inline static std::vector<size_t> file_shape(const T& data) {
         if (std::decay<T>::type::RowsAtCompileTime == 1) {
             return {static_cast<size_t>(data.cols())};
         }
         if (std::decay<T>::type::ColsAtCompileTime == 1) {
             return {static_cast<size_t>(data.rows())};
         }
-        return {static_cast<size_t>(data.rows()), static_cast<size_t>(data.cols())};
+        return inspector<T>::getDimensions(data);
     }
 
-    using EigenIndex = Eigen::DenseIndex;
+    // The shape of an Eigen object as used in HighFive core.
+    inline static std::vector<size_t> mem_shape(const T& data) {
+        return inspector<T>::getDimensions(data);
+    }
 
-    // get the shape of a "DataSet" as size 2 "std::vector<Eigen::Index>"
+    // The shape of an Eigen object as used in HighFive core.
     template <class D>
-    inline static std::vector<EigenIndex> shape(const File& file,
+    inline static std::vector<size_t> mem_shape(const File& file,
                                                 const std::string& path,
-                                                const D& dataset,
-                                                int RowsAtCompileTime) {
+                                                const D& dataset) {
         std::vector<size_t> dims = dataset.getDimensions();
 
-        if (dims.size() == 1 && RowsAtCompileTime == 1) {
-            return std::vector<EigenIndex>{1u, static_cast<EigenIndex>(dims[0])};
+        if (dims.size() == 1 && T::RowsAtCompileTime == 1) {
+            return std::vector<size_t>{1, dims[0]};
         }
-        if (dims.size() == 1) {
-            return std::vector<EigenIndex>{static_cast<EigenIndex>(dims[0]), 1u};
+        if (dims.size() == 1 && T::ColsAtCompileTime == 1) {
+            return std::vector<size_t>{dims[0], 1};
         }
         if (dims.size() == 2) {
-            return std::vector<EigenIndex>{static_cast<EigenIndex>(dims[0]),
-                                           static_cast<EigenIndex>(dims[1])};
+            return dims;
         }
 
         throw detail::error(file, path, "H5Easy::load: Inconsistent rank");
@@ -83,11 +65,12 @@ struct io_impl<T, typename std::enable_if<std::is_base_of<Eigen::DenseBase<T>, T
                                const std::string& path,
                                const T& data,
                                const DumpOptions& options) {
-        using row_major_type = typename types<T>::row_major;
         using value_type = typename std::decay<T>::type::Scalar;
-        row_major_type row_major(data);
-        DataSet dataset = initDataset<value_type>(file, path, shape(data), options);
-        dataset.write_raw(row_major.data());
+
+        std::vector<size_t> file_dims = file_shape(data);
+        std::vector<size_t> mem_dims = mem_shape(data);
+        DataSet dataset = initDataset<value_type>(file, path, file_dims, options);
+        dataset.reshapeMemSpace(mem_dims).write(data);
         if (options.flush()) {
             file.flush();
         }
@@ -96,14 +79,8 @@ struct io_impl<T, typename std::enable_if<std::is_base_of<Eigen::DenseBase<T>, T
 
     inline static T load(const File& file, const std::string& path) {
         DataSet dataset = file.getDataSet(path);
-        std::vector<typename T::Index> dims = shape(file, path, dataset, T::RowsAtCompileTime);
-        T data(dims[0], dims[1]);
-        dataset.read_raw(data.data());
-        if (data.IsVectorAtCompileTime || data.IsRowMajor) {
-            return data;
-        }
-        using col_major = typename types<T>::col_major;
-        return col_major(data.data(), dims[0], dims[1]);
+        std::vector<size_t> dims = mem_shape(file, path, dataset);
+        return dataset.reshapeMemSpace(dims).template read<T>();
     }
 
     inline static Attribute dumpAttribute(File& file,
@@ -111,11 +88,12 @@ struct io_impl<T, typename std::enable_if<std::is_base_of<Eigen::DenseBase<T>, T
                                           const std::string& key,
                                           const T& data,
                                           const DumpOptions& options) {
-        using row_major_type = typename types<T>::row_major;
         using value_type = typename std::decay<T>::type::Scalar;
-        row_major_type row_major(data);
-        Attribute attribute = initAttribute<value_type>(file, path, key, shape(data), options);
-        attribute.write_raw(row_major.data());
+
+        std::vector<size_t> file_dims = file_shape(data);
+        std::vector<size_t> mem_dims = mem_shape(data);
+        Attribute attribute = initAttribute<value_type>(file, path, key, file_dims, options);
+        attribute.reshapeMemSpace(mem_dims).write(data);
         if (options.flush()) {
             file.flush();
         }
@@ -128,14 +106,8 @@ struct io_impl<T, typename std::enable_if<std::is_base_of<Eigen::DenseBase<T>, T
         DataSet dataset = file.getDataSet(path);
         Attribute attribute = dataset.getAttribute(key);
         DataSpace dataspace = attribute.getSpace();
-        std::vector<typename T::Index> dims = shape(file, path, dataspace, T::RowsAtCompileTime);
-        T data(dims[0], dims[1]);
-        attribute.read_raw(data.data());
-        if (data.IsVectorAtCompileTime || data.IsRowMajor) {
-            return data;
-        }
-        using col_major = typename types<T>::col_major;
-        return col_major(data.data(), dims[0], dims[1]);
+        std::vector<size_t> dims = mem_shape(file, path, dataspace);
+        return attribute.reshapeMemSpace(dims).template read<T>();
     }
 };
 
diff --git a/packages/HighFive/include/highfive/h5easy_bits/H5Easy_opencv.hpp b/packages/HighFive/include/highfive/h5easy_bits/H5Easy_opencv.hpp
deleted file mode 100644
index 7be19f1e35484def4a378d7b0e188b3aa40dbedc..0000000000000000000000000000000000000000
--- a/packages/HighFive/include/highfive/h5easy_bits/H5Easy_opencv.hpp
+++ /dev/null
@@ -1,100 +0,0 @@
-/*
- *  Copyright (c), 2017, Adrien Devresse <adrien.devresse@epfl.ch>
- *
- *  Distributed under the Boost Software License, Version 1.0.
- *    (See accompanying file LICENSE_1_0.txt or copy at
- *          http://www.boost.org/LICENSE_1_0.txt)
- *
- */
-#pragma once
-
-#include "../H5Easy.hpp"
-#include "H5Easy_misc.hpp"
-#include "H5Easy_scalar.hpp"
-
-#ifdef H5_USE_OPENCV
-
-namespace H5Easy {
-
-namespace detail {
-
-template <class T>
-struct is_opencv: std::false_type {};
-template <class T>
-struct is_opencv<cv::Mat_<T>>: std::true_type {};
-
-template <typename T>
-struct io_impl<T, typename std::enable_if<is_opencv<T>::value>::type> {
-    inline static std::vector<size_t> shape(const T& data) {
-        return std::vector<size_t>{static_cast<size_t>(data.rows), static_cast<size_t>(data.cols)};
-    }
-
-    inline static std::vector<int> shape(const File& file,
-                                         const std::string& path,
-                                         std::vector<size_t> dims) {
-        if (dims.size() == 1) {
-            return std::vector<int>{static_cast<int>(dims[0]), 1ul};
-        }
-        if (dims.size() == 2) {
-            return std::vector<int>{static_cast<int>(dims[0]), static_cast<int>(dims[1])};
-        }
-
-        throw detail::error(file, path, "H5Easy::load: Inconsistent rank");
-    }
-
-    inline static DataSet dump(File& file,
-                               const std::string& path,
-                               const T& data,
-                               const DumpOptions& options) {
-        using value_type = typename T::value_type;
-        DataSet dataset = initDataset<value_type>(file, path, shape(data), options);
-        std::vector<value_type> v(data.begin(), data.end());
-        dataset.write_raw(v.data());
-        if (options.flush()) {
-            file.flush();
-        }
-        return dataset;
-    }
-
-    inline static T load(const File& file, const std::string& path) {
-        using value_type = typename T::value_type;
-        DataSet dataset = file.getDataSet(path);
-        std::vector<int> dims = shape(file, path, dataset.getDimensions());
-        T data(dims[0], dims[1]);
-        dataset.read_raw(reinterpret_cast<value_type*>(data.data));
-        return data;
-    }
-
-    inline static Attribute dumpAttribute(File& file,
-                                          const std::string& path,
-                                          const std::string& key,
-                                          const T& data,
-                                          const DumpOptions& options) {
-        using value_type = typename T::value_type;
-        Attribute attribute = initAttribute<value_type>(file, path, key, shape(data), options);
-        std::vector<value_type> v(data.begin(), data.end());
-        attribute.write_raw(v.data());
-        if (options.flush()) {
-            file.flush();
-        }
-        return attribute;
-    }
-
-    inline static T loadAttribute(const File& file,
-                                  const std::string& path,
-                                  const std::string& key) {
-        using value_type = typename T::value_type;
-        DataSet dataset = file.getDataSet(path);
-        Attribute attribute = dataset.getAttribute(key);
-        DataSpace dataspace = attribute.getSpace();
-        std::vector<int> dims = shape(file, path, dataspace.getDimensions());
-        T data(dims[0], dims[1]);
-        attribute.read_raw(reinterpret_cast<value_type*>(data.data));
-        return data;
-    }
-};
-
-}  // namespace detail
-}  // namespace H5Easy
-
-#endif  // H5_USE_OPENCV
diff --git a/packages/HighFive/include/highfive/h5easy_bits/H5Easy_scalar.hpp b/packages/HighFive/include/highfive/h5easy_bits/H5Easy_scalar.hpp
index 056d8f2dc5aaa57f9754f94efed25155d7b2341a..0dd42799e798d6d9ddfef0c8816c4fc4d2d600a3 100644
--- a/packages/HighFive/include/highfive/h5easy_bits/H5Easy_scalar.hpp
+++ b/packages/HighFive/include/highfive/h5easy_bits/H5Easy_scalar.hpp
@@ -10,6 +10,7 @@
 
 #include "../H5Easy.hpp"
 #include "H5Easy_misc.hpp"
+#include "default_io_impl.hpp"
 
 namespace H5Easy {
 
@@ -20,49 +21,7 @@ Base template for partial specialization: the fallback if specialized templates
 Used e.g. for scalars.
 */
 template <typename T, typename = void>
-struct io_impl {
-    inline static DataSet dump(File& file,
-                               const std::string& path,
-                               const T& data,
-                               const DumpOptions& options) {
-        DataSet dataset = initScalarDataset(file, path, data, options);
-        dataset.write(data);
-        if (options.flush()) {
-            file.flush();
-        }
-        return dataset;
-    }
-
-    inline static T load(const File& file, const std::string& path) {
-        DataSet dataset = file.getDataSet(path);
-        T data;
-        dataset.read(data);
-        return data;
-    }
-
-    inline static Attribute dumpAttribute(File& file,
-                                          const std::string& path,
-                                          const std::string& key,
-                                          const T& data,
-                                          const DumpOptions& options) {
-        Attribute attribute = initScalarAttribute(file, path, key, data, options);
-        attribute.write(data);
-        if (options.flush()) {
-            file.flush();
-        }
-        return attribute;
-    }
-
-    inline static T loadAttribute(const File& file,
-                                  const std::string& path,
-                                  const std::string& key) {
-        DataSet dataset = file.getDataSet(path);
-        Attribute attribute = dataset.getAttribute(key);
-        T data;
-        attribute.read(data);
-        return data;
-    }
-
+struct io_impl: public default_io_impl<T> {
     inline static DataSet dump_extend(File& file,
                                       const std::string& path,
                                       const T& data,
@@ -93,7 +52,6 @@ struct io_impl {
             return dataset;
         }
 
-        std::vector<size_t> shape = idx;
         const size_t unlim = DataSpace::UNLIMITED;
         std::vector<size_t> unlim_shape(idx.size(), unlim);
         std::vector<hsize_t> chunks(idx.size(), 10);
@@ -103,8 +61,9 @@ struct io_impl {
                 throw error(file, path, "H5Easy::dump: Incorrect dimension ChunkSize");
             }
         }
-        for (size_t& i: shape) {
-            i++;
+        std::vector<size_t> shape(idx.size());
+        for (size_t i = 0; i < idx.size(); ++i) {
+            shape[i] = idx[i] + 1;
         }
         DataSpace dataspace = DataSpace(shape, unlim_shape);
         DataSetCreateProps props;
@@ -121,10 +80,7 @@ struct io_impl {
                               const std::string& path,
                               const std::vector<size_t>& idx) {
         std::vector<size_t> ones(idx.size(), 1);
-        DataSet dataset = file.getDataSet(path);
-        T data;
-        dataset.select(idx, ones).read(data);
-        return data;
+        return file.getDataSet(path).select(idx, ones).read<T>();
     }
 };
 
diff --git a/packages/HighFive/include/highfive/h5easy_bits/H5Easy_xtensor.hpp b/packages/HighFive/include/highfive/h5easy_bits/H5Easy_xtensor.hpp
deleted file mode 100644
index ba27bc84a6b1333279e413bacc7dbdfa4802244b..0000000000000000000000000000000000000000
--- a/packages/HighFive/include/highfive/h5easy_bits/H5Easy_xtensor.hpp
+++ /dev/null
@@ -1,94 +0,0 @@
-/*
- *  Copyright (c), 2017, Adrien Devresse <adrien.devresse@epfl.ch>
- *
- *  Distributed under the Boost Software License, Version 1.0.
- *    (See accompanying file LICENSE_1_0.txt or copy at
- *          http://www.boost.org/LICENSE_1_0.txt)
- *
- */
-#pragma once
-
-#include "../H5Easy.hpp"
-#include "H5Easy_misc.hpp"
-#include "H5Easy_scalar.hpp"
-
-#ifdef H5_USE_XTENSOR
-
-namespace H5Easy {
-
-namespace detail {
-
-template <typename T>
-struct io_impl<T, typename std::enable_if<xt::is_xexpression<T>::value>::type> {
-    inline static void assert_row_major(const File& file, const std::string& path, const T& data) {
-        if (data.layout() != xt::layout_type::row_major) {
-            throw detail::error(file, path, "Only row-major XTensor object are supported.");
-        }
-    }
-
-    inline static std::vector<size_t> shape(const T& data) {
-        return std::vector<size_t>(data.shape().cbegin(), data.shape().cend());
-    }
-
-    inline static DataSet dump(File& file,
-                               const std::string& path,
-                               const T& data,
-                               const DumpOptions& options) {
-        assert_row_major(file, path, data);
-        using value_type = typename std::decay_t<T>::value_type;
-        DataSet dataset = initDataset<value_type>(file, path, shape(data), options);
-        dataset.write_raw(data.data());
-        if (options.flush()) {
-            file.flush();
-        }
-        return dataset;
-    }
-
-    inline static T load(const File& file, const std::string& path) {
-        static_assert(
-            xt::has_data_interface<T>::value,
-            "Cannot load to xt::xfunction or xt::xgenerator, use e.g. xt::xtensor or xt::xarray");
-        DataSet dataset = file.getDataSet(path);
-        std::vector<size_t> dims = dataset.getDimensions();
-        T data = T::from_shape(dims);
-        assert_row_major(file, path, data);
-        dataset.read_raw(data.data());
-        return data;
-    }
-
-    inline static Attribute dumpAttribute(File& file,
-                                          const std::string& path,
-                                          const std::string& key,
-                                          const T& data,
-                                          const DumpOptions& options) {
-        assert_row_major(file, path, data);
-        using value_type = typename std::decay_t<T>::value_type;
-        Attribute attribute = initAttribute<value_type>(file, path, key, shape(data), options);
-        attribute.write_raw(data.data());
-        if (options.flush()) {
-            file.flush();
-        }
-        return attribute;
-    }
-
-    inline static T loadAttribute(const File& file,
-                                  const std::string& path,
-                                  const std::string& key) {
-        static_assert(
-            xt::has_data_interface<T>::value,
-            "Cannot load to xt::xfunction or xt::xgenerator, use e.g. xt::xtensor or xt::xarray");
-        DataSet dataset = file.getDataSet(path);
-        Attribute attribute = dataset.getAttribute(key);
-        DataSpace dataspace = attribute.getSpace();
-        std::vector<size_t> dims = dataspace.getDimensions();
-        T data = T::from_shape(dims);
-        assert_row_major(file, path, data);
-        attribute.read_raw(data.data());
-        return data;
-    }
-};
-
-}  // namespace detail
-}  // namespace H5Easy
-
-#endif  // H5_USE_XTENSOR
diff --git a/packages/HighFive/include/highfive/h5easy_bits/H5Easy_vector.hpp b/packages/HighFive/include/highfive/h5easy_bits/default_io_impl.hpp
similarity index 68%
rename from packages/HighFive/include/highfive/h5easy_bits/H5Easy_vector.hpp
rename to packages/HighFive/include/highfive/h5easy_bits/default_io_impl.hpp
index 4c60f5cfc632d2f77ecd55c219114ad17404241d..ccb6e71504b2590a413314d50738605957c6a899 100644
--- a/packages/HighFive/include/highfive/h5easy_bits/H5Easy_vector.hpp
+++ b/packages/HighFive/include/highfive/h5easy_bits/default_io_impl.hpp
@@ -13,28 +13,22 @@
 #include "H5Easy_scalar.hpp"
 
 namespace H5Easy {
-
 namespace detail {
 
-template <class T>
-struct is_vector: std::false_type {};
-template <class T>
-struct is_vector<std::vector<T>>: std::true_type {};
-
 using HighFive::details::inspector;
 
 template <typename T>
-struct io_impl<T, typename std::enable_if<is_vector<T>::value>::type> {
+struct default_io_impl {
+    inline static std::vector<size_t> shape(const T& data) {
+        return inspector<T>::getDimensions(data);
+    }
+
     inline static DataSet dump(File& file,
                                const std::string& path,
                                const T& data,
                                const DumpOptions& options) {
         using value_type = typename inspector<T>::base_type;
-        auto dims = inspector<T>::getDimensions(data);
-        DataSet dataset = initDataset<value_type>(file,
-                                                  path,
-                                                  std::vector<size_t>(dims.begin(), dims.end()),
-                                                  options);
+        DataSet dataset = initDataset<value_type>(file, path, shape(data), options);
         dataset.write(data);
         if (options.flush()) {
             file.flush();
@@ -43,10 +37,7 @@ struct io_impl<T, typename std::enable_if<is_vector<T>::value>::type> {
     }
 
     inline static T load(const File& file, const std::string& path) {
-        DataSet dataset = file.getDataSet(path);
-        T data;
-        dataset.read(data);
-        return data;
+        return file.getDataSet(path).read<T>();
     }
 
     inline static Attribute dumpAttribute(File& file,
@@ -55,9 +46,7 @@ struct io_impl<T, typename std::enable_if<is_vector<T>::value>::type> {
                                           const T& data,
                                           const DumpOptions& options) {
         using value_type = typename inspector<T>::base_type;
-        auto dims = inspector<T>::getDimensions(data);
-        std::vector<size_t> shape(dims.begin(), dims.end());
-        Attribute attribute = initAttribute<value_type>(file, path, key, shape, options);
+        Attribute attribute = initAttribute<value_type>(file, path, key, shape(data), options);
         attribute.write(data);
         if (options.flush()) {
             file.flush();
@@ -70,9 +59,7 @@ struct io_impl<T, typename std::enable_if<is_vector<T>::value>::type> {
                                   const std::string& key) {
         DataSet dataset = file.getDataSet(path);
         Attribute attribute = dataset.getAttribute(key);
-        T data;
-        attribute.read(data);
-        return data;
+        return attribute.read<T>();
     }
 };
 
diff --git a/packages/HighFive/src/examples/CMakeLists.txt b/packages/HighFive/src/examples/CMakeLists.txt
index 5f4b4edf862dffc6fad9db89e42c6ee0f12cfa8f..fa6f5f739ace8c585765025d2f14a4179a6ff83b 100644
--- a/packages/HighFive/src/examples/CMakeLists.txt
+++ b/packages/HighFive/src/examples/CMakeLists.txt
@@ -71,7 +71,7 @@ foreach(example_source ${core_examples})
 endforeach()
 
 foreach(example_source ${easy_examples})
-  compile_example(${example_source})
+  compile_example(${example_source} HighFiveOptionalDependencies)
 endforeach()
 
 if(HIGHFIVE_TEST_SPAN)
diff --git a/packages/HighFive/src/examples/select_slices.cpp b/packages/HighFive/src/examples/select_slices.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ea428fa3fefa57399c542b5f62fa8860f4010071
--- /dev/null
+++ b/packages/HighFive/src/examples/select_slices.cpp
@@ -0,0 +1,131 @@
+#include <functional>
+#include <iostream>
+#include <string>
+#include <vector>
+
+#include <highfive/highfive.hpp>
+
+void print_mask(const std::vector<std::vector<double>>& values,
+                const std::vector<std::vector<double>>& selected);
+
+void print_values(const std::vector<std::vector<double>>& values);
+
+void pretty_print(const std::vector<std::vector<double>>& values,
+                  const std::vector<std::vector<double>>& selected);
+
+int main(void) {
+    using namespace HighFive;
+    using container_type = std::vector<std::vector<double>>;
+
+    const std::string file_name("select_slices.h5");
+    const std::string dataset_name("dset");
+
+    // Create a new file using the default property lists.
+    File file(file_name, File::Truncate);
+
+    // we have some example values in a 2D vector 2x5
+    // clang-format off
+    container_type values = {
+      {1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7},
+      {2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7},
+      {3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7},
+      {4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7},
+      {5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7},
+      {6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7},
+      {7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7},
+      {8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7},
+      {9.1, 9.2, 9.3, 9.4, 9.5, 9.6, 9.7}
+    };
+    // clang-format on
+
+    auto dset = file.createDataSet(dataset_name, values);
+
+    // Let's start with the selection `values[1:4, 2:4]`.
+    {
+        auto xslice = std::array<size_t, 2>{2, 4};
+        auto yslice = std::array<size_t, 2>{1, 4};
+
+        auto selected = dset.select(ProductSet(yslice, xslice)).read<container_type>();
+        std::cout << " -- values[1:4, 2:4] ------------ \n";
+        pretty_print(values, selected);
+    }
+
+    // Now we'd like the selection `values[[1,2,8], 2:4]`.
+    {
+        auto xslice = std::array<size_t, 2>{2, 4};
+        auto yslice = std::vector<size_t>{1, 2, 8};
+
+        auto selected = dset.select(ProductSet(yslice, xslice)).read<container_type>();
+        std::cout << "\n -- values[[1,2,8], 2:4] -------- \n";
+        pretty_print(values, selected);
+    }
+
+    // ... or the union of multiple slices.
+    {
+        auto xslice = std::array<size_t, 2>{2, 4};
+        auto yslice = std::vector<std::array<size_t, 2>>{{0, 2}, {5, 9}};
+
+        auto selected = dset.select(ProductSet(yslice, xslice)).read<container_type>();
+        std::cout << "\n -- values[[0:2, 5:10], 2:4] -------- \n";
+        pretty_print(values, selected);
+    }
+
+    // Similar for the union of multiple slices in both x- and y-direction.
+    {
+        auto xslice = std::vector<std::array<size_t, 2>>{{0, 1}, {2, 4}, {6, 7}};
+        auto yslice = std::vector<std::array<size_t, 2>>{{0, 2}, {5, 9}};
+
+        auto selected = dset.select(ProductSet(yslice, xslice)).read<container_type>();
+        std::cout << "\n -- values[[0:2, 5:10], [0:1, 2:4, 6:7]] -------- \n";
+        pretty_print(values, selected);
+    }
+
+    // If only a single row/column is need use the following. However,
+    // selecting elements one-by-one in a loop can be a serious performance
+    // issue.
+    {
+        auto xslice = std::vector<std::array<size_t, 2>>{{0, 1}, {2, 4}, {6, 7}};
+        auto row_id = 3;
+
+        auto selected = dset.select(ProductSet(row_id, xslice)).read<container_type>();
+        std::cout << "\n -- values[3, [0:1, 2:4, 6:7]] -------- \n";
+        pretty_print(values, selected);
+    }
+
+    return 0;
+}
+
+void print_mask(const std::vector<std::vector<double>>& values,
+                const std::vector<std::vector<double>>& selected) {
+    std::vector<double> flat_selection;
+    for (const auto& row: selected) {
+        for (const auto& x: row) {
+            flat_selection.push_back(x);
+        }
+    }
+
+    for (const auto& row: values) {
+        for (const auto& x: row) {
+            auto found = std::find(flat_selection.begin(), flat_selection.end(), x) !=
+                         flat_selection.end();
+            std::cout << (found ? "x" : ".") << "  ";
+        }
+        std::cout << "\n";
+    }
+}
+
+void print_values(const std::vector<std::vector<double>>& values) {
+    for (const auto& row: values) {
+        for (const auto& x: row) {
+            std::cout << x << "  ";
+        }
+        std::cout << "\n";
+    }
+}
+
+void pretty_print(const std::vector<std::vector<double>>& values,
+                  const std::vector<std::vector<double>>& selected) {
+    print_mask(values, selected);
+    std::cout << "\n";
+    print_values(selected);
+}
diff --git a/packages/HighFive/tests/unit/CMakeLists.txt b/packages/HighFive/tests/unit/CMakeLists.txt
index c8835ba34747a69e6eb2fcfc97498a8ea75a2507..10cff1f6fe24a5c8c3e4927bc7ec083fd85fdcc5 100644
--- a/packages/HighFive/tests/unit/CMakeLists.txt
+++ b/packages/HighFive/tests/unit/CMakeLists.txt
@@ -6,7 +6,7 @@ if(MSVC)
 endif()
 
 ## Base tests
-foreach(test_name tests_high_five_base tests_high_five_multi_dims tests_high_five_easy test_all_types test_high_five_selection tests_high_five_data_type test_empty_arrays test_legacy test_opencv test_string test_xtensor)
+foreach(test_name tests_high_five_base tests_high_five_easy test_all_types test_high_five_selection tests_high_five_data_type test_boost test_empty_arrays test_legacy test_opencv test_string test_stl test_xtensor)
   add_executable(${test_name} "${test_name}.cpp")
   target_link_libraries(${test_name} HighFive HighFiveWarnings HighFiveFlags Catch2::Catch2WithMain)
   target_link_libraries(${test_name} HighFiveOptionalDependencies)
diff --git a/packages/HighFive/tests/unit/compary_arrays.hpp b/packages/HighFive/tests/unit/compary_arrays.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8a867a55d8607280928a0264af6e1345655c041a
--- /dev/null
+++ b/packages/HighFive/tests/unit/compary_arrays.hpp
@@ -0,0 +1,80 @@
+/*
+ *  Copyright (c), 2023, 2024 Blue Brain Project - EPFL
+ *
+ *  Distributed under the Boost Software License, Version 1.0.
+ *    (See accompanying file LICENSE_1_0.txt or copy at
+ *          http://www.boost.org/LICENSE_1_0.txt)
+ *
+ */
+
+#include <string>
+#include <sstream>
+#include <type_traits>
+
+#include <catch2/catch_template_test_macros.hpp>
+
+#include "data_generator.hpp"
+
+namespace HighFive {
+namespace testing {
+
+template <class T, class = void>
+struct DiffMessageTrait;
+
+template <class T>
+struct DiffMessageTrait<T, typename std::enable_if<std::is_floating_point<T>::value>::type> {
+    static std::string diff(T a, T b) {
+        std::stringstream sstream;
+        sstream << std::scientific << " delta: " << a - b;
+        return sstream.str();
+    }
+};
+
+template <class T>
+struct DiffMessageTrait<T, typename std::enable_if<!std::is_floating_point<T>::value>::type> {
+    static std::string diff(T /* a */, T /* b */) {
+        return "";
+    }
+};
+
+template <class T>
+std::string diff_message(T a, T b) {
+    return DiffMessageTrait<T>::diff(a, b);
+}
+
+template <class Actual, class Expected, class Comp>
+void compare_arrays(const Actual& actual,
+                    const Expected& expected,
+                    const std::vector<size_t>& dims,
+                    Comp comp) {
+    using actual_trait = testing::ContainerTraits<Actual>;
+    using expected_trait = testing::ContainerTraits<Expected>;
+    using base_type = typename actual_trait::base_type;
+
+    auto n = testing::flat_size(dims);
+
+    for (size_t i = 0; i < n; ++i) {
+        auto indices = testing::unravel(i, dims);
+        base_type actual_value = actual_trait::get(actual, indices);
+        base_type expected_value = expected_trait::get(expected, indices);
+        auto c = comp(actual_value, expected_value);
+        if (!c) {
+            std::stringstream sstream;
+            sstream << std::scientific << "i = " << i << ": " << actual_value
+                    << " != " << expected_value << diff_message(actual_value, expected_value);
+            INFO(sstream.str());
+        }
+        REQUIRE(c);
+    }
+}
+
+template <class Actual, class Expected>
+void compare_arrays(const Actual& actual,
+                    const Expected& expected,
+                    const std::vector<size_t>& dims) {
+    using base_type = typename testing::ContainerTraits<Actual>::base_type;
+    compare_arrays(expected, actual, dims, [](base_type a, base_type b) { return a == b; });
+}
+
+}  // namespace testing
+}  // namespace HighFive
diff --git a/packages/HighFive/tests/unit/data_generator.hpp b/packages/HighFive/tests/unit/data_generator.hpp
index d513c3420edd72034263801e8dd56974342eaad5..e4e7dfa85a28f50bdaa53fdd59c5bc9afa5cc4c9 100644
--- a/packages/HighFive/tests/unit/data_generator.hpp
+++ b/packages/HighFive/tests/unit/data_generator.hpp
@@ -283,6 +283,38 @@ struct ContainerTraits<std::span<T, Extent>>: public STLLikeContainerTraits<std:
 #endif
 
 
+template <class T, size_t n>
+struct ContainerTraits<T[n]> {
+    using container_type = T[n];
+    using value_type = T;
+    using base_type = typename ContainerTraits<value_type>::base_type;
+
+    static constexpr bool is_view = ContainerTraits<value_type>::is_view;
+    static constexpr size_t rank = 1 + ContainerTraits<value_type>::rank;
+
+    static void set(container_type& array,
+                    const std::vector<size_t>& indices,
+                    const base_type& value) {
+        return ContainerTraits<value_type>::set(array[indices[0]], lstrip(indices, 1), value);
+    }
+
+    static base_type get(const container_type& array, const std::vector<size_t>& indices) {
+        return ContainerTraits<value_type>::get(array[indices[0]], lstrip(indices, 1));
+    }
+
+    static void assign(container_type& dst, const container_type& src) {
+        for (size_t i = 0; i < n; ++i) {
+            dst[i] = src[i];
+        }
+    }
+
+    static void sanitize_dims(std::vector<size_t>& dims, size_t axis) {
+        dims[axis] = n;
+        ContainerTraits<value_type>::sanitize_dims(dims, axis + 1);
+    }
+};
+
+
 // -- Boost  -------------------------------------------------------------------
 #ifdef HIGHFIVE_TEST_BOOST
 template <class T, size_t n>
@@ -723,6 +755,37 @@ struct MultiDimVector<T, 0> {
     using type = T;
 };
 
+template <class C, class F>
+void initialize_impl(C& array,
+                     const std::vector<size_t>& dims,
+                     std::vector<size_t>& indices,
+                     size_t axis,
+                     F f) {
+    using traits = ContainerTraits<C>;
+    if (axis == indices.size()) {
+        auto value = f(indices);
+        traits::set(array, indices, value);
+    } else {
+        for (size_t i = 0; i < dims[axis]; ++i) {
+            indices[axis] = i;
+            initialize_impl(array, dims, indices, axis + 1, f);
+        }
+    }
+}
+
+template <class C, class F>
+void initialize(C& array, const std::vector<size_t>& dims, F f) {
+    std::vector<size_t> indices(dims.size());
+    initialize_impl(array, dims, indices, 0, f);
+}
+
+template <class C>
+void initialize(C& array, const std::vector<size_t>& dims) {
+    using traits = ContainerTraits<C>;
+    initialize(array, dims, DefaultValues<typename traits::base_type>());
+}
+
+
 template <class Container>
 class DataGenerator {
   public:
@@ -761,30 +824,6 @@ class DataGenerator {
     static void sanitize_dims(std::vector<size_t>& dims) {
         ContainerTraits<Container>::sanitize_dims(dims, /* axis = */ 0);
     }
-
-  private:
-    template <class C, class F>
-    static void initialize(C& array, const std::vector<size_t>& dims, F f) {
-        std::vector<size_t> indices(dims.size());
-        initialize(array, dims, indices, 0, f);
-    }
-
-    template <class C, class F>
-    static void initialize(C& array,
-                           const std::vector<size_t>& dims,
-                           std::vector<size_t>& indices,
-                           size_t axis,
-                           F f) {
-        if (axis == indices.size()) {
-            auto value = f(indices);
-            traits::set(array, indices, value);
-        } else {
-            for (size_t i = 0; i < dims[axis]; ++i) {
-                indices[axis] = i;
-                initialize(array, dims, indices, axis + 1, f);
-            }
-        }
-    }
 };
 
 }  // namespace testing
diff --git a/packages/HighFive/tests/unit/test_all_types.cpp b/packages/HighFive/tests/unit/test_all_types.cpp
index cddc73312fa3a3b7c9f6b7622730749303357d25..0a18e0257a8113195436d056dcdc7c0fa583dc5b 100644
--- a/packages/HighFive/tests/unit/test_all_types.cpp
+++ b/packages/HighFive/tests/unit/test_all_types.cpp
@@ -7,7 +7,6 @@
  *
  */
 #include <string>
-#include <sstream>
 
 #include <catch2/catch_template_test_macros.hpp>
 
@@ -19,6 +18,7 @@
 #include "data_generator.hpp"
 #include "create_traits.hpp"
 #include "supported_types.hpp"
+#include "compary_arrays.hpp"
 
 using namespace HighFive;
 
@@ -52,170 +52,6 @@ TEMPLATE_TEST_CASE("Scalar in DataSet", "[Types]", bool, std::string) {
     }
 }
 
-TEMPLATE_PRODUCT_TEST_CASE("Scalar in std::vector", "[Types]", std::vector, (bool, std::string)) {
-    const std::string file_name("rw_dataset_" + typeNameHelper<TestType>() + ".h5");
-    const std::string dataset_name("dset");
-    TestType t1(5);
-
-    {
-        // Create a new file using the default property lists.
-        File file(file_name, File::ReadWrite | File::Create | File::Truncate);
-
-        // Create the dataset
-        DataSet dataset = file.createDataSet(
-            dataset_name, {5}, create_datatype<typename details::inspector<TestType>::base_type>());
-
-        // Write into the initial part of the dataset
-        dataset.write(t1);
-    }
-
-    // read it back
-    {
-        File file(file_name, File::ReadOnly);
-
-        TestType value;
-        DataSet dataset = file.getDataSet("/" + dataset_name);
-        dataset.read(value);
-        CHECK(t1 == value);
-        CHECK(value.size() == 5);
-    }
-}
-
-TEMPLATE_PRODUCT_TEST_CASE("Scalar in std::vector<std::vector>",
-                           "[Types]",
-                           std::vector,
-                           (bool, std::string)) {
-    const std::string file_name("rw_dataset_vector_" + typeNameHelper<TestType>() + ".h5");
-    const std::string dataset_name("dset");
-    std::vector<TestType> t1(5);
-    for (auto&& e: t1) {
-        e.resize(6);
-    }
-
-    {
-        // Create a new file using the default property lists.
-        File file(file_name, File::ReadWrite | File::Create | File::Truncate);
-
-        // Create the dataset
-        DataSet dataset = file.createDataSet(
-            dataset_name,
-            {5, 6},
-            create_datatype<typename details::inspector<std::vector<TestType>>::base_type>());
-
-        // Write into the initial part of the dataset
-        dataset.write(t1);
-    }
-
-    // read it back
-    {
-        File file(file_name, File::ReadOnly);
-
-        std::vector<TestType> value;
-        DataSet dataset = file.getDataSet("/" + dataset_name);
-        dataset.read(value);
-        CHECK(t1 == value);
-        CHECK(value.size() == 5);
-    }
-}
-
-TEMPLATE_TEST_CASE("Scalar in std::array", "[Types]", bool, std::string) {
-    const std::string file_name("rw_dataset_array_" + typeNameHelper<TestType>() + ".h5");
-    const std::string dataset_name("dset");
-    std::array<TestType, 5> t1{};
-
-    {
-        // Create a new file using the default property lists.
-        File file(file_name, File::ReadWrite | File::Create | File::Truncate);
-
-        // Create the dataset
-        DataSet dataset = file.createDataSet(
-            dataset_name,
-            {5},
-            create_datatype<typename details::inspector<std::array<TestType, 5>>::base_type>());
-
-        // Write into the initial part of the dataset
-        dataset.write(t1);
-    }
-
-    // read it back
-    {
-        File file(file_name, File::ReadOnly);
-
-        std::array<TestType, 5> value;
-        DataSet dataset = file.getDataSet("/" + dataset_name);
-        dataset.read(value);
-        CHECK(t1 == value);
-        CHECK(value.size() == 5);
-    }
-}
-
-TEMPLATE_TEST_CASE("Scalar in std::vector<std::array>", "[Types]", bool, std::string) {
-    const std::string file_name("rw_dataset_vector_array_" + typeNameHelper<TestType>() + ".h5");
-    const std::string dataset_name("dset");
-    std::vector<std::array<TestType, 6>> t1(5);
-
-    {
-        // Create a new file using the default property lists.
-        File file(file_name, File::ReadWrite | File::Create | File::Truncate);
-
-        // Create the dataset
-        DataSet dataset = file.createDataSet(
-            dataset_name,
-            {5, 6},
-            create_datatype<
-                typename details::inspector<std::vector<std::array<TestType, 5>>>::base_type>());
-
-        // Write into the initial part of the dataset
-        dataset.write(t1);
-    }
-
-    // read it back
-    {
-        File file(file_name, File::ReadOnly);
-
-        std::vector<std::array<TestType, 6>> value;
-        DataSet dataset = file.getDataSet("/" + dataset_name);
-        dataset.read(value);
-        CHECK(t1 == value);
-        CHECK(value.size() == 5);
-    }
-}
-
-TEMPLATE_TEST_CASE("Scalar in std::array<std::vector>", "[Types]", bool, std::string) {
-    const std::string file_name("rw_dataset_array_vector" + typeNameHelper<TestType>() + ".h5");
-    const std::string dataset_name("dset");
-    std::array<std::vector<TestType>, 6> t1;
-    for (auto& tt: t1) {
-        tt = std::vector<TestType>(5);
-    }
-
-    {
-        // Create a new file using the default property lists.
-        File file(file_name, File::ReadWrite | File::Create | File::Truncate);
-
-        // Create the dataset
-        DataSet dataset = file.createDataSet(
-            dataset_name,
-            {6, 5},
-            create_datatype<
-                typename details::inspector<std::vector<std::array<TestType, 5>>>::base_type>());
-
-        // Write into the initial part of the dataset
-        dataset.write(t1);
-    }
-
-    // read it back
-    {
-        File file(file_name, File::ReadOnly);
-
-        std::array<std::vector<TestType>, 6> value;
-        DataSet dataset = file.getDataSet("/" + dataset_name);
-        dataset.read(value);
-        CHECK(t1 == value);
-        CHECK(value.size() == 6);
-    }
-}
-
 #if HIGHFIVE_CXX_STD >= 17
 TEMPLATE_PRODUCT_TEST_CASE("Scalar in std::vector<std::byte>", "[Types]", std::vector, std::byte) {
     const std::string file_name("rw_dataset_vector_" + typeNameHelper<TestType>() + ".h5");
@@ -246,68 +82,11 @@ TEMPLATE_PRODUCT_TEST_CASE("Scalar in std::vector<std::byte>", "[Types]", std::v
 }
 #endif
 
-template <class T, class = void>
-struct DiffMessageTrait;
-
-template <class T>
-struct DiffMessageTrait<T, typename std::enable_if<std::is_floating_point<T>::value>::type> {
-    static std::string diff(T a, T b) {
-        std::stringstream sstream;
-        sstream << std::scientific << " delta: " << a - b;
-        return sstream.str();
-    }
-};
-
-template <class T>
-struct DiffMessageTrait<T, typename std::enable_if<!std::is_floating_point<T>::value>::type> {
-    static std::string diff(T /* a */, T /* b */) {
-        return "";
-    }
-};
-
-template <class T>
-std::string diff_message(T a, T b) {
-    return DiffMessageTrait<T>::diff(a, b);
-}
-
-template <class Actual, class Expected, class Comp>
-void compare_arrays(const Actual& actual,
-                    const Expected& expected,
-                    const std::vector<size_t>& dims,
-                    Comp comp) {
-    using actual_trait = testing::ContainerTraits<Actual>;
-    using expected_trait = testing::ContainerTraits<Expected>;
-    using base_type = typename actual_trait::base_type;
-
-    auto n = testing::flat_size(dims);
-
-    for (size_t i = 0; i < n; ++i) {
-        auto indices = testing::unravel(i, dims);
-        base_type actual_value = actual_trait::get(actual, indices);
-        base_type expected_value = expected_trait::get(expected, indices);
-        auto c = comp(actual_value, expected_value);
-        if (!c) {
-            std::stringstream sstream;
-            sstream << std::scientific << "i = " << i << ": " << actual_value
-                    << " != " << expected_value << diff_message(actual_value, expected_value);
-            INFO(sstream.str());
-        }
-        REQUIRE(c);
-    }
-}
-
-template <class Actual, class Expected>
-void compare_arrays(const Actual& actual,
-                    const Expected& expected,
-                    const std::vector<size_t>& dims) {
-    using base_type = typename testing::ContainerTraits<Actual>::base_type;
-    compare_arrays(expected, actual, dims, [](base_type a, base_type b) { return a == b; });
-}
 
 template <class Container, class Expected, class Obj>
 auto check_read_auto(const Expected& expected, const std::vector<size_t>& dims, const Obj& obj) ->
     typename std::enable_if<!testing::ContainerTraits<Container>::is_view>::type {
-    compare_arrays(obj.template read<Container>(), expected, dims);
+    testing::compare_arrays(obj.template read<Container>(), expected, dims);
 }
 
 template <class Container, class Expected, class Obj>
@@ -321,7 +100,7 @@ void check_read_preallocated(const Expected& expected,
     auto actual = testing::DataGenerator<Container>::allocate(dims);
     obj.read(actual);
 
-    compare_arrays(actual, expected, dims);
+    testing::compare_arrays(actual, expected, dims);
 
     testing::ContainerTraits<Container>::deallocate(actual, dims);
 }
@@ -388,7 +167,7 @@ void check_writing(const std::vector<size_t>& dims, Write write) {
     auto actual = testing::DataGenerator<reference_type>::allocate(dims);
     obj.read(actual);
 
-    compare_arrays(actual, expected, dims);
+    testing::compare_arrays(actual, expected, dims);
 
     testing::ContainerTraits<reference_type>::deallocate(actual, dims);
     testing::ContainerTraits<Container>::deallocate(values, dims);
diff --git a/packages/HighFive/tests/unit/test_boost.cpp b/packages/HighFive/tests/unit/test_boost.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..72bca9e7748c1d8c4e994fb7680e18b632aa0efc
--- /dev/null
+++ b/packages/HighFive/tests/unit/test_boost.cpp
@@ -0,0 +1,27 @@
+/*
+ *  Copyright (c), 2024, Blue Brain Project - EPFL
+ *
+ *  Distributed under the Boost Software License, Version 1.0.
+ *    (See accompanying file LICENSE_1_0.txt or copy at
+ *          http://www.boost.org/LICENSE_1_0.txt)
+ *
+ */
+#if HIGHFIVE_TEST_BOOST
+#include <string>
+
+#include <catch2/catch_template_test_macros.hpp>
+
+#include <highfive/highfive.hpp>
+#include <highfive/boost.hpp>
+
+using namespace HighFive;
+
+TEST_CASE("Test boost::multi_array with fortran_storage_order") {
+    const std::string file_name("h5_multi_array_fortran.h5");
+    File file(file_name, File::ReadWrite | File::Create | File::Truncate);
+
+    boost::multi_array<int, 2> ma(boost::extents[2][2], boost::fortran_storage_order());
+    auto dset = file.createDataSet<int>("main_dset", DataSpace::From(ma));
+    CHECK_THROWS_AS(dset.write(ma), DataTypeException);
+}
+#endif
diff --git a/packages/HighFive/tests/unit/test_stl.cpp b/packages/HighFive/tests/unit/test_stl.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..42a76c4b9a71e7ff2cb01f03d22f4ccd02d67a0b
--- /dev/null
+++ b/packages/HighFive/tests/unit/test_stl.cpp
@@ -0,0 +1,58 @@
+/*
+ *  Copyright (c), 2017-2024, Blue Brain Project - EPFL
+ *
+ *  Distributed under the Boost Software License, Version 1.0.
+ *    (See accompanying file LICENSE_1_0.txt or copy at
+ *          http://www.boost.org/LICENSE_1_0.txt)
+ *
+ */
+
+#include <catch2/catch_template_test_macros.hpp>
+
+#include <highfive/highfive.hpp>
+
+
+#include "compary_arrays.hpp"
+#include "create_traits.hpp"
+#include "data_generator.hpp"
+
+using namespace HighFive;
+
+TEST_CASE("std::array undersized", "[stl]") {
+    auto file = File("rw_std_array_undersized.h5", File::Truncate);
+    auto x = std::array<double, 3>{1.0, 2.0, 3.0};
+    auto dset = file.createDataSet("x", x);
+
+    REQUIRE_THROWS(dset.read<std::array<double, 2>>());
+
+    auto xx = std::array<double, 2>();
+    REQUIRE_THROWS(dset.read(xx));
+}
+
+TEST_CASE("T[n][m]") {
+    using reference_container = std::vector<std::vector<double>>;
+    auto file = File("rw_carray.h5", File::Truncate);
+
+    constexpr size_t n = 3;
+    constexpr size_t m = 5;
+
+    double x[n][m];
+
+    SECTION("write") {
+        testing::initialize(x, {n, m});
+
+        auto dset = file.createDataSet("x", x);
+        auto actual = dset.read<reference_container>();
+
+        testing::compare_arrays(x, actual, {n, m});
+    }
+
+    SECTION("read") {
+        auto expected = testing::DataGenerator<reference_container>::create({n, m});
+
+        auto dset = file.createDataSet("x", expected);
+        dset.read(x);
+
+        testing::compare_arrays(expected, x, {n, m});
+    }
+}
diff --git a/packages/HighFive/tests/unit/tests_high_five.hpp b/packages/HighFive/tests/unit/tests_high_five.hpp
index d9a4ed34de07230317584e847131957226d43167..8a4f07bbf9da6eb3c5284860da713a13b1dc3bf5 100644
--- a/packages/HighFive/tests/unit/tests_high_five.hpp
+++ b/packages/HighFive/tests/unit/tests_high_five.hpp
@@ -152,7 +152,7 @@ struct ContentGenerate<std::string> {
         ContentGenerate<char> gen;
         std::string random_string;
         std::mt19937_64 rgen;
-        rgen.seed(88);
+        rgen.seed(42);
         std::uniform_int_distribution<unsigned> int_dist(0, 1000);
         const size_t size_string = int_dist(rgen);
 
diff --git a/packages/HighFive/tests/unit/tests_high_five_base.cpp b/packages/HighFive/tests/unit/tests_high_five_base.cpp
index f52c2dcb893f1aa738101f0116c3b579d72cf335..02ce7a04de746187e096a274b0a6a78f930af450 100644
--- a/packages/HighFive/tests/unit/tests_high_five_base.cpp
+++ b/packages/HighFive/tests/unit/tests_high_five_base.cpp
@@ -1224,6 +1224,207 @@ TEST_CASE("datasetOffset") {
 }
 
 
+// Ensure that "all" combinations of `ProductSet` are being instantiated.
+class CompileProductSet {
+    using Point = size_t;
+    using Points = std::vector<size_t>;
+    using Slice = std::array<size_t, 2>;
+    using Slices = std::vector<Slice>;
+
+    void all() {
+        one_arg();
+        two_args();
+        three_args();
+    }
+
+    template <class... Preamble>
+    void zero_args() {
+        ProductSet(Preamble()...);
+    }
+
+    template <class... Preamble>
+    void one_arg() {
+        zero_args<Preamble..., Point>();
+        zero_args<Preamble..., Points>();
+        zero_args<Preamble..., Slice>();
+        zero_args<Preamble..., Slices>();
+    }
+
+    template <class... Preamble>
+    void two_args() {
+        one_arg<Preamble..., Point>();
+        one_arg<Preamble..., Points>();
+        one_arg<Preamble..., Slice>();
+        one_arg<Preamble..., Slices>();
+    }
+
+    template <class... Preamble>
+    void three_args() {
+        two_args<Preamble..., Point>();
+        two_args<Preamble..., Points>();
+        two_args<Preamble..., Slice>();
+        two_args<Preamble..., Slices>();
+    }
+};
+
+template <class S, class Y, class X>
+void check_product_set_shape(const S& subarray, const Y& yslices, const X& xslices) {
+    std::vector<size_t> subshape{0, 0};
+
+    for (auto yslice: yslices) {
+        subshape[0] += yslice[1] - yslice[0];
+    }
+
+    for (auto xslice: xslices) {
+        subshape[1] += xslice[1] - xslice[0];
+    }
+
+    REQUIRE(subarray.size() == subshape[0]);
+    for (const auto& v: subarray) {
+        REQUIRE(v.size() == subshape[1]);
+    }
+}
+
+template <class S, class Y, class X>
+void check_product_set_values(const S& array,
+                              const S& subarray,
+                              const Y& yslices,
+                              const X& xslices) {
+    size_t il = 0;
+    for (const auto& yslice: yslices) {
+        for (size_t ig = yslice[0]; ig < yslice[1]; ++ig) {
+            size_t jl = 0;
+
+            for (const auto& xslice: xslices) {
+                for (size_t jg = xslice[0]; jg < xslice[1]; ++jg) {
+                    REQUIRE(subarray[il][jl] == array[ig][jg]);
+                    ++jl;
+                }
+            }
+            ++il;
+        }
+    }
+}
+
+template <class S, class Y, class X>
+void check(const S& array, const S& subarray, const Y& yslices, const X& xslices) {
+    check_product_set_shape(subarray, yslices, xslices);
+    check_product_set_values(array, subarray, yslices, xslices);
+}
+
+
+TEST_CASE("productSet") {
+    using Slice = std::array<size_t, 2>;
+    using Slices = std::vector<Slice>;
+    using Point = size_t;
+    using Points = std::vector<size_t>;
+
+    const std::string file_name("h5_test_product_set.h5");
+
+    auto generate = [](size_t n, size_t m, auto f) {
+        auto x = std::vector<std::vector<double>>(n);
+        for (size_t i = 0; i < n; ++i) {
+            x[i] = std::vector<double>(m);
+        }
+
+        for (size_t i = 0; i < n; ++i) {
+            for (size_t j = 0; j < m; ++j) {
+                x[i][j] = f(i, j);
+            }
+        }
+
+        return x;
+    };
+
+    auto array = generate(6, 12, [](size_t i, size_t j) { return double(i) + double(j) * 0.01; });
+
+    auto file = File(file_name, File::Truncate);
+    auto dset = file.createDataSet("dset", array);
+
+    SECTION("rR") {
+        std::vector<std::vector<double>> subarray;
+
+        auto yslice = Slice{1, 3};
+        auto yslices = Slices{yslice};
+        auto xslices = Slices{{0, 1}, {3, 5}};
+
+        dset.select(ProductSet(yslice, xslices)).read(subarray);
+
+        check(array, subarray, yslices, xslices);
+    }
+
+    SECTION("Rr") {
+        std::vector<std::vector<double>> subarray;
+
+        auto yslices = Slices{{0, 1}, {3, 5}};
+        auto xslice = Slice{1, 3};
+        auto xslices = Slices{xslice};
+
+        dset.select(ProductSet(yslices, xslice)).read(subarray);
+
+        check(array, subarray, yslices, xslices);
+    }
+
+    SECTION("RP") {
+        std::vector<std::vector<double>> subarray;
+
+        auto yslices = Slices{{0, 1}, {3, 5}};
+        auto xpoints = Points{2, 4, 5};
+        auto xslices = Slices{{2, 3}, {4, 6}};
+
+        dset.select(ProductSet(yslices, xpoints)).read(subarray);
+
+        check(array, subarray, yslices, xslices);
+    }
+
+    SECTION("pR") {
+        std::vector<std::vector<double>> subarray;
+
+        auto ypoint = Point{2};
+        auto yslices = Slices{{2, 3}};
+        auto xslices = Slices{{0, 1}, {3, 5}};
+
+        dset.select(ProductSet(ypoint, xslices)).read(subarray);
+
+        check(array, subarray, yslices, xslices);
+    }
+
+    SECTION("pp") {
+        std::vector<std::vector<double>> subarray;
+
+        auto xpoint = Point{3};
+        auto ypoint = Point{2};
+        auto yslices = Slices{{2, 3}};
+        auto xslices = Slices{{3, 4}};
+
+        dset.select(ProductSet(ypoint, xpoint)).read(subarray);
+        check(array, subarray, yslices, xslices);
+    }
+
+    SECTION("PP") {
+        std::vector<std::vector<double>> subarray;
+
+        auto xpoints = Points{0, 3, 4};
+        auto ypoints = Points{2, 3};
+        auto yslices = Slices{{2, 4}};
+        auto xslices = Slices{{0, 1}, {3, 5}};
+
+        dset.select(ProductSet(ypoints, xpoints)).read(subarray);
+        check(array, subarray, yslices, xslices);
+    }
+
+    SECTION("RR") {
+        std::vector<std::vector<double>> subarray;
+
+        auto yslices = Slices{{2, 4}};
+        auto xslices = Slices{{0, 1}, {3, 5}};
+
+        dset.select(ProductSet(yslices, xslices)).read(subarray);
+        check(array, subarray, yslices, xslices);
+    }
+}
+
+
 template <typename T>
 void attribute_scalar_rw() {
     std::ostringstream filename;
diff --git a/packages/HighFive/tests/unit/tests_high_five_easy.cpp b/packages/HighFive/tests/unit/tests_high_five_easy.cpp
index d10ef941b6a68fd010380abdd29372855c79c455..8bcce26d095100a8754738db0afe63dff8cebf5a 100644
--- a/packages/HighFive/tests/unit/tests_high_five_easy.cpp
+++ b/packages/HighFive/tests/unit/tests_high_five_easy.cpp
@@ -250,15 +250,17 @@ TEST_CASE("H5Easy_xtensor_column_major") {
 
     xt::xtensor<double, 2> A = 100. * xt::random::randn<double>({20, 5});
 
-    H5Easy::dump(file, "/path/to/A", A);
-
     SECTION("Write column major") {
         column_major_t B = A;
-        REQUIRE_THROWS(H5Easy::dump(file, "path/to/B", B));
+        H5Easy::dump(file, "/path/to/A", B);
+        auto A_r = H5Easy::load<xt::xtensor<double, 2>>(file, "/path/to/A");
+        CHECK(xt::allclose(A, A_r));
     }
 
     SECTION("Read column major") {
-        REQUIRE_THROWS(H5Easy::load<column_major_t>(file, "/path/to/A"));
+        H5Easy::dump(file, "/path/to/A", A);
+        auto A_r = H5Easy::load<column_major_t>(file, "/path/to/A");
+        CHECK(xt::allclose(A, A_r));
     }
 }
 
@@ -269,15 +271,17 @@ TEST_CASE("H5Easy_xarray_column_major") {
 
     xt::xarray<double> A = 100. * xt::random::randn<double>({20, 5});
 
-    H5Easy::dump(file, "/path/to/A", A);
-
     SECTION("Write column major") {
         column_major_t B = A;
-        REQUIRE_THROWS(H5Easy::dump(file, "path/to/B", B));
+        H5Easy::dump(file, "/path/to/A", B);
+        auto A_r = H5Easy::load<xt::xtensor<double, 2>>(file, "/path/to/A");
+        CHECK(xt::allclose(A, A_r));
     }
 
     SECTION("Read column major") {
-        REQUIRE_THROWS(H5Easy::load<column_major_t>(file, "/path/to/A"));
+        H5Easy::dump(file, "/path/to/A", A);
+        auto A_r = H5Easy::load<column_major_t>(file, "/path/to/A");
+        CHECK(xt::allclose(A, A_r));
     }
 }
 
diff --git a/packages/HighFive/tests/unit/tests_high_five_multi_dims.cpp b/packages/HighFive/tests/unit/tests_high_five_multi_dims.cpp
deleted file mode 100644
index a261360e0394fd6e8d8127e50c3e361031393e6f..0000000000000000000000000000000000000000
--- a/packages/HighFive/tests/unit/tests_high_five_multi_dims.cpp
+++ /dev/null
@@ -1,223 +0,0 @@
-/*
- *  Copyright (c), 2017-2019, Blue Brain Project - EPFL
- *
- *  Distributed under the Boost Software License, Version 1.0.
- *    (See accompanying file LICENSE_1_0.txt or copy at
- *          http://www.boost.org/LICENSE_1_0.txt)
- *
- */
-
-#include <string>
-#include <iostream>
-
-#include <highfive/highfive.hpp>
-
-
-#ifdef HIGHFIVE_TEST_BOOST
-#include <boost/multi_array.hpp>
-#include <highfive/boost.hpp>
-#endif
-
-#include <catch2/catch_test_macros.hpp>
-#include <catch2/catch_template_test_macros.hpp>
-
-#include "tests_high_five.hpp"
-
-using namespace HighFive;
-
-/// \brief Test for 2D old-style arrays (T array[x][y])
-template <typename T>
-void readWrite2DArrayTest() {
-    std::ostringstream filename;
-    filename << "h5_rw_2d_array_" << typeNameHelper<T>() << "_test.h5";
-    const std::string DATASET_NAME("dset");
-    const size_t x_size = 100;
-    const size_t y_size = 10;
-
-    // Create a new file using the default property lists.
-    File file(filename.str(), File::ReadWrite | File::Create | File::Truncate);
-
-    // Create the data space for the dataset.
-    std::vector<size_t> dims{x_size, y_size};
-
-    DataSpace dataspace(dims);
-
-    // Create a dataset with arbitrary type
-    DataSet dataset = file.createDataSet<T>(DATASET_NAME, dataspace);
-
-    T array[x_size][y_size];
-
-    ContentGenerate<T> generator;
-    generate2D(array, x_size, y_size, generator);
-
-    dataset.write(array);
-
-    T result[x_size][y_size];
-    dataset.read(result);
-
-    for (size_t i = 0; i < x_size; ++i) {
-        for (size_t j = 0; j < y_size; ++j) {
-            CHECK(result[i][j] == array[i][j]);
-        }
-    }
-}
-
-TEMPLATE_LIST_TEST_CASE("ReadWrite2DArray", "[template]", numerical_test_types) {
-    readWrite2DArrayTest<TestType>();
-}
-
-template <typename T>
-void readWriteArrayTest() {
-    const size_t x_size = 200;
-    typename std::array<T, x_size> vec;
-    ContentGenerate<T> generator;
-    std::generate(vec.begin(), vec.end(), generator);
-
-    typename std::array<T, x_size> result;
-    auto dataset = readWriteDataset<T>(vec, result, 1, "std-array");
-
-    CHECK(result == vec);
-
-    typename std::array<T, 1> tooSmall;
-    CHECK_THROWS_AS(dataset.read(tooSmall), DataSpaceException);
-}
-
-TEMPLATE_LIST_TEST_CASE("readWriteArray", "[template]", numerical_test_types) {
-    readWriteArrayTest<TestType>();
-}
-
-template <typename T, typename VectorSubT>
-void readWriteVectorNDTest(std::vector<VectorSubT>& ndvec, const std::vector<size_t>& dims) {
-    fillVec(ndvec, dims, ContentGenerate<T>());
-
-    std::vector<VectorSubT> result;
-    readWriteDataset<T>(ndvec, result, dims.size(), "vector");
-
-    CHECK(checkLength(result, dims));
-    CHECK(ndvec == result);
-}
-
-TEMPLATE_LIST_TEST_CASE("readWritSimpleVector", "[template]", numerical_test_types) {
-    std::vector<TestType> vec;
-    readWriteVectorNDTest<TestType>(vec, {50});
-}
-
-TEMPLATE_LIST_TEST_CASE("readWrite2DVector", "[template]", numerical_test_types) {
-    std::vector<std::vector<TestType>> _2dvec;
-    readWriteVectorNDTest<TestType>(_2dvec, {10, 8});
-}
-
-TEMPLATE_LIST_TEST_CASE("readWrite3DVector", "[template]", numerical_test_types) {
-    std::vector<std::vector<std::vector<TestType>>> _3dvec;
-    readWriteVectorNDTest<TestType>(_3dvec, {10, 5, 4});
-}
-
-TEMPLATE_LIST_TEST_CASE("readWrite4DVector", "[template]", numerical_test_types) {
-    std::vector<std::vector<std::vector<std::vector<TestType>>>> _4dvec;
-    readWriteVectorNDTest<TestType>(_4dvec, {5, 4, 3, 2});
-}
-
-TEMPLATE_LIST_TEST_CASE("vector of array", "[template]", numerical_test_types) {
-    std::vector<std::array<TestType, 4>> vec{{1, 2, 3, 4}, {1, 2, 3, 4}};
-    std::vector<std::array<TestType, 4>> result;
-    readWriteDataset<TestType>(vec, result, 2, "vector");
-
-    CHECK(vec.size() == result.size());
-    CHECK(vec[0].size() == result[0].size());
-    CHECK(vec == result);
-}
-
-
-#ifdef HIGHFIVE_TEST_BOOST
-
-template <typename T>
-void MultiArray3DTest() {
-    typedef typename boost::multi_array<T, 3> MultiArray;
-
-    std::ostringstream filename;
-    filename << "h5_rw_multiarray_" << typeNameHelper<T>() << "_test.h5";
-
-    const int size_x = 10, size_y = 10, size_z = 10;
-    const std::string DATASET_NAME("dset");
-    MultiArray array(boost::extents[size_x][size_y][size_z]);
-
-    ContentGenerate<T> generator;
-    std::generate(array.data(), array.data() + array.num_elements(), generator);
-
-    // Create a new file using the default property lists.
-    File file(filename.str(), File::ReadWrite | File::Create | File::Truncate);
-
-    DataSet dataset = file.createDataSet<T>(DATASET_NAME, DataSpace::From(array));
-
-    dataset.write(array);
-
-    // read it back
-    MultiArray result;
-
-    dataset.read(result);
-
-    for (long i = 0; i < size_x; ++i) {
-        for (long j = 0; j < size_y; ++j) {
-            for (long k = 0; k < size_z; ++k) {
-                CHECK(array[i][j][k] == result[i][j][k]);
-            }
-        }
-    }
-}
-
-TEMPLATE_LIST_TEST_CASE("MultiArray3D", "[template]", numerical_test_types) {
-    MultiArray3DTest<TestType>();
-}
-
-TEST_CASE("Test boost::multi_array with fortran_storage_order") {
-    const std::string file_name("h5_multi_array_fortran.h5");
-    File file(file_name, File::ReadWrite | File::Create | File::Truncate);
-
-    boost::multi_array<int, 2> ma(boost::extents[2][2], boost::fortran_storage_order());
-    auto dset = file.createDataSet<int>("main_dset", DataSpace::From(ma));
-    CHECK_THROWS_AS(dset.write(ma), DataTypeException);
-}
-
-template <typename T>
-void ublas_matrix_Test() {
-    using Matrix = boost::numeric::ublas::matrix<T>;
-
-    std::ostringstream filename;
-    filename << "h5_rw_multiarray_" << typeNameHelper<T>() << "_test.h5";
-
-    const size_t size_x = 10, size_y = 10;
-    const std::string DATASET_NAME("dset");
-
-    Matrix mat(size_x, size_y);
-
-    ContentGenerate<T> generator;
-    for (std::size_t i = 0; i < mat.size1(); ++i) {
-        for (std::size_t j = 0; j < mat.size2(); ++j) {
-            mat(i, j) = generator();
-        }
-    }
-
-    // Create a new file using the default property lists.
-    File file(filename.str(), File::ReadWrite | File::Create | File::Truncate);
-
-    DataSet dataset = file.createDataSet<T>(DATASET_NAME, DataSpace::From(mat));
-
-    dataset.write(mat);
-
-    // read it back
-    Matrix result;
-
-    dataset.read(result);
-
-    for (size_t i = 0; i < size_x; ++i) {
-        for (size_t j = 0; j < size_y; ++j) {
-            CHECK(mat(i, j) == result(i, j));
-        }
-    }
-}
-
-TEMPLATE_LIST_TEST_CASE("ublas_matrix", "[template]", numerical_test_types) {
-    ublas_matrix_Test<TestType>();
-}
-
-#endif