From 9cf8d26110708a72c875f59baf4c70cab1e8c6ca Mon Sep 17 00:00:00 2001 From: q-posev Date: Tue, 2 Mar 2021 10:55:13 +0100 Subject: [PATCH] hdf5 backend working --- src/trio.c | 14 +- src/trio.org | 390 +++++++++++++++++++++++++++++++++++++----------- src/trio_hdf5.c | 353 +++++++++++++++++++++++++++++++++---------- src/trio_hdf5.h | 6 +- 4 files changed, 591 insertions(+), 172 deletions(-) diff --git a/src/trio.c b/src/trio.c index 9d8d0aa..6ca20f8 100644 --- a/src/trio.c +++ b/src/trio.c @@ -20,7 +20,7 @@ */ trio_t* trio_create(const char* file_name, back_end_t back_end) { - + /* Check that file name is not NULL or empty */ assert (file_name != NULL); assert (file_name[0] != '\0'); @@ -129,11 +129,11 @@ trio_exit_code trio_read_nucleus_num(trio_t* file, uint64_t* num) { case TRIO_TEXT: return trio_text_read_nucleus_num(file, num); break; -/* + case TRIO_HDF5: return trio_hdf5_read_nucleus_num(file, num); break; - +/* case TRIO_JSON: return trio_json_read_nucleus_num(file, num); break; @@ -151,11 +151,11 @@ trio_exit_code trio_write_nucleus_num(trio_t* file, uint64_t num) { case TRIO_TEXT: return trio_text_write_nucleus_num(file, num); break; -/* + case TRIO_HDF5: return trio_hdf5_write_nucleus_num(file, num); break; - +/* case TRIO_JSON: return trio_json_write_nucleus_num(file, num); break; @@ -195,11 +195,11 @@ trio_exit_code trio_write_nucleus_coord(trio_t* file, double* coord) { case TRIO_TEXT: return trio_text_write_nucleus_coord(file, coord); break; -/* + case TRIO_HDF5: return trio_hdf5_write_nucleus_coord(file, coord); break; - +/* case TRIO_JSON: return trio_json_write_nucleus_coord(file, coord); break; diff --git a/src/trio.org b/src/trio.org index e1753c0..378265e 100644 --- a/src/trio.org +++ b/src/trio.org @@ -198,7 +198,7 @@ trio_t* trio_create(const char* file_name, back_end_t back_end); #+begin_src c :tangle trio.c trio_t* trio_create(const char* file_name, back_end_t back_end) { - + /* Check that file name is not NULL or empty */ assert (file_name != NULL); assert (file_name[0] != '\0'); @@ -251,7 +251,7 @@ trio_t* trio_create(const char* file_name, back_end_t back_end) { case TRIO_JSON: rc = trio_json_init(result); break; -*/ +,*/ default: assert (1 == 0); /* Impossible case */ } @@ -324,11 +324,11 @@ trio_exit_code trio_read_nucleus_num(trio_t* file, uint64_t* num) { case TRIO_TEXT: return trio_text_read_nucleus_num(file, num); break; -/* + case TRIO_HDF5: return trio_hdf5_read_nucleus_num(file, num); break; - +/* case TRIO_JSON: return trio_json_read_nucleus_num(file, num); break; @@ -346,11 +346,11 @@ trio_exit_code trio_write_nucleus_num(trio_t* file, uint64_t num) { case TRIO_TEXT: return trio_text_write_nucleus_num(file, num); break; -/* + case TRIO_HDF5: return trio_hdf5_write_nucleus_num(file, num); break; - +/* case TRIO_JSON: return trio_json_write_nucleus_num(file, num); break; @@ -398,11 +398,11 @@ trio_exit_code trio_write_nucleus_coord(trio_t* file, double* coord) { case TRIO_TEXT: return trio_text_write_nucleus_coord(file, coord); break; -/* + case TRIO_HDF5: return trio_hdf5_write_nucleus_coord(file, coord); break; - +/* case TRIO_JSON: return trio_json_write_nucleus_coord(file, coord); break; @@ -839,6 +839,16 @@ trio_exit_code trio_text_write_nucleus_charge(const trio_t* file, const double* #+end_src ** HDF5 Back end +*** HDF5 definitions + +#+begin_src c :tangle trio_hdf5.c + #define NUCLEUS_GROUP_NAME "nucleus" + #define NUCLEUS_NUM_NAME "nucleus_num" + #define NUCLEUS_CHARGE_NAME "nucleus_charge" + #define NUCLEUS_COORD_NAME "nucleus_coord" + +#+end_src + *** HDF5 structures #+begin_src c :tangle trio_hdf5.h @@ -892,6 +902,58 @@ trio_exit_code trio_hdf5_init(trio_t* file); #+begin_src c :tangle trio_hdf5.c + +/* + * Currently H5LTread_dataset_ is used instead of this function + * but keep it for later if we decide to get rid of the H5LT API + */ +dset_t* trio_hdf5_read_dset_low(const trio_hdf5_t* file, const char *dset_name, void *buf) { + + assert (file != NULL); + assert (dset_name != NULL); + assert (buf != NULL); + /* + * Low-level implementation. Involves dealing with all HDF5 handles and dimensions + */ + dset_t* dset = (dset_t*) malloc(sizeof(dset_t)); + assert (dset != NULL); + + dset->dset_id = H5Dopen(file->nucleus_group, + dset_name, + H5P_DEFAULT); + assert (dset->dset_id > 0); + /* + * Get dataspace, datatype and dimensions + * dspace and dtype handles created below have to be closed when not used + */ + dset->dspace_id = H5Dget_space(dset->dset_id); + assert (dset->dspace_id > 0); + + dset->dtype_id = H5Dget_type(dset->dset_id); + assert (dset->dtype_id > 0); + + /* Check dimensions. Usefull, but then additional parameters + * ranks and dims[] have to be passed to the function + int rrank; + const int rank = 1; + hsize_t dims[1] = {0}; + rrank = H5Sget_simple_extent_dims(nucleus->h5_charge->dspace_id, + dims, NULL); + assert (rrank == rank); + for (int i=0; i 0); + } + */ + herr_t status; + status = H5Dread(dset->dset_id, dset->dtype_id, + H5S_ALL, H5S_ALL, H5P_DEFAULT, + buf); + assert (status >= 0); + + return dset; + +} + trio_exit_code trio_hdf5_init(trio_t* file) { trio_hdf5_t* f = (trio_hdf5_t*) file; @@ -911,18 +973,13 @@ trio_exit_code trio_hdf5_init(trio_t* file) { } /* Create groups in the hdf5 file */ - const char* nucleus_group_name = "nucleus"; - //const char* electron_group_name = "electron"; - if (f_ishere == 0){ - f->nucleus_group = H5Gcreate(f->file_id, nucleus_group_name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - //f->electron_group = H5Gcreate(f->file_id, electron_group_name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + f->nucleus_group = H5Gcreate(f->file_id, NUCLEUS_GROUP_NAME, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + //f->electron_group = H5Gcreate(f->file_id, ELECTRON_GROUP_NAME, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); } else { - f->nucleus_group = H5Gopen(f->file_id, nucleus_group_name, H5P_DEFAULT); - //f->electron_group = H5Gopen(f->file_id, electron_group_name, H5P_DEFAULT); + f->nucleus_group = H5Gopen(f->file_id, NUCLEUS_GROUP_NAME, H5P_DEFAULT); + //f->electron_group = H5Gopen(f->file_id, ELECTRON_GROUP_NAME, H5P_DEFAULT); } - /* not sure if assert statement here makes sence - H5Gcreate will raise its own H5 error if somethings is wrong*/ assert (f->nucleus_group > 0L); //assert (f->electron_group > 0L); @@ -948,6 +1005,10 @@ trio_exit_code trio_hdf5_finalize(trio_t* file) { H5Gclose(f->electron_group); f->electron_group = 0; */ + + H5Fclose(f->file_id); + f->file_id = 0; + return TRIO_SUCCESS; } #+end_src @@ -969,14 +1030,16 @@ h5nucleus_t* trio_hdf5_read_nucleus(const trio_hdf5_t* file) { nucleus->h5_coord = NULL; nucleus->h5_charge = NULL; - /* Try to open the file. If HDF5 cannot open, return */ + /* Check that the file was opened/created correctly, return */ if (file->file_id < 0) return nucleus; + /* Quit if the dimensioning attribute is missing in the file */ + if (H5Aexists(file->nucleus_group, NUCLEUS_NUM_NAME) == 0) return nucleus; + herr_t status; /* Read the nucleus_num attribute of nucleus group */ - const char *num_name = "nucleus_num"; hid_t num_id; - num_id = H5Aopen(file->nucleus_group, num_name, H5P_DEFAULT); + num_id = H5Aopen(file->nucleus_group, NUCLEUS_NUM_NAME, H5P_DEFAULT); assert (num_id > 0); status = H5Aread(num_id, H5T_NATIVE_ULLONG, &(nucleus->num)); @@ -986,85 +1049,152 @@ h5nucleus_t* trio_hdf5_read_nucleus(const trio_hdf5_t* file) { nucleus->charge = (double*) calloc(nucleus->num, sizeof(double)); assert (nucleus->charge != NULL); - /* - * High-level H5LT API. No need to deal with dataspaces and datatypes - */ - /* + /* High-level H5LT API. No need to deal with dataspaces and datatypes */ status = H5LTread_dataset_double(file->nucleus_group, - "nucleus_charge", - nucleus->charge) - */ - /* - * Low-level implementation. Involves dealing with all HDF5 handles and dimensions - */ - nucleus->h5_charge = (dset_t*) malloc(sizeof(dset_t)); - assert (nucleus->h5_charge != NULL); - - nucleus->h5_charge->dset_id = H5Dopen(file->nucleus_group, - "nucleus_charge", - H5P_DEFAULT); - assert (nucleus->h5_charge->dset_id > 0); - /* - * Get dataspace, datatype and dimensions - * dspace and dtype handles created below have to be closed when not used - */ - nucleus->h5_charge->dspace_id = H5Dget_space(nucleus->h5_charge->dset_id); - assert (nucleus->h5_charge->dspace_id > 0); - - nucleus->h5_charge->dtype_id = H5Dget_type(nucleus->h5_charge->dset_id); - assert (nucleus->h5_charge->dtype_id > 0); - - int rrank; - const int rank = 1; - hsize_t dims[1] = {0}; - rrank = H5Sget_simple_extent_dims(nucleus->h5_charge->dspace_id, - dims, NULL); - assert (rrank == rank); - for (int i=0; i 0); - } - - status = H5Dread(nucleus->h5_charge->dset_id, nucleus->h5_charge->dtype_id, - H5S_ALL, H5S_ALL, H5P_DEFAULT, - nucleus->charge); - assert (status >= 0); + NUCLEUS_CHARGE_NAME, + nucleus->charge); /* Allocate and read nucleus_coord array */ nucleus->coord = (double*) calloc(3 * nucleus->num, sizeof(double)); assert (nucleus->coord != NULL); + /* High-level H5LT API. No need to deal with dataspaces and datatypes */ status = H5LTread_dataset_double(file->nucleus_group, - "nucleus_coord", + NUCLEUS_COORD_NAME, nucleus->coord); assert (status >= 0); - /* Print arrays */ - /* - for (size_t i=0 ; inum ; i++) { - printf("%lf \n", nucleus->charge[i]); - } - - for (size_t i=0 ; i<3*nucleus->num ; i++) { - printf("%lf \n", nucleus->coord[i]); - } + + /* Low-level read. Do not forget to close the associated IDs (dset,dtype,dspace) + * when not used anymore, see below. Note how this function is similar to H5LTread_dataset_double */ + /* + nucleus->h5_coord = trio_hdf5_read_dset_low(file, NUCLEUS_COORD_NAME, + nucleus->coord); + + H5Sclose(nucleus->h5_coord->dspace_id); + H5Tclose(nucleus->h5_coord->dtype_id); + H5Dclose(nucleus->h5_coord->dset_id); + */ + H5Aclose(num_id); - H5Sclose(nucleus->h5_charge->dspace_id); - H5Tclose(nucleus->h5_charge->dtype_id); - H5Dclose(nucleus->h5_charge->dset_id); - H5Fclose(file->file_id); return nucleus; } trio_exit_code trio_hdf5_write_nucleus(const trio_hdf5_t* file, h5nucleus_t* nucleus) { - assert (nucleus != NULL); assert (file != NULL); + assert (nucleus != NULL); - // TODO - return TRIO_FAILURE; + herr_t status; + hid_t dspace, dtype; + hid_t attr_id; + + dtype = H5Tcopy(H5T_NATIVE_ULLONG); + /* Write the dimensioning variables */ + if (H5Aexists(file->nucleus_group, NUCLEUS_NUM_NAME) == 0) { + dspace = H5Screate(H5S_SCALAR); + attr_id = H5Acreate(file->nucleus_group, NUCLEUS_NUM_NAME, dtype, dspace, + H5P_DEFAULT, H5P_DEFAULT); + assert (attr_id > 0); + + /* High-level routine does not work for some reason + * status = H5LTset_attribute_ulong (file->nucleus_group, "nucleus", NUCLEUS_NUM_NAME, + * &(nucleus->num), 1); + */ + } else { + attr_id = H5Aopen(file->nucleus_group, NUCLEUS_NUM_NAME, H5P_DEFAULT); + assert (attr_id > 0); + } + + status = H5Awrite(attr_id, dtype, &(nucleus->num)); + assert (status >= 0); + + H5Aclose(attr_id); + + /* Write arrays */ + hid_t dset_id; + int charge_rank = 1; + const hsize_t charge_dims[1] = {nucleus->num}; + + if ( H5LTfind_dataset(file->nucleus_group, NUCLEUS_CHARGE_NAME) != 1) { + + status = H5LTmake_dataset_double (file->nucleus_group, NUCLEUS_CHARGE_NAME, + charge_rank, charge_dims, nucleus->charge); + assert (status >= 0); + + } else { + + dset_id = H5Dopen(file->nucleus_group, NUCLEUS_CHARGE_NAME, H5P_DEFAULT); + assert (dset_id > 0); + + dspace = H5Dget_space(dset_id); + assert (dspace > 0); + + dtype = H5Dget_type(dset_id); + assert (dtype > 0); + + + int rrank; + hsize_t dims[1] = {0}; + rrank = H5Sget_simple_extent_dims(dspace, + dims, NULL); + assert (rrank == charge_rank); + // disabling asserts like this allows to overwrite _num variable + for (int i=0; icharge); + assert (status >= 0); + + H5Sclose(dspace); + H5Tclose(dtype); + H5Dclose(dset_id); + + } + + int coord_rank = 2; + const hsize_t coord_dims[2] = {nucleus->num, 3}; + if ( H5LTfind_dataset(file->nucleus_group, NUCLEUS_COORD_NAME) != 1) { + status = H5LTmake_dataset_double (file->nucleus_group, NUCLEUS_COORD_NAME, + coord_rank, coord_dims, nucleus->coord); + assert (status >= 0); + + } else { + + dset_id = H5Dopen(file->nucleus_group, NUCLEUS_COORD_NAME, H5P_DEFAULT); + assert (dset_id > 0); + + dspace = H5Dget_space(dset_id); + assert (dspace > 0); + + dtype = H5Dget_type(dset_id); + assert (dtype > 0); + + + int rrank; + hsize_t dims[2] = {0, 0}; + rrank = H5Sget_simple_extent_dims(dspace, + dims, NULL); + assert (rrank == coord_rank); + for (int i=0; icoord); + assert (status >= 0); + + H5Sclose(dspace); + H5Tclose(dtype); + H5Dclose(dset_id); + + } + + return TRIO_SUCCESS; } + #+end_src *** Free memory @@ -1094,7 +1224,77 @@ trio_exit_code trio_hdf5_free_nucleus(h5nucleus_t* nucleus) { } #+end_src -*** TODO Read/Write the num attribute +*** Read/Write the num attribute + + #+begin_src c :tangle trio_hdf5.h +trio_exit_code trio_hdf5_read_nucleus_num(const trio_t* file, uint64_t* num); +trio_exit_code trio_hdf5_write_nucleus_num(const trio_t* file, const uint64_t num); + #+end_src + + #+begin_src c :tangle trio_hdf5.c +trio_exit_code trio_hdf5_read_nucleus_num(const trio_t* file, uint64_t* num) { + + assert (file != NULL); + assert (num != NULL); + + h5nucleus_t* nucleus = trio_hdf5_read_nucleus((trio_hdf5_t*) file); + + if (nucleus == NULL) { + return TRIO_FAILURE; + } + + /**/ *num = nucleus->num; + + trio_hdf5_free_nucleus(nucleus); + return TRIO_SUCCESS; +} + + +trio_exit_code trio_hdf5_write_nucleus_num(const trio_t* file, const uint64_t num) { + + assert (file != NULL); + assert (num > 0L); + + h5nucleus_t* nucleus = trio_hdf5_read_nucleus((trio_hdf5_t*) file); + + assert (nucleus != NULL); + + if (nucleus->num != num) { + + if (nucleus->num != 0) { + printf("%ld -> %ld %s \n", num, nucleus->num, + "This variable alreasy exists. Overwriting it is not supported"); + trio_hdf5_free_nucleus(nucleus); + return TRIO_FAILURE; + } + + nucleus->num = num; + + if (nucleus->charge != NULL) free(nucleus->charge); + nucleus->charge = NULL; + + nucleus->charge = (double*) calloc(num, sizeof(double)); + assert (nucleus->charge != NULL); + + if (nucleus->coord != NULL) free(nucleus->coord ); + nucleus->coord = NULL; + + nucleus->coord = (double*) calloc(3*num, sizeof(double)); + assert (nucleus->coord != NULL); + + } else { + nucleus->num = num; + } + + trio_exit_code rc = trio_hdf5_write_nucleus((trio_hdf5_t*) file, nucleus); + assert (rc == TRIO_SUCCESS); + + trio_hdf5_free_nucleus(nucleus); + + return TRIO_SUCCESS; +} + #+end_src + *** Read/Write the coord attribute @@ -1102,19 +1302,19 @@ trio_exit_code trio_hdf5_free_nucleus(h5nucleus_t* nucleus) { #+begin_src c :tangle trio_hdf5.h trio_exit_code trio_hdf5_read_nucleus_coord(const trio_t* file, double* coord); -// TODO -//trio_exit_code trio_hdf5_write_nucleus_coord(const trio_t* file, const double* coord); +trio_exit_code trio_hdf5_write_nucleus_coord(const trio_t* file, const double* coord); #+end_src #+begin_src c :tangle trio_hdf5.c trio_exit_code trio_hdf5_read_nucleus_coord(const trio_t* file, double* coord) { assert (file != NULL); - h5nucleus_t* nucleus = trio_hdf5_read_nucleus((trio_hdf5_t*) file); - - if (nucleus == NULL) return TRIO_FAILURE; - assert (coord != NULL); + + h5nucleus_t* nucleus = trio_hdf5_read_nucleus((trio_hdf5_t*) file); + + if (nucleus == NULL) return TRIO_FAILURE; + assert (nucleus->coord != NULL); for (size_t i=0 ; i<3*nucleus->num ; i++) { coord[i] = nucleus->coord[i]; @@ -1126,11 +1326,27 @@ trio_exit_code trio_hdf5_read_nucleus_coord(const trio_t* file, double* coord) { trio_exit_code trio_hdf5_write_nucleus_coord(const trio_t* file, const double* coord) { + assert (file != NULL); assert (coord != NULL); - // TODO - return TRIO_FAILURE; + + h5nucleus_t* nucleus = trio_hdf5_read_nucleus((trio_hdf5_t*) file); + + if (nucleus == NULL) return TRIO_FAILURE; + assert (nucleus->coord != NULL); + + for (size_t i=0 ; i<3*nucleus->num ; i++) { + nucleus->coord[i] = coord[i]; + } + + trio_exit_code rc = trio_hdf5_write_nucleus((trio_hdf5_t*) file, nucleus); + assert (rc == TRIO_SUCCESS); + + trio_hdf5_free_nucleus(nucleus); + + return TRIO_SUCCESS; } + #+end_src *** TODO Read/Write the charge attribute diff --git a/src/trio_hdf5.c b/src/trio_hdf5.c index a189ed9..34df2e1 100644 --- a/src/trio_hdf5.c +++ b/src/trio_hdf5.c @@ -7,6 +7,62 @@ #include "trio_hdf5.h" + #define NUCLEUS_GROUP_NAME "nucleus" + #define NUCLEUS_NUM_NAME "nucleus_num" + #define NUCLEUS_CHARGE_NAME "nucleus_charge" + #define NUCLEUS_COORD_NAME "nucleus_coord" + +/* + * Currently H5LTread_dataset_ is used instead of this function + * but keep it for later if we decide to get rid of the H5LT API + */ +dset_t* trio_hdf5_read_dset_low(const trio_hdf5_t* file, const char *dset_name, void *buf) { + + assert (file != NULL); + assert (dset_name != NULL); + assert (buf != NULL); + /* + * Low-level implementation. Involves dealing with all HDF5 handles and dimensions + */ + dset_t* dset = (dset_t*) malloc(sizeof(dset_t)); + assert (dset != NULL); + + dset->dset_id = H5Dopen(file->nucleus_group, + dset_name, + H5P_DEFAULT); + assert (dset->dset_id > 0); + /* + * Get dataspace, datatype and dimensions + * dspace and dtype handles created below have to be closed when not used + */ + dset->dspace_id = H5Dget_space(dset->dset_id); + assert (dset->dspace_id > 0); + + dset->dtype_id = H5Dget_type(dset->dset_id); + assert (dset->dtype_id > 0); + + /* Check dimensions. Usefull, but then additional parameters + * ranks and dims[] have to be passed to the function + int rrank; + const int rank = 1; + hsize_t dims[1] = {0}; + rrank = H5Sget_simple_extent_dims(nucleus->h5_charge->dspace_id, + dims, NULL); + assert (rrank == rank); + for (int i=0; i 0); + } + */ + herr_t status; + status = H5Dread(dset->dset_id, dset->dtype_id, + H5S_ALL, H5S_ALL, H5P_DEFAULT, + buf); + assert (status >= 0); + + return dset; + +} + trio_exit_code trio_hdf5_init(trio_t* file) { trio_hdf5_t* f = (trio_hdf5_t*) file; @@ -26,18 +82,13 @@ trio_exit_code trio_hdf5_init(trio_t* file) { } /* Create groups in the hdf5 file */ - const char* nucleus_group_name = "nucleus"; - //const char* electron_group_name = "electron"; - if (f_ishere == 0){ - f->nucleus_group = H5Gcreate(f->file_id, nucleus_group_name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - //f->electron_group = H5Gcreate(f->file_id, electron_group_name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + f->nucleus_group = H5Gcreate(f->file_id, NUCLEUS_GROUP_NAME, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + //f->electron_group = H5Gcreate(f->file_id, ELECTRON_GROUP_NAME, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); } else { - f->nucleus_group = H5Gopen(f->file_id, nucleus_group_name, H5P_DEFAULT); - //f->electron_group = H5Gopen(f->file_id, electron_group_name, H5P_DEFAULT); + f->nucleus_group = H5Gopen(f->file_id, NUCLEUS_GROUP_NAME, H5P_DEFAULT); + //f->electron_group = H5Gopen(f->file_id, ELECTRON_GROUP_NAME, H5P_DEFAULT); } - /* not sure if assert statement here makes sence - H5Gcreate will raise its own H5 error if somethings is wrong*/ assert (f->nucleus_group > 0L); //assert (f->electron_group > 0L); @@ -55,6 +106,10 @@ trio_exit_code trio_hdf5_finalize(trio_t* file) { H5Gclose(f->electron_group); f->electron_group = 0; */ + + H5Fclose(f->file_id); + f->file_id = 0; + return TRIO_SUCCESS; } @@ -70,14 +125,16 @@ h5nucleus_t* trio_hdf5_read_nucleus(const trio_hdf5_t* file) { nucleus->h5_coord = NULL; nucleus->h5_charge = NULL; - /* Try to open the file. If HDF5 cannot open, return */ + /* Check that the file was opened/created correctly, return */ if (file->file_id < 0) return nucleus; + /* Quit if the dimensioning attribute is missing in the file */ + if (H5Aexists(file->nucleus_group, NUCLEUS_NUM_NAME) == 0) return nucleus; + herr_t status; /* Read the nucleus_num attribute of nucleus group */ - const char *num_name = "nucleus_num"; hid_t num_id; - num_id = H5Aopen(file->nucleus_group, num_name, H5P_DEFAULT); + num_id = H5Aopen(file->nucleus_group, NUCLEUS_NUM_NAME, H5P_DEFAULT); assert (num_id > 0); status = H5Aread(num_id, H5T_NATIVE_ULLONG, &(nucleus->num)); @@ -87,84 +144,150 @@ h5nucleus_t* trio_hdf5_read_nucleus(const trio_hdf5_t* file) { nucleus->charge = (double*) calloc(nucleus->num, sizeof(double)); assert (nucleus->charge != NULL); - /* - * High-level H5LT API. No need to deal with dataspaces and datatypes - */ - /* + /* High-level H5LT API. No need to deal with dataspaces and datatypes */ status = H5LTread_dataset_double(file->nucleus_group, - "nucleus_charge", - nucleus->charge) - */ - /* - * Low-level implementation. Involves dealing with all HDF5 handles and dimensions - */ - nucleus->h5_charge = (dset_t*) malloc(sizeof(dset_t)); - assert (nucleus->h5_charge != NULL); - - nucleus->h5_charge->dset_id = H5Dopen(file->nucleus_group, - "nucleus_charge", - H5P_DEFAULT); - assert (nucleus->h5_charge->dset_id > 0); - /* - * Get dataspace, datatype and dimensions - * dspace and dtype handles created below have to be closed when not used - */ - nucleus->h5_charge->dspace_id = H5Dget_space(nucleus->h5_charge->dset_id); - assert (nucleus->h5_charge->dspace_id > 0); - - nucleus->h5_charge->dtype_id = H5Dget_type(nucleus->h5_charge->dset_id); - assert (nucleus->h5_charge->dtype_id > 0); - - int rrank; - const int rank = 1; - hsize_t dims[1] = {0}; - rrank = H5Sget_simple_extent_dims(nucleus->h5_charge->dspace_id, - dims, NULL); - assert (rrank == rank); - for (int i=0; i 0); - } - - status = H5Dread(nucleus->h5_charge->dset_id, nucleus->h5_charge->dtype_id, - H5S_ALL, H5S_ALL, H5P_DEFAULT, - nucleus->charge); - assert (status >= 0); + NUCLEUS_CHARGE_NAME, + nucleus->charge); /* Allocate and read nucleus_coord array */ nucleus->coord = (double*) calloc(3 * nucleus->num, sizeof(double)); assert (nucleus->coord != NULL); + /* High-level H5LT API. No need to deal with dataspaces and datatypes */ status = H5LTread_dataset_double(file->nucleus_group, - "nucleus_coord", + NUCLEUS_COORD_NAME, nucleus->coord); assert (status >= 0); - /* Print arrays */ - /* - for (size_t i=0 ; inum ; i++) { - printf("%lf \n", nucleus->charge[i]); - } - - for (size_t i=0 ; i<3*nucleus->num ; i++) { - printf("%lf \n", nucleus->coord[i]); - } + + /* Low-level read. Do not forget to close the associated IDs (dset,dtype,dspace) + * when not used anymore, see below. Note how this function is similar to H5LTread_dataset_double */ + /* + nucleus->h5_coord = trio_hdf5_read_dset_low(file, NUCLEUS_COORD_NAME, + nucleus->coord); + + H5Sclose(nucleus->h5_coord->dspace_id); + H5Tclose(nucleus->h5_coord->dtype_id); + H5Dclose(nucleus->h5_coord->dset_id); + */ + H5Aclose(num_id); - H5Sclose(nucleus->h5_charge->dspace_id); - H5Tclose(nucleus->h5_charge->dtype_id); - H5Dclose(nucleus->h5_charge->dset_id); - H5Fclose(file->file_id); return nucleus; } trio_exit_code trio_hdf5_write_nucleus(const trio_hdf5_t* file, h5nucleus_t* nucleus) { - assert (nucleus != NULL); assert (file != NULL); + assert (nucleus != NULL); - // TODO - return TRIO_FAILURE; + herr_t status; + hid_t dspace, dtype; + hid_t attr_id; + + dtype = H5Tcopy(H5T_NATIVE_ULLONG); + /* Write the dimensioning variables */ + if (H5Aexists(file->nucleus_group, NUCLEUS_NUM_NAME) == 0) { + dspace = H5Screate(H5S_SCALAR); + attr_id = H5Acreate(file->nucleus_group, NUCLEUS_NUM_NAME, dtype, dspace, + H5P_DEFAULT, H5P_DEFAULT); + assert (attr_id > 0); + + /* High-level routine does not work for some reason + * status = H5LTset_attribute_ulong (file->nucleus_group, "nucleus", NUCLEUS_NUM_NAME, + * &(nucleus->num), 1); + */ + } else { + attr_id = H5Aopen(file->nucleus_group, NUCLEUS_NUM_NAME, H5P_DEFAULT); + assert (attr_id > 0); + } + + status = H5Awrite(attr_id, dtype, &(nucleus->num)); + assert (status >= 0); + + H5Aclose(attr_id); + + /* Write arrays */ + hid_t dset_id; + int charge_rank = 1; + const hsize_t charge_dims[1] = {nucleus->num}; + + if ( H5LTfind_dataset(file->nucleus_group, NUCLEUS_CHARGE_NAME) != 1) { + + status = H5LTmake_dataset_double (file->nucleus_group, NUCLEUS_CHARGE_NAME, + charge_rank, charge_dims, nucleus->charge); + assert (status >= 0); + + } else { + + dset_id = H5Dopen(file->nucleus_group, NUCLEUS_CHARGE_NAME, H5P_DEFAULT); + assert (dset_id > 0); + + dspace = H5Dget_space(dset_id); + assert (dspace > 0); + + dtype = H5Dget_type(dset_id); + assert (dtype > 0); + + + int rrank; + hsize_t dims[1] = {0}; + rrank = H5Sget_simple_extent_dims(dspace, + dims, NULL); + assert (rrank == charge_rank); + // disabling asserts like this allows to overwrite _num variable + for (int i=0; icharge); + assert (status >= 0); + + H5Sclose(dspace); + H5Tclose(dtype); + H5Dclose(dset_id); + + } + + int coord_rank = 2; + const hsize_t coord_dims[2] = {nucleus->num, 3}; + if ( H5LTfind_dataset(file->nucleus_group, NUCLEUS_COORD_NAME) != 1) { + status = H5LTmake_dataset_double (file->nucleus_group, NUCLEUS_COORD_NAME, + coord_rank, coord_dims, nucleus->coord); + assert (status >= 0); + + } else { + + dset_id = H5Dopen(file->nucleus_group, NUCLEUS_COORD_NAME, H5P_DEFAULT); + assert (dset_id > 0); + + dspace = H5Dget_space(dset_id); + assert (dspace > 0); + + dtype = H5Dget_type(dset_id); + assert (dtype > 0); + + + int rrank; + hsize_t dims[2] = {0, 0}; + rrank = H5Sget_simple_extent_dims(dspace, + dims, NULL); + assert (rrank == coord_rank); + for (int i=0; icoord); + assert (status >= 0); + + H5Sclose(dspace); + H5Tclose(dtype); + H5Dclose(dset_id); + + } + + return TRIO_SUCCESS; } trio_exit_code trio_hdf5_free_nucleus(h5nucleus_t* nucleus) { @@ -188,14 +311,77 @@ trio_exit_code trio_hdf5_free_nucleus(h5nucleus_t* nucleus) { return TRIO_SUCCESS; } +trio_exit_code trio_hdf5_read_nucleus_num(const trio_t* file, uint64_t* num) { + + assert (file != NULL); + assert (num != NULL); + + h5nucleus_t* nucleus = trio_hdf5_read_nucleus((trio_hdf5_t*) file); + + if (nucleus == NULL) { + return TRIO_FAILURE; + } + + /**/ *num = nucleus->num; + + trio_hdf5_free_nucleus(nucleus); + return TRIO_SUCCESS; +} + + +trio_exit_code trio_hdf5_write_nucleus_num(const trio_t* file, const uint64_t num) { + + assert (file != NULL); + assert (num > 0L); + + h5nucleus_t* nucleus = trio_hdf5_read_nucleus((trio_hdf5_t*) file); + + assert (nucleus != NULL); + + if (nucleus->num != num) { + + if (nucleus->num != 0) { + printf("%ld -> %ld %s \n", num, nucleus->num, + "This variable alreasy exists. Overwriting it is not supported"); + trio_hdf5_free_nucleus(nucleus); + return TRIO_FAILURE; + } + + nucleus->num = num; + + if (nucleus->charge != NULL) free(nucleus->charge); + nucleus->charge = NULL; + + nucleus->charge = (double*) calloc(num, sizeof(double)); + assert (nucleus->charge != NULL); + + if (nucleus->coord != NULL) free(nucleus->coord ); + nucleus->coord = NULL; + + nucleus->coord = (double*) calloc(3*num, sizeof(double)); + assert (nucleus->coord != NULL); + + } else { + nucleus->num = num; + } + + trio_exit_code rc = trio_hdf5_write_nucleus((trio_hdf5_t*) file, nucleus); + assert (rc == TRIO_SUCCESS); + + trio_hdf5_free_nucleus(nucleus); + + return TRIO_SUCCESS; +} + trio_exit_code trio_hdf5_read_nucleus_coord(const trio_t* file, double* coord) { assert (file != NULL); - h5nucleus_t* nucleus = trio_hdf5_read_nucleus((trio_hdf5_t*) file); - - if (nucleus == NULL) return TRIO_FAILURE; - assert (coord != NULL); + + h5nucleus_t* nucleus = trio_hdf5_read_nucleus((trio_hdf5_t*) file); + + if (nucleus == NULL) return TRIO_FAILURE; + assert (nucleus->coord != NULL); for (size_t i=0 ; i<3*nucleus->num ; i++) { coord[i] = nucleus->coord[i]; @@ -207,8 +393,23 @@ trio_exit_code trio_hdf5_read_nucleus_coord(const trio_t* file, double* coord) { trio_exit_code trio_hdf5_write_nucleus_coord(const trio_t* file, const double* coord) { + assert (file != NULL); assert (coord != NULL); - // TODO - return TRIO_FAILURE; + + h5nucleus_t* nucleus = trio_hdf5_read_nucleus((trio_hdf5_t*) file); + + if (nucleus == NULL) return TRIO_FAILURE; + assert (nucleus->coord != NULL); + + for (size_t i=0 ; i<3*nucleus->num ; i++) { + nucleus->coord[i] = coord[i]; + } + + trio_exit_code rc = trio_hdf5_write_nucleus((trio_hdf5_t*) file, nucleus); + assert (rc == TRIO_SUCCESS); + + trio_hdf5_free_nucleus(nucleus); + + return TRIO_SUCCESS; } diff --git a/src/trio_hdf5.h b/src/trio_hdf5.h index ebb815c..aac83dd 100644 --- a/src/trio_hdf5.h +++ b/src/trio_hdf5.h @@ -62,8 +62,10 @@ trio_exit_code trio_hdf5_init(trio_t* file); trio_exit_code trio_hdf5_finalize(trio_t* file); +trio_exit_code trio_hdf5_read_nucleus_num(const trio_t* file, uint64_t* num); +trio_exit_code trio_hdf5_write_nucleus_num(const trio_t* file, const uint64_t num); + trio_exit_code trio_hdf5_read_nucleus_coord(const trio_t* file, double* coord); -// TODO -//trio_exit_code trio_hdf5_write_nucleus_coord(const trio_t* file, const double* coord); +trio_exit_code trio_hdf5_write_nucleus_coord(const trio_t* file, const double* coord); #endif