From 997bc92f52273f0ed88fde647dd70e6aab21ad9e Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 3 Mar 2021 10:09:50 +0100 Subject: [PATCH] Forgot org files --- src/trexio_hdf5.org | 587 +++++++++++++++++++++++++++++++ src/trexio_text.org | 834 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 1421 insertions(+) create mode 100644 src/trexio_hdf5.org create mode 100644 src/trexio_text.org diff --git a/src/trexio_hdf5.org b/src/trexio_hdf5.org new file mode 100644 index 0000000..06e8df2 --- /dev/null +++ b/src/trexio_hdf5.org @@ -0,0 +1,587 @@ +#+Title: HDF5 back end for the TREX Input/Ouput library (TREXIO) + +* File prefixes :noxport: + + #+NAME:header + #+begin_src c +/* This file was generated from the trexio.org org-mode file. + To generate it, open trexio.org in Emacs and execute + M-x org-babel-tangle +*/ + + + + #+end_src + + #+begin_src c :tangle trexio_hdf5.h :noweb yes +<
> +#ifndef _TREXIO_HDF5_H +#define _TREXIO_HDF5_H + +#include "trexio.h" +#include "trexio_s.h" +#include +#include +#include +#include +#include +#include + +#include "hdf5.h" +#include "hdf5_hl.h" // needed for high-level APIs like H5LT, requires additional linking in Makefile + #+end_src + + #+begin_src c :tangle trexio_hdf5.c :noweb yes +<
> +#include "trexio_hdf5.h" + #+end_src + + +* HDF5 Back end + +** HDF5 definitions + +#+begin_src c :tangle trexio_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 trexio_hdf5.h + +typedef struct slab_s { + uint64_t a; + uint64_t b; + uint64_t c; + uint64_t d; +} slab_t; + +typedef struct dset_s { + hid_t dset_id; + hid_t dspace_id; + hid_t dtype_id; + uint64_t* dims; + uint32_t rank; + const char* dset_name; +} dset_t; + +typedef struct h5nucleus_s { + uint64_t num; + double *coord; + double *charge; + dset_t* h5_coord; + dset_t* h5_charge; +} h5nucleus_t; + +typedef struct h5electron_s { + uint64_t alpha_num; + uint64_t beta_num; +} h5electron_t; + +typedef struct trexio_hdf5_s { + trexio_t parent ; + hid_t file_id; + hid_t nucleus_group; + hid_t electron_group; + //... other groups' id + const char* file_name; +} trexio_hdf5_t; + + #+end_src + +** HDF5 basic functions + + + #+begin_src c :tangle trexio_hdf5.h +trexio_exit_code trexio_hdf5_init(trexio_t* file); + #+end_src + + + #+begin_src c :tangle trexio_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* trexio_hdf5_read_dset_low(const trexio_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; + +} + +trexio_exit_code trexio_hdf5_init(trexio_t* file) { + + trexio_hdf5_t* f = (trexio_hdf5_t*) file; + + /* If file doesn't exist, create it */ + int f_ishere = 0; + struct stat st; + + if (stat(file->file_name, &st) == 0) { + printf("%s \n","HDF5 file already exists"); + // RDWR OR RDONLY ??? + f->file_id = H5Fopen(file->file_name, H5F_ACC_RDWR, H5P_DEFAULT); + f_ishere = 1; + } else { + f->file_id = H5Fcreate(file->file_name, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + f_ishere = 0; + } + + /* Create groups in the hdf5 file */ + 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); + } 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); + } + assert (f->nucleus_group > 0L); + //assert (f->electron_group > 0L); + + return TREXIO_SUCCESS; +} + + #+end_src + + + #+begin_src c :tangle trexio_hdf5.h +trexio_exit_code trexio_hdf5_finalize(trexio_t* file); + #+end_src + + #+begin_src c :tangle trexio_hdf5.c +trexio_exit_code trexio_hdf5_finalize(trexio_t* file) { + + trexio_hdf5_t* f = (trexio_hdf5_t*) file; + + H5Gclose(f->nucleus_group); + f->nucleus_group = 0; + +/* + H5Gclose(f->electron_group); + f->electron_group = 0; +*/ + + H5Fclose(f->file_id); + f->file_id = 0; + + return TREXIO_SUCCESS; +} + #+end_src + + + +** Structs for sparse data structures + #+begin_src c :tangle trexio_hdf5.h +typedef struct one_index_s { + double value; + int64_t i; +} one_index_t; + +typedef struct two_index_s { + double value; + int64_t i; + int64_t j; +} two_index_t; + +typedef struct three_index_s { + double value; + int64_t i; + int64_t j; + int64_t k; +} three_index_t; + +typedef struct four_index_s { + double value; + int64_t i; + int64_t j; + int64_t k; + int64_t l; +} four_index_t; + + #+end_src + +** Read/write the nucleus struct + + #+begin_src c :tangle trexio_hdf5.c +h5nucleus_t* trexio_hdf5_read_nucleus(const trexio_hdf5_t* file) { + + /* Allocate the data structure */ + h5nucleus_t* nucleus = (h5nucleus_t*) malloc(sizeof(h5nucleus_t)); + assert (nucleus != NULL); + + nucleus->num = 0; + nucleus->coord = NULL; + nucleus->charge = NULL; + nucleus->h5_coord = NULL; + nucleus->h5_charge = NULL; + + /* 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 */ + hid_t num_id; + num_id = H5Aopen(file->nucleus_group, NUCLEUS_NUM_NAME, H5P_DEFAULT); + assert (num_id > 0); + + status = H5Aread(num_id, H5T_NATIVE_ULLONG, &(nucleus->num)); + assert (status >= 0); + + /* Allocate and read nucleus_charge array */ + nucleus->charge = (double*) calloc(nucleus->num, sizeof(double)); + assert (nucleus->charge != NULL); + + /* High-level H5LT API. No need to deal with dataspaces and datatypes */ + status = H5LTread_dataset_double(file->nucleus_group, + 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_NAME, + nucleus->coord); + assert (status >= 0); + + /* 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 = trexio_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); + + return nucleus; +} + + +trexio_exit_code trexio_hdf5_write_nucleus(const trexio_hdf5_t* file, h5nucleus_t* nucleus) { + + assert (file != NULL); + assert (nucleus != NULL); + + 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 TREXIO_SUCCESS; +} + + #+end_src + +** Free memory + + Memory is allocated when reading. The followig function frees memory. + + #+begin_src c :tangle trexio_hdf5.c +trexio_exit_code trexio_hdf5_free_nucleus(h5nucleus_t* nucleus) { + + if (nucleus == NULL) return TREXIO_FAILURE; + + if (nucleus->coord != NULL) free (nucleus->coord); + nucleus->coord = NULL; + + if (nucleus->charge != NULL) free (nucleus->charge); + nucleus->charge = NULL; + + if (nucleus->h5_coord != NULL) free (nucleus->h5_coord); + nucleus->h5_coord = NULL; + + if (nucleus->h5_charge != NULL) free (nucleus->h5_charge); + nucleus->h5_charge = NULL; + + free (nucleus); + + return TREXIO_SUCCESS; +} + #+end_src + +** Read/Write the num attribute + + #+begin_src c :tangle trexio_hdf5.h +trexio_exit_code trexio_hdf5_read_nucleus_num(const trexio_t* file, uint64_t* num); +trexio_exit_code trexio_hdf5_write_nucleus_num(const trexio_t* file, const uint64_t num); + #+end_src + + #+begin_src c :tangle trexio_hdf5.c +trexio_exit_code trexio_hdf5_read_nucleus_num(const trexio_t* file, uint64_t* num) { + + assert (file != NULL); + assert (num != NULL); + + h5nucleus_t* nucleus = trexio_hdf5_read_nucleus((trexio_hdf5_t*) file); + + if (nucleus == NULL) { + return TREXIO_FAILURE; + } + + /**/ *num = nucleus->num; + + trexio_hdf5_free_nucleus(nucleus); + return TREXIO_SUCCESS; +} + + +trexio_exit_code trexio_hdf5_write_nucleus_num(const trexio_t* file, const uint64_t num) { + + assert (file != NULL); + assert (num > 0L); + + h5nucleus_t* nucleus = trexio_hdf5_read_nucleus((trexio_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"); + trexio_hdf5_free_nucleus(nucleus); + return TREXIO_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; + } + + trexio_exit_code rc = trexio_hdf5_write_nucleus((trexio_hdf5_t*) file, nucleus); + assert (rc == TREXIO_SUCCESS); + + trexio_hdf5_free_nucleus(nucleus); + + return TREXIO_SUCCESS; +} + #+end_src + + +** Read/Write the coord attribute + + The ~coord~ array is assumed allocated with the appropriate size. + + #+begin_src c :tangle trexio_hdf5.h +trexio_exit_code trexio_hdf5_read_nucleus_coord(const trexio_t* file, double* coord); +trexio_exit_code trexio_hdf5_write_nucleus_coord(const trexio_t* file, const double* coord); + #+end_src + + #+begin_src c :tangle trexio_hdf5.c +trexio_exit_code trexio_hdf5_read_nucleus_coord(const trexio_t* file, double* coord) { + + assert (file != NULL); + assert (coord != NULL); + + h5nucleus_t* nucleus = trexio_hdf5_read_nucleus((trexio_hdf5_t*) file); + + if (nucleus == NULL) return TREXIO_FAILURE; + assert (nucleus->coord != NULL); + + for (size_t i=0 ; i<3*nucleus->num ; i++) { + coord[i] = nucleus->coord[i]; + } + + trexio_hdf5_free_nucleus(nucleus); + return TREXIO_SUCCESS; +} + + +trexio_exit_code trexio_hdf5_write_nucleus_coord(const trexio_t* file, const double* coord) { + + assert (file != NULL); + assert (coord != NULL); + + h5nucleus_t* nucleus = trexio_hdf5_read_nucleus((trexio_hdf5_t*) file); + + if (nucleus == NULL) return TREXIO_FAILURE; + assert (nucleus->coord != NULL); + + for (size_t i=0 ; i<3*nucleus->num ; i++) { + nucleus->coord[i] = coord[i]; + } + + trexio_exit_code rc = trexio_hdf5_write_nucleus((trexio_hdf5_t*) file, nucleus); + assert (rc == TREXIO_SUCCESS); + + trexio_hdf5_free_nucleus(nucleus); + + return TREXIO_SUCCESS; +} + + #+end_src +** TODO Read/Write the charge attribute + +* File suffixes :noxport: + + #+begin_src c :tangle trexio_hdf5.h +#endif + #+end_src diff --git a/src/trexio_text.org b/src/trexio_text.org new file mode 100644 index 0000000..208cd5e --- /dev/null +++ b/src/trexio_text.org @@ -0,0 +1,834 @@ +#+Title: TEXT back end of the TREX Input/Ouput library (TREXIO) + +* File prefixes :noxport: + + #+NAME:header + #+begin_src c +/* This file was generated from the trexio.org org-mode file. + To generate it, open trexio.org in Emacs and execute + M-x org-babel-tangle +*/ + + + + #+end_src + + #+begin_src c :tangle trexio_text.h :noweb yes +<
> +#ifndef _TREXIO_TEXT_H +#define _TREXIO_TEXT_H + +#include "trexio.h" +#include "trexio_s.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include + + #+end_src + + #+begin_src c :tangle trexio_text.c :noweb yes +<
> +#include "trexio_text.h" + #+end_src + +* TEXT Back end + + The "file" produced by the text back end is a directory with one + file per group. + + When the file is open, it is locked by the current process. No other + process can read/write the same file. This guarantees that the + representation in memory is consistent with the file and avoid + re-reading the file before writing. + To lock the file, we lock the =.lock= file which is present in the + directory. + + The file is written when closed, or when the flush function is called. + +*** Structs for blocks + #+begin_src c :tangle trexio_text.h +typedef struct nucleus_s { + FILE* file; + uint64_t dim_coord; + uint64_t dim_charge; + double* coord; + double* charge; + uint64_t num; + int to_flush; +} nucleus_t; + +typedef struct electron_s { + FILE* file; + uint64_t alpha_num; + uint64_t beta_num; + int to_flush; +} electron_t; + +typedef struct rdm_s { + FILE* file; + uint64_t dim_one_e; + double* one_e; + char* two_e_file_name; + int to_flush; +} rdm_t; + #+end_src + + +*** TO DO + - to_flush = 1 in write + - to_flush = 0 when flushed + - name + +*** Structs for the text back end + + #+begin_src c :tangle trexio_text.h +typedef struct trexio_text_s { + trexio_t parent ; + int lock_file; + + nucleus_t* nucleus; + electron_t* electron; + rdm_t* rdm; +} trexio_text_t; + + #+end_src + + + #+begin_src c :tangle trexio_text.h +trexio_exit_code trexio_text_init(trexio_t* file); + #+end_src + + #+begin_src c :tangle trexio_text.c +trexio_exit_code trexio_text_init(trexio_t* file) { + if (file == NULL) return TREXIO_INVALID_ARG_1; + + trexio_text_t* f = (trexio_text_t*) file; + + /* If directory doesn't exist, create it in write mode */ + struct stat st; + + if (stat(file->file_name, &st) == 0 && S_ISDIR(st.st_mode)) { + /* Do nothing */ + } else { + if (file->mode == 'r') return TREXIO_READONLY; + + if (mkdir(file->file_name, 0777) != 0) { + return TREXIO_FAILURE; + } + } + + /* Create the lock file in the directory */ + const char* lock_file_name = "/.lock"; + char* file_name = (char*) + calloc( strlen(file->file_name) + strlen(lock_file_name) + 1, + sizeof(char)); + assert (file_name != NULL); + strcpy (file_name, file->file_name); + strcat (file_name, lock_file_name); + + f->lock_file = open(file_name,O_WRONLY|O_CREAT|O_TRUNC, 0644); + assert (f->lock_file > 0); + free(file_name); + + /* Lock the directory for the current process */ + struct flock fl; + + fl.l_type = F_WRLCK; + fl.l_whence = SEEK_SET; + fl.l_start = 0; + fl.l_len = 0; + fl.l_pid = getpid(); + + int rc = fcntl(f->lock_file, F_SETLKW, &fl); + if (rc == -1) return TREXIO_FAILURE; + + trexio_text_t* f_child = (trexio_text_t*) file; + + f_child->nucleus = NULL; + f_child->electron= NULL; + f_child->rdm = NULL; + + return TREXIO_SUCCESS; +} + + #+end_src + + + #+begin_src c :tangle trexio_text.h +trexio_exit_code trexio_text_finalize(trexio_t* file); + #+end_src + + #+begin_src c :tangle trexio_text.c +trexio_exit_code trexio_text_finalize(trexio_t* file) { + if (file == NULL) return TREXIO_INVALID_ARG_1; + + trexio_text_t* f = (trexio_text_t*) file; + + trexio_exit_code rc; + rc = trexio_text_free_nucleus( (trexio_text_t*) file); + assert (rc == TREXIO_SUCCESS); + + rc = trexio_text_free_rdm( (trexio_text_t*) file); + assert (rc == TREXIO_SUCCESS); + + /* Release lock */ + struct flock fl; + + fl.l_type = F_WRLCK; + fl.l_whence = SEEK_SET; + fl.l_start = 0; + fl.l_len = 0; + fl.l_pid = getpid(); + fl.l_type = F_UNLCK; + fcntl(f->lock_file, F_SETLK, &fl); + + close(f->lock_file); + return TREXIO_SUCCESS; +} + #+end_src + +*** Nucleus struct + +**** Read the struct + + #+begin_src c :tangle trexio_text.h +nucleus_t* trexio_text_read_nucleus(trexio_text_t* file); + #+end_src + + #+begin_src c :tangle trexio_text.c +nucleus_t* trexio_text_read_nucleus(trexio_text_t* file) { + if (file == NULL) return NULL; + + if (file->nucleus != NULL) return file->nucleus; + + /* Allocate the data structure */ + nucleus_t* nucleus = (nucleus_t*) malloc(sizeof(nucleus_t)); + assert (nucleus != NULL); + + nucleus->file = NULL; + nucleus->num = 0; + nucleus->coord = NULL; + nucleus->charge = NULL; + nucleus->to_flush = 0; + + /* Try to open the file. If the file does not exist, return */ + const char* nucleus_file_name = "/nucleus.txt"; + char * file_name = (char*) + calloc( strlen(file->parent.file_name) + strlen(nucleus_file_name) + 1, + sizeof(char)); + assert (file_name != NULL); + strcpy (file_name, file->parent.file_name); + strcat (file_name, nucleus_file_name); + + /* If the file exists, read it */ + FILE* f = fopen(file_name,"r"); + if (f != NULL) { + + /* Find size of file to allocate the max size of the string buffer */ + fseek(f, 0L, SEEK_END); + size_t sz = ftell(f); + fseek(f, 0L, SEEK_SET); + char* buffer = (char*) malloc(sz*sizeof(char)); + + /* Read the dimensioning variables */ + int rc; + rc = fscanf(f, "%s", buffer); + assert (rc == 1); + assert (strcmp(buffer, "dim_charge") == 0); + + rc = fscanf(f, "%lu", &(nucleus->dim_charge)); + assert (rc == 1); + + rc = fscanf(f, "%s", buffer); + assert (rc == 1); + assert (strcmp(buffer, "dim_coord") == 0); + + rc = fscanf(f, "%lu", &(nucleus->dim_coord)); + assert (rc == 1); + + /* Allocate arrays */ + nucleus->charge = (double*) calloc(nucleus->dim_charge, sizeof(double)); + assert (nucleus->charge != NULL); + + nucleus->coord = (double*) calloc(nucleus->dim_coord, sizeof(double)); + assert (nucleus->coord != NULL); + + /* Read data */ + rc = fscanf(f, "%s", buffer); + assert (rc == 1); + assert (strcmp(buffer, "num") == 0); + + rc = fscanf(f, "%lu", &(nucleus->num)); + assert (rc == 1); + + rc = fscanf(f, "%s", buffer); + assert (rc == 1); + assert (strcmp(buffer, "charge") == 0); + + for (uint64_t i=0 ; idim_charge ; i++) { + rc = fscanf(f, "%lf", &(nucleus->charge[i])); + assert (rc == 1); + } + + rc = fscanf(f, "%s", buffer); + assert (rc == 1); + assert (strcmp(buffer, "coord") == 0); + + for (uint64_t i=0 ; idim_coord ; i++) { + rc = fscanf(f, "%lf", &(nucleus->coord[i])); + assert (rc == 1); + } + free(buffer); + fclose(f); + f = NULL; + } + if (file->parent.mode == 'w') { + nucleus->file = fopen(file_name,"a"); + } else { + nucleus->file = fopen(file_name,"r"); + } + free(file_name); + file->nucleus = nucleus; + return nucleus; +} + #+end_src + +**** Flush the struct + + #+begin_src c :tangle trexio_text.h +trexio_exit_code trexio_text_flush_nucleus(const trexio_text_t* file); + #+end_src + + #+begin_src c :tangle trexio_text.c +trexio_exit_code trexio_text_flush_nucleus(const trexio_text_t* file) { + if (file == NULL) return TREXIO_INVALID_ARG_1; + + if (file->parent.mode == 'r') return TREXIO_READONLY; + + nucleus_t* nucleus = file->nucleus; + + if (nucleus == NULL) return TREXIO_SUCCESS; + + if (nucleus->to_flush == 0) return TREXIO_SUCCESS; + + FILE* f = nucleus->file; + assert (f != NULL); + rewind(f); + + /* Write the dimensioning variables */ + fprintf(f, "dim_charge %ld\n", nucleus->dim_charge); + fprintf(f, "dim_coord %ld\n", nucleus->dim_coord ); + + /* Write arrays */ + fprintf(f, "num %ld\n", nucleus->num); + fprintf(f, "charge\n"); + for (uint64_t i=0 ; idim_charge ; i++) { + fprintf(f, "%lf\n", nucleus->charge[i]); + } + + fprintf(f, "coord\n"); + for (uint64_t i=0 ; idim_coord ; i++) { + fprintf(f, "%lf\n", nucleus->coord[i]); + } + fflush(f); + file->nucleus->to_flush = 0; + return TREXIO_SUCCESS; +} + #+end_src + +**** Free memory + + Memory is allocated when reading. The following function frees memory. + + #+begin_src c :tangle trexio_text.h +trexio_exit_code trexio_text_free_nucleus(trexio_text_t* file); + #+end_src + + #+begin_src c :tangle trexio_text.c +trexio_exit_code trexio_text_free_nucleus(trexio_text_t* file) { + if (file == NULL) return TREXIO_INVALID_ARG_1; + + trexio_exit_code rc; + + if (file->parent.mode != 'r') { + rc = trexio_text_flush_nucleus(file); + if (rc != TREXIO_SUCCESS) return TREXIO_FAILURE; + } + + nucleus_t* nucleus = file->nucleus; + if (nucleus == NULL) return TREXIO_SUCCESS; + + if (nucleus->file != NULL) { + fclose(nucleus->file); + nucleus->file = NULL; + } + + if (nucleus->coord != NULL) { + free (nucleus->coord); + nucleus->coord = NULL; + } + + if (nucleus->charge != NULL) { + free (nucleus->charge); + nucleus->charge = NULL; + } + + free (nucleus); + file->nucleus = NULL; + return TREXIO_SUCCESS; +} + #+end_src + +**** Read/Write the num attribute + + #+begin_src c :tangle trexio_text.h +trexio_exit_code trexio_text_read_nucleus_num(const trexio_t* file, uint64_t* num); +trexio_exit_code trexio_text_write_nucleus_num(const trexio_t* file, const uint64_t num); + #+end_src + + #+begin_src c :tangle trexio_text.c +trexio_exit_code trexio_text_read_nucleus_num(const trexio_t* file, uint64_t* num) { + + if (file == NULL) return TREXIO_INVALID_ARG_1; + if (num == NULL) return TREXIO_INVALID_ARG_2; + + nucleus_t* nucleus = trexio_text_read_nucleus((trexio_text_t*) file); + if (nucleus == NULL) return TREXIO_FAILURE; + + /**/ *num = nucleus->num; + + return TREXIO_SUCCESS; +} + + +trexio_exit_code trexio_text_write_nucleus_num(const trexio_t* file, const uint64_t num) { + + if (file == NULL) return TREXIO_INVALID_ARG_1; + + if (file->mode == 'r') return TREXIO_READONLY; + + nucleus_t* nucleus = trexio_text_read_nucleus((trexio_text_t*) file); + if (nucleus == NULL) return TREXIO_FAILURE; + + nucleus->num = num; + nucleus->to_flush = 1; + + return TREXIO_SUCCESS; +} + #+end_src + +**** Read/Write the coord attribute + + The ~coord~ array is assumed allocated with the appropriate size. + + #+begin_src c :tangle trexio_text.h +trexio_exit_code trexio_text_read_nucleus_coord(const trexio_t* file, double* coord, const uint64_t dim_coord); +trexio_exit_code trexio_text_write_nucleus_coord(const trexio_t* file, const double* coord, const uint64_t dim_coord); + #+end_src + + #+begin_src c :tangle trexio_text.c +trexio_exit_code trexio_text_read_nucleus_coord(const trexio_t* file, double* coord, const uint64_t dim_coord) { + + if (file == NULL) return TREXIO_INVALID_ARG_1; + if (coord == NULL) return TREXIO_INVALID_ARG_2; + + nucleus_t* nucleus = trexio_text_read_nucleus((trexio_text_t*) file); + if (nucleus == NULL) return TREXIO_FAILURE; + + if (dim_coord != nucleus->dim_coord) return TREXIO_INVALID_ARG_3; + + for (uint64_t i=0 ; icoord[i]; + } + + return TREXIO_SUCCESS; +} + + +trexio_exit_code trexio_text_write_nucleus_coord(const trexio_t* file, const double* coord, const uint64_t dim_coord) { + if (file == NULL) return TREXIO_INVALID_ARG_1; + if (coord == NULL) return TREXIO_INVALID_ARG_2; + + if (file->mode == 'r') return TREXIO_READONLY; + + nucleus_t* nucleus = trexio_text_read_nucleus((trexio_text_t*) file); + if (nucleus == NULL) return TREXIO_FAILURE; + + if (nucleus->coord != NULL) { + free(nucleus->coord); + nucleus->coord = NULL; + } + + nucleus->dim_coord = dim_coord; + nucleus->coord = (double*) calloc(dim_coord, sizeof(double)); + + for (uint64_t i=0 ; icoord[i] = coord[i]; + } + + nucleus->to_flush = 1; + return TREXIO_SUCCESS; +} + #+end_src +**** Read/Write the charge attribute + + The ~charge~ array is assumed allocated with the appropriate size. + + #+begin_src c :tangle trexio_text.h +trexio_exit_code trexio_text_read_nucleus_charge(const trexio_t* file, double* charge, const uint64_t dim_charge); +trexio_exit_code trexio_text_write_nucleus_charge(const trexio_t* file, const double* charge, const uint64_t dim_charge); + #+end_src + + #+begin_src c :tangle trexio_text.c +trexio_exit_code trexio_text_read_nucleus_charge(const trexio_t* file, double* charge, const uint64_t dim_charge) { + + if (file == NULL) return TREXIO_INVALID_ARG_1; + if (charge == NULL) return TREXIO_INVALID_ARG_2; + + nucleus_t* nucleus = trexio_text_read_nucleus((trexio_text_t*) file); + if (nucleus == NULL) return TREXIO_FAILURE; + + if (dim_charge != nucleus->dim_charge) return TREXIO_INVALID_ARG_3; + + for (uint64_t i=0 ; icharge[i]; + } + + return TREXIO_SUCCESS; +} + + +trexio_exit_code trexio_text_write_nucleus_charge(const trexio_t* file, const double* charge, const uint64_t dim_charge) { + if (file == NULL) return TREXIO_INVALID_ARG_1; + if (charge == NULL) return TREXIO_INVALID_ARG_2; + + if (file->mode == 'r') return TREXIO_READONLY; + + nucleus_t* nucleus = trexio_text_read_nucleus((trexio_text_t*) file); + if (nucleus == NULL) return TREXIO_FAILURE; + + if (nucleus->charge != NULL) { + free(nucleus->charge); + nucleus->charge = NULL; + } + + nucleus->dim_charge = dim_charge; + nucleus->charge = (double*) calloc(dim_charge, sizeof(double)); + + for (uint64_t i=0 ; icharge[i] = charge[i]; + } + + nucleus->to_flush = 1; + return TREXIO_SUCCESS; +} + #+end_src + +*** RDM struct +**** Read the complete struct + + #+begin_src c :tangle trexio_text.h +rdm_t* trexio_text_read_rdm(trexio_text_t* file); + #+end_src + + #+begin_src c :tangle trexio_text.c +rdm_t* trexio_text_read_rdm(trexio_text_t* file) { + if (file == NULL) return NULL; + + if (file->rdm != NULL) return file->rdm; + + /* Allocate the data structure */ + rdm_t* rdm = (rdm_t*) malloc(sizeof(rdm_t)); + assert (rdm != NULL); + + rdm->one_e = NULL; + rdm->two_e_file_name = NULL; + rdm->file = NULL; + rdm->to_flush = 0; + + /* Try to open the file. If the file does not exist, return */ + const char* rdm_file_name = "/rdm.txt"; + char * file_name = (char*) + calloc( strlen(file->parent.file_name) + strlen(rdm_file_name) + 1, + sizeof(char)); + assert (file_name != NULL); + strcpy (file_name, file->parent.file_name); + strcat (file_name, rdm_file_name); + + /* If the file exists, read it */ + FILE* f = fopen(file_name,"r"); + if (f != NULL) { + + /* Find size of file to allocate the max size of the string buffer */ + fseek(f, 0L, SEEK_END); + size_t sz = ftell(f); + fseek(f, 0L, SEEK_SET); + char* buffer = (char*) malloc(sz*sizeof(char)); + + /* Read the dimensioning variables */ + int rc; + rc = fscanf(f, "%s", buffer); + assert (rc == 1); + assert (strcmp(buffer, "dim_one_e") == 0); + + rc = fscanf(f, "%lu", &(rdm->dim_one_e)); + assert (rc == 1); + + /* Allocate arrays */ + rdm->one_e = (double*) calloc(rdm->dim_one_e, sizeof(double)); + assert (rdm->one_e != NULL); + + /* Read one_e */ + rc = fscanf(f, "%s", buffer); + assert (rc == 1); + assert (strcmp(buffer, "one_e") == 0); + + for (uint64_t i=0 ; idim_one_e; i++) { + rc = fscanf(f, "%lf", &(rdm->one_e[i])); + assert (rc == 1); + } + + /* Read two_e */ + rc = fscanf(f, "%s", buffer); + assert (rc == 1); + assert (strcmp(buffer, "two_e_file_name") == 0); + + rc = fscanf(f, "%s", buffer); + assert (rc == 1); + rdm->two_e_file_name = (char*) malloc (strlen(buffer)*sizeof(char)); + strcpy(rdm->two_e_file_name, buffer); + + free(buffer); + fclose(f); + f = NULL; + } + if (file->parent.mode == 'w') { + rdm->file = fopen(file_name,"a"); + } else { + rdm->file = fopen(file_name,""); + } + free(file_name); + file->rdm = rdm ; + return rdm; +} + #+end_src + +**** Flush the complete struct + + #+begin_src c :tangle trexio_text.h +trexio_exit_code trexio_text_flush_rdm(const trexio_text_t* file); + #+end_src + + #+begin_src c :tangle trexio_text.c +trexio_exit_code trexio_text_flush_rdm(const trexio_text_t* file) { + if (file == NULL) return TREXIO_INVALID_ARG_1; + + if (file->parent.mode == 'r') return TREXIO_READONLY; + + rdm_t* rdm = file->rdm; + if (rdm == NULL) return TREXIO_SUCCESS; + + if (rdm->to_flush == 0) return TREXIO_SUCCESS; + + FILE* f = rdm->file; + assert (f != NULL); + rewind(f); + + /* Write the dimensioning variables */ + fprintf(f, "num %ld\n", rdm->dim_one_e); + + /* Write arrays */ + fprintf(f, "one_e\n"); + for (uint64_t i=0 ; i< rdm->dim_one_e; i++) { + fprintf(f, "%lf\n", rdm->one_e[i]); + } + + fprintf(f, "two_e_file_name\n"); + fprintf(f, "%s\n", rdm->two_e_file_name); + + fflush(f); + file->rdm->to_flush = 0; + return TREXIO_SUCCESS; +} + #+end_src + +**** Free memory + + Memory is allocated when reading. The followig function frees memory. + + #+begin_src c :tangle trexio_text.h +trexio_exit_code trexio_text_free_rdm(trexio_text_t* file); + #+end_src + + #+begin_src c :tangle trexio_text.c +trexio_exit_code trexio_text_free_rdm(trexio_text_t* file) { + if (file == NULL) return TREXIO_INVALID_ARG_1; + + trexio_exit_code rc; + if (file->parent.mode != 'r') { + rc = trexio_text_flush_rdm(file); + if (rc != TREXIO_SUCCESS) return TREXIO_FAILURE; + } + + rdm_t* rdm = file->rdm; + if (rdm == NULL) return TREXIO_SUCCESS; + + if (rdm->file != NULL) { + fclose(rdm->file); + rdm->file = NULL; + } + + if (rdm->one_e != NULL) { + free (rdm->one_e); + rdm->one_e = NULL; + } + + if (rdm->two_e_file_name != NULL) { + free (rdm->two_e_file_name); + rdm->two_e_file_name = NULL; + } + + free (rdm); + file->rdm = NULL; + return TREXIO_SUCCESS; +} + #+end_src + +**** Read/Write the one_e attribute + + The ~one_e~ array is assumed allocated with the appropriate size. + + #+begin_src c :tangle trexio_text.h +trexio_exit_code trexio_text_read_rdm_one_e(const trexio_t* file, double* one_e, const uint64_t dim_one_e); +trexio_exit_code trexio_text_write_rdm_one_e(const trexio_t* file, const double* one_e, const uint64_t dim_one_e); + #+end_src + + #+begin_src c :tangle trexio_text.c +trexio_exit_code trexio_text_read_rdm_one_e(const trexio_t* file, double* one_e, const uint64_t dim_one_e) { + + if (file == NULL) return TREXIO_INVALID_ARG_1; + if (one_e == NULL) return TREXIO_INVALID_ARG_2; + + rdm_t* rdm = trexio_text_read_rdm((trexio_text_t*) file); + if (rdm == NULL) return TREXIO_FAILURE; + + if (dim_one_e != rdm->dim_one_e) return TREXIO_INVALID_ARG_3; + + for (uint64_t i=0 ; ione_e[i]; + } + + return TREXIO_SUCCESS; +} + + +trexio_exit_code trexio_text_write_rdm_one_e(const trexio_t* file, const double* one_e, const uint64_t dim_one_e) { + if (file == NULL) return TREXIO_INVALID_ARG_1; + if (one_e == NULL) return TREXIO_INVALID_ARG_2; + if (file->mode != 'r') return TREXIO_READONLY; + + rdm_t* rdm = trexio_text_read_rdm((trexio_text_t*) file); + if (rdm == NULL) return TREXIO_FAILURE; + + rdm->dim_one_e = dim_one_e; + for (uint64_t i=0 ; ione_e[i] = one_e[i]; + } + + rdm->to_flush = 1; + return TREXIO_SUCCESS; +} + #+end_src + +**** Read/Write the two_e attribute + + ~two_e~ is a sparse data structure, which can be too large to fit + in memory. So we provide functions to read and write it by + chunks. + In the text back end, the easiest way to do it is to create a + file for each sparse float structure. + + #+begin_src c :tangle trexio_text.h +trexio_exit_code trexio_text_buffered_read_rdm_two_e(const trexio_t* file, const uint64_t offset, const uint64_t size, int64_t* index, double* value); +trexio_exit_code trexio_text_buffered_write_rdm_two_e(const trexio_t* file, const uint64_t offset, const uint64_t size, const int64_t* index, const double* value); + #+end_src + + #+begin_src c :tangle trexio_text.c +trexio_exit_code trexio_text_buffered_read_rdm_two_e(const trexio_t* file, const uint64_t offset, const uint64_t size, int64_t* index, double* value) { + + if (file == NULL) return TREXIO_INVALID_ARG_1; + if (index == NULL) return TREXIO_INVALID_ARG_4; + if (value == NULL) return TREXIO_INVALID_ARG_5; + + rdm_t* rdm = trexio_text_read_rdm((trexio_text_t*) file); + if (rdm == NULL) return TREXIO_FAILURE; + + FILE* f = fopen(rdm->two_e_file_name, "r"); + if (f == NULL) return TREXIO_END; + + const uint64_t line_length = 64; + fseek(f, (long) offset * line_length, SEEK_SET); + + int rc; + for (uint64_t i=0 ; imode != 'r') return TREXIO_READONLY; + + rdm_t* rdm = trexio_text_read_rdm((trexio_text_t*) file); + if (rdm == NULL) return TREXIO_FAILURE; + + FILE* f = fopen(rdm->two_e_file_name, "w"); + if (f == NULL) return TREXIO_FAILURE; + + const uint64_t line_length = 64; + fseek(f, (long) offset * line_length, SEEK_SET); + + int rc; + for (uint64_t i=0 ; i