From 8adb4ac81354be0b5fe8dc736e82f495f07937e6 Mon Sep 17 00:00:00 2001 From: Yohai Meiron Date: Wed, 28 Oct 2020 22:46:27 -0400 Subject: [PATCH] Fixed memory issue with reading of HDF5 input file and improved style --- Makefile | 1 - TODO.md | 2 -- io.cpp | 76 +++++++++++++++++++++++++++++----------------------- io.h | 12 +++++++-- phigrape.cpp | 16 ++++++----- 5 files changed, 63 insertions(+), 44 deletions(-) diff --git a/Makefile b/Makefile index 06637d0..5b0a0b6 100644 --- a/Makefile +++ b/Makefile @@ -11,7 +11,6 @@ yebisu: GRAPEHOME = ../yebisu yebisu: GRAPELIB = -L$(GRAPEHOME) -lyebisug6 GRAPEINC = -I$(GRAPEHOME) -CXXFLAGS ?= -mcmodel=medium CXXFLAGS += -O$(OPTIMIZATION) INC = $(GRAPEINC) $(CUDAINC) LIB = $(GRAPELIB) $(CUDALIB) -lm diff --git a/TODO.md b/TODO.md index 7573942..22d82ff 100644 --- a/TODO.md +++ b/TODO.md @@ -1,8 +1,6 @@ TODO ==== -* Memory bug when reading HDF5? x and v not allocated. - * Break main() into smaller chunks; operations that are timed should become independent functions. * Const everything \ No newline at end of file diff --git a/io.cpp b/io.cpp index d0a187d..34c4bbf 100644 --- a/io.cpp +++ b/io.cpp @@ -3,6 +3,7 @@ #endif #include "double3.h" +#include "io.h" #include #include #include @@ -24,8 +25,9 @@ bool is_hdf5(std::string file_name) return result; } -void ascii_read(const std::string file_name, int &step_num, int &N, double& t, std::vector &m, std::vector &x, std::vector &v) +Input_data ascii_read(const std::string &file_name) { + Input_data input_data; char rest[512]; int result; @@ -34,30 +36,31 @@ void ascii_read(const std::string file_name, int &step_num, int &N, double& t, s std::string str; std::getline(file, str); - result = sscanf(str.c_str(), "%d%s", &step_num, rest); + result = sscanf(str.c_str(), "%d%s", &input_data.step_num, rest); if (result!=1) throw std::runtime_error("Error parsing line 1: expected one integer"); std::getline(file, str); - result = sscanf(str.c_str(), "%d%s", &N, rest); + result = sscanf(str.c_str(), "%d%s", &input_data.N, rest); if (result!=1) throw std::runtime_error("Error parsing line 2: expected one integer"); std::getline(file, str); - result = sscanf(str.c_str(), "%lf%s", &t, rest); + result = sscanf(str.c_str(), "%lf%s", &input_data.t, rest); if (result!=1) throw std::runtime_error("Error parsing line 3: expected one real number"); - m.resize(N); - x.resize(N); - v.resize(N); + input_data.m.resize(input_data.N); + input_data.x.resize(input_data.N); + input_data.v.resize(input_data.N); int i = -1; while (std::getline(file, str)) { - if (++i > N) throw std::runtime_error("Error parsing line " + std::to_string(i+4) + ": particle out of range"); - result = sscanf(str.c_str(), "%*s %lf %lf %lf %lf %lf %lf %lf%s", &m[i], &(x[i].x), &(x[i].y), &(x[i].z), &(v[i].x), &(v[i].y), &(v[i].z), rest); + if (++i > input_data.N) throw std::runtime_error("Error parsing line " + std::to_string(i+4) + ": particle out of range"); + result = sscanf(str.c_str(), "%*s %lf %lf %lf %lf %lf %lf %lf%s", &input_data.m[i], &(input_data.x[i].x), &(input_data.x[i].y), &(input_data.x[i].z), &(input_data.v[i].x), &(input_data.v[i].y), &(input_data.v[i].z), rest); } file.close(); + return input_data; } -void ascii_write(const std::string file_name, const int step_num, const int N, const double t, const std::vector &m, const std::vector &x, const std::vector &v, int precision=10) +void ascii_write(const std::string file_name, const int step_num, const int N, const double t, const std::vector &m, const std::vector &x, const std::vector &v, int precision) { auto file = std::ofstream(file_name); if (!file.is_open()) throw std::runtime_error("Cannot open file for output"); @@ -75,9 +78,11 @@ void ascii_write(const std::string file_name, const int step_num, const int N, c file.close(); } -void h5_read(const std::string file_name, int *step_num, int *N, double *t, double m[], double3 x[], double3 v[]) +Input_data h5_read(const std::string &file_name) { #ifdef HAS_HDF5 + Input_data input_data; + // Open file and root group; count number of top level objects hid_t file_id = H5Fopen(file_name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); hid_t group_id = H5Gopen2(file_id, "/", H5P_DEFAULT); @@ -92,7 +97,7 @@ void h5_read(const std::string file_name, int *step_num, int *N, double *t, doub if (name.substr(0, 5) == "Step#") step_num_arr.push_back(std::stoi(name.substr(5, std::string::npos))); } // Find the highest step in the file - *step_num = *max_element(step_num_arr.begin(), step_num_arr.end()); + input_data.step_num = *max_element(step_num_arr.begin(), step_num_arr.end()); // Prepare to read the data sets dimensionality data. char path[256]; @@ -100,53 +105,58 @@ void h5_read(const std::string file_name, int *step_num, int *N, double *t, doub int ndims; hsize_t dims[2]; - sprintf(path, "/Step#%d", *step_num); + sprintf(path, "/Step#%d", input_data.step_num); attr_id = H5Aopen_by_name(file_id, path, "Time", H5P_DEFAULT, H5P_DEFAULT); - H5Aread(attr_id, H5T_NATIVE_DOUBLE, t); + H5Aread(attr_id, H5T_NATIVE_DOUBLE, &input_data.t); H5Aclose(attr_id); - sprintf(path, "/Step#%d/MASS", *step_num); + sprintf(path, "/Step#%d/MASS", input_data.step_num); dataset_id = H5Dopen2(file_id, path, H5P_DEFAULT); dataspace_id = H5Dget_space(dataset_id); ndims = H5Sget_simple_extent_ndims(dataspace_id); if (ndims != 1) - throw std::runtime_error("Dataset MASS in Step#" + std::to_string(*step_num) + " of file " + file_name + " is " + std::to_string(ndims) + "-dimensional (expected 1-dimensional)"); + throw std::runtime_error("Dataset MASS in Step#" + std::to_string(input_data.step_num) + " of file " + file_name + " is " + std::to_string(ndims) + "-dimensional (expected 1-dimensional)"); H5Sget_simple_extent_dims(dataspace_id, dims, NULL); - H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, m); - H5Sclose(dataspace_id); - H5Dclose(dataset_id); - *N = dims[0]; + input_data.N = dims[0]; - sprintf(path, "/Step#%d/X", *step_num); - dataset_id = H5Dopen2(file_id, path, H5P_DEFAULT); - dataspace_id = H5Dget_space(dataset_id); - ndims = H5Sget_simple_extent_ndims(dataspace_id); - if (ndims != 2) throw std::runtime_error("Bad dimensionality"); - H5Sget_simple_extent_dims(dataspace_id, dims, NULL); - if ((dims[0] != *N) || (dims[1] != 3)) throw std::runtime_error("Bad dimensionality"); - H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, (double*)x); + input_data.m.resize(input_data.N); + H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, input_data.m.data()); H5Sclose(dataspace_id); H5Dclose(dataset_id); - sprintf(path, "/Step#%d/V", *step_num); + sprintf(path, "/Step#%d/X", input_data.step_num); dataset_id = H5Dopen2(file_id, path, H5P_DEFAULT); dataspace_id = H5Dget_space(dataset_id); ndims = H5Sget_simple_extent_ndims(dataspace_id); - if (ndims != 2) throw std::runtime_error("Bad dimensionality"); + if (ndims != 2) throw std::runtime_error("Bad dimensionality in " + std::string(path)); H5Sget_simple_extent_dims(dataspace_id, dims, NULL); - if ((dims[0] != *N) || (dims[1] != 3)) throw std::runtime_error("Bad dimensionality"); - H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, (double*)v); + if ((dims[0] != input_data.N) || (dims[1] != 3)) throw std::runtime_error("Bad dimensionality"); + input_data.x.resize(input_data.N); + H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, (double*)input_data.x.data()); + H5Sclose(dataspace_id); + H5Dclose(dataset_id); + + sprintf(path, "/Step#%d/V", input_data.step_num); + dataset_id = H5Dopen2(file_id, path, H5P_DEFAULT); + dataspace_id = H5Dget_space(dataset_id); + ndims = H5Sget_simple_extent_ndims(dataspace_id); + if (ndims != 2) throw std::runtime_error("Bad dimensionality in" + std::string(path)); + H5Sget_simple_extent_dims(dataspace_id, dims, NULL); + if ((dims[0] != input_data.N) || (dims[1] != 3)) throw std::runtime_error("Bad dimensionality"); + input_data.v.resize(input_data.N); + H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, (double*)input_data.v.data()); H5Sclose(dataspace_id); H5Dclose(dataset_id); H5Gclose(group_id); H5Fclose(file_id); + return input_data; #else throw std::runtime_error("h5_read was called but compiled without HDF5 support"); #endif } -void h5_write(const std::string file_name, const int step_num, const int N, const double t, const std::vector &m, const std::vector &x, const std::vector &v, const std::vector &pot, const std::vector &acc, const std::vector &jrk, const int extra_mode=0, const bool use_double_precision=true) +void h5_write(const std::string file_name, const int step_num, const int N, const double t, const std::vector &m, const std::vector &x, const std::vector &v, const std::vector &pot, const std::vector &acc, const std::vector &jrk, const int extra_mode, const bool use_double_precision) { #ifdef HAS_HDF5 hid_t file_id, group_id, attribute_id, dataspace_id; diff --git a/io.h b/io.h index 7571592..097ebb9 100644 --- a/io.h +++ b/io.h @@ -1,15 +1,23 @@ #pragma once #include +#include #include "double3.h" +struct Input_data { + int N, step_num; + double t; + std::vector m; + std::vector x, v; +}; + bool is_hdf5(std::string file_name); // This function is implemented independently of the HDF5 library -void ascii_read(const std::string file_name, int &step_num, int &N, double &t, std::vector &m, std::vector &x, std::vector &v); +Input_data ascii_read(const std::string &file_name); void ascii_write(const std::string file_name, const int step_num, const int N, const double t, const std::vector &m, const std::vector &x, const std::vector &v, int precision=10); -void h5_read(const std::string file_name, int *step_num, int *N, double *t, double m[], double3 x[], double3 v[]); +Input_data h5_read(const std::string &file_name); // In case the code is compiled without HDF5 support, the implementation of this function just throws an error void h5_write(const std::string file_name, const int step_num, const int N, const double t, const std::vector &m, const std::vector &x, const std::vector &v, const std::vector &pot, const std::vector &acc, const std::vector &jrk, const int extra_mode=0, const bool use_double_precision=true); diff --git a/phigrape.cpp b/phigrape.cpp index 845557f..1cab74e 100644 --- a/phigrape.cpp +++ b/phigrape.cpp @@ -276,19 +276,23 @@ int main(int argc, char *argv[]) /* Print the Rank and the names of processors */ printf("Rank of the processor %03d on %s \n", myRank, processor_name); - int diskstep, N; - double time_cur; - std::vector m; - std::vector x, v; + Input_data input_data; if (is_hdf5(config.input_file_name)) { #ifndef HAS_HDF5 fprintf(stderr, "ERROR: input file is in HDF5 format, but the code was compiled without HDF5 support\n"); return -1; #endif - h5_read(config.input_file_name, &diskstep, &N, &time_cur, m.data(), x.data(), v.data()); + input_data = h5_read(config.input_file_name); } else - ascii_read(config.input_file_name, diskstep, N, time_cur, m, x, v); + input_data = ascii_read(config.input_file_name); + + int N = input_data.N; + int diskstep = input_data.step_num; + double time_cur = input_data.t; + auto &m = input_data.m; + auto &x = input_data.x; + auto &v = input_data.v; if (myRank == rootRank) { printf("\n");