diff --git a/ops/c/include/ops_hdf5_common.h b/ops/c/include/ops_hdf5_common.h index 1d38fa30d4..a4e15da9fd 100644 --- a/ops/c/include/ops_hdf5_common.h +++ b/ops/c/include/ops_hdf5_common.h @@ -49,6 +49,18 @@ hid_t h5_type(const char *type); +hid_t h5_storage_type(const char *data_type, REAL_PRECISION real_precision); + +/* H5Dwrite after converting buf to the dataset type when it differs from + * mem_type. Needed so parallel HDF5 can use collective MPI-IO. */ +herr_t H5Dwrite_matching_types(hid_t dataset_id, hid_t mem_type, hid_t mem_space, + hid_t file_space, hid_t xfer_plist, + const void *buf, size_t nelem); + +herr_t H5Dread_matching_types(hid_t dataset_id, hid_t mem_type, hid_t mem_space, + hid_t file_space, hid_t xfer_plist, void *buf, + size_t nelem); + void split_h5_name(const char *data_name, std::vector &h5_name_list); diff --git a/ops/c/src/externlib/ops_hdf5.cpp b/ops/c/src/externlib/ops_hdf5.cpp index 79215f1100..9990bf4622 100755 --- a/ops/c/src/externlib/ops_hdf5.cpp +++ b/ops/c/src/externlib/ops_hdf5.cpp @@ -608,7 +608,8 @@ void ops_fetch_dat_hdf5_file(ops_dat dat, char const *file_name) { hid_t dset_id = H5Dopen(group_id, dat->name, H5P_DEFAULT); hid_t dataspace = H5Dget_space(dset_id); hid_t datatype = h5_type(dat->type); - H5Dwrite(dset_id, datatype, H5S_ALL, dataspace, H5P_DEFAULT, data); + H5Dwrite_matching_types(dset_id, datatype, H5S_ALL, dataspace, H5P_DEFAULT, + data, (size_t)t_size * (size_t)dat->dim); } free(data); H5Gclose(group_id); @@ -1226,7 +1227,8 @@ void ops_get_const_hdf5(char const *name, int dim, char const *type, size_t bytesize = H5Tget_size(datatype); char *data = (char *)xmalloc(bytesize * const_dim); - H5Dread(dset_id, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, data); + H5Dread_matching_types(dset_id, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, data, + (size_t)const_dim); memcpy((void *)const_data, (void *)data, bytesize * const_dim); ops_free(typ); @@ -1307,10 +1309,10 @@ void ops_write_const_hdf5(char const *name, int dim, char const *type, OPS_instance::getOPSInstance()->ostream() << " overwriting" << "\n"; dataspace = H5Dget_space(dset_id); - hid_t datatype = h5_type(typ); + hid_t datatype = h5_type(type); // write to the existing dataset with default properties - H5Dwrite(dset_id, datatype, H5S_ALL, dataspace, H5P_DEFAULT, - const_data); + H5Dwrite_matching_types(dset_id, datatype, H5S_ALL, dataspace, H5P_DEFAULT, + const_data, (size_t)dim); ops_free(typ); H5Dclose(dset_id); H5Fclose(file_id); @@ -1329,8 +1331,8 @@ void ops_write_const_hdf5(char const *name, int dim, char const *type, dset_id = H5Dcreate(file_id, name, datatype, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); // write data - H5Dwrite(dset_id, datatype, H5S_ALL, dataspace, H5P_DEFAULT, - const_data); + H5Dwrite_matching_types(dset_id, datatype, H5S_ALL, dataspace, H5P_DEFAULT, + const_data, (size_t)dim); H5Dclose(dset_id); H5Sclose(dataspace); @@ -1488,8 +1490,11 @@ void write_buf_hdf5(char const *file_name, const char *data_name, H5_dataset_space(file_id, dims, size_f, h5_name_list, dat->type, groupid_list, dataset_id, file_space); - H5Dwrite(dataset_id, h5_type(dat->type), H5S_ALL, file_space, H5P_DEFAULT, - buf); + size_t nelem = 1; + for (int d = 0; d < dims; d++) + nelem *= (size_t)size[d]; + H5Dwrite_matching_types(dataset_id, h5_type(dat->type), H5S_ALL, file_space, + H5P_DEFAULT, buf, nelem); H5Sclose(file_space); H5Dclose(dataset_id); for (int grp = groupid_list.size() - 1; grp >= 0; grp--) { @@ -1517,8 +1522,11 @@ void write_buf_hdf5(char const *file_name, const char *data_name, H5_dataset_space(file_id, dims, size_f, h5_name_list, dat->type, float_precision, groupid_list, dataset_id, file_space); - H5Dwrite(dataset_id, h5_type(dat->type), H5S_ALL, file_space, H5P_DEFAULT, - buf); + size_t nelem = 1; + for (int d = 0; d < dims; d++) + nelem *= (size_t)size[d]; + H5Dwrite_matching_types(dataset_id, h5_type(dat->type), H5S_ALL, file_space, + H5P_DEFAULT, buf, nelem); H5Sclose(file_space); H5Dclose(dataset_id); for (int grp = groupid_list.size() - 1; grp >= 0; grp--) { diff --git a/ops/c/src/externlib/ops_hdf5_common.cpp b/ops/c/src/externlib/ops_hdf5_common.cpp index db5edf5b2b..9ba0ac390e 100644 --- a/ops/c/src/externlib/ops_hdf5_common.cpp +++ b/ops/c/src/externlib/ops_hdf5_common.cpp @@ -46,6 +46,8 @@ // hdf5 header #include "ops_hdf5_common.h" #include +#include +#include int half_type_init = 0; hid_t H5T_IEEE_FP16; @@ -161,6 +163,136 @@ hid_t h5_type(const char *type) { return h5t; } +hid_t h5_storage_type(const char *data_type, REAL_PRECISION real_precision) { + hid_t type = h5_type(data_type); + if (real_precision == REAL_PRECISION::Undefined) + return type; + + if (type == H5T_NATIVE_DOUBLE) { + if (real_precision == REAL_PRECISION::Single) + return H5T_NATIVE_FLOAT; + if (real_precision == REAL_PRECISION::Half) + return H5T_IEEE_FP16; + return H5T_NATIVE_DOUBLE; + } + if (type == H5T_NATIVE_FLOAT) { + if (real_precision == REAL_PRECISION::Half) + return H5T_IEEE_FP16; + return H5T_NATIVE_FLOAT; + } + return type; +} + +/* Convert buf to the dataset's type before H5Dwrite. Parallel HDF5 disables + * collective MPI-IO when H5Dwrite itself has to convert + * (H5D_MPIO_DATATYPE_CONVERSION). */ +herr_t H5Dwrite_matching_types(hid_t dataset_id, hid_t mem_type, hid_t mem_space, + hid_t file_space, hid_t xfer_plist, + const void *buf, size_t nelem) { + hid_t dset_type = H5Dget_type(dataset_id); + if (dset_type < 0) { + OPSException ex(OPS_HDF5_ERROR); + ex << "Error: H5Dget_type failed"; + throw ex; + } + + char *converted = nullptr; + const void *write_buf = buf; + hid_t write_type = mem_type; + htri_t equal = H5Tequal(mem_type, dset_type); + if (equal < 0) { + H5Tclose(dset_type); + OPSException ex(OPS_HDF5_ERROR); + ex << "Error: H5Tequal failed comparing memory and dataset types"; + throw ex; + } + if (equal == 0 && nelem > 0 && buf != nullptr) { + size_t src_sz = H5Tget_size(mem_type); + size_t dst_sz = H5Tget_size(dset_type); + size_t nbytes = nelem * (src_sz > dst_sz ? src_sz : dst_sz); + converted = (char *)malloc(nbytes); + if (converted == nullptr) { + H5Tclose(dset_type); + OPSException ex(OPS_HDF5_ERROR); + ex << "Error: out of memory converting HDF5 write buffer"; + throw ex; + } + memcpy(converted, buf, nelem * src_sz); + if (H5Tconvert(mem_type, dset_type, nelem, converted, NULL, H5P_DEFAULT) < + 0) { + free(converted); + H5Tclose(dset_type); + OPSException ex(OPS_HDF5_ERROR); + ex << "Error: H5Tconvert failed for HDF5 write buffer"; + throw ex; + } + write_buf = converted; + write_type = dset_type; + } + + herr_t err = + H5Dwrite(dataset_id, write_type, mem_space, file_space, xfer_plist, + write_buf); + free(converted); + H5Tclose(dset_type); + return err; +} + +/* Read with the dataset's type, then convert into mem_type. Same collective + * MPI-IO constraint as the write path. */ +herr_t H5Dread_matching_types(hid_t dataset_id, hid_t mem_type, hid_t mem_space, + hid_t file_space, hid_t xfer_plist, void *buf, + size_t nelem) { + hid_t dset_type = H5Dget_type(dataset_id); + if (dset_type < 0) { + OPSException ex(OPS_HDF5_ERROR); + ex << "Error: H5Dget_type failed"; + throw ex; + } + + htri_t equal = H5Tequal(mem_type, dset_type); + if (equal < 0) { + H5Tclose(dset_type); + OPSException ex(OPS_HDF5_ERROR); + ex << "Error: H5Tequal failed comparing memory and dataset types"; + throw ex; + } + if (equal > 0) { + herr_t err = + H5Dread(dataset_id, mem_type, mem_space, file_space, xfer_plist, buf); + H5Tclose(dset_type); + return err; + } + + size_t src_sz = H5Tget_size(dset_type); + size_t dst_sz = H5Tget_size(mem_type); + size_t nbytes = nelem * (src_sz > dst_sz ? src_sz : dst_sz); + if (nbytes == 0) + nbytes = 1; + char *tmp = (char *)malloc(nbytes); + if (tmp == nullptr) { + H5Tclose(dset_type); + OPSException ex(OPS_HDF5_ERROR); + ex << "Error: out of memory converting HDF5 read buffer"; + throw ex; + } + herr_t err = + H5Dread(dataset_id, dset_type, mem_space, file_space, xfer_plist, tmp); + if (err >= 0 && nelem > 0 && buf != nullptr) { + if (H5Tconvert(dset_type, mem_type, nelem, tmp, NULL, H5P_DEFAULT) < 0) { + free(tmp); + H5Tclose(dset_type); + OPSException ex(OPS_HDF5_ERROR); + ex << "Error: H5Tconvert failed for HDF5 read buffer"; + throw ex; + } + memcpy(buf, tmp, nelem * dst_sz); + } + free(tmp); + H5Tclose(dset_type); + return err; +} + void split_h5_name(const char *data_name, std::vector &h5_name_list) { std::stringstream name_stream(data_name); @@ -197,32 +329,7 @@ void H5_dataset_space(const hid_t file_id, const int data_dims, hid_t type = h5_type(data_type); if (type == H5T_NATIVE_FLOAT || type == H5T_NATIVE_DOUBLE || type == H5T_IEEE_FP16) { - hid_t real_type; - if (type == H5T_NATIVE_DOUBLE) { - if (real_precision == REAL_PRECISION::Double) { - real_type = H5T_NATIVE_DOUBLE; - } - if (real_precision == REAL_PRECISION::Single) { - real_type = H5T_NATIVE_FLOAT; - } - if (real_precision == REAL_PRECISION::Half) { - real_type = H5T_IEEE_FP16; - } - } - - if (type == H5T_NATIVE_FLOAT) { - if (real_precision == REAL_PRECISION::Single) { - real_type = H5T_NATIVE_FLOAT; - } - if (real_precision == REAL_PRECISION::Half) { - real_type = H5T_IEEE_FP16; - } - } - - if (type == H5T_IEEE_FP16) { - real_type = H5T_IEEE_FP16; - } - + hid_t real_type = h5_storage_type(data_type, real_precision); dataset_id = H5Dcreate(parent_group, data_name, real_type, file_space, H5P_DEFAULT, data_plist_id, H5P_DEFAULT); } else { diff --git a/ops/c/src/mpi/ops_mpi_hdf5.cpp b/ops/c/src/mpi/ops_mpi_hdf5.cpp index 29fe317115..8e48942c6e 100755 --- a/ops/c/src/mpi/ops_mpi_hdf5.cpp +++ b/ops/c/src/mpi/ops_mpi_hdf5.cpp @@ -37,6 +37,7 @@ */ #include +#include #include #include #include @@ -62,6 +63,23 @@ static char *copy_str(char const *src) { return strncpy(dest, src, len); } +/* Report whether HDF5 honoured the collective MPI-IO request. Datatype + * conversion (mem type != file type) is a documented reason for falling + * back to independent I/O (H5D_MPIO_DATATYPE_CONVERSION = 0x02). */ +static void report_mpio_io_mode(hid_t xfer_plist, const char *tag) { + if (OPS_instance::getOPSInstance()->OPS_diags <= 3) + return; + H5D_mpio_actual_io_mode_t mode = H5D_MPIO_NO_COLLECTIVE; + herr_t e1 = H5Pget_mpio_actual_io_mode(xfer_plist, &mode); + uint32_t local_cause = 0, global_cause = 0; + herr_t e2 = + H5Pget_mpio_no_collective_cause(xfer_plist, &local_cause, &global_cause); + ops_printf("HDF5 %s actual_io_mode=%d (4=contiguous collective, 0=none) " + "no_collective_cause local=0x%x global=0x%x " + "(0x2=DATATYPE_CONVERSION) get_status=%d/%d\n", + tag, (int)mode, local_cause, global_cause, (int)e1, (int)e2); +} + // // MPI Communicator for parallel I/O // @@ -860,7 +878,11 @@ void ops_fetch_dat_hdf5_file(ops_dat dat, char const *file_name) { // write data hid_t datatype = h5_type(dat->type); - H5Dwrite(dset_id, datatype, memspace, filespace, plist_id, data); + size_t nelem = 1; + for (int d = 0; d < block->dims; d++) + nelem *= (size_t)SIZE[d]; + H5Dwrite_matching_types(dset_id, datatype, memspace, filespace, plist_id, + data, nelem); MPI_Barrier(OPS_MPI_HDF5_BLOCK_WORLD); free(data); @@ -1567,7 +1589,11 @@ void ops_read_dat_hdf5(ops_dat dat) { // read data hid_t datatype = h5_type(dat->type); - H5Dread(dset_id, datatype, memspace, filespace, plist_id, data); + size_t nelem = 1; + for (int d = 0; d < block->dims; d++) + nelem *= (size_t)size[d]; + H5Dread_matching_types(dset_id, datatype, memspace, filespace, plist_id, + data, nelem); ops_dat_set_data(dat, 0, data); @@ -1789,7 +1815,6 @@ void ops_get_const_hdf5(char const *name, int dim, char const *type, "type of constant %s in file %s and requested type %s do not match, " "performing automatic type conversion\n", typ, file_name, type); - typ = type; } H5Dclose(dset_id); @@ -1797,10 +1822,11 @@ void ops_get_const_hdf5(char const *name, int dim, char const *type, dset_id = H5Dopen(file_id, name, H5P_DEFAULT); char *data = nullptr; - hid_t datatype = h5_type(typ); + hid_t datatype = h5_type(type); size_t type_size = H5Tget_size(datatype); data = (char *)xmalloc(type_size * const_dim); - H5Dread(dset_id, datatype, H5S_ALL, H5S_ALL, plist_id, data); + H5Dread_matching_types(dset_id, datatype, H5S_ALL, H5S_ALL, plist_id, data, + (size_t)const_dim); memcpy((void *)const_data, (void *)data, type_size * const_dim); free(data); @@ -1895,7 +1921,6 @@ void ops_write_const_hdf5(char const *name, int dim, char const *type, ops_printf("type of constant %s in file %s and requested type %s do " "not match, performing automatic type conversion\n", typ, file_name, type); - typ = type; } // existing const attributes matches with const to be written .. overwriting @@ -1904,8 +1929,8 @@ void ops_write_const_hdf5(char const *name, int dim, char const *type, dataspace = H5Dget_space(dset_id); hid_t datatype = h5_type(type); - H5Dwrite(dset_id, datatype, H5S_ALL, dataspace, plist_id, - const_data); + H5Dwrite_matching_types(dset_id, datatype, H5S_ALL, dataspace, plist_id, + const_data, (size_t)dim); H5Pclose(plist_id); H5Sclose(dataspace); @@ -1932,8 +1957,8 @@ void ops_write_const_hdf5(char const *name, int dim, char const *type, dset_id = H5Dcreate(file_id, name, datatype, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); // write data - H5Dwrite(dset_id, datatype, H5S_ALL, dataspace, plist_id, - const_data); + H5Dwrite_matching_types(dset_id, datatype, H5S_ALL, dataspace, plist_id, + const_data, (size_t)dim); H5Dclose(dset_id); H5Pclose(plist_id); @@ -2127,8 +2152,12 @@ void write_plane_buf_hdf5(const char *file_name, const char *data_name, hid_t xfer_data_plist_id = H5Pcreate(H5P_DATASET_XFER); H5Pset_dxpl_mpio(xfer_data_plist_id, H5FD_MPIO_COLLECTIVE); - H5Dwrite(dataset_id, h5_type(dat->type), memspace, file_space, - xfer_data_plist_id, buf); + size_t nelem = 1; + for (int d = 0; d < data_dims; d++) + nelem *= (size_t)local_data_size_c[d]; + H5Dwrite_matching_types(dataset_id, h5_type(dat->type), memspace, file_space, + xfer_data_plist_id, buf, nelem); + report_mpio_io_mode(xfer_data_plist_id, "plane write"); H5Pclose(xfer_data_plist_id); H5Sclose(file_space); H5Sclose(memspace); @@ -2250,8 +2279,12 @@ void write_slab_buf_hdf5(const char *file_name, const char *data_name, hid_t xfer_data_plist_id = H5Pcreate(H5P_DATASET_XFER); H5Pset_dxpl_mpio(xfer_data_plist_id, H5FD_MPIO_COLLECTIVE); - H5Dwrite(dataset_id, h5_type(dat->type), memspace, file_space, - xfer_data_plist_id, buf); + size_t nelem = 1; + for (int d = 0; d < space_dim; d++) + nelem *= (size_t)local_data_size_c[d]; + H5Dwrite_matching_types(dataset_id, h5_type(dat->type), memspace, file_space, + xfer_data_plist_id, buf, nelem); + report_mpio_io_mode(xfer_data_plist_id, "slab write"); H5Pclose(xfer_data_plist_id); H5Sclose(file_space); H5Sclose(memspace);