#ifdef DEAL_II_WITH_HDF5
+# include <deal.II/base/array_view.h>
+
# include <deal.II/lac/full_matrix.h>
# include <hdf5.h>
-# include <vector>
+# include <numeric>
DEAL_II_NAMESPACE_OPEN
const bool mpi,
const MPI_Comm mpi_communicator);
};
+
+
+
+ // definitions
+
+ namespace internal
+ {
+ // This function gives the HDF5 datatype corresponding to the C++ type. In
+ // the case of std::complex types the HDF5 handlers are automatically freed
+ // using the destructor of std::shared_ptr.
+ // std::shared_ptr is used instead of std::unique_ptr because the destructor
+ // of std::shared_ptr doesn't have to be defined in the template argument.
+ // In the other hand, the destructor of std::unique has to be defined in the
+ // template argument. Native types such as H5T_NATIVE_DOUBLE do not
+ // require a destructor, but compound types such as std::complex<double>
+ // require a destructor to free the HDF5 resources.
+ template <typename number>
+ std::shared_ptr<hid_t>
+ get_hdf5_datatype()
+ {
+ std::shared_ptr<hid_t> t_type;
+ if (std::is_same<number, float>::value)
+ {
+ return std::make_shared<hid_t>(H5T_NATIVE_FLOAT);
+ }
+ else if (std::is_same<number, double>::value)
+ {
+ return std::make_shared<hid_t>(H5T_NATIVE_DOUBLE);
+ }
+ else if (std::is_same<number, int>::value)
+ {
+ return std::make_shared<hid_t>(H5T_NATIVE_INT);
+ }
+ else if (std::is_same<number, unsigned int>::value)
+ {
+ return std::make_shared<hid_t>(H5T_NATIVE_UINT);
+ }
+ else if (std::is_same<number, std::complex<float>>::value)
+ {
+ 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<number, std::complex<double>>::value)
+ {
+ 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
+ Assert(false, ExcInternalError());
+ return t_type;
+ }
+
+
+
+ // Several HDF5 functions such as H5Screate_simple() require a
+ // one-dimensional array that specifies the size of each dimension of the
+ // container, see:
+ // https://support.hdfgroup.org/HDF5/doc1.8/RM/RM_H5S.html#Dataspace-CreateSimple
+ // The function get_container_dimensions returns a vector with the
+ // dimensions of the container.
+ // For a std::vector the function returns std::vector<hsize_t>{vector_size}
+ // For a Vector the function returns std::vector<hsize_t>{vector_size}
+ // For a FullMatrix the function returns std::vector<hsize_t>{rows, columns}
+ 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;
+ }
+
+
+
+ // This function returns the total size of the container.
+ // For a std::vector the function returns int(vector_size)
+ // For a Vector the function returns int(vector_size)
+ // For a FullMatrix the function returns int(rows*columns)
+ 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());
+ }
+
+
+
+ // This function initializes and returns a container of type std::vector,
+ // Vector or FullMatrix. The function does not set the values of the
+ // elements of the container. The container can store data of a HDF5 dataset
+ // or a HDF5 selection. The dimensions parameter holds the dimensions of the
+ // HDF5 dataset or selection.
+ //
+ // In the case of a std::vector, the size of the vector will be the total
+ // size given by dimensions. For example in the case of a dataset of rank 3,
+ // the dimensions are std::vector<hsize_t>{dim_0,dim_1,dim_2}. The size of
+ // the returned std::vector will be dim_0*dim_1*dim_2
+ //
+ // In the case of a dealii::Vector, the size of the returned dealii::Vector
+ // will be as well dim_0*dim_1*dim_2
+ //
+ // A FullMatrix can store only data of HDF5 datasets with rank 2. The size
+ // of the FullMatrix will be FullMatrix(dim_0,dim_2)
+ template <typename Container>
+ typename std::enable_if<
+ std::is_same<Container,
+ std::vector<typename Container::value_type>>::value,
+ Container>::type
+ initialize_container(const std::vector<hsize_t> &dimensions)
+ {
+ return Container(std::accumulate(
+ dimensions.begin(), dimensions.end(), 1, std::multiplies<int>()));
+ }
+
+
+
+ template <typename Container>
+ typename std::enable_if<
+ std::is_same<Container, Vector<typename Container::value_type>>::value,
+ Container>::type
+ initialize_container(const std::vector<hsize_t> &dimensions)
+ {
+ return Container(std::accumulate(
+ dimensions.begin(), dimensions.end(), 1, std::multiplies<int>()));
+ }
+
+
+
+ template <typename Container>
+ typename std::enable_if<
+ std::is_same<Container,
+ FullMatrix<typename Container::value_type>>::value,
+ Container>::type
+ 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]);
+ }
+
+
+ // This helper function sets the property list of the read and write
+ // operations of DataSet. A property list has to be created for the MPI
+ // driver. For the serial driver the default H5P_DEFAULT can be used. In
+ // addition H5Pset_dxpl_mpio is used to set the MPI mode to collective.
+ 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;
+ }
+
+
+ // This helper function release the property list handler of the read and
+ // write operations of DataSet. For the serial there is no need to release
+ // the property list handler because H5P_DEFAULT has been used. If
+ // query_io_mode is True then H5Pget_mpio_actual_io_mode and
+ // H5Pget_mpio_no_collective_cause are used to check if the operation has
+ // been collective.
+ void
+ release_plist(hid_t & plist,
+ H5D_mpio_actual_io_mode_t &io_mode,
+ uint32_t & local_no_collective_cause,
+ 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;
+ }
+
+
+ // Convert a HDF5 no_collective_cause code to a human readable string
+ std::string
+ no_collective_cause_to_string(const uint32_t no_collective_cause)
+ {
+ std::string message;
+
+ auto append_to_message = [&message](const char *p) {
+ if (message.length() > 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 <>
+ bool
+ HDF5Object::get_attribute(const std::string &attr_name) const
+ {
+ // The enum field generated by h5py can be casted to int
+ int int_value;
+ hid_t attr;
+ herr_t ret;
+
+ attr = H5Aopen(*hdf5_reference, attr_name.data(), H5P_DEFAULT);
+ Assert(attr >= 0, ExcMessage("Error at H5Aopen"));
+ ret = H5Aread(attr, H5T_NATIVE_INT, &int_value);
+ (void)ret;
+ Assert(ret >= 0, ExcMessage("Error at H5Aread"));
+ ret = H5Aclose(attr);
+ (void)ret;
+ Assert(ret >= 0, ExcMessage("Error at H5Aclose"));
+
+ return (int_value != 0);
+ }
+
+
+
+ template <>
+ 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 <>
+ 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
#ifdef DEAL_II_WITH_HDF5
-# include <deal.II/base/array_view.h>
# include <deal.II/base/hdf5.h>
-# include <deal.II/lac/full_matrix.h>
-# include <deal.II/lac/vector.h>
-
# include <hdf5.h>
-# include <memory>
-# include <numeric>
-# include <vector>
-
DEAL_II_NAMESPACE_OPEN
namespace HDF5
{
- namespace internal
- {
- // This function gives the HDF5 datatype corresponding to the C++ type. In
- // the case of std::complex types the HDF5 handlers are automatically freed
- // using the destructor of std::shared_ptr.
- // std::shared_ptr is used instead of std::unique_ptr because the destructor
- // of std::shared_ptr doesn't have to be defined in the template argument.
- // In the other hand, the destructor of std::unique has to be defined in the
- // template argument. Native types such as H5T_NATIVE_DOUBLE do not
- // require a destructor, but compound types such as std::complex<double>
- // require a destructor to free the HDF5 resources.
- template <typename number>
- std::shared_ptr<hid_t>
- get_hdf5_datatype()
- {
- std::shared_ptr<hid_t> t_type;
- if (std::is_same<number, float>::value)
- {
- return std::make_shared<hid_t>(H5T_NATIVE_FLOAT);
- }
- else if (std::is_same<number, double>::value)
- {
- return std::make_shared<hid_t>(H5T_NATIVE_DOUBLE);
- }
- else if (std::is_same<number, int>::value)
- {
- return std::make_shared<hid_t>(H5T_NATIVE_INT);
- }
- else if (std::is_same<number, unsigned int>::value)
- {
- return std::make_shared<hid_t>(H5T_NATIVE_UINT);
- }
- else if (std::is_same<number, std::complex<float>>::value)
- {
- 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<number, std::complex<double>>::value)
- {
- 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
- Assert(false, ExcInternalError());
- return t_type;
- }
-
-
-
- // Several HDF5 functions such as H5Screate_simple() require a
- // one-dimensional array that specifies the size of each dimension of the
- // container, see:
- // https://support.hdfgroup.org/HDF5/doc1.8/RM/RM_H5S.html#Dataspace-CreateSimple
- // The function get_container_dimensions returns a vector with the
- // dimensions of the container.
- // For a std::vector the function returns std::vector<hsize_t>{vector_size}
- // For a Vector the function returns std::vector<hsize_t>{vector_size}
- // For a FullMatrix the function returns std::vector<hsize_t>{rows, columns}
- 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;
- }
-
-
-
- // This function returns the total size of the container.
- // For a std::vector the function returns int(vector_size)
- // For a Vector the function returns int(vector_size)
- // For a FullMatrix the function returns int(rows*columns)
- 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());
- }
-
-
-
- // This function initializes and returns a container of type std::vector,
- // Vector or FullMatrix. The function does not set the values of the
- // elements of the container. The container can store data of a HDF5 dataset
- // or a HDF5 selection. The dimensions parameter holds the dimensions of the
- // HDF5 dataset or selection.
- //
- // In the case of a std::vector, the size of the vector will be the total
- // size given by dimensions. For example in the case of a dataset of rank 3,
- // the dimensions are std::vector<hsize_t>{dim_0,dim_1,dim_2}. The size of
- // the returned std::vector will be dim_0*dim_1*dim_2
- //
- // In the case of a dealii::Vector, the size of the returned dealii::Vector
- // will be as well dim_0*dim_1*dim_2
- //
- // A FullMatrix can store only data of HDF5 datasets with rank 2. The size
- // of the FullMatrix will be FullMatrix(dim_0,dim_2)
- template <typename Container>
- typename std::enable_if<
- std::is_same<Container,
- std::vector<typename Container::value_type>>::value,
- Container>::type
- initialize_container(const std::vector<hsize_t> &dimensions)
- {
- return Container(std::accumulate(
- dimensions.begin(), dimensions.end(), 1, std::multiplies<int>()));
- }
-
-
-
- template <typename Container>
- typename std::enable_if<
- std::is_same<Container, Vector<typename Container::value_type>>::value,
- Container>::type
- initialize_container(const std::vector<hsize_t> &dimensions)
- {
- return Container(std::accumulate(
- dimensions.begin(), dimensions.end(), 1, std::multiplies<int>()));
- }
-
-
-
- template <typename Container>
- typename std::enable_if<
- std::is_same<Container,
- FullMatrix<typename Container::value_type>>::value,
- Container>::type
- 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]);
- }
-
-
- // This helper function sets the property list of the read and write
- // operations of DataSet. A property list has to be created for the MPI
- // driver. For the serial driver the default H5P_DEFAULT can be used. In
- // addition H5Pset_dxpl_mpio is used to set the MPI mode to collective.
- 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;
- }
-
-
- // This helper function release the property list handler of the read and
- // write operations of DataSet. For the serial there is no need to release
- // the property list handler because H5P_DEFAULT has been used. If
- // query_io_mode is True then H5Pget_mpio_actual_io_mode and
- // H5Pget_mpio_no_collective_cause are used to check if the operation has
- // been collective.
- void
- release_plist(hid_t & plist,
- H5D_mpio_actual_io_mode_t &io_mode,
- uint32_t & local_no_collective_cause,
- 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;
- }
-
-
- // Convert a HDF5 no_collective_cause code to a human readable string
- std::string
- no_collective_cause_to_string(const uint32_t no_collective_cause)
- {
- std::string message;
-
- auto append_to_message = [&message](const char *p) {
- if (message.length() > 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
-
-
-
HDF5Object::HDF5Object(const std::string &name, const bool mpi)
: name(name)
, mpi(mpi)
- 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 <>
- bool
- HDF5Object::get_attribute(const std::string &attr_name) const
- {
- // The enum field generated by h5py can be casted to int
- int int_value;
- hid_t attr;
- herr_t ret;
-
- attr = H5Aopen(*hdf5_reference, attr_name.data(), H5P_DEFAULT);
- Assert(attr >= 0, ExcMessage("Error at H5Aopen"));
- ret = H5Aread(attr, H5T_NATIVE_INT, &int_value);
- (void)ret;
- Assert(ret >= 0, ExcMessage("Error at H5Aread"));
- ret = H5Aclose(attr);
- (void)ret;
- Assert(ret >= 0, ExcMessage("Error at H5Aclose"));
-
- return (int_value != 0);
- }
-
-
-
- template <>
- 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 <>
- 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;
- }
-
-
-
std::string
HDF5Object::get_name() const
{
- 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;
- }
-
-
-
void
DataSet::set_query_io_mode(const bool new_query_io_mode)
{
- 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);
- }
-
-
-
File::File(const std::string &name, const FileAccessMode mode)
: File(name,
mode,
(void)ret;
(void)mpi_communicator;
}
-}
+} // namespace HDF5
DEAL_II_NAMESPACE_CLOSE
projects / dealii.git / commitdiff