hdf5.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2019 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
 
#ifndef dealii_hdf5_h
#define dealii_hdf5_h
 
#include <deal.II/base/config.h>
 
#ifdef DEAL_II_WITH_HDF5
 
#  include <deal.II/base/array_view.h>
 
#  include <deal.II/lac/full_matrix.h>
 
#  include <hdf5.h>
 
#  include <numeric>
 
DEAL_II_NAMESPACE_OPEN
 
// It is necessary to turn clang-format off in order to maintain the Doxygen
// links because they are longer than 80 characters
// clang-format off
// clang-format on
namespace HDF5
{
  class HDF5Object
  {
  protected:
    HDF5Object(const std::string &name, const bool mpi);
 
  public:
    template <typename T>
    T
    get_attribute(const std::string &attr_name) const;
 
    template <typename T>
    void
    set_attribute(const std::string &attr_name, const T value);
 
    std::string
    get_name() const;
 
  protected:
    const std::string name;
 
    std::shared_ptr<hid_t> hdf5_reference;
 
    const bool mpi;
  };
 
  class DataSet : public HDF5Object
  {
    friend class Group;
 
  protected:
    DataSet(const std::string &name, const hid_t &parent_group_id, bool mpi);
 
    DataSet(const std::string            &name,
            const hid_t                  &parent_group_id,
            const std::vector<hsize_t>   &dimensions,
            const std::shared_ptr<hid_t> &t_type,
            const bool                    mpi);
 
  public:
    template <typename Container>
    Container
    read();
 
    template <typename Container>
    Container
    read_selection(const std::vector<hsize_t> &coordinates);
 
    // clang-format off
    // clang-format on
    template <typename Container>
    Container
    read_hyperslab(const std::vector<hsize_t> &offset,
                   const std::vector<hsize_t> &count);
 
    template <typename Container>
    Container
    read_hyperslab(const std::vector<hsize_t> &data_dimensions,
                   const std::vector<hsize_t> &offset,
                   const std::vector<hsize_t> &stride,
                   const std::vector<hsize_t> &count,
                   const std::vector<hsize_t> &block);
 
    template <typename number>
    void
    read_none();
 
    template <typename Container>
    void
    write(const Container &data);
 
    template <typename Container>
    void
    write_selection(const Container            &data,
                    const std::vector<hsize_t> &coordinates);
 
    // clang-format off
    // clang-format on
    template <typename Container>
    void
    write_hyperslab(const Container            &data,
                    const std::vector<hsize_t> &offset,
                    const std::vector<hsize_t> &count);
 
    template <typename Container>
    void
    write_hyperslab(const Container            &data,
                    const std::vector<hsize_t> &data_dimensions,
                    const std::vector<hsize_t> &offset,
                    const std::vector<hsize_t> &stride,
                    const std::vector<hsize_t> &count,
                    const std::vector<hsize_t> &block);
 
    template <typename number>
    void
    write_none();
 
    bool
    get_query_io_mode() const;
 
    void
    set_query_io_mode(const bool new_query_io_mode);
 
    std::string
    get_io_mode();
 
    H5D_mpio_actual_io_mode_t
    get_io_mode_as_hdf5_type();
 
    std::string
    get_local_no_collective_cause();
 
    std::uint32_t
    get_local_no_collective_cause_as_hdf5_type();
 
    std::string
    get_global_no_collective_cause();
 
    std::uint32_t
    get_global_no_collective_cause_as_hdf5_type();
 
    std::vector<hsize_t>
    get_dimensions() const;
 
    unsigned int
    get_size() const;
 
    unsigned int
    get_rank() const;
 
  private:
    unsigned int rank;
 
    std::vector<hsize_t> dimensions;
 
    std::shared_ptr<hid_t> dataspace;
 
    unsigned int size;
 
    bool query_io_mode;
 
    H5D_mpio_actual_io_mode_t io_mode;
 
    std::uint32_t local_no_collective_cause;
 
    std::uint32_t global_no_collective_cause;
  };
 
  class Group : public HDF5Object
  {
  protected:
    enum class GroupAccessMode
    {
      open,
      create
    };
    Group(const std::string    &name,
          const Group          &parent_group,
          const bool            mpi,
          const GroupAccessMode mode);
 
    Group(const std::string &name, const bool mpi);
 
  public:
    Group
    open_group(const std::string &name) const;
 
    Group
    create_group(const std::string &name) const;
 
    DataSet
    open_dataset(const std::string &name) const;
 
    template <typename number>
    DataSet
    create_dataset(const std::string          &name,
                   const std::vector<hsize_t> &dimensions) const;
 
    template <typename Container>
    void
    write_dataset(const std::string &name, const Container &data) const;
  };
 
  class File : public Group
  {
  public:
    enum class FileAccessMode
    {
      open,
      create
    };
 
    File(const std::string &name, const FileAccessMode mode);
 
    File(const std::string   &name,
         const FileAccessMode mode,
         const MPI_Comm       mpi_communicator);
 
  private:
    File(const std::string   &name,
         const FileAccessMode mode,
         const bool           mpi,
         const MPI_Comm       mpi_communicator);
  };
 
  namespace internal
  {
    template <typename number>
    std::shared_ptr<hid_t>
    get_hdf5_datatype();
 
    template <typename number>
    std::vector<hsize_t>
    get_container_dimensions(const std::vector<number> &data);
 
    template <typename number>
    std::vector<hsize_t>
    get_container_dimensions(const Vector<number> &data);
 
    template <typename number>
    std::vector<hsize_t>
    get_container_dimensions(const FullMatrix<number> &data);
 
    template <typename number>
    unsigned int
    get_container_size(const std::vector<number> &data);
 
    template <typename number>
    unsigned int
    get_container_size(const Vector<number> &data);
 
    template <typename number>
    unsigned int
    get_container_size(const FullMatrix<number> &data);
 
    template <typename Container>
    std::enable_if_t<
      std::is_same_v<Container, std::vector<typename Container::value_type>>,
      Container>
    initialize_container(const std::vector<hsize_t> &dimensions);
 
    template <typename Container>
    std::enable_if_t<
      std::is_same_v<Container, Vector<typename Container::value_type>>,
      Container>
    initialize_container(const std::vector<hsize_t> &dimensions);
 
    template <typename Container>
    std::enable_if_t<
      std::is_same_v<Container, FullMatrix<typename Container::value_type>>,
      Container>
    initialize_container(const std::vector<hsize_t> &dimensions);
 
    inline void
    set_plist(hid_t &plist, const bool mpi);
 
    inline void
    release_plist(hid_t                     &plist,
                  H5D_mpio_actual_io_mode_t &io_mode,
                  std::uint32_t             &local_no_collective_cause,
                  std::uint32_t             &global_no_collective_cause,
                  const bool                 mpi,
                  const bool                 query_io_mode);
 
    inline std::string
    no_collective_cause_to_string(const std::uint32_t no_collective_cause);
  } // namespace internal
 
 
 
  // definitions
 
  namespace internal
  {
    template <typename number>
    std::shared_ptr<hid_t>
    get_hdf5_datatype()
    {
      static_assert(std::is_same_v<number, float> ||
                      std::is_same_v<number, double> ||
                      std::is_same_v<number, int> ||
                      std::is_same_v<number, bool> ||
                      std::is_same_v<number, unsigned int> ||
                      std::is_same_v<number, std::complex<float>> ||
                      std::is_same_v<number, std::complex<double>>,
                    "The data type you are trying to get the HDF5 tag for "
                    "is not supported by this function.");
 
      if (std::is_same_v<number, float>)
        {
          return std::make_shared<hid_t>(H5T_NATIVE_FLOAT);
        }
      else if (std::is_same_v<number, double>)
        {
          return std::make_shared<hid_t>(H5T_NATIVE_DOUBLE);
        }
      else if (std::is_same_v<number, int>)
        {
          return std::make_shared<hid_t>(H5T_NATIVE_INT);
        }
      else if (std::is_same_v<number, unsigned int>)
        {
          return std::make_shared<hid_t>(H5T_NATIVE_UINT);
        }
      else if (std::is_same_v<number, bool>)
        {
          return std::make_shared<hid_t>(H5T_NATIVE_HBOOL);
        }
      else if (std::is_same_v<number, std::complex<float>>)
        {
          std::shared_ptr<hid_t> t_type =
            std::shared_ptr<hid_t>(new hid_t, [](hid_t *pointer) {
              // Release the HDF5 resource
              const herr_t ret = H5Tclose(*pointer);
              AssertNothrow(ret >= 0, ExcInternalError());
              (void)ret;
              delete pointer;
            });
 
          *t_type = H5Tcreate(H5T_COMPOUND, sizeof(std::complex<float>));
          //  The C++ standards committee agreed to mandate that the storage
          //  format used for the std::complex type be binary-compatible with
          //  the C99 type, i.e. an array T[2] with consecutive real [0] and
          //  imaginary [1] parts.
          herr_t ret = H5Tinsert(*t_type, "r", 0, H5T_NATIVE_FLOAT);
          Assert(ret >= 0, ExcInternalError());
          ret = H5Tinsert(*t_type, "i", sizeof(float), H5T_NATIVE_FLOAT);
          Assert(ret >= 0, ExcInternalError());
          (void)ret;
          return t_type;
        }
      else if (std::is_same_v<number, std::complex<double>>)
        {
          std::shared_ptr<hid_t> t_type =
            std::shared_ptr<hid_t>(new hid_t, [](hid_t *pointer) {
              // Release the HDF5 resource
              const herr_t ret = H5Tclose(*pointer);
              AssertNothrow(ret >= 0, ExcInternalError());
              (void)ret;
              delete pointer;
            });
          *t_type = H5Tcreate(H5T_COMPOUND, sizeof(std::complex<double>));
          //  The C++ standards committee agreed to mandate that the storage
          //  format used for the std::complex type be binary-compatible with
          //  the C99 type, i.e. an array T[2] with consecutive real [0] and
          //  imaginary [1] parts.
          herr_t ret = H5Tinsert(*t_type, "r", 0, H5T_NATIVE_DOUBLE);
          Assert(ret >= 0, ExcInternalError());
          ret = H5Tinsert(*t_type, "i", sizeof(double), H5T_NATIVE_DOUBLE);
          Assert(ret >= 0, ExcInternalError());
          (void)ret;
          return t_type;
        }
 
      // The function should not reach this point
      DEAL_II_ASSERT_UNREACHABLE();
      return {};
    }
 
 
 
    template <typename number>
    std::vector<hsize_t>
    get_container_dimensions(const std::vector<number> &data)
    {
      std::vector<hsize_t> dimensions = {data.size()};
      return dimensions;
    }
 
 
 
    template <typename number>
    std::vector<hsize_t>
    get_container_dimensions(const Vector<number> &data)
    {
      std::vector<hsize_t> dimensions = {data.size()};
      return dimensions;
    }
 
 
 
    template <typename number>
    std::vector<hsize_t>
    get_container_dimensions(const FullMatrix<number> &data)
    {
      std::vector<hsize_t> dimensions = {data.m(), data.n()};
      return dimensions;
    }
 
 
 
    template <typename number>
    unsigned int
    get_container_size(const std::vector<number> &data)
    {
      return static_cast<unsigned int>(data.size());
    }
 
 
 
    template <typename number>
    unsigned int
    get_container_size(const Vector<number> &data)
    {
      return static_cast<unsigned int>(data.size());
    }
 
 
 
    template <typename number>
    unsigned int
    get_container_size(const FullMatrix<number> &data)
    {
      return static_cast<unsigned int>(data.m() * data.n());
    }
 
 
 
    template <typename Container>
    std::enable_if_t<
      std::is_same_v<Container, std::vector<typename Container::value_type>>,
      Container>
    initialize_container(const std::vector<hsize_t> &dimensions)
    {
      return Container(std::accumulate(
        dimensions.begin(), dimensions.end(), 1, std::multiplies<int>()));
    }
 
 
 
    template <typename Container>
    std::enable_if_t<
      std::is_same_v<Container, Vector<typename Container::value_type>>,
      Container>
    initialize_container(const std::vector<hsize_t> &dimensions)
    {
      return Container(std::accumulate(
        dimensions.begin(), dimensions.end(), 1, std::multiplies<int>()));
    }
 
 
 
    template <typename Container>
    std::enable_if_t<
      std::is_same_v<Container, FullMatrix<typename Container::value_type>>,
      Container>
    initialize_container(const std::vector<hsize_t> &dimensions)
    {
      // If the rank is higher than 2, then remove single-dimensional entries
      // from the shape defined by dimensions. This is equivalent to the squeeze
      // function of python/numpy. For example the following code would convert
      // the vector {1,3,1,2} to {3,2}
      std::vector<hsize_t> squeezed_dimensions;
 
      if (dimensions.size() > 2)
        {
          for (const auto &dimension : dimensions)
            {
              if (dimension > 1)
                squeezed_dimensions.push_back(dimension);
            }
        }
      else
        {
          squeezed_dimensions = dimensions;
        }
 
      AssertDimension(squeezed_dimensions.size(), 2);
      return Container(squeezed_dimensions[0], squeezed_dimensions[1]);
    }
 
 
    inline void
    set_plist(hid_t &plist, const bool mpi)
    {
      if (mpi)
        {
#  ifdef DEAL_II_WITH_MPI
          plist = H5Pcreate(H5P_DATASET_XFER);
          Assert(plist >= 0, ExcInternalError());
          const herr_t ret = H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
          (void)ret;
          Assert(ret >= 0, ExcInternalError());
#  else
          AssertThrow(false, ExcNotImplemented());
#  endif
        }
      else
        {
          plist = H5P_DEFAULT;
        }
 
      (void)plist;
      (void)mpi;
    }
 
 
    inline void
    release_plist(hid_t                     &plist,
                  H5D_mpio_actual_io_mode_t &io_mode,
                  std::uint32_t             &local_no_collective_cause,
                  std::uint32_t             &global_no_collective_cause,
                  const bool                 mpi,
                  const bool                 query_io_mode)
    {
      if (mpi)
        {
#  ifdef DEAL_II_WITH_MPI
          herr_t ret;
          (void)ret;
          if (query_io_mode)
            {
              ret = H5Pget_mpio_actual_io_mode(plist, &io_mode);
              Assert(ret >= 0, ExcInternalError());
              ret =
                H5Pget_mpio_no_collective_cause(plist,
                                                &local_no_collective_cause,
                                                &global_no_collective_cause);
              Assert(ret >= 0, ExcInternalError());
            }
          ret = H5Pclose(plist);
          Assert(ret >= 0, ExcInternalError());
#  else
          AssertThrow(false, ExcNotImplemented());
#  endif
        }
 
      (void)plist;
      (void)io_mode;
      (void)local_no_collective_cause;
      (void)global_no_collective_cause;
      (void)mpi;
      (void)query_io_mode;
    }
 
 
    inline std::string
    no_collective_cause_to_string(const std::uint32_t no_collective_cause)
    {
      std::string message;
 
      auto append_to_message = [&message](const char *p) {
        if (message.size() > 0)
          message += ", ";
        message += p;
      };
 
      // The first is not a bitmask comparison, the rest are bitmask
      // comparisons.
      // https://support.hdfgroup.org/HDF5/doc/RM/RM_H5P.html#Property-GetMpioNoCollectiveCause
      // See H5Ppublic.h
      // Hex codes are used because the HDF5 Group can deprecate some of the
      // enum codes. For example the enum code H5D_MPIO_FILTERS is not defined
      // in 1.10.2 because it is possible to use compressed datasets with the
      // MPI/IO driver.
 
      // H5D_MPIO_COLLECTIVE
      if (no_collective_cause == 0x00)
        {
          append_to_message("H5D_MPIO_COLLECTIVE");
        }
      // H5D_MPIO_SET_INDEPENDENT
      if ((no_collective_cause & 0x01) == 0x01)
        {
          append_to_message("H5D_MPIO_SET_INDEPENDENT");
        }
      // H5D_MPIO_DATATYPE_CONVERSION
      if ((no_collective_cause & 0x02) == 0x02)
        {
          append_to_message("H5D_MPIO_DATATYPE_CONVERSION");
        }
      // H5D_MPIO_DATA_TRANSFORMS
      if ((no_collective_cause & 0x04) == 0x04)
        {
          append_to_message("H5D_MPIO_DATA_TRANSFORMS");
        }
      // H5D_MPIO_NOT_SIMPLE_OR_SCALAR_DATASPACES
      if ((no_collective_cause & 0x10) == 0x10)
        {
          append_to_message("H5D_MPIO_NOT_SIMPLE_OR_SCALAR_DATASPACES");
        }
      // H5D_MPIO_NOT_CONTIGUOUS_OR_CHUNKED_DATASET
      if ((no_collective_cause & 0x20) == 0x20)
        {
          append_to_message("H5D_MPIO_NOT_CONTIGUOUS_OR_CHUNKED_DATASET");
        }
      // H5D_MPIO_FILTERS
      if ((no_collective_cause & 0x40) == 0x40)
        {
          append_to_message("H5D_MPIO_FILTERS");
        }
      return message;
    }
  } // namespace internal
 
 
  template <typename T>
  T
  HDF5Object::get_attribute(const std::string &attr_name) const
  {
    const std::shared_ptr<hid_t> t_type = internal::get_hdf5_datatype<T>();
    T                            value;
    hid_t                        attr;
    herr_t                       ret;
 
    attr = H5Aopen(*hdf5_reference, attr_name.data(), H5P_DEFAULT);
    Assert(attr >= 0, ExcMessage("Error at H5Aopen"));
    (void)ret;
    ret = H5Aread(attr, *t_type, &value);
    Assert(ret >= 0, ExcMessage("Error at H5Aread"));
    (void)ret;
    ret = H5Aclose(attr);
    Assert(ret >= 0, ExcMessage("Error at H5Aclose"));
 
    return value;
  }
 
 
 
  template <>
  inline std::string
  HDF5Object::get_attribute(const std::string &attr_name) const
  {
    // Reads a UTF8 variable string
    //
    // code inspired from
    // https://support.hdfgroup.org/ftp/HDF5/examples/misc-examples/vlstratt.c
    //
    // In the case of a variable length string the user does not have to reserve
    // memory for string_out. H5Aread will reserve the memory and the
    // user has to free the memory.
    //
    // Todo:
    // - Use H5Dvlen_reclaim instead of free
 
    char  *string_out;
    hid_t  attr;
    hid_t  type;
    herr_t ret;
 
    /* Create a datatype to refer to. */
    type = H5Tcopy(H5T_C_S1);
    Assert(type >= 0, ExcInternalError());
 
    // Python strings are encoded in UTF8
    ret = H5Tset_cset(type, H5T_CSET_UTF8);
    Assert(type >= 0, ExcInternalError());
 
    ret = H5Tset_size(type, H5T_VARIABLE);
    Assert(ret >= 0, ExcInternalError());
 
    attr = H5Aopen(*hdf5_reference, attr_name.data(), H5P_DEFAULT);
    Assert(attr >= 0, ExcInternalError());
 
    ret = H5Aread(attr, type, &string_out);
    Assert(ret >= 0, ExcInternalError());
 
    std::string string_value(string_out);
    // The memory of the variable length string has to be freed.
    // H5Dvlen_reclaim could be also used
    free(string_out);
    ret = H5Tclose(type);
    Assert(ret >= 0, ExcInternalError());
 
    ret = H5Aclose(attr);
    Assert(ret >= 0, ExcInternalError());
 
    (void)ret;
    return string_value;
  }
 
 
 
  template <typename T>
  void
  HDF5Object::set_attribute(const std::string &attr_name, const T value)
  {
    hid_t  attr;
    hid_t  aid;
    herr_t ret;
 
    const std::shared_ptr<hid_t> t_type = internal::get_hdf5_datatype<T>();
 
 
    /*
     * Create scalar attribute.
     */
    aid = H5Screate(H5S_SCALAR);
    Assert(aid >= 0, ExcMessage("Error at H5Screate"));
    attr = H5Acreate2(*hdf5_reference,
                      attr_name.data(),
                      *t_type,
                      aid,
                      H5P_DEFAULT,
                      H5P_DEFAULT);
    Assert(attr >= 0, ExcMessage("Error at H5Acreate2"));
 
    /*
     * Write scalar attribute.
     */
    ret = H5Awrite(attr, *t_type, &value);
    Assert(ret >= 0, ExcMessage("Error at H5Awrite"));
 
    ret = H5Sclose(aid);
    Assert(ret >= 0, ExcMessage("Error at H5Sclose"));
    ret = H5Aclose(attr);
    Assert(ret >= 0, ExcMessage("Error at H5Aclose"));
 
    (void)ret;
  }
 
 
 
  template <>
  inline void
  HDF5Object::set_attribute(const std::string &attr_name,
                            const std::string  value) // NOLINT
  {
    // Writes a UTF8 variable string
    //
    // code inspired from
    // https://support.hdfgroup.org/ftp/HDF5/examples/misc-examples/vlstratt.c
 
    hid_t  attr;
    hid_t  aid;
    hid_t  t_type;
    herr_t ret;
 
    /* Create a datatype to refer to. */
    t_type = H5Tcopy(H5T_C_S1);
    Assert(t_type >= 0, ExcInternalError());
 
    // Python strings are encoded in UTF8
    ret = H5Tset_cset(t_type, H5T_CSET_UTF8);
    Assert(t_type >= 0, ExcInternalError());
 
    ret = H5Tset_size(t_type, H5T_VARIABLE);
    Assert(ret >= 0, ExcInternalError());
 
    /*
     * Create scalar attribute.
     */
    aid = H5Screate(H5S_SCALAR);
    Assert(aid >= 0, ExcMessage("Error at H5Screate"));
    attr = H5Acreate2(
      *hdf5_reference, attr_name.data(), t_type, aid, H5P_DEFAULT, H5P_DEFAULT);
    Assert(attr >= 0, ExcMessage("Error at H5Acreate2"));
 
    /*
     * Write scalar attribute.
     * In most of the cases H5Awrite and H5Dwrite take a pointer to the data.
     * But in the particular case of a variable length string, H5Awrite takes
     * the address of the pointer of the string.
     */
    const char *c_string_value = value.c_str();
    ret                        = H5Awrite(attr, t_type, &c_string_value);
    Assert(ret >= 0, ExcInternalError());
 
    ret = H5Sclose(aid);
    Assert(ret >= 0, ExcMessage("Error at H5Sclose"));
    ret = H5Aclose(attr);
    Assert(ret >= 0, ExcMessage("Error at H5Aclose"));
 
    (void)ret;
  }
 
 
 
  template <typename Container>
  Container
  DataSet::read()
  {
    const std::shared_ptr<hid_t> t_type =
      internal::get_hdf5_datatype<typename Container::value_type>();
    hid_t  plist;
    herr_t ret;
 
    Container data = internal::initialize_container<Container>(dimensions);
 
    internal::set_plist(plist, mpi);
 
    ret = H5Dread(*hdf5_reference,
                  *t_type,
                  H5S_ALL,
                  H5S_ALL,
                  plist,
                  make_array_view(data).data());
    Assert(ret >= 0, ExcInternalError());
 
 
    internal::release_plist(plist,
                            io_mode,
                            local_no_collective_cause,
                            global_no_collective_cause,
                            mpi,
                            query_io_mode);
 
    (void)ret;
    return data;
  }
 
 
 
  template <typename Container>
  Container
  DataSet::read_selection(const std::vector<hsize_t> &coordinates)
  {
    Assert(coordinates.size() % rank == 0,
           ExcMessage(
             "The dimension of coordinates has to be divisible by the rank"));
    const std::shared_ptr<hid_t> t_type =
      internal::get_hdf5_datatype<typename Container::value_type>();
    hid_t  plist;
    hid_t  memory_dataspace;
    herr_t ret;
 
    std::vector<hsize_t> data_dimensions{
      static_cast<hsize_t>(coordinates.size() / rank)};
 
    Container data = internal::initialize_container<Container>(data_dimensions);
 
    memory_dataspace = H5Screate_simple(1, data_dimensions.data(), nullptr);
    Assert(memory_dataspace >= 0, ExcMessage("Error at H5Screate_simple"));
    ret = H5Sselect_elements(*dataspace,
                             H5S_SELECT_SET,
                             data.size(),
                             coordinates.data());
    Assert(ret >= 0, ExcMessage("Error at H5Sselect_elements"));
 
    internal::set_plist(plist, mpi);
 
    ret = H5Dread(*hdf5_reference,
                  *t_type,
                  memory_dataspace,
                  *dataspace,
                  plist,
                  make_array_view(data).data());
    Assert(ret >= 0, ExcMessage("Error at H5Dread"));
 
    internal::release_plist(plist,
                            io_mode,
                            local_no_collective_cause,
                            global_no_collective_cause,
                            mpi,
                            query_io_mode);
 
    ret = H5Sclose(memory_dataspace);
    Assert(ret >= 0, ExcMessage("Error at H5SClose"));
 
    (void)ret;
    return data;
  }
 
 
 
  template <typename Container>
  Container
  DataSet::read_hyperslab(const std::vector<hsize_t> &offset,
                          const std::vector<hsize_t> &count)
  {
    const std::shared_ptr<hid_t> t_type =
      internal::get_hdf5_datatype<typename Container::value_type>();
    hid_t  plist;
    hid_t  memory_dataspace;
    herr_t ret;
 
    // In this particular overload of read_hyperslab the data_dimensions are
    // the same as count
    std::vector<hsize_t> data_dimensions = count;
 
    Container data = internal::initialize_container<Container>(data_dimensions);
 
    memory_dataspace =
      H5Screate_simple(data_dimensions.size(), data_dimensions.data(), nullptr);
    Assert(memory_dataspace >= 0, ExcMessage("Error at H5Screate_simple"));
    ret = H5Sselect_hyperslab(*dataspace,
                              H5S_SELECT_SET,
                              offset.data(),
                              nullptr,
                              count.data(),
                              nullptr);
    Assert(ret >= 0, ExcMessage("Error at H5Sselect_hyperslab"));
 
    internal::set_plist(plist, mpi);
 
    ret = H5Dread(*hdf5_reference,
                  *t_type,
                  memory_dataspace,
                  *dataspace,
                  plist,
                  make_array_view(data).data());
    Assert(ret >= 0, ExcMessage("Error at H5Dread"));
 
    internal::release_plist(plist,
                            io_mode,
                            local_no_collective_cause,
                            global_no_collective_cause,
                            mpi,
                            query_io_mode);
 
    ret = H5Sclose(memory_dataspace);
    Assert(ret >= 0, ExcMessage("Error at H5SClose"));
 
    (void)ret;
    return data;
  }
 
 
 
  template <typename Container>
  Container
  DataSet::read_hyperslab(const std::vector<hsize_t> &data_dimensions,
                          const std::vector<hsize_t> &offset,
                          const std::vector<hsize_t> &stride,
                          const std::vector<hsize_t> &count,
                          const std::vector<hsize_t> &block)
  {
    const std::shared_ptr<hid_t> t_type =
      internal::get_hdf5_datatype<typename Container::value_type>();
    hid_t  plist;
    hid_t  memory_dataspace;
    herr_t ret;
 
    Container data = internal::initialize_container<Container>(data_dimensions);
 
    memory_dataspace =
      H5Screate_simple(data_dimensions.size(), data_dimensions.data(), nullptr);
    Assert(memory_dataspace >= 0, ExcMessage("Error at H5Screate_simple"));
    ret = H5Sselect_hyperslab(*dataspace,
                              H5S_SELECT_SET,
                              offset.data(),
                              stride.data(),
                              count.data(),
                              block.data());
    Assert(ret >= 0, ExcMessage("Error at H5Sselect_hyperslab"));
 
    internal::set_plist(plist, mpi);
 
    ret = H5Dread(*hdf5_reference,
                  *t_type,
                  memory_dataspace,
                  *dataspace,
                  plist,
                  make_array_view(data).data());
    Assert(ret >= 0, ExcMessage("Error at H5Dread"));
 
    internal::release_plist(plist,
                            io_mode,
                            local_no_collective_cause,
                            global_no_collective_cause,
                            mpi,
                            query_io_mode);
 
    ret = H5Sclose(memory_dataspace);
    Assert(ret >= 0, ExcMessage("Error at H5SClose"));
 
    (void)ret;
    return data;
  }
 
 
 
  template <typename number>
  void
  DataSet::read_none()
  {
    const std::shared_ptr<hid_t> t_type = internal::get_hdf5_datatype<number>();
    const std::vector<hsize_t>   data_dimensions = {0};
 
    hid_t  memory_dataspace;
    hid_t  plist;
    herr_t ret;
 
    memory_dataspace = H5Screate_simple(1, data_dimensions.data(), nullptr);
    Assert(memory_dataspace >= 0, ExcMessage("Error at H5Screate_simple"));
    ret = H5Sselect_none(*dataspace);
    Assert(ret >= 0, ExcMessage("H5Sselect_none"));
 
    internal::set_plist(plist, mpi);
 
    // The pointer of data can safely be nullptr, see the discussion at the HDF5
    // forum:
    // https://forum.hdfgroup.org/t/parallel-i-o-does-not-support-filters-yet/884/17
    ret = H5Dread(
      *hdf5_reference, *t_type, memory_dataspace, *dataspace, plist, nullptr);
    Assert(ret >= 0, ExcMessage("Error at H5Dread"));
 
    internal::release_plist(plist,
                            io_mode,
                            local_no_collective_cause,
                            global_no_collective_cause,
                            mpi,
                            query_io_mode);
 
    ret = H5Sclose(memory_dataspace);
    Assert(ret >= 0, ExcMessage("Error at H5SClose"));
 
    (void)ret;
  }
 
 
 
  template <typename Container>
  void
  DataSet::write(const Container &data)
  {
    AssertDimension(size, internal::get_container_size(data));
    const std::shared_ptr<hid_t> t_type =
      internal::get_hdf5_datatype<typename Container::value_type>();
    hid_t  plist;
    herr_t ret;
 
    internal::set_plist(plist, mpi);
 
    ret = H5Dwrite(*hdf5_reference,
                   *t_type,
                   H5S_ALL,
                   H5S_ALL,
                   plist,
                   make_array_view(data).data());
    Assert(ret >= 0, ExcMessage("Error at H5Dwrite"));
 
    internal::release_plist(plist,
                            io_mode,
                            local_no_collective_cause,
                            global_no_collective_cause,
                            mpi,
                            query_io_mode);
 
    (void)ret;
  }
 
 
 
  template <typename Container>
  void
  DataSet::write_selection(const Container            &data,
                           const std::vector<hsize_t> &coordinates)
  {
    AssertDimension(coordinates.size(), data.size() * rank);
    const std::shared_ptr<hid_t> t_type =
      internal::get_hdf5_datatype<typename Container::value_type>();
    const std::vector<hsize_t> data_dimensions =
      internal::get_container_dimensions(data);
 
    hid_t  memory_dataspace;
    hid_t  plist;
    herr_t ret;
 
 
    memory_dataspace = H5Screate_simple(1, data_dimensions.data(), nullptr);
    Assert(memory_dataspace >= 0, ExcMessage("Error at H5Screate_simple"));
    ret = H5Sselect_elements(*dataspace,
                             H5S_SELECT_SET,
                             data.size(),
                             coordinates.data());
    Assert(ret >= 0, ExcMessage("Error at H5Sselect_elements"));
 
    internal::set_plist(plist, mpi);
 
    ret = H5Dwrite(*hdf5_reference,
                   *t_type,
                   memory_dataspace,
                   *dataspace,
                   plist,
                   make_array_view(data).data());
    Assert(ret >= 0, ExcMessage("Error at H5Dwrite"));
 
    internal::release_plist(plist,
                            io_mode,
                            local_no_collective_cause,
                            global_no_collective_cause,
                            mpi,
                            query_io_mode);
 
    ret = H5Sclose(memory_dataspace);
    Assert(ret >= 0, ExcMessage("Error at H5SClose"));
 
    (void)ret;
  }
 
 
 
  template <typename Container>
  void
  DataSet::write_hyperslab(const Container            &data,
                           const std::vector<hsize_t> &offset,
                           const std::vector<hsize_t> &count)
  {
    AssertDimension(std::accumulate(count.begin(),
                                    count.end(),
                                    1,
                                    std::multiplies<unsigned int>()),
                    internal::get_container_size(data));
    const std::shared_ptr<hid_t> t_type =
      internal::get_hdf5_datatype<typename Container::value_type>();
    // In this particular overload of write_hyperslab the data_dimensions are
    // the same as count
    const std::vector<hsize_t> &data_dimensions = count;
 
    hid_t  memory_dataspace;
    hid_t  plist;
    herr_t ret;
 
    memory_dataspace =
      H5Screate_simple(data_dimensions.size(), data_dimensions.data(), nullptr);
    Assert(memory_dataspace >= 0, ExcMessage("Error at H5Screate_simple"));
    ret = H5Sselect_hyperslab(*dataspace,
                              H5S_SELECT_SET,
                              offset.data(),
                              nullptr,
                              count.data(),
                              nullptr);
    Assert(ret >= 0, ExcMessage("Error at H5Sselect_hyperslab"));
 
    internal::set_plist(plist, mpi);
 
    ret = H5Dwrite(*hdf5_reference,
                   *t_type,
                   memory_dataspace,
                   *dataspace,
                   plist,
                   make_array_view(data).data());
    Assert(ret >= 0, ExcMessage("Error at H5Dwrite"));
 
    internal::release_plist(plist,
                            io_mode,
                            local_no_collective_cause,
                            global_no_collective_cause,
                            mpi,
                            query_io_mode);
 
    ret = H5Sclose(memory_dataspace);
    Assert(ret >= 0, ExcMessage("Error at H5Sclose"));
 
    (void)ret;
  }
 
 
 
  template <typename Container>
  void
  DataSet::write_hyperslab(const Container            &data,
                           const std::vector<hsize_t> &data_dimensions,
                           const std::vector<hsize_t> &offset,
                           const std::vector<hsize_t> &stride,
                           const std::vector<hsize_t> &count,
                           const std::vector<hsize_t> &block)
  {
    const std::shared_ptr<hid_t> t_type =
      internal::get_hdf5_datatype<typename Container::value_type>();
 
    hid_t  memory_dataspace;
    hid_t  plist;
    herr_t ret;
 
    memory_dataspace =
      H5Screate_simple(data_dimensions.size(), data_dimensions.data(), nullptr);
    Assert(memory_dataspace >= 0, ExcMessage("Error at H5Screate_simple"));
    ret = H5Sselect_hyperslab(*dataspace,
                              H5S_SELECT_SET,
                              offset.data(),
                              stride.data(),
                              count.data(),
                              block.data());
    Assert(ret >= 0, ExcMessage("Error at H5Sselect_hyperslab"));
 
    internal::set_plist(plist, mpi);
 
    ret = H5Dwrite(*hdf5_reference,
                   *t_type,
                   memory_dataspace,
                   *dataspace,
                   plist,
                   make_array_view(data).data());
    Assert(ret >= 0, ExcMessage("Error at H5Dwrite"));
 
    internal::release_plist(plist,
                            io_mode,
                            local_no_collective_cause,
                            global_no_collective_cause,
                            mpi,
                            query_io_mode);
 
    ret = H5Sclose(memory_dataspace);
    Assert(ret >= 0, ExcMessage("Error at H5Sclose"));
 
    (void)ret;
  }
 
 
 
  template <typename number>
  void
  DataSet::write_none()
  {
    std::shared_ptr<hid_t> t_type = internal::get_hdf5_datatype<number>();
    std::vector<hsize_t>   data_dimensions = {0};
 
    hid_t  memory_dataspace;
    hid_t  plist;
    herr_t ret;
 
    memory_dataspace = H5Screate_simple(1, data_dimensions.data(), nullptr);
    Assert(memory_dataspace >= 0, ExcMessage("Error at H5Screate_simple"));
    ret = H5Sselect_none(*dataspace);
    Assert(ret >= 0, ExcMessage("Error at H5PSselect_none"));
 
    internal::set_plist(plist, mpi);
 
    // The pointer of data can safely be nullptr, see the discussion at the HDF5
    // forum:
    // https://forum.hdfgroup.org/t/parallel-i-o-does-not-support-filters-yet/884/17
    ret = H5Dwrite(
      *hdf5_reference, *t_type, memory_dataspace, *dataspace, plist, nullptr);
    Assert(ret >= 0, ExcMessage("Error at H5Dwrite"));
 
    internal::release_plist(plist,
                            io_mode,
                            local_no_collective_cause,
                            global_no_collective_cause,
                            mpi,
                            query_io_mode);
 
    ret = H5Sclose(memory_dataspace);
    Assert(ret >= 0, ExcMessage("Error at H5Sclose"));
 
    (void)ret;
  }
 
 
 
  template <typename number>
  DataSet
  Group::create_dataset(const std::string          &name,
                        const std::vector<hsize_t> &dimensions) const
  {
    std::shared_ptr<hid_t> t_type = internal::get_hdf5_datatype<number>();
    return {name, *hdf5_reference, dimensions, t_type, mpi};
  }
 
 
 
  template <typename Container>
  void
  Group::write_dataset(const std::string &name, const Container &data) const
  {
    std::vector<hsize_t> dimensions = internal::get_container_dimensions(data);
    auto                 dataset =
      create_dataset<typename Container::value_type>(name, dimensions);
    dataset.write(data);
  }
} // namespace HDF5
 
DEAL_II_NAMESPACE_CLOSE
 
 
#endif // DEAL_II_WITH_HDF5
 
#endif // dealii_hdf5_h