diff --git a/.github/workflows/run_tests/action.yml b/.github/workflows/run_tests/action.yml index 928991ce3..8413a880b 100644 --- a/.github/workflows/run_tests/action.yml +++ b/.github/workflows/run_tests/action.yml @@ -23,6 +23,8 @@ runs: mpirun -n 3 --oversubscribe tests/test_StaticSVD ./tests/test_IncrementalSVDBrand mpirun -n 3 --oversubscribe tests/test_IncrementalSVDBrand + ./tests/test_HDFDatabase + mpirun -n 3 --oversubscribe tests/test_HDFDatabase ./tests/test_NNLS mpirun -n 3 --oversubscribe tests/test_NNLS shell: bash diff --git a/CMakeLists.txt b/CMakeLists.txt index f3f1eaf15..4acea3b82 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -265,6 +265,7 @@ if(GTEST_FOUND) set(unit_test_stems Vector Matrix + HDFDatabase DEIM DMD GNAT diff --git a/docker/Dockerfile b/docker/Dockerfile index 1cfdfc866..b7161cbfc 100644 --- a/docker/Dockerfile +++ b/docker/Dockerfile @@ -10,7 +10,7 @@ WORKDIR /$ENVDIR # install packages RUN sudo apt-get install -yq git RUN sudo apt-get install --no-install-recommends -yq make gcc gfortran libssl-dev cmake -RUN sudo apt-get install -yq libopenblas-dev libmpich-dev libblas-dev liblapack-dev libscalapack-mpi-dev libhdf5-serial-dev +RUN sudo apt-get install -yq libopenblas-dev libmpich-dev libblas-dev liblapack-dev libscalapack-mpi-dev libhdf5-mpi-dev hdf5-tools RUN sudo apt-get install -yq vim RUN sudo apt-get install -yq git-lfs RUN sudo apt-get install -yq valgrind diff --git a/lib/CAROM_config.h.in b/lib/CAROM_config.h.in index 9a0246992..569332ebd 100644 --- a/lib/CAROM_config.h.in +++ b/lib/CAROM_config.h.in @@ -57,6 +57,9 @@ /* Defined if you have HDF5 support */ #cmakedefine01 CAROM_HAVE_HDF5 +/* Defined if you have parallel HDF5 support */ +#cmakedefine01 HDF5_IS_PARALLEL + /* Define to 1 if you have the header file. */ #cmakedefine01 CAROM_HAVE_INTTYPES_H diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt index c21ae8169..0e546061d 100644 --- a/lib/CMakeLists.txt +++ b/lib/CMakeLists.txt @@ -54,6 +54,7 @@ set(module_list hyperreduction/Hyperreduction utils/Database utils/HDFDatabase + utils/HDFDatabaseMPIO utils/CSVDatabase utils/Utilities utils/ParallelBuffer diff --git a/lib/linalg/BasisGenerator.cpp b/lib/linalg/BasisGenerator.cpp index 46b9bda9f..9e2cec4ae 100644 --- a/lib/linalg/BasisGenerator.cpp +++ b/lib/linalg/BasisGenerator.cpp @@ -31,6 +31,7 @@ BasisGenerator::BasisGenerator( bool incremental, const std::string& basis_file_name, Database::formats file_format) : + d_dim(options.dim), d_incremental(incremental), d_basis_writer(0), d_basis_reader(0), @@ -157,7 +158,7 @@ BasisGenerator::loadSampleRange(const std::string& base_file_name, if (d_basis_reader) delete d_basis_reader; - d_basis_reader = new BasisReader(base_file_name, db_format); + d_basis_reader = new BasisReader(base_file_name, db_format, d_dim); const Matrix* mat; const Vector* singular_vals; diff --git a/lib/linalg/BasisGenerator.h b/lib/linalg/BasisGenerator.h index ec6f7898d..87604447e 100644 --- a/lib/linalg/BasisGenerator.h +++ b/lib/linalg/BasisGenerator.h @@ -64,7 +64,7 @@ class BasisGenerator Options options, bool incremental, const std::string& basis_file_name = "", - Database::formats file_format = Database::HDF5); + Database::formats file_format = Database::formats::HDF5); /** * @brief Destructor. @@ -155,7 +155,7 @@ class BasisGenerator const std::string& kind = "basis", int col_min = 0, int col_max = 1e9, - Database::formats db_format = Database::HDF5); + Database::formats db_format = Database::formats::HDF5); /** * @brief Load previously saved sample (basis or state). @@ -170,8 +170,8 @@ class BasisGenerator void loadSamples(const std::string& base_file_name, const std::string& kind = "basis", - int cutoff = 1e9, - Database::formats db_format = Database::HDF5); + int cut_off = 1e9, + Database::formats db_format = Database::formats::HDF5); /** * @brief Computes next time an svd sample is needed. @@ -322,7 +322,7 @@ class BasisGenerator int getDim() { - return d_svd->getDim(); + return d_dim; } /** @@ -376,6 +376,13 @@ class BasisGenerator * @brief The number of processors being run on. */ int d_num_procs; + + /** + * @brief Dimension of the system on this processor. + * + * Equivalent to d_svd->getDim(). + */ + const int d_dim; }; } diff --git a/lib/linalg/BasisReader.cpp b/lib/linalg/BasisReader.cpp index b974b3cd1..1fdb11e65 100644 --- a/lib/linalg/BasisReader.cpp +++ b/lib/linalg/BasisReader.cpp @@ -12,18 +12,22 @@ #include "BasisReader.h" #include "utils/HDFDatabase.h" -#include "utils/CSVDatabase.h" +#include "utils/HDFDatabaseMPIO.h" #include "Matrix.h" #include "Vector.h" #include "mpi.h" +#include "utils/mpi_utils.h" namespace CAROM { BasisReader::BasisReader( const std::string& base_file_name, - Database::formats db_format) : + Database::formats db_format, + const int dim) : + d_dim(dim), full_file_name(""), - base_file_name_(base_file_name) + base_file_name_(base_file_name), + d_format(db_format) { CAROM_ASSERT(!base_file_name.empty()); @@ -37,17 +41,29 @@ BasisReader::BasisReader( rank = 0; } - char tmp[100]; - sprintf(tmp, ".%06d", rank); - full_file_name = base_file_name + tmp; - if (db_format == Database::HDF5) { + full_file_name = base_file_name; + + // Enforce hdf data format. + if (d_format == Database::formats::HDF5) + { d_database = new HDFDatabase(); } - else if (db_format == Database::CSV) { - d_database = new CSVDatabase(); + else if (d_format == Database::formats::HDF5_MPIO) + { + /* + For MPIO case, local dimension needs to be specified. + We allow 0 local dimension. (global dimension still needs to be positive) + */ + std::vector tmp; + d_global_dim = get_global_offsets(d_dim, tmp, MPI_COMM_WORLD); + CAROM_VERIFY(d_dim >= 0); + CAROM_VERIFY(d_global_dim > 0); + d_database = new HDFDatabaseMPIO(); } + else + CAROM_ERROR("BasisWriter only supports HDF5/HDF5_MPIO data format!\n"); - d_database->open(full_file_name, "r"); + d_database->open(full_file_name, "r", MPI_COMM_WORLD); } BasisReader::~BasisReader() @@ -66,7 +82,8 @@ BasisReader::getSpatialBasis() d_database->getDoubleArray("spatial_basis", &spatial_basis_vectors->item(0, 0), - num_rows*num_cols); + num_rows*num_cols, + true); return spatial_basis_vectors; } @@ -97,7 +114,8 @@ BasisReader::getSpatialBasis( num_rows*num_cols_to_read, start_col - 1, num_cols_to_read, - num_cols); + num_cols, + true); return spatial_basis_vectors; } @@ -135,7 +153,7 @@ BasisReader::getTemporalBasis() int num_cols = getNumSamples("temporal_basis"); char tmp[100]; - Matrix* temporal_basis_vectors = new Matrix(num_rows, num_cols, true); + Matrix* temporal_basis_vectors = new Matrix(num_rows, num_cols, false); sprintf(tmp, "temporal_basis"); d_database->getDoubleArray(tmp, &temporal_basis_vectors->item(0, 0), @@ -163,7 +181,7 @@ BasisReader::getTemporalBasis( CAROM_VERIFY(start_col <= end_col && end_col <= num_cols); int num_cols_to_read = end_col - start_col + 1; - Matrix* temporal_basis_vectors = new Matrix(num_rows, num_cols_to_read, true); + Matrix* temporal_basis_vectors = new Matrix(num_rows, num_cols_to_read, false); sprintf(tmp, "temporal_basis"); d_database->getDoubleArray(tmp, &temporal_basis_vectors->item(0, 0), @@ -268,7 +286,18 @@ BasisReader::getDim( "temporal_basis_num_rows"); d_database->getInteger(tmp, num_rows); - return num_rows; + /* only basis and snapshot are stored as distributed matrices */ + if ((kind != "temporal_basis") && (d_format == Database::formats::HDF5_MPIO)) + { + /* + for MPIO database, return specified local dimension. + only checks the global dimension match. + */ + CAROM_VERIFY(d_global_dim == num_rows); + return d_dim; + } + else + return num_rows; } int @@ -298,11 +327,12 @@ BasisReader::getSnapshotMatrix() int num_cols = getNumSamples("snapshot"); char tmp[100]; - Matrix* snapshots = new Matrix(num_rows, num_cols, false); + Matrix* snapshots = new Matrix(num_rows, num_cols, true); sprintf(tmp, "snapshot_matrix"); d_database->getDoubleArray(tmp, &snapshots->item(0, 0), - num_rows*num_cols); + num_rows*num_cols, + true); return snapshots; } @@ -326,14 +356,15 @@ BasisReader::getSnapshotMatrix( int num_cols_to_read = end_col - start_col + 1; char tmp[100]; - Matrix* snapshots = new Matrix(num_rows, num_cols_to_read, false); + Matrix* snapshots = new Matrix(num_rows, num_cols_to_read, true); sprintf(tmp, "snapshot_matrix"); d_database->getDoubleArray(tmp, &snapshots->item(0, 0), num_rows*num_cols_to_read, start_col - 1, num_cols_to_read, - num_cols); + num_cols, + true); return snapshots; } } diff --git a/lib/linalg/BasisReader.h b/lib/linalg/BasisReader.h index eb17998a2..fee556e8f 100644 --- a/lib/linalg/BasisReader.h +++ b/lib/linalg/BasisReader.h @@ -44,10 +44,13 @@ class BasisReader { * @param[in] db_format Format of the file to read. * One of the implemented file formats defined in * Database. + * @param[in] dim Number of rows of basis that will be read from a file. + * If negative, will use the dimension from the rank-specific local file. */ BasisReader( const std::string& base_file_name, - Database::formats db_format = Database::HDF5); + Database::formats db_format = Database::formats::HDF5, + const int dim = -1); /** * @brief Destructor. @@ -197,6 +200,11 @@ class BasisReader { * * @brief Returns the dimension of the system on this processor. * + * @param[in] kind Type of matrix whose row-dimension is returned. + * "basis" - local dimension of spatial basis + * "snapshot" - local dimension of snapshot matrix + * "temporal_basis" - global dimension of temporal basis + * * @return The dimension of the system on this processor. */ int @@ -277,6 +285,11 @@ class BasisReader { */ Database* d_database; + /** + * @brief The database being read from. + */ + Database::formats d_format; + /** * @brief Base file name stored for consistency between reading and writing. */ @@ -286,6 +299,20 @@ class BasisReader { * @brief Full file name of database including rank. */ std::string full_file_name; + + /** + * @brief Dimension of the basis on this processor. + * + * If negative, use the dimension from the rank-specific local file. + */ + const int d_dim; + + /** + * @brief Dimension of the basis on this processor. + * + * If negative, use the dimension from the rank-specific local file. + */ + int d_global_dim; }; } diff --git a/lib/linalg/BasisWriter.cpp b/lib/linalg/BasisWriter.cpp index b2bfd0ffd..1257006b8 100644 --- a/lib/linalg/BasisWriter.cpp +++ b/lib/linalg/BasisWriter.cpp @@ -12,6 +12,7 @@ #include "BasisWriter.h" #include "utils/HDFDatabase.h" +#include "utils/HDFDatabaseMPIO.h" #include "Matrix.h" #include "Vector.h" #include "BasisGenerator.h" @@ -45,18 +46,22 @@ BasisWriter::BasisWriter( rank = 0; } - char tmp[100]; - sprintf(tmp, ".%06d", rank); - full_file_name = base_file_name + tmp; - - char tmp2[100]; - sprintf(tmp2, "_snapshot.%06d", rank); - snap_file_name = base_file_name + tmp2; + full_file_name = base_file_name; + snap_file_name = base_file_name + "_snapshot"; // create and open snapshot/basis database - CAROM_VERIFY(db_format_ == Database::HDF5); - d_snap_database = new HDFDatabase(); - d_database = new HDFDatabase(); + if (db_format_ == Database::formats::HDF5) + { + d_snap_database = new HDFDatabase(); + d_database = new HDFDatabase(); + } + else if (db_format_ == Database::formats::HDF5_MPIO) + { + d_snap_database = new HDFDatabaseMPIO(); + d_database = new HDFDatabaseMPIO(); + } + else + CAROM_ERROR("BasisWriter only supports HDF5/HDF5_MPIO data format!\n"); } BasisWriter::~BasisWriter() @@ -73,20 +78,27 @@ BasisWriter::writeBasis(const std::string& kind) char tmp[100]; if (kind == "basis") { - d_database->create(full_file_name); + d_database->create(full_file_name, MPI_COMM_WORLD); const Matrix* basis = d_basis_generator->getSpatialBasis(); + /* spatial basis is always distributed */ + CAROM_VERIFY(basis->distributed()); int num_rows = basis->numRows(); + int nrows_infile = num_rows; + if (db_format_ == Database::formats::HDF5_MPIO) + MPI_Allreduce(MPI_IN_PLACE, &nrows_infile, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); sprintf(tmp, "spatial_basis_num_rows"); - d_database->putInteger(tmp, num_rows); + d_database->putInteger(tmp, nrows_infile); int num_cols = basis->numColumns(); sprintf(tmp, "spatial_basis_num_cols"); d_database->putInteger(tmp, num_cols); sprintf(tmp, "spatial_basis"); - d_database->putDoubleArray(tmp, &basis->item(0, 0), num_rows*num_cols); + d_database->putDoubleArray(tmp, &basis->item(0, 0), num_rows*num_cols, true); if(d_basis_generator->updateRightSV()) { const Matrix* tbasis = d_basis_generator->getTemporalBasis(); + /* temporal basis is always not distributed */ + CAROM_VERIFY(!tbasis->distributed()); num_rows = tbasis->numRows(); sprintf(tmp, "temporal_basis_num_rows"); d_database->putInteger(tmp, num_rows); @@ -98,6 +110,8 @@ BasisWriter::writeBasis(const std::string& kind) } const Vector* sv = d_basis_generator->getSingularValues(); + /* singular values are always not distributed */ + CAROM_VERIFY(!sv->distributed()); int sv_dim = sv->dim(); sprintf(tmp, "singular_value_size"); d_database->putInteger(tmp, sv_dim); @@ -108,17 +122,23 @@ BasisWriter::writeBasis(const std::string& kind) } if (kind == "snapshot") { - d_snap_database->create(snap_file_name); + d_snap_database->create(snap_file_name, MPI_COMM_WORLD); const Matrix* snapshots = d_basis_generator->getSnapshotMatrix(); + /* snapshot matrix is always distributed */ + CAROM_VERIFY(snapshots->distributed()); int num_rows = snapshots->numRows(); // d_dim + int nrows_infile = num_rows; + if (db_format_ == Database::formats::HDF5_MPIO) + MPI_Allreduce(MPI_IN_PLACE, &nrows_infile, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); sprintf(tmp, "snapshot_matrix_num_rows"); - d_snap_database->putInteger(tmp, num_rows); + d_snap_database->putInteger(tmp, nrows_infile); int num_cols = snapshots->numColumns(); // d_num_samples sprintf(tmp, "snapshot_matrix_num_cols"); d_snap_database->putInteger(tmp, num_cols); sprintf(tmp, "snapshot_matrix"); - d_snap_database->putDoubleArray(tmp, &snapshots->item(0,0), num_rows*num_cols); + d_snap_database->putDoubleArray(tmp, &snapshots->item(0,0), num_rows*num_cols, + true); d_snap_database->close(); } diff --git a/lib/linalg/BasisWriter.h b/lib/linalg/BasisWriter.h index d495b6ef0..9684460d3 100644 --- a/lib/linalg/BasisWriter.h +++ b/lib/linalg/BasisWriter.h @@ -42,7 +42,7 @@ class BasisWriter { BasisWriter( BasisGenerator* basis_generator, const std::string& base_file_name, - Database::formats db_format = Database::HDF5); + Database::formats db_format = Database::formats::HDF5); /** * @brief Destructor. diff --git a/lib/linalg/svd/SVD.h b/lib/linalg/svd/SVD.h index 7ddbc21b8..c8fcddc70 100644 --- a/lib/linalg/svd/SVD.h +++ b/lib/linalg/svd/SVD.h @@ -143,7 +143,7 @@ class SVD /** * @brief Dimension of the system. */ - int d_dim; + const int d_dim; /** * @brief Number of samples stored for the current time interval. diff --git a/lib/linalg/svd/StaticSVD.cpp b/lib/linalg/svd/StaticSVD.cpp index 6eff9b9a8..a926a2e1a 100644 --- a/lib/linalg/svd/StaticSVD.cpp +++ b/lib/linalg/svd/StaticSVD.cpp @@ -220,7 +220,7 @@ StaticSVD::getSnapshotMatrix() " To preserve the snapshots, set Options::static_svd_preserve_snapshot to be true!\n"); if (d_snapshots) delete d_snapshots; - d_snapshots = new Matrix(d_dim, d_num_samples, false); + d_snapshots = new Matrix(d_dim, d_num_samples, true); for (int rank = 0; rank < d_num_procs; ++rank) { int nrows = d_dims[static_cast(rank)]; diff --git a/lib/utils/CSVDatabase.cpp b/lib/utils/CSVDatabase.cpp index 99547d6a7..263447594 100644 --- a/lib/utils/CSVDatabase.cpp +++ b/lib/utils/CSVDatabase.cpp @@ -28,18 +28,20 @@ CSVDatabase::~CSVDatabase() bool CSVDatabase::create( - const std::string& file_name) + const std::string& file_name, + const MPI_Comm comm) { - Database::create(file_name); + Database::create(file_name, comm); return true; } bool CSVDatabase::open( const std::string& file_name, - const std::string& type) + const std::string& type, + const MPI_Comm comm) { - Database::open(file_name, type); + Database::open(file_name, type, comm); return true; } @@ -53,7 +55,8 @@ void CSVDatabase::putIntegerArray( const std::string& file_name, const int* const data, - int nelements) + int nelements, + const bool distributed) { CAROM_VERIFY(!file_name.empty()); CAROM_VERIFY(data != nullptr); @@ -71,7 +74,8 @@ void CSVDatabase::putDoubleArray( const std::string& file_name, const double* const data, - int nelements) + int nelements, + const bool distributed) { CAROM_VERIFY(!file_name.empty()); CAROM_VERIFY(data != nullptr); @@ -90,7 +94,8 @@ void CSVDatabase::putDoubleVector( const std::string& file_name, const std::vector& data, - int nelements) + int nelements, + const bool distributed) { CAROM_VERIFY(!file_name.empty()); CAROM_VERIFY(nelements > 0); @@ -146,7 +151,8 @@ void CSVDatabase::getIntegerArray( const std::string& file_name, int* data, - int nelements) + int nelements, + const bool distributed) { CAROM_VERIFY(!file_name.empty()); #ifndef DEBUG_CHECK_ASSERTIONS @@ -188,7 +194,8 @@ void CSVDatabase::getDoubleArray( const std::string& file_name, double* data, - int nelements) + int nelements, + const bool distributed) { CAROM_VERIFY(!file_name.empty()); #ifndef DEBUG_CHECK_ASSERTIONS @@ -211,7 +218,8 @@ CSVDatabase::getDoubleArray( const std::string& file_name, double* data, int nelements, - const std::vector& idx) + const std::vector& idx, + const bool distributed) { CAROM_VERIFY(!file_name.empty()); #ifndef DEBUG_CHECK_ASSERTIONS @@ -252,7 +260,8 @@ CSVDatabase::getDoubleArray( int nelements, int offset, int block_size, - int stride) + int stride, + const bool distributed) { CAROM_VERIFY(!file_name.empty()); #ifndef DEBUG_CHECK_ASSERTIONS diff --git a/lib/utils/CSVDatabase.h b/lib/utils/CSVDatabase.h index 9eefa1f8a..c9132a15f 100644 --- a/lib/utils/CSVDatabase.h +++ b/lib/utils/CSVDatabase.h @@ -39,36 +39,43 @@ class CSVDatabase : public Database /** * @brief Creates a new CSV database file with the supplied name. + * NOTE: CSVDatabase does not actually create a file with this. + * This function will only print out the file_name. * * @param[in] file_name Name of CSV database file to create. + * @param[in] comm MPI communicator for distributed data I/O. + * CSVDatabase I/O is always performed serially, regardless. * * @return True if file create was successful. */ - virtual bool create( - const std::string& file_name) override; + const std::string& file_name, + const MPI_Comm comm=MPI_COMM_NULL) override; /** * @brief Opens an existing CSV database file with the supplied name. + * NOTE: CSVDatabase does not actually open a file with this. + * This function will only print out the file_name. * * @param[in] file_name Name of existing CSV database file to open. * @param[in] type Read/write type ("r"/"wr") + * @param[in] comm MPI communicator for distributed data I/O. + * CSVDatabase I/O is always performed serially, regardless. * * @return True if file open was successful. */ - virtual bool open( const std::string& file_name, - const std::string& type) override; + const std::string& type, + const MPI_Comm comm=MPI_COMM_NULL) override; /** * @brief Closes the currently open CSV database file. * * @return True if the file close was successful. */ - virtual bool close(); @@ -83,13 +90,15 @@ class CSVDatabase : public Database * written. * @param[in] data The array of integer values to be written. * @param[in] nelements The number of integers in the array. + * @param[in] distributed True if data is a distributed integer array. + * CSVDatabase writes the array serially whether or not distributed. */ - virtual void putIntegerArray( const std::string& file_name, const int* const data, - int nelements); + int nelements, + const bool distributed=false) override; /** * @brief Writes an array of doubles associated with the supplied filename. @@ -102,13 +111,15 @@ class CSVDatabase : public Database * written. * @param[in] data The array of double values to be written. * @param[in] nelements The number of doubles in the array. + * @param[in] distributed True if data is a distributed double array. + * CSVDatabase writes the array serially whether or not distributed. */ - virtual void putDoubleArray( const std::string& file_name, const double* const data, - int nelements); + int nelements, + const bool distributed=false) override; /** * @brief Writes a vector of doubles associated with the supplied filename to @@ -122,13 +133,15 @@ class CSVDatabase : public Database * written. * @param[in] data The vector of double values to be written. * @param[in] nelements The number of doubles in the vector. + * @param[in] distributed True if data is a distributed double vector. + * CSVDatabase writes the vector serially whether or not distributed. */ - virtual void putDoubleVector( const std::string& file_name, const std::vector& data, - int nelements); + int nelements, + const bool distributed=false) override; /** * @brief Writes a vector of complex doubles associated with the supplied filename. @@ -142,7 +155,6 @@ class CSVDatabase : public Database * @param[in] data The vector of complex double values to be written. * @param[in] nelements The number of complex doubles in the vector. */ - virtual void putComplexVector( const std::string& file_name, @@ -161,7 +173,6 @@ class CSVDatabase : public Database * @param[in] data The vector of strings to be written. * @param[in] nelements The number of strings in the vector. */ - virtual void putStringVector( const std::string& file_name, @@ -178,13 +189,15 @@ class CSVDatabase : public Database * read. * @param[out] data The allocated array of integer values to be read. * @param[in] nelements The number of integers in the array. + * @param[in] distributed True if data is a distributed integer array. + * CSVDatabase reads the array serially whether or not distributed. */ - virtual void getIntegerArray( const std::string& file_name, int* data, - int nelements); + int nelements, + const bool distributed=false) override; /** * @brief Reads a vector of integers associated with the supplied filename. @@ -211,9 +224,8 @@ class CSVDatabase : public Database * @param[in] file_name The filename associated with the array of values to be * read. */ - virtual int - getDoubleArraySize(const std::string& file_name) + getDoubleArraySize(const std::string& file_name) override { std::vector tmp; getDoubleVector(file_name, tmp); @@ -230,13 +242,15 @@ class CSVDatabase : public Database * read. * @param[out] data The allocated array of double values to be read. * @param[in] nelements The number of doubles in the array. + * @param[in] distributed True if data is a distributed double array. + * CSVDatabase reads the array serially whether or not distributed. */ - virtual void getDoubleArray( const std::string& file_name, double* data, - int nelements); + int nelements, + const bool distributed=false) override; /** * @brief Reads a sub-array of doubles associated with the supplied filename. @@ -249,14 +263,16 @@ class CSVDatabase : public Database * @param[out] data The allocated sub-array of double values to be read. * @param[in] nelements The number of doubles in the full array. * @param[in] idx The set of indices in the sub-array. + * @param[in] distributed True if data is a distributed integer array. + * CSVDatabase reads the array serially whether or not distributed. */ - virtual void getDoubleArray( const std::string& file_name, double* data, int nelements, - const std::vector& idx); + const std::vector& idx, + const bool distributed=false) override; /** * @brief Reads an array of doubles associated with the supplied filename. @@ -271,8 +287,9 @@ class CSVDatabase : public Database * @param[in] offset The initial offset in the array. * @param[in] block_size The block size to read from the CSV dataset. * @param[in] stride The stride to read from the CSV dataset. + * @param[in] distributed True if data is a distributed integer array. + * CSVDatabase reads the array serially whether or not distributed. */ - virtual void getDoubleArray( const std::string& file_name, @@ -280,7 +297,8 @@ class CSVDatabase : public Database int nelements, int offset, int block_size, - int stride); + int stride, + const bool distributed=false) override; /** * @brief Reads a vector of doubles associated with the supplied filename. diff --git a/lib/utils/Database.cpp b/lib/utils/Database.cpp index 49c6c9eb7..df2cb01b6 100644 --- a/lib/utils/Database.cpp +++ b/lib/utils/Database.cpp @@ -32,17 +32,22 @@ Database::~Database() } bool -Database::create(const std::string& file_name) +Database::create( + const std::string& file_name, + const MPI_Comm comm) { std::cout << "Creating file: " << file_name << std::endl; + return true; } bool Database::open( const std::string& file_name, - const std::string& type) + const std::string& type, + const MPI_Comm comm) { std::cout << "Opening file: " << file_name << std::endl; + return true; } } diff --git a/lib/utils/Database.h b/lib/utils/Database.h index e839da863..1008f391d 100644 --- a/lib/utils/Database.h +++ b/lib/utils/Database.h @@ -13,8 +13,10 @@ #ifndef included_Database_h #define included_Database_h +#include "Utilities.h" #include #include +#include "mpi.h" namespace CAROM { @@ -43,19 +45,26 @@ class Database * @brief Creates a new database file with the supplied name. * * @param[in] file_name Name of database file to create. + * @param[in] comm MPI communicator for distributed data I/O. + * If MPI_COMM_NULL, the file will be created serially. + * The behavior for distributed I/O will vary by classes. * * @return True if file create was successful. */ virtual bool create( - const std::string& file_name); + const std::string& file_name, + const MPI_Comm comm=MPI_COMM_NULL); /** * @brief Opens an existing database file with the supplied name. * * @param[in] file_name Name of existing database file to open. * @param[in] type Read/write type ("r"/"wr") + * @param[in] comm MPI communicator for distributed data I/O. + * If MPI_COMM_NULL, the file will be opened serially. + * The behavior for distributed I/O will vary by classes. * * @return True if file open was successful. */ @@ -63,7 +72,8 @@ class Database bool open( const std::string& file_name, - const std::string& type); + const std::string& type, + const MPI_Comm comm=MPI_COMM_NULL); /** * @brief Closes the currently open database file. @@ -97,13 +107,16 @@ class Database * written. * @param[in] data The array of integer values to be written. * @param[in] nelements The number of integers in the array. + * @param[in] distributed If true, distributed integer array will be written. + * the distributed I/O behavior varies with classes. */ virtual void putIntegerArray( const std::string& key, const int* const data, - int nelements) = 0; + int nelements, + const bool distributed=false) = 0; /** * @brief Writes a double associated with the supplied key to currently @@ -128,13 +141,16 @@ class Database * written. * @param[in] data The array of double values to be written. * @param[in] nelements The number of doubles in the array. + * @param[in] distributed If true, distributed double array will be written. + * the distributed I/O behavior varies with classes. */ virtual void putDoubleArray( const std::string& key, const double* const data, - int nelements) = 0; + int nelements, + const bool distributed=false) = 0; /** * @brief Writes a vector of doubles associated with the supplied key to @@ -144,13 +160,16 @@ class Database * written. * @param[in] data The vector of double values to be written. * @param[in] nelements The number of doubles in the vector. + * @param[in] distributed If true, distributed double vector will be written. + * the distributed I/O behavior varies with classes. */ virtual void putDoubleVector( const std::string& key, const std::vector& data, - int nelements) = 0; + int nelements, + const bool distributed=false) = 0; /** * @brief Reads an integer associated with the supplied key from the @@ -175,13 +194,16 @@ class Database * read. * @param[out] data The allocated array of integer values to be read. * @param[in] nelements The number of integers in the array. + * @param[in] distributed If true, distributed integer array will be read. + * the distributed I/O behavior varies with classes. */ virtual void getIntegerArray( const std::string& key, int* data, - int nelements) = 0; + int nelements, + const bool distributed=false) = 0; /** * @brief Reads a double associated with the supplied key from the @@ -220,13 +242,16 @@ class Database * read. * @param[out] data The allocated array of double values to be read. * @param[in] nelements The number of doubles in the array. + * @param[in] distributed If true, distributed double array will be read. + * the distributed I/O behavior varies with classes. */ virtual void getDoubleArray( const std::string& key, double* data, - int nelements) = 0; + int nelements, + const bool distributed=false) = 0; /** * @brief Reads a sub-array of doubles associated with the supplied key @@ -240,6 +265,8 @@ class Database * @param[out] data The allocated sub-array of double values to be read. * @param[in] nelements The number of doubles in the full array. * @param[in] idx The set of indices in the sub-array. + * @param[in] distributed If true, distributed double array will be read. + * the distributed I/O behavior varies with classes. */ virtual void @@ -247,7 +274,8 @@ class Database const std::string& key, double* data, int nelements, - const std::vector& idx) = 0; + const std::vector& idx, + const bool distributed=false) = 0; /** * @brief Reads an array of doubles associated with the supplied key @@ -258,8 +286,13 @@ class Database * @param[out] data The allocated array of double values to be read. * @param[in] nelements The number of doubles in the array. * @param[in] offset The initial offset in the array. + * Typically, this is a zero-based column index of the matrix data. * @param[in] block_size The block size to read from the HDF5 dataset. + * Typically, this is a number of columns of the matrix data. * @param[in] stride The stride to read from the HDF5 dataset. + * Typically, this is the total number of columns of the matrix data. + * @param[in] distributed If true, distributed double array will be read. + * the distributed I/O behavior varies with classes. */ virtual void @@ -269,15 +302,17 @@ class Database int nelements, int offset, int block_size, - int stride) = 0; + int stride, + const bool distributed=false) = 0; /** * @brief Implemented database file formats. Add to this enum each time a * new database format is implemented. */ - enum formats { + enum class formats { HDF5, - CSV + CSV, + HDF5_MPIO }; private: diff --git a/lib/utils/HDFDatabase.cpp b/lib/utils/HDFDatabase.cpp index 60bd0ad69..49e9136a8 100644 --- a/lib/utils/HDFDatabase.cpp +++ b/lib/utils/HDFDatabase.cpp @@ -43,12 +43,23 @@ HDFDatabase::~HDFDatabase() bool HDFDatabase::create( - const std::string& file_name) + const std::string& file_name, + const MPI_Comm comm) { - Database::create(file_name); - CAROM_VERIFY(!file_name.empty()); - hid_t file_id = H5Fcreate(file_name.c_str(), + std::string ext_filename(file_name); + + if (comm != MPI_COMM_NULL) + { + MPI_Comm_rank(comm, &d_rank); + char tmp[10]; + sprintf(tmp, ".%06d", d_rank); + ext_filename += tmp; + } + + Database::create(ext_filename, comm); + + hid_t file_id = H5Fcreate(ext_filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); @@ -64,22 +75,33 @@ HDFDatabase::create( bool HDFDatabase::open( const std::string& file_name, - const std::string& type) + const std::string& type, + const MPI_Comm comm) { - Database::open(file_name, type); - CAROM_VERIFY(!file_name.empty()); + std::string ext_filename(file_name); + + if (comm != MPI_COMM_NULL) + { + MPI_Comm_rank(comm, &d_rank); + char tmp[10]; + sprintf(tmp, ".%06d", d_rank); + ext_filename += tmp; + } + + Database::open(ext_filename, type, comm); + CAROM_VERIFY(type == "r" || type == "wr"); hid_t file_id; if (type == "r") { - file_id = H5Fopen(file_name.c_str(), + file_id = H5Fopen(ext_filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); } else if (type == "wr") { - file_id = H5Fopen(file_name.c_str(), + file_id = H5Fopen(ext_filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT); } @@ -114,7 +136,8 @@ void HDFDatabase::putIntegerArray( const std::string& key, const int* const data, - int nelements) + int nelements, + const bool distributed) { CAROM_VERIFY(!key.empty()); CAROM_VERIFY(data != nullptr); @@ -166,7 +189,8 @@ void HDFDatabase::putDoubleArray( const std::string& key, const double* const data, - int nelements) + int nelements, + const bool distributed) { CAROM_VERIFY(!key.empty()); CAROM_VERIFY(data != nullptr); @@ -218,7 +242,8 @@ void HDFDatabase::putDoubleVector( const std::string& key, const std::vector& data, - int nelements) + int nelements, + const bool distributed) { putDoubleArray(key, data.data(), nelements); } @@ -227,7 +252,8 @@ void HDFDatabase::getIntegerArray( const std::string& key, int* data, - int nelements) + int nelements, + const bool distributed) { if (nelements == 0) return; @@ -297,7 +323,8 @@ void HDFDatabase::getDoubleArray( const std::string& key, double* data, - int nelements) + int nelements, + const bool distributed) { if (nelements == 0) return; @@ -340,7 +367,8 @@ HDFDatabase::getDoubleArray( const std::string& key, double* data, int nelements, - const std::vector& idx) + const std::vector& idx, + const bool distributed) { if (idx.size() == 0) { @@ -373,7 +401,8 @@ HDFDatabase::getDoubleArray( int nelements, int offset, int block_size, - int stride) + int stride, + const bool distributed) { CAROM_VERIFY(!key.empty()); #ifndef DEBUG_CHECK_ASSERTIONS diff --git a/lib/utils/HDFDatabase.h b/lib/utils/HDFDatabase.h index d91e707d6..3faadcf60 100644 --- a/lib/utils/HDFDatabase.h +++ b/lib/utils/HDFDatabase.h @@ -14,6 +14,7 @@ #define included_HDFDatabase_h #include "Database.h" +#include "Utilities.h" #include "hdf5.h" #include @@ -40,19 +41,25 @@ class HDFDatabase : public Database * @brief Creates a new HDF5 database file with the supplied name. * * @param[in] file_name Name of HDF5 database file to create. + * @param[in] comm MPI communicator for distributed data I/O. + * If not MPI_COMM_NULL, each process creates a file + * whose name is file_name extended with process rank in 6 digits. * * @return True if file create was successful. */ virtual bool create( - const std::string& file_name) override; + const std::string& file_name, + const MPI_Comm comm=MPI_COMM_NULL) override; /** * @brief Opens an existing HDF5 database file with the supplied name. * * @param[in] file_name Name of existing HDF5 database file to open. * @param[in] type Read/write type ("r"/"wr") + * If not MPI_COMM_NULL, each process opens a file + * whose name is file_name extended with process rank in 6 digits. * * @return True if file open was successful. */ @@ -60,7 +67,8 @@ class HDFDatabase : public Database bool open( const std::string& file_name, - const std::string& type) override; + const std::string& type, + const MPI_Comm comm=MPI_COMM_NULL) override; /** * @brief Closes the currently open HDF5 database file. @@ -83,13 +91,17 @@ class HDFDatabase : public Database * written. * @param[in] data The array of integer values to be written. * @param[in] nelements The number of integers in the array. + * @param[in] distributed True if data is a distributed integer array. + * HDFDatabase writes the array in file-per-process, + * where each file is written serially by one process. */ virtual void putIntegerArray( const std::string& key, const int* const data, - int nelements); + int nelements, + const bool distributed=false); /** * @brief Writes an array of doubles associated with the supplied key to @@ -103,13 +115,17 @@ class HDFDatabase : public Database * written. * @param[in] data The array of double values to be written. * @param[in] nelements The number of doubles in the array. + * @param[in] distributed True if data is a distributed double array. + * HDFDatabase writes the array in file-per-process, + * where each file is written serially by one process. */ virtual void putDoubleArray( const std::string& key, const double* const data, - int nelements); + int nelements, + const bool distributed=false); /** * @brief Writes a vector of doubles associated with the supplied key to @@ -123,13 +139,17 @@ class HDFDatabase : public Database * written. * @param[in] data The vector of double values to be written. * @param[in] nelements The number of doubles in the vector. + * @param[in] distributed True if data is a distributed double array. + * HDFDatabase writes the array in file-per-process, + * where each file is written serially by one process. */ virtual void putDoubleVector( const std::string& key, const std::vector& data, - int nelements); + int nelements, + const bool distributed=false); /** * @brief Reads an array of integers associated with the supplied key @@ -142,13 +162,17 @@ class HDFDatabase : public Database * read. * @param[out] data The allocated array of integer values to be read. * @param[in] nelements The number of integers in the array. + * @param[in] distributed True if data is a distributed integer array. + * HDFDatabase reads the array in file-per-process, + * where each file is read serially by one process. */ virtual void getIntegerArray( const std::string& key, int* data, - int nelements); + int nelements, + const bool distributed=false); /** * @brief Count the number of elements in an array of doubles associated @@ -173,13 +197,17 @@ class HDFDatabase : public Database * read. * @param[out] data The allocated array of double values to be read. * @param[in] nelements The number of doubles in the array. + * @param[in] distributed True if data is a distributed double array. + * HDFDatabase reads the array in file-per-process, + * where each file is read serially by one process. */ virtual void getDoubleArray( const std::string& key, double* data, - int nelements); + int nelements, + const bool distributed=false); /** * @brief Reads a sub-array of doubles associated with the supplied key @@ -193,6 +221,9 @@ class HDFDatabase : public Database * @param[out] data The allocated sub-array of double values to be read. * @param[in] nelements The number of doubles in the full array. * @param[in] idx The set of indices in the sub-array. + * @param[in] distributed True if data is a distributed double array. + * HDFDatabase reads the array in file-per-process, + * where each file is read serially by one process. */ virtual void @@ -200,7 +231,8 @@ class HDFDatabase : public Database const std::string& key, double* data, int nelements, - const std::vector& idx); + const std::vector& idx, + const bool distributed=false); /** * @brief Reads an array of doubles associated with the supplied key @@ -214,8 +246,14 @@ class HDFDatabase : public Database * @param[out] data The allocated array of double values to be read. * @param[in] nelements The number of doubles in the array. * @param[in] offset The initial offset in the array. + * Typically, this is a column index of the matrix data. * @param[in] block_size The block size to read from the HDF5 dataset. + * Typically, this is a number of columns of the matrix data. * @param[in] stride The stride to read from the HDF5 dataset. + * Typically, this is the total number of columns of the matrix data. + * @param[in] distributed True if data is a distributed double array. + * HDFDatabase reads the array in file-per-process, + * where each file is read serially by one process. */ virtual void @@ -225,21 +263,10 @@ class HDFDatabase : public Database int nelements, int offset, int block_size, - int stride); + int stride, + const bool distributed=false); -private: - /** - * @brief Unimplemented copy constructor. - */ - HDFDatabase( - const HDFDatabase& other); - - /** - * @brief Unimplemented assignment operator. - */ - HDFDatabase& - operator = ( - const HDFDatabase& rhs); +protected: /** * @brief Returns true if the specified key represents an integer entry. @@ -273,6 +300,7 @@ class HDFDatabase : public Database * @param[in] type_key The attribute to be written. * @param[in] dataset_id ID of the dataset key will be written to. */ + virtual void writeAttribute( int type_key, @@ -313,6 +341,25 @@ class HDFDatabase : public Database * @brief The key representing an integer array. * */ static const int KEY_INT_ARRAY; + + /** + * @brief MPI process rank. + */ + int d_rank; + +private: + /** + * @brief Unimplemented copy constructor. + */ + HDFDatabase( + const HDFDatabase& other); + + /** + * @brief Unimplemented assignment operator. + */ + HDFDatabase& + operator = ( + const HDFDatabase& rhs); }; } diff --git a/lib/utils/HDFDatabaseMPIO.cpp b/lib/utils/HDFDatabaseMPIO.cpp new file mode 100644 index 000000000..5e2d9bc40 --- /dev/null +++ b/lib/utils/HDFDatabaseMPIO.cpp @@ -0,0 +1,638 @@ +/****************************************************************************** + * + * Copyright (c) 2013-2024, Lawrence Livermore National Security, LLC + * and other libROM project developers. See the top-level COPYRIGHT + * file for details. + * + * SPDX-License-Identifier: (Apache-2.0 OR MIT) + * + *****************************************************************************/ + +// Description: The concrete database implementation using HDF5. + +#include "HDFDatabaseMPIO.h" +#include "utils/mpi_utils.h" +#include + +namespace CAROM { + +HDFDatabaseMPIO::HDFDatabaseMPIO() : +#if HDF5_IS_PARALLEL + HDFDatabase(), + d_rank(-1), + d_comm(MPI_COMM_NULL) +#else + HDFDatabase() +#endif +{} + +#if HDF5_IS_PARALLEL + +bool +HDFDatabaseMPIO::create( + const std::string& file_name, + const MPI_Comm comm) +{ + int mpi_init; + MPI_Initialized(&mpi_init); + if (mpi_init) { + CAROM_VERIFY(comm != MPI_COMM_NULL); + d_comm = comm; + MPI_Comm_rank(d_comm, &d_rank); + } + else + d_rank = 0; + + /* Create a single file with a name compatible to HDFDatabase. */ + std::string file_name_ext(file_name + ".000000"); + if (d_rank == 0) + Database::create(file_name_ext); + CAROM_VERIFY(!file_name.empty()); + + hid_t plist_id; + /* + * Set up file access property list with parallel I/O access + */ + plist_id = H5Pcreate(H5P_FILE_ACCESS); + H5Pset_fapl_mpio(plist_id, d_comm, MPI_INFO_NULL); + /* + * OPTIONAL: It is generally recommended to set collective + * metadata reads/writes on FAPL to perform metadata reads + * collectively, which usually allows datasets + * to perform better at scale, although it is not + * strictly necessary. + */ + H5Pset_all_coll_metadata_ops(plist_id, true); + H5Pset_coll_metadata_write(plist_id, true); + + hid_t file_id = H5Fcreate(file_name_ext.c_str(), + H5F_ACC_TRUNC, + H5P_DEFAULT, + plist_id); + bool result = file_id >= 0; + CAROM_VERIFY(result); + d_is_file = true; + d_file_id = file_id; + d_group_id = file_id; + + H5Pclose(plist_id); + + return result; +} + +bool +HDFDatabaseMPIO::open( + const std::string& file_name, + const std::string& type, + const MPI_Comm comm) +{ + int mpi_init; + MPI_Initialized(&mpi_init); + if (mpi_init) { + CAROM_VERIFY(comm != MPI_COMM_NULL); + d_comm = comm; + MPI_Comm_rank(d_comm, &d_rank); + } + else + d_rank = 0; + + std::string file_name_ext(file_name + ".000000"); + if (d_rank == 0) + Database::open(file_name_ext, type); + CAROM_VERIFY(!file_name.empty()); + CAROM_VERIFY(type == "r" || type == "wr"); + + hid_t plist_id; + /* + * Set up file access property list with parallel I/O access + */ + plist_id = H5Pcreate(H5P_FILE_ACCESS); + H5Pset_fapl_mpio(plist_id, d_comm, MPI_INFO_NULL); + /* + * OPTIONAL: It is generally recommended to set collective + * metadata reads/writes on FAPL to perform metadata reads + * collectively, which usually allows datasets + * to perform better at scale, although it is not + * strictly necessary. + */ + H5Pset_all_coll_metadata_ops(plist_id, true); + H5Pset_coll_metadata_write(plist_id, true); + + hid_t file_id; + if (type == "r") + { + file_id = H5Fopen(file_name_ext.c_str(), + H5F_ACC_RDONLY, + plist_id); + } + else if (type == "wr") + { + file_id = H5Fopen(file_name_ext.c_str(), + H5F_ACC_RDWR, + plist_id); + } + bool result = file_id >= 0; + CAROM_VERIFY(result); + d_is_file = true; + d_file_id = file_id; + d_group_id = file_id; + + H5Pclose(plist_id); + + return result; +} + +void +HDFDatabaseMPIO::putIntegerArray_parallel( + const std::string& key, + const int* const data, + int nelem_local) +{ + CAROM_VERIFY(!key.empty()); + CAROM_VERIFY(data != nullptr); + CAROM_VERIFY(nelem_local >= 0); + + /* determine global nelements and offsets */ + std::vector offsets; + int nelements = CAROM::get_global_offsets(nelem_local, offsets, d_comm); + CAROM_VERIFY(nelements > 0); + + const int dim_rank = 1; + hsize_t dim[dim_rank] = { static_cast(nelements) }; + hid_t filespace = H5Screate_simple(dim_rank, dim, NULL); + CAROM_VERIFY(filespace >= 0); + +#if (H5_VERS_MAJOR > 1) || ((H5_VERS_MAJOR == 1) && (H5_VERS_MINOR > 6)) + hid_t dataset = H5Dcreate(d_group_id, + key.c_str(), + H5T_NATIVE_INT, + filespace, + H5P_DEFAULT, + H5P_DEFAULT, + H5P_DEFAULT); +#else + CAROM_ERROR("HDFDatabaseMPIO is not compatible with current HDF5 version!\n"); +#endif + CAROM_VERIFY(dataset >= 0); + H5Sclose(filespace); + + /* + * Each process defines dataset in memory and writes it to the hyperslab + * in the file. + */ + /* hyperslab selection parameters */ + hsize_t count[dim_rank]; + hsize_t offset[dim_rank]; + count[0] = nelem_local; + offset[0] = offsets[d_rank]; + hid_t memspace = H5Screate_simple(dim_rank, count, NULL); + + /* + * Select hyperslab in the file. + */ + filespace = H5Dget_space(dataset); + if (nelem_local > 0) + H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL); + else + H5Sselect_none(filespace); + + /* + * Create property list for collective dataset write. + */ + hid_t plist_id = H5Pcreate(H5P_DATASET_XFER); + H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE); + + herr_t errf = H5Dwrite(dataset, + H5T_NATIVE_INT, + memspace, + filespace, + plist_id, + data); + CAROM_VERIFY(errf >= 0); + + // Write attribute so we know what kind of data this is. + writeAttribute(KEY_INT_ARRAY, dataset); + + errf = H5Sclose(filespace); + CAROM_VERIFY(errf >= 0); + + errf = H5Sclose(memspace); + CAROM_VERIFY(errf >= 0); + + errf = H5Pclose(plist_id); + CAROM_VERIFY(errf >= 0); + + errf = H5Dclose(dataset); + CAROM_VERIFY(errf >= 0); +#ifndef DEBUG_CHECK_ASSERTIONS + CAROM_NULL_USE(errf); +#endif +} + +void +HDFDatabaseMPIO::putDoubleArray_parallel( + const std::string& key, + const double* const data, + int nelem_local) +{ + CAROM_VERIFY(!key.empty()); + CAROM_VERIFY(data != nullptr); + CAROM_VERIFY(nelem_local >= 0); + + /* determine global nelements and offsets */ + std::vector offsets; + int nelements = CAROM::get_global_offsets(nelem_local, offsets, d_comm); + CAROM_VERIFY(nelements > 0); + + const int dim_rank = 1; + hsize_t dim[dim_rank] = { static_cast(nelements) }; + hid_t filespace = H5Screate_simple(dim_rank, dim, 0); + CAROM_VERIFY(filespace >= 0); + +#if (H5_VERS_MAJOR > 1) || ((H5_VERS_MAJOR == 1) && (H5_VERS_MINOR > 6)) + hid_t dataset = H5Dcreate(d_group_id, + key.c_str(), + H5T_NATIVE_DOUBLE, + filespace, + H5P_DEFAULT, + H5P_DEFAULT, + H5P_DEFAULT); +#else + CAROM_ERROR("HDFDatabaseMPIO is not compatible with current HDF5 version!\n"); +#endif + CAROM_VERIFY(dataset >= 0); + H5Sclose(filespace); + + /* + * Each process defines dataset in memory and writes it to the hyperslab + * in the file. + */ + /* hyperslab selection parameters */ + hsize_t count[dim_rank]; + hsize_t offset[dim_rank]; + count[0] = nelem_local; + offset[0] = offsets[d_rank]; + hid_t memspace = H5Screate_simple(dim_rank, count, NULL); + + /* + * Select hyperslab in the file. + */ + filespace = H5Dget_space(dataset); + if (nelem_local > 0) + H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL); + else + H5Sselect_none(filespace); + + /* + * Create property list for collective dataset write. + */ + hid_t plist_id = H5Pcreate(H5P_DATASET_XFER); + H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE); + + herr_t errf = H5Dwrite(dataset, + H5T_NATIVE_DOUBLE, + memspace, + filespace, + plist_id, + data); + CAROM_VERIFY(errf >= 0); + + // Write attribute so we know what kind of data this is. + writeAttribute(KEY_DOUBLE_ARRAY, dataset); + + errf = H5Sclose(filespace); + CAROM_VERIFY(errf >= 0); + + errf = H5Sclose(memspace); + CAROM_VERIFY(errf >= 0); + + errf = H5Pclose(plist_id); + CAROM_VERIFY(errf >= 0); + + errf = H5Dclose(dataset); + CAROM_VERIFY(errf >= 0); +#ifndef DEBUG_CHECK_ASSERTIONS + CAROM_NULL_USE(errf); +#endif +} + +void +HDFDatabaseMPIO::getIntegerArray_parallel( + const std::string& key, + int* data, + int nelem_local) +{ + CAROM_VERIFY(nelem_local >= 0); + /* determine global nelements and offsets */ + std::vector offsets; + int nelements = CAROM::get_global_offsets(nelem_local, offsets, d_comm); + if (nelements == 0) return; + + CAROM_VERIFY(!key.empty()); +#ifndef DEBUG_CHECK_ASSERTIONS + CAROM_NULL_USE(nelem_local); +#endif + +#if (H5_VERS_MAJOR > 1) || ((H5_VERS_MAJOR == 1) && (H5_VERS_MINOR > 6)) + hid_t dset = H5Dopen(d_group_id, key.c_str(), H5P_DEFAULT); +#else + CAROM_ERROR("HDFDatabaseMPIO is not compatible with current HDF5 version!\n"); +#endif + CAROM_VERIFY(dset >= 0); + + hid_t filespace = H5Dget_space(dset); + CAROM_VERIFY(filespace >= 0); + + hsize_t nsel = H5Sget_select_npoints(filespace); + CAROM_VERIFY(static_cast(nsel) == nelements); + + const int dim_rank = 1; + H5Sclose(filespace); + + /* + * Each process defines dataset in memory and writes it to the hyperslab + * in the file. + */ + /* hyperslab selection parameters */ + hsize_t count[dim_rank]; + hsize_t offset[dim_rank]; + count[0] = nelem_local; + offset[0] = offsets[d_rank]; + hid_t memspace = H5Screate_simple(dim_rank, count, NULL); + + /* + * Select hyperslab in the file. + */ + filespace = H5Dget_space(dset); + if (nelem_local > 0) + H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL); + else + H5Sselect_none(filespace); + + /* + * Create property list for collective dataset write. + */ + hid_t plist_id = H5Pcreate(H5P_DATASET_XFER); + H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE); + + herr_t errf; + errf = H5Dread(dset, H5T_NATIVE_INT, memspace, filespace, plist_id, data); + CAROM_VERIFY(errf >= 0); + + errf = H5Sclose(filespace); + CAROM_VERIFY(errf >= 0); + + errf = H5Sclose(memspace); + CAROM_VERIFY(errf >= 0); + + errf = H5Pclose(plist_id); + CAROM_VERIFY(errf >= 0); + + errf = H5Dclose(dset); + CAROM_VERIFY(errf >= 0); +#ifndef DEBUG_CHECK_ASSERTIONS + CAROM_NULL_USE(errf); +#endif +} + +void +HDFDatabaseMPIO::getDoubleArray_parallel( + const std::string& key, + double* data, + int nelem_local) +{ + CAROM_VERIFY(nelem_local >= 0); + /* determine global nelements and offsets */ + std::vector offsets; + int nelements = CAROM::get_global_offsets(nelem_local, offsets, d_comm); + if (nelements == 0) return; + + CAROM_VERIFY(!key.empty()); +#ifndef DEBUG_CHECK_ASSERTIONS + CAROM_NULL_USE(nelements); +#endif + +#if (H5_VERS_MAJOR > 1) || ((H5_VERS_MAJOR == 1) && (H5_VERS_MINOR > 6)) + hid_t dset = H5Dopen(d_group_id, key.c_str(), H5P_DEFAULT); +#else + hid_t dset = H5Dopen(d_group_id, key.c_str()); +#endif + CAROM_VERIFY(dset >= 0); + + hid_t filespace = H5Dget_space(dset); + CAROM_VERIFY(filespace >= 0); + + hsize_t nsel = H5Sget_select_npoints(filespace); + CAROM_VERIFY(static_cast(nsel) == nelements); + + const int dim_rank = 1; + H5Sclose(filespace); + + /* + * Each process defines dataset in memory and writes it to the hyperslab + * in the file. + */ + /* hyperslab selection parameters */ + hsize_t count[dim_rank]; + hsize_t offset[dim_rank]; + count[0] = nelem_local; + offset[0] = offsets[d_rank]; + hid_t memspace = H5Screate_simple(dim_rank, count, NULL); + + /* + * Select hyperslab in the file. + */ + filespace = H5Dget_space(dset); + if (nelem_local > 0) + H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL); + else + H5Sselect_none(filespace); + + /* + * Create property list for collective dataset write. + */ + hid_t plist_id = H5Pcreate(H5P_DATASET_XFER); + H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE); + + herr_t errf; + errf = H5Dread(dset, H5T_NATIVE_DOUBLE, memspace, filespace, plist_id, data); + CAROM_VERIFY(errf >= 0); + + errf = H5Sclose(filespace); + CAROM_VERIFY(errf >= 0); + + errf = H5Sclose(memspace); + CAROM_VERIFY(errf >= 0); + + errf = H5Pclose(plist_id); + CAROM_VERIFY(errf >= 0); + + errf = H5Dclose(dset); + CAROM_VERIFY(errf >= 0); +#ifndef DEBUG_CHECK_ASSERTIONS + CAROM_NULL_USE(errf); +#endif +} + +void +HDFDatabaseMPIO::getDoubleArray_parallel( + const std::string& key, + double* data, + int nelem_local, + const std::vector& idx_local) +{ + std::vector idx_tmp = idx_local; + if (idx_local.size() == 0) + { + for (int i = 0; i < nelem_local; i++) + idx_tmp.push_back(i); + } + + std::vector alldata(nelem_local); + getDoubleArray_parallel(key, alldata.data(), nelem_local); + int k = 0; + for (int i = 0; i < nelem_local; ++i) + { + if (idx_tmp[k] == i) + { + data[k++] = alldata[i]; + } + if (k == idx_tmp.size()) + { + break; + } + } + CAROM_VERIFY(k == idx_tmp.size()); +} + +void +HDFDatabaseMPIO::getDoubleArray_parallel( + const std::string& key, + double* data, + int nelem_local, + int block_offset_global, + int block_size_global, + int stride_global) +{ + CAROM_VERIFY(nelem_local >= 0); + CAROM_VERIFY(CAROM::is_same(stride_global, d_comm)); + CAROM_VERIFY(CAROM::is_same(block_size_global, d_comm)); + CAROM_VERIFY(CAROM::is_same(block_offset_global, d_comm)); + /* We assume stride sets the maximum block size */ + CAROM_VERIFY(block_offset_global + block_size_global <= stride_global); + /* determine global nelements and offsets */ + hsize_t num_local_blocks = nelem_local / block_size_global; + std::vector global_offsets; + int nelements = CAROM::get_global_offsets(num_local_blocks * stride_global, + global_offsets, d_comm); + + CAROM_VERIFY(!key.empty()); +#ifndef DEBUG_CHECK_ASSERTIONS + CAROM_NULL_USE(nelements); +#endif + +#if (H5_VERS_MAJOR > 1) || ((H5_VERS_MAJOR == 1) && (H5_VERS_MINOR > 6)) + hid_t dset = H5Dopen(d_group_id, key.c_str(), H5P_DEFAULT); +#else + hid_t dset = H5Dopen(d_group_id, key.c_str()); +#endif + CAROM_VERIFY(dset >= 0); + + hid_t filespace = H5Dget_space(dset); + CAROM_VERIFY(filespace >= 0); + + hsize_t nsel = H5Sget_select_npoints(filespace); + CAROM_VERIFY((nsel == 0) || (nsel == nelements)); + + const int dim_rank = 1; + H5Sclose(filespace); + + herr_t errf; + if (nsel > 0) { + /* + * Each process defines dataset in memory and writes it to the hyperslab + * in the file. + */ + /* hyperslab selection parameters */ + hsize_t num_blocks[dim_rank] = {num_local_blocks}; + hsize_t buffer_array_size[dim_rank] = {(hsize_t) nelem_local}; + hsize_t offset[dim_rank] = {(hsize_t) global_offsets[d_rank] + block_offset_global}; + hsize_t strides[dim_rank] = {(hsize_t) stride_global}; + hsize_t block_sizes[1] = {(hsize_t) block_size_global}; + hid_t memspace = H5Screate_simple(dim_rank, buffer_array_size, NULL); + + /* + * Select hyperslab in the file. + */ + filespace = H5Dget_space(dset); + if (nelem_local > 0) + H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, + strides, num_blocks, block_sizes); + else + H5Sselect_none(filespace); + + /* + * Create property list for collective dataset write. + */ + hid_t plist_id = H5Pcreate(H5P_DATASET_XFER); + H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE); + + errf = H5Dread(dset, H5T_NATIVE_DOUBLE, memspace, filespace, plist_id, data); + CAROM_VERIFY(errf >= 0); + + errf = H5Sclose(memspace); + CAROM_VERIFY(errf >= 0); + + errf = H5Pclose(plist_id); + CAROM_VERIFY(errf >= 0); + } + + errf = H5Sclose(filespace); + CAROM_VERIFY(errf >= 0); + + errf = H5Dclose(dset); + CAROM_VERIFY(errf >= 0); +#ifndef DEBUG_CHECK_ASSERTIONS + CAROM_NULL_USE(errf); +#endif +} + +void +HDFDatabaseMPIO::writeAttribute( + int type_key, + hid_t dataset_id) +{ + CAROM_VERIFY(CAROM::is_same(type_key, d_comm)); + + hid_t attr_id = H5Screate(H5S_SCALAR); + CAROM_VERIFY(attr_id >= 0); + +#if (H5_VERS_MAJOR > 1) || ((H5_VERS_MAJOR == 1) && (H5_VERS_MINOR > 6)) + hid_t attr = H5Acreate(dataset_id, + "Type", + H5T_NATIVE_INT, + attr_id, + H5P_DEFAULT, + H5P_DEFAULT); +#else + CAROM_ERROR("HDFDatabaseMPIO is not compatible with current HDF5 version!\n"); +#endif + CAROM_VERIFY(attr >= 0); + + /* only the last process writes the attribute */ + herr_t errf = H5Awrite(attr, H5T_NATIVE_INT, &type_key); + CAROM_VERIFY(errf >= 0); + + errf = H5Aclose(attr); + CAROM_VERIFY(errf >= 0); + + errf = H5Sclose(attr_id); + CAROM_VERIFY(errf >= 0); +#ifndef DEBUG_CHECK_ASSERTIONS + CAROM_NULL_USE(errf); +#endif +} + +#endif + +} diff --git a/lib/utils/HDFDatabaseMPIO.h b/lib/utils/HDFDatabaseMPIO.h new file mode 100644 index 000000000..6686a95fe --- /dev/null +++ b/lib/utils/HDFDatabaseMPIO.h @@ -0,0 +1,425 @@ +/****************************************************************************** + * + * Copyright (c) 2013-2024, Lawrence Livermore National Security, LLC + * and other libROM project developers. See the top-level COPYRIGHT + * file for details. + * + * SPDX-License-Identifier: (Apache-2.0 OR MIT) + * + *****************************************************************************/ + +// Description: The concrete database implementation using HDF5. + +#ifndef included_HDFDatabaseMPI_h +#define included_HDFDatabaseMPI_h + +#include "HDFDatabase.h" + +namespace CAROM { + +/** + * HDFDatabaseMPIO implements Multi-Path Input/Output (MPIO) of HDF5. + */ +class HDFDatabaseMPIO : public HDFDatabase +{ +public: + /** + * @brief Default constructor. + */ + HDFDatabaseMPIO(); + + /** + * @brief Destructor. + */ + virtual + ~HDFDatabaseMPIO() {} + +#if HDF5_IS_PARALLEL + /** + * @brief Creates a new HDF5 database file for distributed data + * with the supplied name. + * + * @param[in] file_name Name of HDF5 database file to create. + * @param[in] comm MPI communicator for distributed data I/O. + * By default, HDFDatabaseMPIO uses MPI_COMM_WORLD + * and does not allow MPI_COMM_NULL. + * + * @return True if file create was successful. + */ + bool + create( + const std::string& file_name, + const MPI_Comm comm=MPI_COMM_WORLD) override; + + /** + * @brief Opens an existing HDF5 database file for distributed data + * with the supplied name. + * + * @param[in] file_name Name of existing HDF5 database file to open. + * @param[in] type Read/write type ("r"/"wr") + * @param[in] comm MPI communicator for distributed data I/O. + * By default, HDFDatabaseMPIO uses MPI_COMM_WORLD + * and does not allow MPI_COMM_NULL. + * + * @return True if file open was successful. + */ + bool + open( + const std::string& file_name, + const std::string& type, + const MPI_Comm comm=MPI_COMM_WORLD) override; + + /** + * @brief Writes a local array of integers only for root rank, + * with associated supplied key to + * the currently open HDF5 database file. + * + * @pre !key.empty() + * @pre data != nullptr + * @pre nelements > 0 + * + * @param[in] key The key associated with the array of values to be + * written. + * @param[in] data The array of integer values to be written. + * @param[in] nelements The number of integers in the array. + * @param[in] distributed If true, distributed integer array will be written. + * If not, only the root process writes its integer array. + */ + void + putIntegerArray( + const std::string& key, + const int* const data, + int nelements, + const bool distributed=false) override + { + if ((!distributed) && (d_rank != 0)) + nelements = 0; + putIntegerArray_parallel(key, data, nelements); + } + + /** + * @brief Writes an array of doubles in the root rank + * associated with the supplied key to + * the currently open HDF5 database file. + * + * @pre !key.empty() + * @pre data != nullptr + * @pre nelements > 0 + * + * @param[in] key The key associated with the array of values to be + * written. + * @param[in] data The array of double values to be written. + * @param[in] nelements The number of doubles in the array. + * @param[in] distributed If true, distributed double array will be written. + * If not, only the root process writes its double array. + */ + void + putDoubleArray( + const std::string& key, + const double* const data, + int nelements, + const bool distributed=false) override + { + if ((!distributed) && (d_rank != 0)) + nelements = 0; + putDoubleArray_parallel(key, data, nelements); + } + + /** + * @brief Reads an array of integers associated with the supplied key + * from the currently open HDF5 database file. + * All processes share the same non-distributed integer array. + * + * @pre !key.empty() + * @pre data != nullptr || nelements == 0 + * + * @param[in] key The key associated with the array of values to be + * read. + * @param[out] data The allocated array of integer values to be read. + * @param[in] nelements The number of integers in the array. + * @param[in] distributed If true, the integer array will be read in a distributed way. + * If not, the root process reads the entire array and broadcast to all processes. + */ + void + getIntegerArray( + const std::string& key, + int* data, + int nelements, + const bool distributed=false) override + { + if (distributed) + { + getIntegerArray_parallel(key, data, nelements); + return; + } + + int read_size = (d_rank == 0) ? nelements : 0; + getIntegerArray_parallel(key, data, read_size); + + CAROM_VERIFY(d_comm != MPI_COMM_NULL); + MPI_Bcast(data, nelements, MPI_INT, 0, d_comm); + } + + /** + * @brief Reads an array of doubles associated with the supplied key + * from the currently open HDF5 database file. + * All processes share the same non-distributed double array. + * + * @pre !key.empty() + * @pre data != nullptr || nelements == 0 + * + * @param[in] key The key associated with the array of values to be + * read. + * @param[out] data The allocated array of double values to be read. + * @param[in] nelements The number of doubles in the array. + * @param[in] distributed If true, the double array will be read in a distributed way. + * If not, the root process reads the entire array and broadcast to all processes. + */ + void + getDoubleArray( + const std::string& key, + double* data, + int nelements, + const bool distributed=false) override + { + if (distributed) + { + getDoubleArray_parallel(key, data, nelements); + return; + } + + int read_size = (d_rank == 0) ? nelements : 0; + getDoubleArray_parallel(key, data, read_size); + + CAROM_VERIFY(d_comm != MPI_COMM_NULL); + MPI_Bcast(data, nelements, MPI_DOUBLE, 0, d_comm); + } + + /** + * @brief Reads a sub-array of doubles associated with the supplied key + * from the currently open HDF5 database file. + * All processes share the same non-distributed double array. + * + * @pre !key.empty() + * @pre data != nullptr || nelements == 0 + * + * @param[in] key The key associated with the array of values to be + * read. + * @param[out] data The allocated sub-array of double values to be read. + * @param[in] nelements The number of doubles in the full array. + * @param[in] idx The set of indices in the sub-array. + * @param[in] distributed If true, the double array will be read in a distributed way. + * If not, the root process reads the entire array and broadcast to all processes. + */ + void + getDoubleArray( + const std::string& key, + double* data, + int nelements, + const std::vector& idx, + const bool distributed=false) override + { + if (distributed) + { + getDoubleArray_parallel(key, data, nelements, idx); + return; + } + + int read_size = (d_rank == 0) ? nelements : 0; + getDoubleArray_parallel(key, data, read_size, idx); + + CAROM_VERIFY(d_comm != MPI_COMM_NULL); + MPI_Bcast(data, nelements, MPI_DOUBLE, 0, d_comm); + } + + /** + * @brief Reads an array of doubles associated with the supplied key + * from the currently open HDF5 database file. + * All processes share the same non-distributed double array. + * + * @pre !key.empty() + * @pre data != nullptr || nelements == 0 + * + * @param[in] key The key associated with the array of values to be + * read. + * @param[out] data The allocated array of double values to be read. + * @param[in] nelements The number of doubles in the array. + * @param[in] offset The initial offset in the array. + * @param[in] block_size The block size to read from the HDF5 dataset. + * @param[in] stride The stride to read from the HDF5 dataset. + * @param[in] distributed If true, the double array will be read in a distributed way. + * If not, the root process reads the entire array and broadcast to all processes. + */ + void + getDoubleArray( + const std::string& key, + double* data, + int nelements, + int offset, + int block_size, + int stride, + const bool distributed=false) override + { + if (distributed) + { + getDoubleArray_parallel(key, data, nelements, offset, block_size, stride); + return; + } + + int read_size = (d_rank == 0) ? nelements : 0; + getDoubleArray_parallel(key, data, read_size, offset, block_size, stride); + + CAROM_VERIFY(d_comm != MPI_COMM_NULL); + MPI_Bcast(data, nelements, MPI_DOUBLE, 0, d_comm); + } + + void + writeAttribute( + int type_key, + hid_t dataset_id) override; + +private: + MPI_Comm d_comm; + int d_rank; + + /** + * @brief Writes a distributed array of integers + * associated with the supplied key + * to the currently open HDF5 database file. + * + * @pre !key.empty() + * @pre data != nullptr + * @pre nelem_local >= 0 + * @pre nelements > 0 + * + * @param[in] key The key associated with the array of values to be + * written. + * @param[in] data The array of integer values to be written. + * @param[in] nelem_local The local number of integers in the array. + */ + virtual + void + putIntegerArray_parallel( + const std::string& key, + const int* const data, + int nelem_local); + + /** + * @brief Writes a distributed array of doubles + * associated with the supplied key to + * the currently open HDF5 database file. + * + * @pre !key.empty() + * @pre data != nullptr + * @pre nelem_local >= 0 + * @pre nelements > 0 + * + * @param[in] key The key associated with the array of values to be + * written. + * @param[in] data The array of double values to be written. + * @param[in] nelem_local The number of doubles in the array. + */ + virtual + void + putDoubleArray_parallel( + const std::string& key, + const double* const data, + int nelem_local); + + /** + * @brief Reads a distributed array of integers + * associated with the supplied key + * from the currently open HDF5 database file. + * + * @pre !key.empty() + * @pre data != nullptr || nelements == 0 + * + * @param[in] key The key associated with the array of values to be + * read. + * @param[out] data The allocated array of integer values to be read. + * @param[in] nelem_local The local number of integers in the array. + */ + virtual + void + getIntegerArray_parallel( + const std::string& key, + int* data, + int nelem_local); + + /** + * @brief Reads a distributed array of doubles + * associated with the supplied key + * from the currently open HDF5 database file. + * + * @pre !key.empty() + * @pre data != nullptr || nelements == 0 + * + * @param[in] key The key associated with the array of values to be + * read. + * @param[out] data The allocated array of double values to be read. + * @param[in] nelem_local The number of doubles in the array. + */ + virtual + void + getDoubleArray_parallel( + const std::string& key, + double* data, + int nelem_local); + + /** + * @brief Reads a distributed sub-array of doubles + * associated with the supplied key + * from the currently open HDF5 database file. + * + * @pre !key.empty() + * @pre data != nullptr || nelements == 0 + * + * @param[in] key The key associated with the array of values to be + * read. + * @param[out] data The allocated sub-array of double values to be read. + * @param[in] nelem_local The number of doubles in the full array. + * @param[in] idx_local The set of indices in the sub-array. + */ + virtual + void + getDoubleArray_parallel( + const std::string& key, + double* data, + int nelem_local, + const std::vector& idx_local); + + /** + * @brief Reads a distributed array of doubles + * associated with the supplied key + * from the currently open HDF5 database file. + * + * @pre !key.empty() + * @pre data != nullptr || nelements == 0 + * @pre block_offset_global + block_size_global <= stride_global + * + * @param[in] key The key associated with the array of values to be + * read. + * @param[out] data The allocated array of double values to be read. + * @param[in] nelem_local The number of doubles in the array at the local process. + * @param[in] block_offset_global The initial offset in the global array. + * Typically, this is a global column index of the matrix data. + * @param[in] block_size_global The total block size to read from the HDF5 dataset. + * Typically, this is a number of columns (in global) of the matrix data. + * @param[in] stride_global The global stride to read from the HDF5 dataset. + * Typically, this is the total number of columns of the matrix data. + */ + virtual + void + getDoubleArray_parallel( + const std::string& key, + double* data, + int nelem_local, + int block_offset_global, + int block_size_global, + int stride_global); +#endif +}; + +} + +#endif diff --git a/lib/utils/mpi_utils.cpp b/lib/utils/mpi_utils.cpp index a92d23094..7670bd5ac 100644 --- a/lib/utils/mpi_utils.cpp +++ b/lib/utils/mpi_utils.cpp @@ -68,4 +68,11 @@ get_global_offsets(const int local_dim, std::vector &offsets, return dim; } +bool +is_same(int x, const MPI_Comm &comm) { + int p[2] = {-x, x}; + MPI_Allreduce(MPI_IN_PLACE, p, 2, MPI_INT, MPI_MIN, comm); + return (p[0] == -p[1]); +} + } diff --git a/lib/utils/mpi_utils.h b/lib/utils/mpi_utils.h index bb2a02a3b..fdfd5e142 100644 --- a/lib/utils/mpi_utils.h +++ b/lib/utils/mpi_utils.h @@ -43,6 +43,15 @@ int get_global_offsets(const int local_dim, std::vector &offsets, const MPI_Comm &comm=MPI_COMM_WORLD); +/** + * @brief Check if an integer is equal over all MPI processes. + * + * @param[in] x Input integer value to test equality. + * @param[in] comm MPI communicator. default value MPI_COMM_WORLD. + */ +bool +is_same(int x, const MPI_Comm &comm=MPI_COMM_WORLD); + } #endif diff --git a/unit_tests/test_HDFDatabase.cpp b/unit_tests/test_HDFDatabase.cpp new file mode 100644 index 000000000..68b43beb7 --- /dev/null +++ b/unit_tests/test_HDFDatabase.cpp @@ -0,0 +1,871 @@ +/****************************************************************************** + * + * Copyright (c) 2013-2024, Lawrence Livermore National Security, LLC + * and other libROM project developers. See the top-level COPYRIGHT + * file for details. + * + * SPDX-License-Identifier: (Apache-2.0 OR MIT) + * + *****************************************************************************/ + +#ifdef CAROM_HAS_GTEST + +#include +#include "linalg/BasisGenerator.h" +#include "utils/HDFDatabase.h" +#include +#include +#include +#include // for memcpy +#include +#include "mpi.h" +#include "utils/mpi_utils.h" +#include +using namespace std::chrono; + +static const int nrow = 123, ncol = 21; +static const double threshold = 1.0e-13; + +/** + * Simple smoke test to make sure Google Test is properly linked + */ +TEST(GoogleTestFramework, GoogleTestFrameworkFound) { + SUCCEED(); +} + +#if HDF5_IS_PARALLEL + +TEST(HDF5, Test_parallel_writing) +{ + int nproc, rank; + MPI_Comm_size(MPI_COMM_WORLD, &nproc); + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + + const int dim_rank = 2; + const int nrow_local = CAROM::split_dimension(nrow, MPI_COMM_WORLD); + std::vector offsets; + int dummy = CAROM::get_global_offsets(nrow_local, offsets, MPI_COMM_WORLD); + + /* + * Initialize data buffer + */ + int data[nrow_local * ncol]; + for (int d = 0; d < nrow_local * ncol; d++) + data[d] = d + offsets[rank] * ncol; + + MPI_Barrier(MPI_COMM_WORLD); + const auto start = steady_clock::now(); + + hid_t plist_id; + hid_t file_id; + herr_t errf = 0; + + /* + * Set up file access property list with parallel I/O access + */ + plist_id = H5Pcreate(H5P_FILE_ACCESS); + H5Pset_fapl_mpio(plist_id, MPI_COMM_WORLD, MPI_INFO_NULL); + + /* + * OPTIONAL: It is generally recommended to set collective + * metadata reads/writes on FAPL to perform metadata reads + * collectively, which usually allows datasets + * to perform better at scale, although it is not + * strictly necessary. + */ + H5Pset_all_coll_metadata_ops(plist_id, true); + H5Pset_coll_metadata_write(plist_id, true); + + /* + * Create a new file collectively and release property list identifier. + */ + file_id = H5Fcreate("test1.h5", H5F_ACC_TRUNC, H5P_DEFAULT, plist_id); + H5Pclose(plist_id); + + /* + * Create the dataspace for the dataset. + */ + hsize_t dim_global[dim_rank] = { static_cast(nrow), static_cast(ncol) }; + hid_t filespace = H5Screate_simple(dim_rank, dim_global, NULL); + + /* + * Create the dataset with default properties and close filespace. + */ + hid_t dset_id = + H5Dcreate(file_id, "test-int", H5T_NATIVE_INT, filespace, H5P_DEFAULT, + H5P_DEFAULT, H5P_DEFAULT); + H5Sclose(filespace); + + /* + * Each process defines dataset in memory and writes it to the hyperslab + * in the file. + */ + /* hyperslab selection parameters */ + hsize_t count[dim_rank]; + hsize_t offset[dim_rank]; + count[0] = nrow_local; + count[1] = ncol; + offset[0] = offsets[rank]; + offset[1] = 0; + hid_t memspace = H5Screate_simple(dim_rank, count, NULL); + + /* + * Select hyperslab in the file. + */ + filespace = H5Dget_space(dset_id); + H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL); + + /* + * Create property list for collective dataset write. + */ + plist_id = H5Pcreate(H5P_DATASET_XFER); + H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE); + + /* + * Write the data collectively. + */ + errf = H5Dwrite(dset_id, H5T_NATIVE_INT, memspace, filespace, plist_id, data); + CAROM_VERIFY(errf >= 0); + + /* + * Write a integer scalar attribute + */ + // Only one process attribute writes the attribute. + // Make sure all processes have the same value. + const int test_attr = 1234; + hid_t attr_id = H5Screate(H5S_SCALAR); + hid_t attr = H5Acreate(dset_id, "test_attr", H5T_NATIVE_INT, + attr_id, H5P_DEFAULT, H5P_DEFAULT); + + errf = H5Awrite(attr, H5T_NATIVE_INT, &test_attr); + errf = H5Aclose(attr); + errf = H5Sclose(attr_id); + MPI_Barrier(MPI_COMM_WORLD); + + /* + * Close/release resources. + */ + H5Dclose(dset_id); + H5Sclose(filespace); + H5Sclose(memspace); + H5Pclose(plist_id); + H5Fclose(file_id); + + MPI_Barrier(MPI_COMM_WORLD); + const auto stop = steady_clock::now(); + const auto duration = duration_cast(stop-start); + printf("rank: %d, duration: %dms\n", rank, duration); +} + +TEST(HDF5, Test_parallel_reading) +{ + int nproc, rank; + MPI_Comm_size(MPI_COMM_WORLD, &nproc); + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + + const int dim_rank = 2; + const int nrow_local = CAROM::split_dimension(nrow, MPI_COMM_WORLD); + std::vector offsets; + int dummy = CAROM::get_global_offsets(nrow_local, offsets, MPI_COMM_WORLD); + + /* + * Initialize data buffer + */ + int data[nrow_local * ncol]; + + MPI_Barrier(MPI_COMM_WORLD); + const auto start = steady_clock::now(); + + hid_t plist_id; + hid_t file_id; + herr_t errf = 0; + + /* + * Set up file access property list with parallel I/O access + */ + plist_id = H5Pcreate(H5P_FILE_ACCESS); + H5Pset_fapl_mpio(plist_id, MPI_COMM_WORLD, MPI_INFO_NULL); + + /* + * OPTIONAL: It is generally recommended to set collective + * metadata reads/writes on FAPL to perform metadata reads + * collectively, which usually allows datasets + * to perform better at scale, although it is not + * strictly necessary. + */ + H5Pset_all_coll_metadata_ops(plist_id, true); + H5Pset_coll_metadata_write(plist_id, true); + + /* + * Open a new file collectively and release property list identifier. + */ + file_id = H5Fopen("test1.h5", H5F_ACC_RDONLY, plist_id); + H5Pclose(plist_id); + + /* + * Open the dataset with default properties. + */ + hid_t dset_id = H5Dopen(file_id, "test-int", H5P_DEFAULT); + + /* + * Get filespace and read the global size. + */ + hid_t filespace = H5Dget_space(dset_id); + CAROM_VERIFY(filespace >= 0); + int ndims = H5Sget_simple_extent_ndims(filespace); + CAROM_VERIFY(ndims == dim_rank); + hsize_t dim_global[dim_rank]; + errf = H5Sget_simple_extent_dims(filespace, dim_global, NULL); + CAROM_VERIFY(errf >= 0); + CAROM_VERIFY(dim_global[0] == nrow); + CAROM_VERIFY(dim_global[1] == ncol); + + H5Sclose(filespace); + + /* + * Each process defines dataset in memory and writes it to the hyperslab + * in the file. + */ + /* hyperslab selection parameters */ + hsize_t count[dim_rank]; + hsize_t offset[dim_rank]; + count[0] = nrow_local; + count[1] = ncol; + offset[0] = offsets[rank]; + offset[1] = 0; + hid_t memspace = H5Screate_simple(dim_rank, count, NULL); + + /* + * Select hyperslab in the file. + */ + filespace = H5Dget_space(dset_id); + H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL); + + /* + * Create property list for collective dataset write. + */ + plist_id = H5Pcreate(H5P_DATASET_XFER); + H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE); + + /* + * Read the data collectively. + */ + errf = H5Dread(dset_id, H5T_NATIVE_INT, memspace, filespace, plist_id, data); + CAROM_VERIFY(errf >= 0); + + /* + * Write a integer scalar attribute + */ + int test_attr = -1; + hid_t attr = H5Aopen_name(dset_id, "test_attr"); + errf = H5Aread(attr, H5T_NATIVE_INT, &test_attr); + // MPI_Barrier(MPI_COMM_WORLD); + errf = H5Aclose(attr); + if (rank == 0) + EXPECT_EQ(test_attr, 1234); + + /* + * Close/release resources. + */ + H5Dclose(dset_id); + H5Sclose(filespace); + H5Sclose(memspace); + H5Pclose(plist_id); + H5Fclose(file_id); + + MPI_Barrier(MPI_COMM_WORLD); + const auto stop = steady_clock::now(); + const auto duration = duration_cast(stop-start); + printf("rank: %d, duration: %dms\n", rank, duration); + + for (int d = 0; d < nrow_local * ncol; d++) + EXPECT_TRUE(data[d] == d + offsets[rank] * ncol); +} + +TEST(HDF5, Test_selective_parallel_writing) +{ + int nproc, rank; + MPI_Comm_size(MPI_COMM_WORLD, &nproc); + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + + /* Only the first 2 processes will participate in I/O. */ + const int ioproc = std::min(2, nproc); + int iorank[ioproc]; + for (int r = 0; r < ioproc; r++) + iorank[r] = r; + MPI_Group world_grp, io_grp; + MPI_Comm_group(MPI_COMM_WORLD, &world_grp); + MPI_Group_incl(world_grp, ioproc, iorank, &io_grp); + MPI_Comm io_comm; + MPI_Comm_create(MPI_COMM_WORLD, io_grp, &io_comm); + + const int dim_rank = 2; + int nrow_local = 0; + if (rank < ioproc) + nrow_local = CAROM::split_dimension(nrow, io_comm); + printf("rank %d, I/O row size: %d\n", rank, nrow_local); + std::vector offsets; + int dummy = CAROM::get_global_offsets(nrow_local, offsets, MPI_COMM_WORLD); + for (int k = offsets.size(); k < nproc; k++) + offsets.push_back(offsets.back()); + + /* + * Initialize data buffer + */ + int data[nrow_local * ncol]; + for (int d = 0; d < nrow_local * ncol; d++) + data[d] = d + offsets[rank] * ncol; + + MPI_Barrier(MPI_COMM_WORLD); + const auto start = steady_clock::now(); + + hid_t plist_id; + hid_t file_id; + herr_t errf = 0; + + /* + * Set up file access property list with parallel I/O access + */ + plist_id = H5Pcreate(H5P_FILE_ACCESS); + H5Pset_fapl_mpio(plist_id, MPI_COMM_WORLD, MPI_INFO_NULL); + + /* + * OPTIONAL: It is generally recommended to set collective + * metadata reads/writes on FAPL to perform metadata reads + * collectively, which usually allows datasets + * to perform better at scale, although it is not + * strictly necessary. + */ + H5Pset_all_coll_metadata_ops(plist_id, true); + H5Pset_coll_metadata_write(plist_id, true); + + /* + * Create a new file collectively and release property list identifier. + */ + file_id = H5Fcreate("test.h5", H5F_ACC_TRUNC, H5P_DEFAULT, plist_id); + H5Pclose(plist_id); + + /* + * Create the dataspace for the dataset. + */ + hsize_t dim_global[dim_rank] = { static_cast(nrow), static_cast(ncol) }; + hid_t filespace = H5Screate_simple(dim_rank, dim_global, NULL); + + /* + * Create the dataset with default properties and close filespace. + */ + hid_t dset_id = + H5Dcreate(file_id, "test-int", H5T_NATIVE_INT, filespace, H5P_DEFAULT, + H5P_DEFAULT, H5P_DEFAULT); + H5Sclose(filespace); + + /* + * Each process defines dataset in memory and writes it to the hyperslab + * in the file. + */ + /* hyperslab selection parameters */ + hsize_t count[dim_rank]; + hsize_t offset[dim_rank]; + count[0] = nrow_local; + count[1] = ncol; + offset[0] = offsets[rank]; + offset[1] = 0; + hid_t memspace = H5Screate_simple(dim_rank, count, NULL); + + /* + * Select hyperslab in the file. + */ + filespace = H5Dget_space(dset_id); + if (rank < ioproc) + H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL); + else + H5Sselect_none(filespace); + + /* + * Create property list for collective dataset write. + */ + plist_id = H5Pcreate(H5P_DATASET_XFER); + H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE); + + /* + * Write the data collectively. + */ + errf = H5Dwrite(dset_id, H5T_NATIVE_INT, memspace, filespace, plist_id, data); + CAROM_VERIFY(errf >= 0); + + /* + * Close/release resources. + */ + H5Dclose(dset_id); + H5Sclose(filespace); + H5Sclose(memspace); + H5Pclose(plist_id); + H5Fclose(file_id); + + MPI_Barrier(MPI_COMM_WORLD); + const auto stop = steady_clock::now(); + const auto duration = duration_cast(stop-start); + printf("rank: %d, duration: %dms\n", rank, duration); +} + +TEST(HDF5, Test_selective_parallel_reading) +{ + int nproc, rank; + MPI_Comm_size(MPI_COMM_WORLD, &nproc); + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + + /* Only the first 2 processes will participate in I/O. */ + const int ioproc = std::min(2, nproc); + int iorank[ioproc]; + for (int r = 0; r < ioproc; r++) + iorank[r] = r; + MPI_Group world_grp, io_grp; + MPI_Comm_group(MPI_COMM_WORLD, &world_grp); + MPI_Group_incl(world_grp, ioproc, iorank, &io_grp); + MPI_Comm io_comm; + MPI_Comm_create(MPI_COMM_WORLD, io_grp, &io_comm); + + const int dim_rank = 2; + int nrow_local = 0; + if (rank < ioproc) + nrow_local = CAROM::split_dimension(nrow, io_comm); + printf("rank %d, I/O row size: %d\n", rank, nrow_local); + std::vector offsets; + int dummy = CAROM::get_global_offsets(nrow_local, offsets, MPI_COMM_WORLD); + + /* + * Initialize data buffer + */ + int data[nrow_local * ncol]; + + MPI_Barrier(MPI_COMM_WORLD); + const auto start = steady_clock::now(); + + hid_t plist_id; + hid_t file_id; + herr_t errf = 0; + + /* + * Set up file access property list with parallel I/O access + */ + plist_id = H5Pcreate(H5P_FILE_ACCESS); + H5Pset_fapl_mpio(plist_id, MPI_COMM_WORLD, MPI_INFO_NULL); + + /* + * OPTIONAL: It is generally recommended to set collective + * metadata reads/writes on FAPL to perform metadata reads + * collectively, which usually allows datasets + * to perform better at scale, although it is not + * strictly necessary. + */ + H5Pset_all_coll_metadata_ops(plist_id, true); + H5Pset_coll_metadata_write(plist_id, true); + + /* + * Open a new file collectively and release property list identifier. + */ + file_id = H5Fopen("test.h5", H5F_ACC_RDONLY, plist_id); + H5Pclose(plist_id); + + /* + * Open the dataset with default properties. + */ + hid_t dset_id = H5Dopen(file_id, "test-int", H5P_DEFAULT); + + /* + * Get filespace and read the global size. + */ + hid_t filespace = H5Dget_space(dset_id); + CAROM_VERIFY(filespace >= 0); + int ndims = H5Sget_simple_extent_ndims(filespace); + CAROM_VERIFY(ndims == dim_rank); + hsize_t dim_global[dim_rank]; + errf = H5Sget_simple_extent_dims(filespace, dim_global, NULL); + CAROM_VERIFY(errf >= 0); + CAROM_VERIFY(dim_global[0] == nrow); + CAROM_VERIFY(dim_global[1] == ncol); + + H5Sclose(filespace); + + /* + * Each process defines dataset in memory and writes it to the hyperslab + * in the file. + */ + /* hyperslab selection parameters */ + hsize_t count[dim_rank]; + hsize_t offset[dim_rank]; + count[0] = nrow_local; + count[1] = ncol; + offset[0] = offsets[rank]; + offset[1] = 0; + hid_t memspace = H5Screate_simple(dim_rank, count, NULL); + + /* + * Select hyperslab in the file. + */ + filespace = H5Dget_space(dset_id); + if (rank < ioproc) + H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL); + else + H5Sselect_none(filespace); + + /* + * Create property list for collective dataset write. + */ + plist_id = H5Pcreate(H5P_DATASET_XFER); + H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE); + + /* + * Read the data collectively. + */ + errf = H5Dread(dset_id, H5T_NATIVE_INT, memspace, filespace, plist_id, data); + CAROM_VERIFY(errf >= 0); + + /* + * Close/release resources. + */ + H5Dclose(dset_id); + H5Sclose(filespace); + H5Sclose(memspace); + H5Pclose(plist_id); + H5Fclose(file_id); + + MPI_Barrier(MPI_COMM_WORLD); + const auto stop = steady_clock::now(); + const auto duration = duration_cast(stop-start); + printf("rank: %d, duration: %dms\n", rank, duration); + + for (int d = 0; d < nrow_local * ncol; d++) + EXPECT_TRUE(data[d] == d + offsets[rank] * ncol); +} + +#endif + +TEST(BasisGeneratorIO, HDFDatabase) +{ + // Get the rank of this process, and the number of processors. + int mpi_init, d_rank, d_num_procs; + MPI_Initialized(&mpi_init); + if (mpi_init == 0) { + MPI_Init(nullptr, nullptr); + } + + MPI_Comm_rank(MPI_COMM_WORLD, &d_rank); + MPI_Comm_size(MPI_COMM_WORLD, &d_num_procs); + + int nrow_local = CAROM::split_dimension(nrow, MPI_COMM_WORLD); + std::vector row_offset(d_num_procs + 1); + const int dummy = CAROM::get_global_offsets(nrow_local, row_offset, + MPI_COMM_WORLD); + EXPECT_EQ(nrow, dummy); + + std::default_random_engine generator; + generator.seed( + 1234); // fix the seed to keep the same result for different nproc. + std::uniform_real_distribution<> uniform_distribution(0.0, 1.0); + std::normal_distribution normal_distribution(0.0, 1.0); + + // distribute from a global matrix to keep the same system for different nproc. + CAROM::Matrix snapshots(nrow, ncol, false); + for (int i = 0; i < nrow; i++) + for (int j = 0; j < ncol; j++) + snapshots(i, j) = normal_distribution(generator); + snapshots.distribute(nrow_local); + + CAROM::Options svd_options = CAROM::Options(nrow_local, ncol, 1); + svd_options.setMaxBasisDimension(nrow); + svd_options.setRandomizedSVD(false); + svd_options.setDebugMode(true); + CAROM::BasisGenerator sampler(svd_options, false, "test_basis"); + CAROM::Vector sample(nrow_local, true); + for (int s = 0; s < ncol; s++) + { + for (int d = 0; d < nrow_local; d++) + sample(d) = snapshots(d, s); + + sampler.takeSample(sample.getData()); + } + + sampler.writeSnapshot(); + const CAROM::Matrix *snapshot = sampler.getSnapshotMatrix(); + + sampler.endSamples(); + const CAROM::Matrix *spatial_basis = sampler.getSpatialBasis(); + + CAROM::BasisReader basis_reader("test_basis"); + const CAROM::Matrix *spatial_basis1 = basis_reader.getSpatialBasis(); + + EXPECT_EQ(spatial_basis->numRows(), spatial_basis1->numRows()); + EXPECT_EQ(spatial_basis->numColumns(), spatial_basis1->numColumns()); + for (int i = 0; i < spatial_basis->numRows(); i++) + for (int j = 0; j < spatial_basis->numColumns(); j++) + EXPECT_NEAR((*spatial_basis)(i, j), (*spatial_basis1)(i, j), threshold); + + CAROM::BasisReader snapshot_reader("test_basis_snapshot"); + const CAROM::Matrix *snapshot1 = snapshot_reader.getSnapshotMatrix(); + + EXPECT_EQ(snapshot->numRows(), snapshot1->numRows()); + EXPECT_EQ(snapshot->numColumns(), snapshot1->numColumns()); + for (int i = 0; i < snapshot->numRows(); i++) + for (int j = 0; j < snapshot->numColumns(); j++) + EXPECT_NEAR((*snapshot)(i, j), (*snapshot1)(i, j), threshold); +} + +TEST(BasisGeneratorIO, Base_MPIO_combination) +{ + // Get the rank of this process, and the number of processors. + int mpi_init, d_rank, d_num_procs; + MPI_Initialized(&mpi_init); + if (mpi_init == 0) { + MPI_Init(nullptr, nullptr); + } + + MPI_Comm_rank(MPI_COMM_WORLD, &d_rank); + MPI_Comm_size(MPI_COMM_WORLD, &d_num_procs); + + int nrow_local = CAROM::split_dimension(nrow, MPI_COMM_WORLD); + std::vector row_offset(d_num_procs + 1); + const int dummy = CAROM::get_global_offsets(nrow_local, row_offset, + MPI_COMM_WORLD); + EXPECT_EQ(nrow, dummy); + + std::string base_name = "test_basis"; + std::string mpio_name = "test_mpio"; + CAROM::Options svd_options = CAROM::Options(nrow_local, ncol, 1); + svd_options.setMaxBasisDimension(nrow); + svd_options.setRandomizedSVD(false); + svd_options.setDebugMode(true); + CAROM::BasisGenerator sampler(svd_options, false, mpio_name, + CAROM::Database::formats::HDF5_MPIO); + + sampler.loadSamples(base_name + "_snapshot", "snapshot", 1e9, + CAROM::Database::formats::HDF5); + sampler.writeSnapshot(); + const CAROM::Matrix *snapshot = sampler.getSnapshotMatrix(); + + CAROM::BasisReader snapshot_reader("test_basis_snapshot"); + const CAROM::Matrix *snapshot1 = snapshot_reader.getSnapshotMatrix(); + + EXPECT_EQ(snapshot->numRows(), snapshot1->numRows()); + EXPECT_EQ(snapshot->numColumns(), snapshot1->numColumns()); + for (int i = 0; i < snapshot->numRows(); i++) + for (int j = 0; j < snapshot->numColumns(); j++) + EXPECT_NEAR((*snapshot)(i, j), (*snapshot1)(i, j), threshold); + + sampler.endSamples(); + const CAROM::Matrix *spatial_basis = sampler.getSpatialBasis(); + + CAROM::BasisReader basis_reader("test_basis"); + const CAROM::Matrix *spatial_basis1 = basis_reader.getSpatialBasis(); + + EXPECT_EQ(spatial_basis->numRows(), spatial_basis1->numRows()); + EXPECT_EQ(spatial_basis->numColumns(), spatial_basis1->numColumns()); + for (int i = 0; i < spatial_basis->numRows(); i++) + for (int j = 0; j < spatial_basis->numColumns(); j++) + EXPECT_NEAR((*spatial_basis)(i, j), (*spatial_basis1)(i, j), threshold); +} + +TEST(BasisGeneratorIO, MPIO_Base_combination) +{ + // Get the rank of this process, and the number of processors. + int mpi_init, d_rank, d_num_procs; + MPI_Initialized(&mpi_init); + if (mpi_init == 0) { + MPI_Init(nullptr, nullptr); + } + + MPI_Comm_rank(MPI_COMM_WORLD, &d_rank); + MPI_Comm_size(MPI_COMM_WORLD, &d_num_procs); + + int nrow_local = CAROM::split_dimension(nrow, MPI_COMM_WORLD); + std::vector row_offset(d_num_procs + 1); + const int dummy = CAROM::get_global_offsets(nrow_local, row_offset, + MPI_COMM_WORLD); + EXPECT_EQ(nrow, dummy); + + std::string test_name = "test_basis2"; + std::string mpio_name = "test_mpio"; + CAROM::Options svd_options = CAROM::Options(nrow_local, ncol, 1); + svd_options.setMaxBasisDimension(nrow); + svd_options.setRandomizedSVD(false); + svd_options.setDebugMode(true); + CAROM::BasisGenerator sampler(svd_options, false, test_name, + CAROM::Database::formats::HDF5); + + sampler.loadSamples(mpio_name + "_snapshot", "snapshot", 1e9, + CAROM::Database::formats::HDF5_MPIO); + const CAROM::Matrix *snapshot = sampler.getSnapshotMatrix(); + + CAROM::BasisReader snapshot_reader("test_basis_snapshot"); + const CAROM::Matrix *snapshot1 = snapshot_reader.getSnapshotMatrix(); + + EXPECT_EQ(snapshot->numRows(), snapshot1->numRows()); + EXPECT_EQ(snapshot->numColumns(), snapshot1->numColumns()); + for (int i = 0; i < snapshot->numRows(); i++) + for (int j = 0; j < snapshot->numColumns(); j++) + EXPECT_NEAR((*snapshot)(i, j), (*snapshot1)(i, j), threshold); + + sampler.endSamples(); + const CAROM::Matrix *spatial_basis = sampler.getSpatialBasis(); + + CAROM::BasisReader basis_reader("test_basis"); + const CAROM::Matrix *spatial_basis1 = basis_reader.getSpatialBasis(); + + EXPECT_EQ(spatial_basis->numRows(), spatial_basis1->numRows()); + EXPECT_EQ(spatial_basis->numColumns(), spatial_basis1->numColumns()); + for (int i = 0; i < spatial_basis->numRows(); i++) + for (int j = 0; j < spatial_basis->numColumns(); j++) + EXPECT_NEAR((*spatial_basis)(i, j), (*spatial_basis1)(i, j), threshold); +} + +TEST(BasisReaderIO, partial_getSpatialBasis) +{ + // Get the rank of this process, and the number of processors. + int mpi_init, d_rank, d_num_procs; + MPI_Initialized(&mpi_init); + if (mpi_init == 0) { + MPI_Init(nullptr, nullptr); + } + + MPI_Comm_rank(MPI_COMM_WORLD, &d_rank); + MPI_Comm_size(MPI_COMM_WORLD, &d_num_procs); + + int nrow_local = CAROM::split_dimension(nrow, MPI_COMM_WORLD); + std::vector row_offset(d_num_procs + 1); + const int dummy = CAROM::get_global_offsets(nrow_local, row_offset, + MPI_COMM_WORLD); + EXPECT_EQ(nrow, dummy); + + std::string base_name = "test_basis"; + std::string mpio_name = "test_mpio"; + + CAROM::BasisReader basis_reader("test_basis"); + const CAROM::Matrix *spatial_basis = basis_reader.getSpatialBasis(); + + CAROM::BasisReader basis_reader1("test_mpio", + CAROM::Database::formats::HDF5_MPIO, + nrow_local); + const CAROM::Matrix *spatial_basis1 = basis_reader1.getSpatialBasis(); + + EXPECT_EQ(spatial_basis->numRows(), spatial_basis1->numRows()); + EXPECT_EQ(spatial_basis->numColumns(), spatial_basis1->numColumns()); + for (int i = 0; i < spatial_basis->numRows(); i++) + for (int j = 0; j < spatial_basis->numColumns(); j++) + EXPECT_NEAR((*spatial_basis)(i, j), (*spatial_basis1)(i, j), threshold); + + delete spatial_basis; + delete spatial_basis1; + const int col1 = 3, col2 = 7; + spatial_basis = basis_reader.getSpatialBasis(col1, col2); + spatial_basis1 = basis_reader1.getSpatialBasis(col1, col2); + + EXPECT_EQ(spatial_basis->numRows(), spatial_basis1->numRows()); + EXPECT_EQ(spatial_basis->numColumns(), spatial_basis1->numColumns()); + for (int i = 0; i < spatial_basis->numRows(); i++) + for (int j = 0; j < spatial_basis->numColumns(); j++) + EXPECT_NEAR((*spatial_basis)(i, j), (*spatial_basis1)(i, j), threshold); + + delete spatial_basis; + delete spatial_basis1; +} + +TEST(BasisGeneratorIO, Scaling_test) +{ + int nproc, rank; + MPI_Comm_size(MPI_COMM_WORLD, &nproc); + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + + const int scale_nrow = 500, scale_ncol = 100; + int nrow_local = CAROM::split_dimension(scale_nrow, MPI_COMM_WORLD); + std::vector row_offset(nproc + 1); + const int dummy = CAROM::get_global_offsets(nrow_local, row_offset, + MPI_COMM_WORLD); + EXPECT_EQ(scale_nrow, dummy); + + std::default_random_engine generator; + generator.seed( + 1234); // fix the seed to keep the same result for different nproc. + std::uniform_real_distribution<> uniform_distribution(0.0, 1.0); + std::normal_distribution normal_distribution(0.0, 1.0); + + // distribute from a global matrix to keep the same system for different nproc. + CAROM::Matrix snapshots(scale_nrow, scale_ncol, false); + for (int i = 0; i < scale_nrow; i++) + for (int j = 0; j < scale_ncol; j++) + snapshots(i, j) = normal_distribution(generator); + snapshots.distribute(nrow_local); + + CAROM::Options svd_options = CAROM::Options(nrow_local, scale_ncol, 1); + svd_options.setMaxBasisDimension(scale_nrow); + svd_options.setRandomizedSVD(false); + + CAROM::BasisGenerator base_sampler(svd_options, false, "base"); + CAROM::BasisGenerator mpio_sampler(svd_options, false, "mpio", + CAROM::Database::formats::HDF5_MPIO); + + CAROM::Vector sample(nrow_local, true); + for (int s = 0; s < scale_ncol; s++) + { + for (int d = 0; d < nrow_local; d++) + sample(d) = snapshots(d, s); + + base_sampler.takeSample(sample.getData()); + mpio_sampler.takeSample(sample.getData()); + } + + MPI_Barrier(MPI_COMM_WORLD); + auto start = steady_clock::now(); + base_sampler.writeSnapshot(); + MPI_Barrier(MPI_COMM_WORLD); + auto stop = steady_clock::now(); + auto duration = duration_cast(stop-start); + if (rank == 0) + printf("Base writeSnapshot- duration: %dms\n", duration); + + MPI_Barrier(MPI_COMM_WORLD); + start = steady_clock::now(); + mpio_sampler.writeSnapshot(); + MPI_Barrier(MPI_COMM_WORLD); + stop = steady_clock::now(); + duration = duration_cast(stop-start); + if (rank == 0) + printf("MPIO writeSnapshot- duration: %dms\n", duration); + + MPI_Barrier(MPI_COMM_WORLD); + start = steady_clock::now(); + base_sampler.endSamples(); + MPI_Barrier(MPI_COMM_WORLD); + stop = steady_clock::now(); + duration = duration_cast(stop-start); + if (rank == 0) + printf("Base endSamples- duration: %dms\n", duration); + + MPI_Barrier(MPI_COMM_WORLD); + start = steady_clock::now(); + mpio_sampler.endSamples(); + MPI_Barrier(MPI_COMM_WORLD); + stop = steady_clock::now(); + duration = duration_cast(stop-start); + if (rank == 0) + printf("MPIO endSamples- duration: %dms\n", duration); +} + +int main(int argc, char* argv[]) +{ + ::testing::InitGoogleTest(&argc, argv); + MPI_Init(&argc, &argv); + int result = RUN_ALL_TESTS(); + MPI_Finalize(); + return result; +} + +#else // #ifndef CAROM_HAS_GTEST + +int main() +{ + std::cout << "libROM was compiled without Google Test support, so unit " + << "tests have been disabled. To enable unit tests, compile " + << "libROM with Google Test support." << std::endl; +} + +#endif // #endif CAROM_HAS_GTEST