From 16143ccd17fbb95c4f4bef7de10ebb110ff2dd5c Mon Sep 17 00:00:00 2001
From: Daniel Garcia-Sanchez <daniel.garcia-sanchez@insp.upmc.fr>
Date: Sun, 19 Aug 2018 21:16:44 +0200
Subject: [PATCH] Add functions read_*()

Add read_selection()

Add read_hyperslab()

Add read_none()
---
 include/deal.II/base/hdf5.h | 104 ++++++++++--
 source/base/hdf5.cc         | 307 ++++++++++++++++++++++++++++++++++--
 2 files changed, 383 insertions(+), 28 deletions(-)

diff --git a/include/deal.II/base/hdf5.h b/include/deal.II/base/hdf5.h
index ebdd9b4af2..bce1896103 100644
--- a/include/deal.II/base/hdf5.h
+++ b/include/deal.II/base/hdf5.h
@@ -128,8 +128,9 @@ namespace HDF5
     };
 
     /**
-     * Reads an attribute. T can be double, int, unsigned int, bool,
-     * std::complex<double> or std::string. Note that the encoding of
+     * Reads an attribute. T can be float, double, long double,
+     * std::complex<float>, std::complex<double>, std::complex<long double>,
+     * int, unsigned int, bool or std::string. Note that the encoding of
      * std::string is UTF8 in order to be compatible with python3.
      *
      * Datatype conversion takes place at the time of a read or write and is
@@ -143,8 +144,9 @@ namespace HDF5
     attr(const std::string attr_name) const;
 
     /**
-     * Writes an attribute. T can be double, int, unsigned int, bool,
-     * std::complex<double> or std::string. Note that the encoding of
+     * Reads an attribute. T can be float, double, long double,
+     * std::complex<float>, std::complex<double>, std::complex<long double>,
+     * int, unsigned int, bool or std::string. Note that the encoding of
      * std::string is UTF8 in order to be compatible with python3.
      *
      * Datatype conversion takes place at the time of a read or write and is
@@ -188,8 +190,9 @@ namespace HDF5
 
   public:
     /**
-     * Reads data of the dataset. T can be double, int, unsigned int, bool
-     * or std::complex<double>.
+     * Reads data of the dataset. T can be float, double, long double,
+     * std::complex<float>, std::complex<double>, std::complex<long double>, int
+     * or unsigned int.
      *
      * Datatype conversion takes place at the time of a read or write and is
      * automatic. See the <a
@@ -202,8 +205,81 @@ namespace HDF5
     read();
 
     /**
-     * Writes data in the dataset. T can be double, int, unsigned int,
-     * or std::complex<double>.
+     * Reads data of a subset of the dataset. T can be float, double, long
+     * double, std::complex<float>, std::complex<double>,
+     * std::complex<long double>, int or unsigned int.
+     *
+     * The selected elements can be scattered and take any shape in the dataset.
+     * For examples, in the 4D case, we will be selecting three points and a 4D
+     * dataspace has rank 4, so the selection will be described in a 3-by-4
+     * array. To select the points (1,1,1,1), (14,6,12,18), and (8,22,30,22),
+     * the point selection array would be as follows:
+     *
+     *    0  0  0  0
+     *
+     *   13  5 11 17
+     *
+     *    7 21 29 21
+     *
+     * <a
+     * href="https://support.hdfgroup.org/newsletters/newsletter140.html">Parallel
+     * HDF5 supports collective I/O on point selections.</a>
+     *
+     * Datatype conversion takes place at the time of a read or write and is
+     * automatic. See the <a
+     * href="https://support.hdfgroup.org/HDF5/doc/UG/HDF5_Users_Guide-Responsive%20HTML5/index.html#t=HDF5_Users_Guide%2FDatatypes%2FHDF5_Datatypes.htm%23TOC_6_10_Data_Transferbc-26&rhtocid=6.5_2">Data
+     * Transfer: Datatype Conversion and Selection</a>  section in the HDF5
+     * User's Guide.
+     */
+    template <typename T>
+    std::vector<T>
+    read_selection(const std::vector<hsize_t> coordinates);
+
+    /**
+     * Reads a hyperslab from the dataset. The parameters are summarized
+     * below:
+     *  - Offset: The starting location for the hyperslab.
+     *  - Count: The number of elements to select along each dimension.
+     *
+     * Stride and block are set to NULL.
+     *
+     * See the <a
+     * href="https://support.hdfgroup.org/HDF5/doc/UG/HDF5_Users_Guide-Responsive%20HTML5/HDF5_Users_Guide/Dataspaces/HDF5_Dataspaces_and_Partial_I_O.htm?rhtocid=7.2#TOC_7_4_Dataspaces_and_Databc-6">Dataspaces
+     * and Data Transfer</a>  section in the HDF5 User's Guide. See as well the
+     * <a
+     * href="https://support.hdfgroup.org/HDF5/doc1.8/RM/RM_H5S.html#Dataspace-SelectHyperslab">H5Sselect_hyperslab
+     * definition</a>.
+     *
+     * Datatype conversion takes place at the time of a read or read and is
+     * automatic. See the <a
+     * href="https://support.hdfgroup.org/HDF5/doc/UG/HDF5_Users_Guide-Responsive%20HTML5/index.html#t=HDF5_Users_Guide%2FDatatypes%2FHDF5_Datatypes.htm%23TOC_6_10_Data_Transferbc-26&rhtocid=6.5_2">Data
+     * Transfer: Datatype Conversion and Selection</a>  section in the HDF5
+     * User's Guide.
+     */
+    template <template <class...> class Container, typename T>
+    Container<T>
+    read_hyperslab(const std::vector<hsize_t> &offset,
+                   const std::vector<hsize_t> &count);
+
+    /**
+     * This function does not read any data, but can contribute to a collective
+     * read call. T can be float, double, long double, std::complex<float>,
+     * std::complex<double>, std::complex<long double>, int or unsigned int.
+     *
+     * Datatype conversion takes place at the time of a read or write and is
+     * automatic. See the <a
+     * href="https://support.hdfgroup.org/HDF5/doc/UG/HDF5_Users_Guide-Responsive%20HTML5/index.html#t=HDF5_Users_Guide%2FDatatypes%2FHDF5_Datatypes.htm%23TOC_6_10_Data_Transferbc-26&rhtocid=6.5_2">Data
+     * Transfer: Datatype Conversion and Selection</a>  section in the HDF5
+     * User's Guide.
+     */
+    template <typename T>
+    void
+    read_none();
+
+    /**
+     * Writes data in the dataset. T can be float, double, long double,
+     * std::complex<float>, std::complex<double>, std::complex<long double>, int
+     * or unsigned int.
      *
      * Datatype conversion takes place at the time of a read or write and is
      * automatic. See the <a
@@ -216,8 +292,9 @@ namespace HDF5
     write(const Container<T> &data);
 
     /**
-     * Writes data to a subset of the dataset. T can be double, int, unsigned
-     * int, bool or std::complex<double>.
+     * Writes data to a subset of the dataset. T can be float, double, long
+     * double, std::complex<float>, std::complex<double>,
+     * std::complex<long double>, int or unsigned int.
      *
      * The selected elements can be scattered and take any shape in the dataset.
      * For examples, in the 4D case, we will be selecting three points and a 4D
@@ -255,7 +332,7 @@ namespace HDF5
     // clang-format off
     /**
      *
-     * stride and block are set to NULL. For more complex hyperslabs see
+     * Stride and block are set to NULL. For more complex hyperslabs see
      * write_hyperslab(const Container<T> &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).
      */
     // clang-format on
@@ -313,11 +390,10 @@ namespace HDF5
                     const std::vector<hsize_t> &count,
                     const std::vector<hsize_t> &block);
 
-
     /**
      * This function does not write any data, but can contribute to a collective
-     * write call. T can be double, int, unsigned int, bool or
-     * std::complex<double>.
+     * write call. T can be float, double, long double, std::complex<float>,
+     * std::complex<double>, std::complex<long double>, int or unsigned int.
      *
      * Datatype conversion takes place at the time of a read or write and is
      * automatic. See the <a
diff --git a/source/base/hdf5.cc b/source/base/hdf5.cc
index 519166e14e..f637700c60 100644
--- a/source/base/hdf5.cc
+++ b/source/base/hdf5.cc
@@ -502,6 +502,194 @@ namespace HDF5
     return data;
   }
 
+  template <typename T>
+  std::vector<T>
+  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"));
+    std::shared_ptr<hid_t> t_type = internal::get_hdf5_datatype<T>();
+    hid_t                  plist;
+    hid_t                  memory_dataspace;
+    herr_t                 ret;
+
+    std::vector<hsize_t> data_dimensions{
+      static_cast<hsize_t>(coordinates.size() / _rank)};
+
+    std::vector<T> data(static_cast<unsigned int>(coordinates.size() / _rank));
+
+    memory_dataspace = H5Screate_simple(1, data_dimensions.data(), NULL);
+    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"));
+
+    if (mpi)
+      {
+        plist = H5Pcreate(H5P_DATASET_XFER);
+        Assert(plist >= 0, ExcMessage("Error at H5Pcreate"));
+        ret = H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
+        Assert(ret >= 0, ExcMessage("Error at H5Pset_dxpl_mpio"));
+      }
+    else
+      {
+        plist = H5P_DEFAULT;
+      }
+
+    ret = H5Dread(*hdf5_reference,
+                  *t_type,
+                  memory_dataspace,
+                  *dataspace,
+                  plist,
+                  data.data());
+    Assert(ret >= 0, ExcMessage("Error at H5Dread"));
+
+    if (mpi)
+      {
+        if (_check_io_mode)
+          {
+            ret = H5Pget_mpio_actual_io_mode(plist, &_io_mode);
+            Assert(ret >= 0, ExcMessage("Error at H5Pget_mpio_actual_io_mode"));
+            ret = H5Pget_mpio_no_collective_cause(plist,
+                                                  &_local_no_collective_cause,
+                                                  &_global_no_collective_cause);
+            Assert(ret >= 0,
+                   ExcMessage("Error at H5Pget_mpio_no_collective_cause"));
+          }
+        ret = H5Pclose(plist);
+        Assert(ret >= 0, ExcMessage("Error at H5Pclose"));
+      }
+
+    ret = H5Sclose(memory_dataspace);
+    Assert(ret >= 0, ExcMessage("Error at H5SClose"));
+
+    return data;
+  }
+
+  template <template <class...> class Container, typename T>
+  Container<T>
+  DataSet::read_hyperslab(const std::vector<hsize_t> &offset,
+                          const std::vector<hsize_t> &count)
+  {
+    std::shared_ptr<hid_t> t_type = internal::get_hdf5_datatype<T>();
+    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<T> data =
+      internal::initialize_container<Container, T>(data_dimensions);
+
+    memory_dataspace =
+      H5Screate_simple(data_dimensions.size(), data_dimensions.data(), NULL);
+    Assert(memory_dataspace >= 0, ExcMessage("Error at H5Screate_simple"));
+    ret = H5Sselect_hyperslab(
+      *dataspace, H5S_SELECT_SET, offset.data(), NULL, count.data(), NULL);
+    Assert(ret >= 0, ExcMessage("Error at H5Sselect_hyperslab"));
+
+    if (mpi)
+      {
+        plist = H5Pcreate(H5P_DATASET_XFER);
+        Assert(plist >= 0, ExcMessage("Error at H5Pcreate"));
+        ret = H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
+        Assert(ret >= 0, ExcMessage("Error at H5Pset_dxpl_mpio"));
+      }
+    else
+      {
+        plist = H5P_DEFAULT;
+      }
+
+    ret = H5Dread(*hdf5_reference,
+                  *t_type,
+                  memory_dataspace,
+                  *dataspace,
+                  plist,
+                  internal::get_container_pointer<Container, T>(data));
+    Assert(ret >= 0, ExcMessage("Error at H5Dread"));
+
+    if (mpi)
+      {
+        if (_check_io_mode)
+          {
+            ret = H5Pget_mpio_actual_io_mode(plist, &_io_mode);
+            Assert(ret >= 0, ExcMessage("Error at H5Pget_mpio_actual_io_mode"));
+            ret = H5Pget_mpio_no_collective_cause(plist,
+                                                  &_local_no_collective_cause,
+                                                  &_global_no_collective_cause);
+            Assert(ret >= 0,
+                   ExcMessage("Error at H5Pget_mpio_no_collective_cause"));
+          }
+        ret = H5Pclose(plist);
+        Assert(ret >= 0, ExcMessage("Error at H5Pclose"));
+      }
+
+    ret = H5Sclose(memory_dataspace);
+    Assert(ret >= 0, ExcMessage("Error at H5SClose"));
+
+    return data;
+  }
+
+  template <typename T>
+  void
+  DataSet::read_none()
+  {
+    std::shared_ptr<hid_t> t_type          = internal::get_hdf5_datatype<T>();
+    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(), NULL);
+    Assert(memory_dataspace >= 0, ExcMessage("Error at H5Screate_simple"));
+    ret = H5Sselect_none(*dataspace);
+    Assert(ret >= 0, ExcMessage("H5Sselect_none"));
+
+    if (mpi)
+      {
+        plist = H5Pcreate(H5P_DATASET_XFER);
+        Assert(plist >= 0, ExcMessage("Error at H5Pcreate"));
+        ret = H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
+        Assert(ret >= 0, ExcMessage("Error at H5Pset_dxpl_mpio"));
+      }
+    else
+      {
+        plist = H5P_DEFAULT;
+      }
+
+    // The pointer of data can safely be NULL, 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, NULL);
+    Assert(ret >= 0, ExcMessage("Error at H5Dread"));
+
+    if (mpi)
+      {
+        if (_check_io_mode)
+          {
+            ret = H5Pget_mpio_actual_io_mode(plist, &_io_mode);
+            Assert(ret >= 0, ExcMessage("Error at H5Pget_mpio_actual_io_mode"));
+            ret = H5Pget_mpio_no_collective_cause(plist,
+                                                  &_local_no_collective_cause,
+                                                  &_global_no_collective_cause);
+            Assert(ret >= 0,
+                   ExcMessage("Error at H5Pget_mpio_no_collective_cause"));
+          }
+        ret = H5Pclose(plist);
+        Assert(ret >= 0, ExcMessage("Error at H5Pclose"));
+      }
+
+    ret = H5Sclose(memory_dataspace);
+    Assert(ret >= 0, ExcMessage("Error at H5SClose"));
+  }
+
   template <template <class...> class Container, typename T>
   void
   DataSet::write(const Container<T> &data)
@@ -558,42 +746,49 @@ namespace HDF5
 
 
     memory_dataspace = H5Screate_simple(1, data_dimensions.data(), NULL);
-    H5Sselect_elements(*dataspace,
-                       H5S_SELECT_SET,
-                       data.size(),
-                       coordinates.data());
+    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"));
 
     if (mpi)
       {
         plist = H5Pcreate(H5P_DATASET_XFER);
-        H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
+        Assert(plist >= 0, ExcMessage("Error at H5Pcreate"));
+        ret = H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
+        Assert(ret >= 0, ExcMessage("Error at H5Pset_dxpl_mpio"));
       }
     else
       {
         plist = H5P_DEFAULT;
       }
 
-    H5Dwrite(*hdf5_reference,
-             *t_type,
-             memory_dataspace,
-             *dataspace,
-             plist,
-             data.data());
+    ret = H5Dwrite(*hdf5_reference,
+                   *t_type,
+                   memory_dataspace,
+                   *dataspace,
+                   plist,
+                   data.data());
+    Assert(ret >= 0, ExcMessage("Error at H5Dwrite"));
 
     if (mpi)
       {
         if (_check_io_mode)
           {
             ret = H5Pget_mpio_actual_io_mode(plist, &_io_mode);
-            Assert(ret >= 0, ExcInternalError());
+            Assert(ret >= 0, ExcMessage("Error at H5Pget_mpio_actual_io_mode"));
             ret = H5Pget_mpio_no_collective_cause(plist,
                                                   &_local_no_collective_cause,
                                                   &_global_no_collective_cause);
-            Assert(ret >= 0, ExcInternalError());
+            Assert(ret >= 0,
+                   ExcMessage("Error at H5Pget_mpio_no_collective_cause"));
           }
         H5Pclose(plist);
       }
-    H5Sclose(memory_dataspace);
+    ret = H5Sclose(memory_dataspace);
+    Assert(ret >= 0, ExcMessage("Error at H5SClose"));
   }
 
   template <template <class...> class Container, typename T>
@@ -1222,6 +1417,90 @@ namespace HDF5
   template FullMatrix<std::complex<double>>
   DataSet::read<FullMatrix, std::complex<double>>();
 
+  template std::vector<float>
+  DataSet::read_selection<float>(const std::vector<hsize_t> coordinates);
+  template std::vector<double>
+  DataSet::read_selection<double>(const std::vector<hsize_t> coordinates);
+  template std::vector<long double>
+  DataSet::read_selection<long double>(const std::vector<hsize_t> coordinates);
+  template std::vector<std::complex<float>>
+  DataSet::read_selection<std::complex<float>>(
+    const std::vector<hsize_t> coordinates);
+  template std::vector<std::complex<double>>
+  DataSet::read_selection<std::complex<double>>(
+    const std::vector<hsize_t> coordinates);
+  template std::vector<std::complex<long double>>
+  DataSet::read_selection<std::complex<long double>>(
+    const std::vector<hsize_t> coordinates);
+  template std::vector<int>
+  DataSet::read_selection<int>(const std::vector<hsize_t> coordinates);
+  template std::vector<unsigned int>
+  DataSet::read_selection<unsigned int>(const std::vector<hsize_t> coordinates);
+
+  template std::vector<float>
+  DataSet::read_hyperslab<std::vector, float>(
+    const std::vector<hsize_t> &offset,
+    const std::vector<hsize_t> &count);
+  template std::vector<double>
+  DataSet::read_hyperslab<std::vector, double>(
+    const std::vector<hsize_t> &offset,
+    const std::vector<hsize_t> &count);
+  template std::vector<long double>
+  DataSet::read_hyperslab<std::vector, long double>(
+    const std::vector<hsize_t> &offset,
+    const std::vector<hsize_t> &count);
+  template std::vector<std::complex<float>>
+  DataSet::read_hyperslab<std::vector, std::complex<float>>(
+    const std::vector<hsize_t> &offset,
+    const std::vector<hsize_t> &count);
+  template std::vector<std::complex<double>>
+  DataSet::read_hyperslab<std::vector, std::complex<double>>(
+    const std::vector<hsize_t> &offset,
+    const std::vector<hsize_t> &count);
+  template std::vector<std::complex<long double>>
+  DataSet::read_hyperslab<std::vector, std::complex<long double>>(
+    const std::vector<hsize_t> &offset,
+    const std::vector<hsize_t> &count);
+  template std::vector<int>
+  DataSet::read_hyperslab<std::vector, int>(const std::vector<hsize_t> &offset,
+                                            const std::vector<hsize_t> &count);
+  template std::vector<unsigned int>
+  DataSet::read_hyperslab<std::vector, unsigned int>(
+    const std::vector<hsize_t> &offset,
+    const std::vector<hsize_t> &count);
+  template FullMatrix<float>
+  DataSet::read_hyperslab<FullMatrix, float>(const std::vector<hsize_t> &offset,
+                                             const std::vector<hsize_t> &count);
+  template FullMatrix<double>
+  DataSet::read_hyperslab<FullMatrix, double>(
+    const std::vector<hsize_t> &offset,
+    const std::vector<hsize_t> &count);
+  template FullMatrix<std::complex<float>>
+  DataSet::read_hyperslab<FullMatrix, std::complex<float>>(
+    const std::vector<hsize_t> &offset,
+    const std::vector<hsize_t> &count);
+  template FullMatrix<std::complex<double>>
+  DataSet::read_hyperslab<FullMatrix, std::complex<double>>(
+    const std::vector<hsize_t> &offset,
+    const std::vector<hsize_t> &count);
+
+  template void
+  DataSet::read_none<float>();
+  template void
+  DataSet::read_none<double>();
+  template void
+  DataSet::read_none<long double>();
+  template void
+  DataSet::read_none<std::complex<float>>();
+  template void
+  DataSet::read_none<std::complex<double>>();
+  template void
+  DataSet::read_none<std::complex<long double>>();
+  template void
+  DataSet::read_none<int>();
+  template void
+  DataSet::read_none<unsigned int>();
+
   template void
   DataSet::write<std::vector, float>(const std::vector<float> &data);
   template void
-- 
2.39.5