From 463bd216ddaf897df32eb0b2ef35086224039c43 Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Tue, 16 Jan 2024 18:09:42 -0800 Subject: [PATCH 01/14] BasisWriter::writeBasis create/open a file and closes the file at the end. --- lib/linalg/BasisReader.cpp | 1 - lib/linalg/BasisWriter.cpp | 45 +++++++++++++++++++------------------- lib/utils/CSVDatabase.cpp | 2 ++ lib/utils/CSVDatabase.h | 4 ++-- lib/utils/Database.cpp | 23 +++++++++++++++++++ lib/utils/Database.h | 6 +++-- lib/utils/HDFDatabase.cpp | 4 ++++ lib/utils/HDFDatabase.h | 4 ++-- 8 files changed, 59 insertions(+), 30 deletions(-) diff --git a/lib/linalg/BasisReader.cpp b/lib/linalg/BasisReader.cpp index 900f85f09..3b72c0185 100644 --- a/lib/linalg/BasisReader.cpp +++ b/lib/linalg/BasisReader.cpp @@ -48,7 +48,6 @@ BasisReader::BasisReader( d_database = new CSVDatabase(); } - std::cout << "Opening file: " << full_file_name << std::endl; d_database->open(full_file_name, "r"); int num_time_intervals; diff --git a/lib/linalg/BasisWriter.cpp b/lib/linalg/BasisWriter.cpp index 15924526e..e3e06ff6d 100644 --- a/lib/linalg/BasisWriter.cpp +++ b/lib/linalg/BasisWriter.cpp @@ -53,20 +53,19 @@ BasisWriter::BasisWriter( char tmp2[100]; sprintf(tmp2, "_snapshot.%06d", rank); snap_file_name = base_file_name + tmp2; + + // create and open snapshot/basis database + // TODO(kevin): can it be CSV at all? we might want to remove if statement here. + if (db_format_ == Database::HDF5) { + d_snap_database = new HDFDatabase(); + d_database = new HDFDatabase(); + } } BasisWriter::~BasisWriter() { - if (d_database) { - d_database->putInteger("num_time_intervals", d_num_intervals_written); - d_database->close(); - delete d_database; - } - if (d_snap_database) { - d_snap_database->putInteger("num_time_intervals", d_num_intervals_written); - d_snap_database->close(); - delete d_snap_database; - } + if (d_database) delete d_database; + if (d_snap_database) delete d_snap_database; } void @@ -81,13 +80,10 @@ BasisWriter::writeBasis(const std::string& kind) sprintf(tmp, "time_%06d", d_num_intervals_written); if (kind == "basis") { - - // create and open basis database - if (db_format_ == Database::HDF5) { - d_database = new HDFDatabase(); - } - std::cout << "Creating file: " << full_file_name << std::endl; - d_database->create(full_file_name); + if (fileExists(full_file_name)) + d_database->open(full_file_name, "wr"); + else + d_database->create(full_file_name); d_database->putDouble(tmp, time_interval_start_time); @@ -122,15 +118,15 @@ BasisWriter::writeBasis(const std::string& kind) ++d_num_intervals_written; + d_database->putInteger("num_time_intervals", d_num_intervals_written); + d_database->close(); } if (kind == "snapshot") { - // create and open snapshot database - if (db_format_ == Database::HDF5) { - d_snap_database = new HDFDatabase(); - } - std::cout << "Creating file: " << snap_file_name << std::endl; - d_snap_database->create(snap_file_name); + if (fileExists(snap_file_name)) + d_snap_database->open(snap_file_name, "wr"); + else + d_snap_database->create(snap_file_name); d_snap_database->putDouble(tmp, time_interval_start_time); @@ -143,6 +139,9 @@ BasisWriter::writeBasis(const std::string& kind) d_snap_database->putInteger(tmp, num_cols); sprintf(tmp, "snapshot_matrix_%06d", d_num_intervals_written); d_snap_database->putDoubleArray(tmp, &snapshots->item(0,0), num_rows*num_cols); + + d_snap_database->putInteger("num_time_intervals", d_num_intervals_written); + d_snap_database->close(); } } diff --git a/lib/utils/CSVDatabase.cpp b/lib/utils/CSVDatabase.cpp index 6b3a9752f..01eb9a525 100644 --- a/lib/utils/CSVDatabase.cpp +++ b/lib/utils/CSVDatabase.cpp @@ -30,6 +30,7 @@ bool CSVDatabase::create( const std::string& file_name) { + Database::create(file_name); return true; } @@ -38,6 +39,7 @@ CSVDatabase::open( const std::string& file_name, const std::string& type) { + Database::open(file_name, type); return true; } diff --git a/lib/utils/CSVDatabase.h b/lib/utils/CSVDatabase.h index aa3969a70..359d02cbe 100644 --- a/lib/utils/CSVDatabase.h +++ b/lib/utils/CSVDatabase.h @@ -47,7 +47,7 @@ class CSVDatabase : public Database virtual bool create( - const std::string& file_name); + const std::string& file_name) override; /** * @brief Opens an existing CSV database file with the supplied name. @@ -61,7 +61,7 @@ class CSVDatabase : public Database bool open( const std::string& file_name, - const std::string& type); + const std::string& type) override; /** * @brief Closes the currently open CSV database file. diff --git a/lib/utils/Database.cpp b/lib/utils/Database.cpp index ca9f2b9ca..4533afd38 100644 --- a/lib/utils/Database.cpp +++ b/lib/utils/Database.cpp @@ -11,9 +11,18 @@ // Description: The abstract database class defines interface for databases. #include "Database.h" +#include +#include namespace CAROM { +bool fileExists(const std::string& name) +{ + std::ifstream f(name.c_str()); + return f.good(); + // ifstream f will be closed upon the end of the function. +} + Database::Database() { } @@ -22,4 +31,18 @@ Database::~Database() { } +bool +Database::create(const std::string& file_name) +{ + std::cout << "Creating file: " << file_name << std::endl; +} + +bool +Database::open( + const std::string& file_name, + const std::string& type) +{ + std::cout << "Opening file: " << file_name << std::endl; +} + } diff --git a/lib/utils/Database.h b/lib/utils/Database.h index 9e7c39588..8ef6cbb06 100644 --- a/lib/utils/Database.h +++ b/lib/utils/Database.h @@ -18,6 +18,8 @@ namespace CAROM { +bool fileExists(const std::string& name); + /** * Class Database is an abstract base class that provides basic ability to * write to and read from a file. It's capabilities are limited to what the @@ -47,7 +49,7 @@ class Database virtual bool create( - const std::string& file_name) = 0; + const std::string& file_name); /** * @brief Opens an existing database file with the supplied name. @@ -61,7 +63,7 @@ class Database bool open( const std::string& file_name, - const std::string& type) = 0; + const std::string& type); /** * @brief Closes the currently open database file. diff --git a/lib/utils/HDFDatabase.cpp b/lib/utils/HDFDatabase.cpp index aa1b1a22e..4c00f7e40 100644 --- a/lib/utils/HDFDatabase.cpp +++ b/lib/utils/HDFDatabase.cpp @@ -44,6 +44,8 @@ bool HDFDatabase::create( const std::string& file_name) { + Database::create(file_name); + CAROM_VERIFY(!file_name.empty()); hid_t file_id = H5Fcreate(file_name.c_str(), H5F_ACC_TRUNC, @@ -63,6 +65,8 @@ HDFDatabase::open( const std::string& file_name, const std::string& type) { + Database::open(file_name, type); + CAROM_VERIFY(!file_name.empty()); CAROM_VERIFY(type == "r" || type == "wr"); hid_t file_id; diff --git a/lib/utils/HDFDatabase.h b/lib/utils/HDFDatabase.h index 328d4cacb..b65178d68 100644 --- a/lib/utils/HDFDatabase.h +++ b/lib/utils/HDFDatabase.h @@ -46,7 +46,7 @@ class HDFDatabase : public Database virtual bool create( - const std::string& file_name); + const std::string& file_name) override; /** * @brief Opens an existing HDF5 database file with the supplied name. @@ -60,7 +60,7 @@ class HDFDatabase : public Database bool open( const std::string& file_name, - const std::string& type); + const std::string& type) override; /** * @brief Closes the currently open HDF5 database file. From 25d40a3e7772402da3046cbf242e78d46d6d1716 Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Tue, 16 Jan 2024 18:23:41 -0800 Subject: [PATCH 02/14] stylization. --- lib/utils/Database.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/lib/utils/Database.cpp b/lib/utils/Database.cpp index 4533afd38..1ee734a18 100644 --- a/lib/utils/Database.cpp +++ b/lib/utils/Database.cpp @@ -18,9 +18,9 @@ namespace CAROM { bool fileExists(const std::string& name) { - std::ifstream f(name.c_str()); - return f.good(); - // ifstream f will be closed upon the end of the function. + std::ifstream f(name.c_str()); + return f.good(); + // ifstream f will be closed upon the end of the function. } Database::Database() From 68d9ed11123ca937ce087b31c28cabc8d34ca771 Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Tue, 16 Jan 2024 18:48:53 -0800 Subject: [PATCH 03/14] HDFDatabase::putIntegerArray - overwrites if the dataset exists. --- lib/linalg/BasisWriter.cpp | 3 ++- lib/utils/HDFDatabase.cpp | 35 +++++++++++++++++++++++------------ 2 files changed, 25 insertions(+), 13 deletions(-) diff --git a/lib/linalg/BasisWriter.cpp b/lib/linalg/BasisWriter.cpp index e3e06ff6d..4ff3dd700 100644 --- a/lib/linalg/BasisWriter.cpp +++ b/lib/linalg/BasisWriter.cpp @@ -80,7 +80,8 @@ BasisWriter::writeBasis(const std::string& kind) sprintf(tmp, "time_%06d", d_num_intervals_written); if (kind == "basis") { - if (fileExists(full_file_name)) + bool file_exists = fileExists(full_file_name); + if (file_exists) d_database->open(full_file_name, "wr"); else d_database->create(full_file_name); diff --git a/lib/utils/HDFDatabase.cpp b/lib/utils/HDFDatabase.cpp index 4c00f7e40..e99af4ef6 100644 --- a/lib/utils/HDFDatabase.cpp +++ b/lib/utils/HDFDatabase.cpp @@ -12,6 +12,7 @@ #include "HDFDatabase.h" #include "Utilities.h" +#include namespace CAROM { @@ -123,21 +124,31 @@ HDFDatabase::putIntegerArray( hid_t space = H5Screate_simple(1, dim, 0); CAROM_VERIFY(space >= 0); + const bool key_exists = (H5Lexists(d_group_id, key.c_str(), H5P_DEFAULT) > 0); + hid_t dataset; + if (key_exists) + { + std::cout << "HDF5Database: overwriting dataset " << key << std::endl; + dataset = H5Dopen(d_group_id, key.c_str(), H5P_DEFAULT); + } + else + { #if (H5_VERS_MAJOR > 1) || ((H5_VERS_MAJOR == 1) && (H5_VERS_MINOR > 6)) - hid_t dataset = H5Dcreate(d_group_id, - key.c_str(), - H5T_STD_I32BE, - space, - H5P_DEFAULT, - H5P_DEFAULT, - H5P_DEFAULT); + dataset = H5Dcreate(d_group_id, + key.c_str(), + H5T_STD_I32BE, + space, + H5P_DEFAULT, + H5P_DEFAULT, + H5P_DEFAULT); #else - hid_t dataset = H5Dcreate(d_group_id, - key.c_str(), - H5T_STD_I32BE, - space, - H5P_DEFAULT); + dataset = H5Dcreate(d_group_id, + key.c_str(), + H5T_STD_I32BE, + space, + H5P_DEFAULT); #endif + } CAROM_VERIFY(dataset >= 0); herr_t errf = H5Dwrite(dataset, From d0e89ee1f6d8201ee30099532c40d488fdb3f733 Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Thu, 8 Feb 2024 16:36:27 -0800 Subject: [PATCH 04/14] enforce single time interval in Options. --- lib/linalg/Options.h | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/lib/linalg/Options.h b/lib/linalg/Options.h index 082d4a2b2..4d1751dbe 100644 --- a/lib/linalg/Options.h +++ b/lib/linalg/Options.h @@ -50,7 +50,15 @@ class Options samples_per_time_interval(samples_per_time_interval_), max_time_intervals(max_time_intervals_), update_right_SV(update_right_SV_), - write_snapshots(write_snapshots_) {}; + write_snapshots(write_snapshots_) + { + if (max_time_intervals > 1) + { + printf("time interval is obsolete and will be removed in the future. Set max_time_intervals=%d to 1 or -1!\n", + max_time_intervals); + } + CAROM_VERIFY(max_time_intervals <= 1); + }; /** * @brief Sets the maximum basis dimension of the SVD algorithm. From 19d78e36ab03670e5dc3879399eb158ad78a4cc9 Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Thu, 8 Feb 2024 16:47:43 -0800 Subject: [PATCH 05/14] HDFDatabase::putIntegerArray does not allow overwrite. --- lib/utils/HDFDatabase.cpp | 34 ++++++++++++---------------------- 1 file changed, 12 insertions(+), 22 deletions(-) diff --git a/lib/utils/HDFDatabase.cpp b/lib/utils/HDFDatabase.cpp index e99af4ef6..aeaba1970 100644 --- a/lib/utils/HDFDatabase.cpp +++ b/lib/utils/HDFDatabase.cpp @@ -124,31 +124,21 @@ HDFDatabase::putIntegerArray( hid_t space = H5Screate_simple(1, dim, 0); CAROM_VERIFY(space >= 0); - const bool key_exists = (H5Lexists(d_group_id, key.c_str(), H5P_DEFAULT) > 0); - hid_t dataset; - if (key_exists) - { - std::cout << "HDF5Database: overwriting dataset " << key << std::endl; - dataset = H5Dopen(d_group_id, key.c_str(), H5P_DEFAULT); - } - else - { #if (H5_VERS_MAJOR > 1) || ((H5_VERS_MAJOR == 1) && (H5_VERS_MINOR > 6)) - dataset = H5Dcreate(d_group_id, - key.c_str(), - H5T_STD_I32BE, - space, - H5P_DEFAULT, - H5P_DEFAULT, - H5P_DEFAULT); + hid_t dataset = H5Dcreate(d_group_id, + key.c_str(), + H5T_STD_I32BE, + space, + H5P_DEFAULT, + H5P_DEFAULT, + H5P_DEFAULT); #else - dataset = H5Dcreate(d_group_id, - key.c_str(), - H5T_STD_I32BE, - space, - H5P_DEFAULT); + hid_t dataset = H5Dcreate(d_group_id, + key.c_str(), + H5T_STD_I32BE, + space, + H5P_DEFAULT); #endif - } CAROM_VERIFY(dataset >= 0); herr_t errf = H5Dwrite(dataset, From 4ecec6da5a0c96a2f31fe947737d9c00aa0ddc71 Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Thu, 8 Feb 2024 16:53:11 -0800 Subject: [PATCH 06/14] BasisWriter::writeBasis always create the file, which will overwrite the exisiting file. --- lib/linalg/BasisWriter.cpp | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/lib/linalg/BasisWriter.cpp b/lib/linalg/BasisWriter.cpp index 4ff3dd700..37ffeceb7 100644 --- a/lib/linalg/BasisWriter.cpp +++ b/lib/linalg/BasisWriter.cpp @@ -81,10 +81,7 @@ BasisWriter::writeBasis(const std::string& kind) if (kind == "basis") { bool file_exists = fileExists(full_file_name); - if (file_exists) - d_database->open(full_file_name, "wr"); - else - d_database->create(full_file_name); + d_database->create(full_file_name); d_database->putDouble(tmp, time_interval_start_time); @@ -124,10 +121,7 @@ BasisWriter::writeBasis(const std::string& kind) } if (kind == "snapshot") { - if (fileExists(snap_file_name)) - d_snap_database->open(snap_file_name, "wr"); - else - d_snap_database->create(snap_file_name); + d_snap_database->create(snap_file_name); d_snap_database->putDouble(tmp, time_interval_start_time); From bbaa999318f440583629cfb8161aa12fa3197949 Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Thu, 8 Feb 2024 17:03:38 -0800 Subject: [PATCH 07/14] add a header and stylization. --- lib/linalg/Options.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/lib/linalg/Options.h b/lib/linalg/Options.h index 4d1751dbe..803d1ce0b 100644 --- a/lib/linalg/Options.h +++ b/lib/linalg/Options.h @@ -14,6 +14,8 @@ #ifndef included_Options_h #define included_Options_h +#include "utils/Utilities.h" + namespace CAROM { /** @@ -55,7 +57,7 @@ class Options if (max_time_intervals > 1) { printf("time interval is obsolete and will be removed in the future. Set max_time_intervals=%d to 1 or -1!\n", - max_time_intervals); + max_time_intervals); } CAROM_VERIFY(max_time_intervals <= 1); }; From 31ad436b8ce4c4c274f251835f58ae8c849c5223 Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Thu, 8 Feb 2024 17:13:33 -0800 Subject: [PATCH 08/14] remove increase time interval test, as time interval is fixed to 1. --- unit_tests/test_SVD.cpp | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/unit_tests/test_SVD.cpp b/unit_tests/test_SVD.cpp index 52288d144..ad165bf23 100644 --- a/unit_tests/test_SVD.cpp +++ b/unit_tests/test_SVD.cpp @@ -212,20 +212,6 @@ TEST(SVDSerialTest, Test_getBasisIntervalStartTime) EXPECT_DOUBLE_EQ(svd.getBasisIntervalStartTime(2), 2); } - -TEST(SVDSerialTest, Test_increaseTimeInterval) -{ - FakeSVD svd(CAROM::Options(5, 2, 2)); - - ASSERT_NO_THROW(svd.takeSample(NULL, 0, true)); - ASSERT_NO_THROW(svd.takeSample(NULL, 0.5, true)); - ASSERT_NO_THROW(svd.takeSample(NULL, 1, true)); - ASSERT_NO_THROW(svd.takeSample(NULL, 1.5, true)); - - /* The maximum number of time intervals is surpassed */ - EXPECT_DEATH(svd.takeSample(NULL, 2, true), ".*"); -} - int main(int argc, char* argv[]) { ::testing::InitGoogleTest(&argc, argv); From 9c0f8b615587c233752be6246be7e76b13ea408f Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Sun, 11 Feb 2024 15:43:16 -0800 Subject: [PATCH 09/14] StaticSVD::ComputeSVD - copies the sample snapshot matrix for SVD operation and deletes it after SVD. --- lib/linalg/svd/StaticSVD.cpp | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/lib/linalg/svd/StaticSVD.cpp b/lib/linalg/svd/StaticSVD.cpp index 9684b6032..c8f0ad493 100644 --- a/lib/linalg/svd/StaticSVD.cpp +++ b/lib/linalg/svd/StaticSVD.cpp @@ -262,10 +262,9 @@ StaticSVD::computeSVD() // transpose matrix if sample size > dimension. bool transpose = d_total_dim < d_num_samples; + snapshot_matrix = new SLPK_Matrix; if (transpose) { // create a transposed matrix if sample size > dimension. - snapshot_matrix = new SLPK_Matrix; - int d_blocksize_tr = d_num_samples / d_nprow; if (d_num_samples % d_nprow != 0) { d_blocksize_tr += 1; @@ -282,8 +281,12 @@ StaticSVD::computeSVD() } } else { - // use d_samples if sample size <= dimension. - snapshot_matrix = d_samples.get(); + // copy d_samples if sample size <= dimension. + // SVD operation does not preserve the original snapshot matrix. + initialize_matrix(snapshot_matrix, d_total_dim, d_num_samples, + d_nprow, d_npcol, d_blocksize, d_blocksize); + copy_matrix(snapshot_matrix, 1, 1, d_samples.get(), 1, 1, d_total_dim, + d_num_samples); } // This block does the actual ScaLAPACK call to do the factorization. @@ -372,11 +375,9 @@ StaticSVD::computeSVD() } } - if (transpose) { - // Delete the transposed snapshot matrix. - free_matrix_data(snapshot_matrix); - delete snapshot_matrix; - } + // Delete the snapshot matrix. + free_matrix_data(snapshot_matrix); + delete snapshot_matrix; } void From dd9d84e610a0724d84068af5878d274a5d2a897f Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Tue, 13 Feb 2024 14:43:19 -0800 Subject: [PATCH 10/14] add communicator for Vector. --- lib/linalg/Vector.cpp | 39 ++++++++++++++++++++++----------------- lib/linalg/Vector.h | 16 ++++++++++++++-- 2 files changed, 36 insertions(+), 19 deletions(-) diff --git a/lib/linalg/Vector.cpp b/lib/linalg/Vector.cpp index dccf72995..b62da2bdd 100644 --- a/lib/linalg/Vector.cpp +++ b/lib/linalg/Vector.cpp @@ -15,8 +15,6 @@ #include "Vector.h" #include "utils/HDFDatabase.h" -#include "mpi.h" - #include #include @@ -26,23 +24,26 @@ Vector::Vector() : d_vec(NULL), d_alloc_size(0), d_distributed(false), - d_owns_data(true) + d_owns_data(true), + d_comm(MPI_COMM_WORLD) {} Vector::Vector( int dim, - bool distributed) : + bool distributed, + MPI_Comm comm) : d_vec(NULL), d_alloc_size(0), d_distributed(distributed), - d_owns_data(true) + d_owns_data(true), + d_comm(comm) { CAROM_VERIFY(dim > 0); setSize(dim); int mpi_init; MPI_Initialized(&mpi_init); if (mpi_init) { - MPI_Comm_size(MPI_COMM_WORLD, &d_num_procs); + MPI_Comm_size(d_comm, &d_num_procs); } else { d_num_procs = 1; @@ -53,11 +54,13 @@ Vector::Vector( double* vec, int dim, bool distributed, - bool copy_data) : + bool copy_data, + MPI_Comm comm) : d_vec(NULL), d_alloc_size(0), d_distributed(distributed), - d_owns_data(copy_data) + d_owns_data(copy_data), + d_comm(comm) { CAROM_VERIFY(vec != 0); CAROM_VERIFY(dim > 0); @@ -73,7 +76,7 @@ Vector::Vector( int mpi_init; MPI_Initialized(&mpi_init); if (mpi_init) { - MPI_Comm_size(MPI_COMM_WORLD, &d_num_procs); + MPI_Comm_size(d_comm, &d_num_procs); } else { d_num_procs = 1; @@ -85,13 +88,14 @@ Vector::Vector( d_vec(NULL), d_alloc_size(0), d_distributed(other.d_distributed), - d_owns_data(true) + d_owns_data(true), + d_comm(other.getComm()) { setSize(other.d_dim); int mpi_init; MPI_Initialized(&mpi_init); if (mpi_init) { - MPI_Comm_size(MPI_COMM_WORLD, &d_num_procs); + MPI_Comm_size(d_comm, &d_num_procs); } else { d_num_procs = 1; @@ -224,13 +228,14 @@ Vector::inner_product( { CAROM_VERIFY(dim() == other.dim()); CAROM_ASSERT(distributed() == other.distributed()); + CAROM_ASSERT(d_comm == other.getComm()); double ip; double local_ip = 0.0; for (int i = 0; i < d_dim; ++i) { local_ip += d_vec[i]*other.d_vec[i]; } if (d_num_procs > 1 && d_distributed) { - MPI_Allreduce(&local_ip, &ip, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&local_ip, &ip, 1, MPI_DOUBLE, MPI_SUM, d_comm); } else { ip = local_ip; @@ -451,7 +456,7 @@ Vector::write(const std::string& base_file_name) MPI_Initialized(&mpi_init); int rank; if (mpi_init) { - MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_rank(d_comm, &rank); } else { rank = 0; @@ -476,7 +481,7 @@ void Vector::print(const char * prefix) const { int my_rank; - const bool success = MPI_Comm_rank(MPI_COMM_WORLD, &my_rank); + const bool success = MPI_Comm_rank(d_comm, &my_rank); CAROM_ASSERT(success); std::string filename_str = prefix + std::to_string(my_rank); @@ -497,7 +502,7 @@ Vector::read(const std::string& base_file_name) MPI_Initialized(&mpi_init); int rank; if (mpi_init) { - MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_rank(d_comm, &rank); } else { rank = 0; @@ -521,7 +526,7 @@ Vector::read(const std::string& base_file_name) database.getDoubleArray(tmp, d_vec, d_alloc_size); d_owns_data = true; if (mpi_init) { - MPI_Comm_size(MPI_COMM_WORLD, &d_num_procs); + MPI_Comm_size(d_comm, &d_num_procs); } else { d_num_procs = 1; @@ -555,7 +560,7 @@ Vector::local_read(const std::string& base_file_name, int rank) database.getDoubleArray(tmp, d_vec, d_alloc_size); d_owns_data = true; if (mpi_init) { - MPI_Comm_size(MPI_COMM_WORLD, &d_num_procs); + MPI_Comm_size(d_comm, &d_num_procs); } else { d_num_procs = 1; diff --git a/lib/linalg/Vector.h b/lib/linalg/Vector.h index 0daea6237..d7c070544 100644 --- a/lib/linalg/Vector.h +++ b/lib/linalg/Vector.h @@ -18,6 +18,7 @@ #include "utils/Utilities.h" #include #include +#include "mpi.h" namespace CAROM { @@ -41,10 +42,12 @@ class Vector * the Vector on this processor. * @param[in] distributed If true the dimensions of the Vector are spread * over all processors. + * @param[in] comm MPI communicator of the processors that own this Vector. */ Vector( int dim, - bool distributed); + bool distributed, + MPI_Comm comm = MPI_COMM_WORLD); /** * @brief Constructor in which the Vector is given its initial values. @@ -61,12 +64,14 @@ class Vector * @param[in] copy_data If true the vector allocates its own storage and * copies the contents of vec into its own storage. * Otherwise it uses vec as its storage. + * @param[in] comm MPI communicator of the processors that own this Vector. */ Vector( double* vec, int dim, bool distributed, - bool copy_data = true); + bool copy_data = true, + MPI_Comm comm = MPI_COMM_WORLD); /** * @brief Copy constructor. @@ -759,6 +764,8 @@ class Vector */ double localMin(int nmax = 0); + const MPI_Comm getComm() const { return d_comm; } + private: /** * @brief The storage for the Vector's values on this processor. @@ -797,6 +804,11 @@ class Vector * If d_owns_data is false, then the object may not reallocate d_vec. */ bool d_owns_data; + + /** + * @brief MPI communicator of the processors that own this Vector. + */ + MPI_Comm d_comm; }; /** From 962b76c55115119629f1e379b04cff0ec7a7d93e Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Tue, 13 Feb 2024 15:29:05 -0800 Subject: [PATCH 11/14] stylization. --- lib/linalg/Vector.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/lib/linalg/Vector.h b/lib/linalg/Vector.h index d7c070544..965fe671b 100644 --- a/lib/linalg/Vector.h +++ b/lib/linalg/Vector.h @@ -764,7 +764,9 @@ class Vector */ double localMin(int nmax = 0); - const MPI_Comm getComm() const { return d_comm; } + const MPI_Comm getComm() const { + return d_comm; + } private: /** From 1ceb57b3e4f3e88b0b687aa4bb4e71e6b9353127 Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Tue, 13 Feb 2024 16:28:28 -0800 Subject: [PATCH 12/14] enforce all operations are executed within the same MPI communicator. --- lib/linalg/Vector.cpp | 35 ++++++++++++++++++++++++++++++----- lib/linalg/Vector.h | 10 +++++++++- 2 files changed, 39 insertions(+), 6 deletions(-) diff --git a/lib/linalg/Vector.cpp b/lib/linalg/Vector.cpp index b62da2bdd..653f4ab6d 100644 --- a/lib/linalg/Vector.cpp +++ b/lib/linalg/Vector.cpp @@ -20,12 +20,12 @@ namespace CAROM { -Vector::Vector() : +Vector::Vector(MPI_Comm comm) : d_vec(NULL), d_alloc_size(0), d_distributed(false), d_owns_data(true), - d_comm(MPI_COMM_WORLD) + d_comm(comm) {} Vector::Vector( @@ -116,6 +116,7 @@ Vector::operator = ( { d_distributed = rhs.d_distributed; d_num_procs = rhs.d_num_procs; + d_comm = rhs.getComm(); setSize(rhs.d_dim); memcpy(d_vec, rhs.d_vec, d_dim*sizeof(double)); return *this; @@ -126,6 +127,7 @@ Vector::operator += ( const Vector& rhs) { CAROM_VERIFY(d_dim == rhs.d_dim); + CAROM_VERIFY(d_comm == rhs.getComm()); for(int i=0; i void Vector::transform(Vector& result, std::function transformer) const { + CAROM_VERIFY(d_comm == result.getComm()); result.setSize(d_dim); result.d_distributed = d_distributed; transformer(d_dim, result.d_vec); @@ -171,6 +175,8 @@ Vector::transform(Vector& result, void Vector::transform(Vector*& result, std::function transformer) const { + CAROM_ASSERT(result == 0 || result->getComm() == d_comm); + // If the result has not been allocated then do so. Otherwise size it // correctly. if (result == 0) { @@ -197,6 +203,8 @@ void Vector::transform(Vector& result, std::function transformer) const { + CAROM_VERIFY(d_comm == result.getComm()); + Vector* origVector = new Vector(*this); result.setSize(d_dim); result.d_distributed = d_distributed; @@ -208,6 +216,8 @@ void Vector::transform(Vector*& result, std::function transformer) const { + CAROM_ASSERT(result == 0 || result->getComm() == d_comm); + // If the result has not been allocated then do so. Otherwise size it // correctly. Vector* origVector = new Vector(*this); @@ -273,7 +283,9 @@ Vector::plus( Vector*& result) const { CAROM_ASSERT(result == 0 || result->distributed() == distributed()); + CAROM_ASSERT(result == 0 || result->getComm() == d_comm); CAROM_ASSERT(distributed() == other.distributed()); + CAROM_ASSERT(d_comm == other.getComm()); CAROM_VERIFY(dim() == other.dim()); // If the result has not been allocated then do so. Otherwise size it @@ -297,7 +309,9 @@ Vector::plus( Vector& result) const { CAROM_ASSERT(result.distributed() == distributed()); + CAROM_ASSERT(result.getComm() == d_comm); CAROM_ASSERT(distributed() == other.distributed()); + CAROM_ASSERT(d_comm == other.getComm()); CAROM_VERIFY(dim() == other.dim()); // Size result correctly. @@ -316,7 +330,9 @@ Vector::plusAx( Vector*& result) const { CAROM_ASSERT(result == 0 || result->distributed() == distributed()); + CAROM_ASSERT(result == 0 || result->getComm() == d_comm); CAROM_ASSERT(distributed() == other.distributed()); + CAROM_ASSERT(d_comm == other.getComm()); CAROM_VERIFY(dim() == other.dim()); // If the result has not been allocated then do so. Otherwise size it @@ -341,7 +357,9 @@ Vector::plusAx( Vector& result) const { CAROM_ASSERT(result.distributed() == distributed()); + CAROM_ASSERT(result.getComm() == d_comm); CAROM_ASSERT(distributed() == other.distributed()); + CAROM_ASSERT(d_comm == other.getComm()); CAROM_VERIFY(dim() == other.dim()); // Size result correctly. @@ -359,6 +377,7 @@ Vector::plusEqAx( const Vector& other) { CAROM_ASSERT(distributed() == other.distributed()); + CAROM_ASSERT(d_comm == other.getComm()); CAROM_VERIFY(dim() == other.dim()); // Do the addition. @@ -373,7 +392,9 @@ Vector::minus( Vector*& result) const { CAROM_ASSERT(result == 0 || result->distributed() == distributed()); + CAROM_ASSERT(result == 0 || result->getComm() == d_comm); CAROM_ASSERT(distributed() == other.distributed()); + CAROM_ASSERT(d_comm == other.getComm()); CAROM_VERIFY(dim() == other.dim()); // If the result has not been allocated then do so. Otherwise size it @@ -397,7 +418,9 @@ Vector::minus( Vector& result) const { CAROM_ASSERT(result.distributed() == distributed()); + CAROM_ASSERT(result.getComm() == d_comm); CAROM_ASSERT(distributed() == other.distributed()); + CAROM_ASSERT(d_comm == other.getComm()); CAROM_VERIFY(dim() == other.dim()); // Size result correctly. @@ -415,6 +438,7 @@ Vector::mult( Vector*& result) const { CAROM_ASSERT(result == 0 || result->distributed() == distributed()); + CAROM_ASSERT(result == 0 || result->getComm() == d_comm); // If the result has not been allocated then do so. Otherwise size it // correctly. @@ -437,6 +461,7 @@ Vector::mult( Vector& result) const { CAROM_ASSERT(result.distributed() == distributed()); + CAROM_ASSERT(result.getComm() == d_comm); // Size result correctly. result.setSize(d_dim); @@ -599,7 +624,7 @@ int getCenterPoint(std::vector& points, double min_dist_to_centroid; for (int i = 0; i < points.size(); i++) { - Vector diff; + Vector diff(points[0]->getComm()); centroid.minus(*points[i], diff); double dist_to_centroid = diff.norm(); if (i == 0 || dist_to_centroid < min_dist_to_centroid) @@ -617,7 +642,7 @@ int getCenterPoint(std::vector& points, std::vector dist_to_points(points.size(), 0); for (int i = 0; i < points.size(); i++) { for (int j = i + 1; j < points.size(); j++) { - Vector diff; + Vector diff(points[0]->getComm()); points[i]->minus(*points[j], diff); double dist = diff.norm(); dist_to_points[i] += dist; @@ -655,7 +680,7 @@ int getClosestPoint(std::vector& points, int closest_point = 0; double closest_dist_to_test_point = INT_MAX; for (int i = 0; i < points.size(); i++) { - Vector diff; + Vector diff(points[0]->getComm()); test_point->minus(*points[i], diff); double dist = diff.norm(); if (dist < closest_dist_to_test_point) diff --git a/lib/linalg/Vector.h b/lib/linalg/Vector.h index 965fe671b..1aed9b837 100644 --- a/lib/linalg/Vector.h +++ b/lib/linalg/Vector.h @@ -30,7 +30,12 @@ namespace CAROM { class Vector { public: - Vector(); + /** + * @brief Constructor creating a Vector with uninitialized values and size. + * + * @param[in] comm MPI communicator of the processors that own this Vector. + */ + Vector(MPI_Comm comm = MPI_COMM_WORLD); /** * @brief Constructor creating a Vector with uninitialized values. @@ -764,6 +769,9 @@ class Vector */ double localMin(int nmax = 0); + /** + * @brief Get MPI communicator of this Vector. + */ const MPI_Comm getComm() const { return d_comm; } From bdaf073424163152446cd46c0a4115ef34ec990b Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Tue, 13 Feb 2024 16:36:56 -0800 Subject: [PATCH 13/14] communicator cannot be changed once initialized. So MPI scope of Vector never changes. --- lib/linalg/Vector.cpp | 2 +- lib/linalg/Vector.h | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/lib/linalg/Vector.cpp b/lib/linalg/Vector.cpp index 653f4ab6d..040ef2d89 100644 --- a/lib/linalg/Vector.cpp +++ b/lib/linalg/Vector.cpp @@ -114,9 +114,9 @@ Vector& Vector::operator = ( const Vector& rhs) { + CAROM_VERIFY(d_comm == rhs.getComm()); d_distributed = rhs.d_distributed; d_num_procs = rhs.d_num_procs; - d_comm = rhs.getComm(); setSize(rhs.d_dim); memcpy(d_vec, rhs.d_vec, d_dim*sizeof(double)); return *this; diff --git a/lib/linalg/Vector.h b/lib/linalg/Vector.h index 1aed9b837..01494f0f8 100644 --- a/lib/linalg/Vector.h +++ b/lib/linalg/Vector.h @@ -32,7 +32,7 @@ class Vector public: /** * @brief Constructor creating a Vector with uninitialized values and size. - * + * * @param[in] comm MPI communicator of the processors that own this Vector. */ Vector(MPI_Comm comm = MPI_COMM_WORLD); @@ -817,8 +817,9 @@ class Vector /** * @brief MPI communicator of the processors that own this Vector. + * Cannot change once initialized. */ - MPI_Comm d_comm; + const MPI_Comm d_comm; }; /** From 2007c300628e036fc512d458a475c0d258598bfa Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Tue, 13 Feb 2024 17:29:11 -0800 Subject: [PATCH 14/14] mpi communicator for matrix. if scalapack routines are used, communicator is enforced to be MPI_COMM_WORLD. --- lib/linalg/Matrix.cpp | 200 +++++++++++++++++++++++++++++------------- lib/linalg/Matrix.h | 19 +++- 2 files changed, 155 insertions(+), 64 deletions(-) diff --git a/lib/linalg/Matrix.cpp b/lib/linalg/Matrix.cpp index 4918b1ba1..956e50b10 100644 --- a/lib/linalg/Matrix.cpp +++ b/lib/linalg/Matrix.cpp @@ -16,7 +16,6 @@ #include "utils/HDFDatabase.h" #include "utils/mpi_utils.h" -#include "mpi.h" #include #include #include @@ -64,26 +63,29 @@ Matrix::Matrix() : d_mat(NULL), d_alloc_size(0), d_distributed(false), - d_owns_data(true) + d_owns_data(true), + d_comm(MPI_COMM_WORLD) {} Matrix::Matrix( int num_rows, int num_cols, bool distributed, - bool randomized) : + bool randomized, + MPI_Comm comm) : d_mat(0), d_alloc_size(0), d_distributed(distributed), - d_owns_data(true) + d_owns_data(true), + d_comm(comm) { CAROM_VERIFY(num_rows > 0); CAROM_VERIFY(num_cols > 0); int mpi_init; MPI_Initialized(&mpi_init); if (mpi_init) { - MPI_Comm_size(MPI_COMM_WORLD, &d_num_procs); - MPI_Comm_rank(MPI_COMM_WORLD, &d_rank); + MPI_Comm_size(d_comm, &d_num_procs); + MPI_Comm_rank(d_comm, &d_rank); } else { d_num_procs = 1; @@ -106,11 +108,13 @@ Matrix::Matrix( int num_rows, int num_cols, bool distributed, - bool copy_data) : + bool copy_data, + MPI_Comm comm) : d_mat(0), d_alloc_size(0), d_distributed(distributed), - d_owns_data(copy_data) + d_owns_data(copy_data), + d_comm(comm) { CAROM_VERIFY(mat != 0); CAROM_VERIFY(num_rows > 0); @@ -118,8 +122,8 @@ Matrix::Matrix( int mpi_init; MPI_Initialized(&mpi_init); if (mpi_init) { - MPI_Comm_size(MPI_COMM_WORLD, &d_num_procs); - MPI_Comm_rank(MPI_COMM_WORLD, &d_rank); + MPI_Comm_size(d_comm, &d_num_procs); + MPI_Comm_rank(d_comm, &d_rank); } else { d_num_procs = 1; @@ -145,13 +149,14 @@ Matrix::Matrix( d_mat(0), d_alloc_size(0), d_distributed(other.d_distributed), - d_owns_data(true) + d_owns_data(true), + d_comm(other.getComm()) { int mpi_init; MPI_Initialized(&mpi_init); if (mpi_init) { - MPI_Comm_size(MPI_COMM_WORLD, &d_num_procs); - MPI_Comm_rank(MPI_COMM_WORLD, &d_rank); + MPI_Comm_size(d_comm, &d_num_procs); + MPI_Comm_rank(d_comm, &d_rank); } else { d_num_procs = 1; @@ -172,6 +177,7 @@ Matrix& Matrix::operator = ( const Matrix& rhs) { + CAROM_VERIFY(d_comm == rhs.getComm()); d_distributed = rhs.d_distributed; d_num_procs = rhs.d_num_procs; setSize(rhs.d_num_rows, rhs.d_num_cols); @@ -183,6 +189,7 @@ Matrix& Matrix::operator += ( const Matrix& rhs) { + CAROM_VERIFY(d_comm == rhs.getComm()); CAROM_VERIFY(rhs.d_num_rows == d_num_rows); CAROM_VERIFY(rhs.d_num_cols == d_num_cols); for(int i=0; idistributed() == distributed()); + CAROM_VERIFY(result == 0 || result->getComm() == d_comm); CAROM_VERIFY(n > 0 && n <= d_num_cols); // If the result has not been allocated then do so. Otherwise size it // correctly. if (result == 0) { - result = new Matrix(d_num_rows, n, d_distributed); + result = new Matrix(d_num_rows, n, d_distributed, d_comm); } else { @@ -296,6 +305,7 @@ Matrix::getFirstNColumns( Matrix& result) const { CAROM_VERIFY(result.distributed() == distributed()); + CAROM_VERIFY(result.getComm() == d_comm); CAROM_VERIFY(n > 0 && n <= d_num_cols); // Size result correctly. @@ -316,13 +326,15 @@ Matrix::mult( Matrix*& result) const { CAROM_VERIFY(result == 0 || result->distributed() == distributed()); + CAROM_VERIFY(result == 0 || result->getComm() == d_comm); CAROM_VERIFY(!other.distributed()); + CAROM_VERIFY(d_comm == other.getComm()); CAROM_VERIFY(numColumns() == other.numRows()); // If the result has not been allocated then do so. Otherwise size it // correctly. if (result == 0) { - result = new Matrix(d_num_rows, other.d_num_cols, d_distributed); + result = new Matrix(d_num_rows, other.d_num_cols, d_distributed, d_comm); } else { result->setSize(d_num_rows, other.d_num_cols); @@ -346,7 +358,9 @@ Matrix::mult( Matrix& result) const { CAROM_VERIFY(result.distributed() == distributed()); + CAROM_VERIFY(result.getComm() == d_comm); CAROM_VERIFY(!other.distributed()); + CAROM_VERIFY(d_comm == other.getComm()); CAROM_VERIFY(numColumns() == other.numRows()); // Size result correctly. @@ -370,7 +384,9 @@ Matrix::mult( Vector*& result) const { CAROM_VERIFY(result == 0 || result->distributed() == distributed()); + CAROM_VERIFY(result == 0 || result->getComm() == d_comm); CAROM_VERIFY(!other.distributed()); + CAROM_VERIFY(d_comm == other.getComm()); CAROM_VERIFY(numColumns() == other.dim()); // If the result has not been allocated then do so. Otherwise size it @@ -398,7 +414,9 @@ Matrix::mult( Vector& result) const { CAROM_VERIFY(result.distributed() == distributed()); + CAROM_VERIFY(result.getComm() == d_comm); CAROM_VERIFY(!other.distributed()); + CAROM_VERIFY(d_comm == other.getComm()); CAROM_VERIFY(numColumns() == other.dim()); // Size result correctly. @@ -421,6 +439,8 @@ Matrix::pointwise_mult( Vector& result) const { // TODO: change the CAROM_ASSERTs to CAROM_VERIFYs or generalize and eliminate the checks. + CAROM_VERIFY(d_comm == other.getComm()); + CAROM_VERIFY(d_comm == result.getComm()); CAROM_ASSERT(!result.distributed()); CAROM_ASSERT(!distributed()); CAROM_VERIFY(!other.distributed()); @@ -437,6 +457,7 @@ Matrix::pointwise_mult( int this_row, Vector& other) const { + CAROM_VERIFY(d_comm == other.getComm()); CAROM_ASSERT(!distributed()); CAROM_VERIFY(!other.distributed()); CAROM_VERIFY(numColumns() == other.dim()); @@ -453,14 +474,16 @@ Matrix::elementwise_mult( Matrix*& result) const { CAROM_VERIFY(result == 0 || result->distributed() == distributed()); + CAROM_VERIFY(result == 0 || result->getComm() == d_comm); CAROM_VERIFY(distributed() == other.distributed()); + CAROM_VERIFY(other.getComm() == d_comm); CAROM_VERIFY(numRows() == other.numRows()); CAROM_VERIFY(numColumns() == other.numColumns()); // If the result has not been allocated then do so. Otherwise size it // correctly. if (result == 0) { - result = new Matrix(d_num_rows, d_num_cols, d_distributed); + result = new Matrix(d_num_rows, d_num_cols, d_distributed, d_comm); } else { result->setSize(d_num_rows, d_num_cols); @@ -481,7 +504,9 @@ Matrix::elementwise_mult( Matrix& result) const { CAROM_VERIFY(result.distributed() == distributed()); + CAROM_VERIFY(result.getComm() == d_comm); CAROM_VERIFY(distributed() == other.distributed()); + CAROM_VERIFY(other.getComm() == d_comm); CAROM_VERIFY(numRows() == other.numRows()); CAROM_VERIFY(numColumns() == other.numColumns()); @@ -501,10 +526,13 @@ void Matrix::elementwise_square( Matrix*& result) const { + CAROM_VERIFY(result == 0 || result->distributed() == distributed()); + CAROM_VERIFY(result == 0 || result->getComm() == d_comm); + // If the result has not been allocated then do so. Otherwise size it // correctly. if (result == 0) { - result = new Matrix(d_num_rows, d_num_cols, d_distributed); + result = new Matrix(d_num_rows, d_num_cols, d_distributed, d_comm); } else { result->setSize(d_num_rows, d_num_cols); @@ -523,6 +551,9 @@ void Matrix::elementwise_square( Matrix& result) const { + CAROM_VERIFY(result.distributed() == distributed()); + CAROM_VERIFY(result.getComm() == d_comm); + // Size result correctly. result.setSize(d_num_rows, d_num_cols); @@ -542,6 +573,7 @@ Matrix::multPlus( double c) const { CAROM_VERIFY(a.distributed() == distributed()); + CAROM_VERIFY(a.getComm() == d_comm); CAROM_VERIFY(!b.distributed()); CAROM_VERIFY(numColumns() == b.dim()); CAROM_VERIFY(numRows() == a.dim()); @@ -561,13 +593,15 @@ Matrix::transposeMult( Matrix*& result) const { CAROM_VERIFY(result == 0 || !result->distributed()); + CAROM_VERIFY(result == 0 || result->getComm() == d_comm); CAROM_VERIFY(distributed() == other.distributed()); + CAROM_VERIFY(d_comm == other.getComm()); CAROM_VERIFY(numRows() == other.numRows()); // If the result has not been allocated then do so. Otherwise size it // correctly. if (result == 0) { - result = new Matrix(d_num_cols, other.d_num_cols, false); + result = new Matrix(d_num_cols, other.d_num_cols, false, d_comm); } else { result->setSize(d_num_cols, other.d_num_cols); @@ -590,7 +624,7 @@ Matrix::transposeMult( new_mat_size, MPI_DOUBLE, MPI_SUM, - MPI_COMM_WORLD); + d_comm); } } @@ -600,7 +634,9 @@ Matrix::transposeMult( Matrix& result) const { CAROM_VERIFY(!result.distributed()); + CAROM_VERIFY(d_comm == result.getComm()); CAROM_VERIFY(distributed() == other.distributed()); + CAROM_VERIFY(d_comm == other.getComm()); CAROM_VERIFY(numRows() == other.numRows()); // Size result correctly. @@ -623,7 +659,7 @@ Matrix::transposeMult( new_mat_size, MPI_DOUBLE, MPI_SUM, - MPI_COMM_WORLD); + d_comm); } } @@ -633,13 +669,15 @@ Matrix::transposeMult( Vector*& result) const { CAROM_VERIFY(result == 0 || !result->distributed()); + CAROM_VERIFY(result == 0 || result->getComm() == d_comm); CAROM_VERIFY(distributed() == other.distributed()); + CAROM_VERIFY(d_comm == other.getComm()); CAROM_VERIFY(numRows() == other.dim()); // If the result has not been allocated then do so. Otherwise size it // correctly. if (result == 0) { - result = new Vector(d_num_cols, false); + result = new Vector(d_num_cols, false, d_comm); } else { result->setSize(d_num_cols); @@ -659,7 +697,7 @@ Matrix::transposeMult( d_num_cols, MPI_DOUBLE, MPI_SUM, - MPI_COMM_WORLD); + d_comm); } } @@ -669,7 +707,9 @@ Matrix::transposeMult( Vector& result) const { CAROM_VERIFY(!result.distributed()); + CAROM_VERIFY(d_comm == result.getComm()); CAROM_VERIFY(distributed() == other.distributed()); + CAROM_VERIFY(d_comm == other.getComm()); CAROM_VERIFY(numRows() == other.dim()); // If the result has not been allocated then do so. Otherwise size it @@ -690,7 +730,7 @@ Matrix::transposeMult( d_num_cols, MPI_DOUBLE, MPI_SUM, - MPI_COMM_WORLD); + d_comm); } } @@ -698,12 +738,14 @@ void Matrix::getColumn(int column, Vector*& result) const { + CAROM_VERIFY(result == 0 || result->getComm() == d_comm); + if (result == 0) { if (d_distributed) { - result = new Vector(d_num_rows, true); + result = new Vector(d_num_rows, true, d_comm); } else { - result = new Vector(d_num_rows, false); + result = new Vector(d_num_rows, false, d_comm); } } getColumn(column, *result); @@ -713,6 +755,8 @@ void Matrix::getColumn(int column, Vector& result) const { + CAROM_VERIFY(d_comm == result.getComm()); + result.setSize(d_num_rows); for (int i = 0; i < d_num_rows; i++) { result.item(i) = item(i, column); @@ -727,6 +771,10 @@ Matrix::inverse( (!result->distributed() && result->numRows() == numRows() && result->numColumns() == numColumns())); + // enforce the communicator to be MPI_COMM_WORLD if scalapack routines are used. + // For using local communicator with scalapack routines, we need to figure out how to set up blacs context. + CAROM_VERIFY(d_comm == MPI_COMM_WORLD); + CAROM_VERIFY(result == 0 || result->getComm() == MPI_COMM_WORLD); CAROM_VERIFY(!distributed()); CAROM_VERIFY(numRows() == numColumns()); @@ -776,6 +824,8 @@ Matrix::inverse( { CAROM_VERIFY(!result.distributed() && result.numRows() == numRows() && result.numColumns() == numColumns()); + CAROM_VERIFY(d_comm == MPI_COMM_WORLD); + CAROM_VERIFY(result.getComm() == MPI_COMM_WORLD); CAROM_VERIFY(!distributed()); CAROM_VERIFY(numRows() == numColumns()); @@ -832,6 +882,7 @@ void Matrix::inverse() { CAROM_VERIFY(!distributed()); + CAROM_VERIFY(d_comm == MPI_COMM_WORLD); CAROM_VERIFY(numRows() == numColumns()); // Call lapack routines to do the inversion. @@ -884,8 +935,8 @@ void Matrix::transposePseudoinverse() AtA->inverse(); // Pseudoinverse is (AtA)^{-1}*this^T, but we store the transpose of the result in this, namely this*(AtA)^{-T}. - Vector row(numColumns(), false); - Vector res(numColumns(), false); + Vector row(numColumns(), false, d_comm); + Vector res(numColumns(), false, d_comm); for (int i=0; i row_offsets; int num_total_rows = get_global_offsets(local_num_rows, row_offsets, - MPI_COMM_WORLD); + d_comm); CAROM_VERIFY(num_total_rows == d_num_rows); int local_offset = row_offsets[d_rank] * d_num_cols; const int new_size = local_num_rows * d_num_cols; @@ -1065,7 +1116,7 @@ Matrix::gather() std::vector row_offsets; const int num_total_rows = get_global_offsets(d_num_rows, row_offsets, - MPI_COMM_WORLD); + d_comm); CAROM_VERIFY(num_total_rows == d_num_distributed_rows); const int new_size = d_num_distributed_rows * d_num_cols; @@ -1080,7 +1131,7 @@ Matrix::gather() double *d_new_mat = new double [new_size] {0.0}; CAROM_VERIFY(MPI_Allgatherv(d_mat, d_alloc_size, MPI_DOUBLE, d_new_mat, data_cnts, data_offsets, MPI_DOUBLE, - MPI_COMM_WORLD) == MPI_SUCCESS); + d_comm) == MPI_SUCCESS); delete [] d_mat; delete [] data_offsets, data_cnts; @@ -1101,7 +1152,7 @@ Matrix::calculateNumDistributedRows() { 1, MPI_INT, MPI_SUM, - MPI_COMM_WORLD) == MPI_SUCCESS); + d_comm) == MPI_SUCCESS); d_num_distributed_rows = num_total_rows; } else { @@ -1112,8 +1163,10 @@ Matrix::calculateNumDistributedRows() { Matrix* Matrix::qr_factorize() const { + CAROM_VERIFY(d_comm == MPI_COMM_WORLD); + int myid; - MPI_Comm_rank(MPI_COMM_WORLD, &myid); + MPI_Comm_rank(d_comm, &myid); std::vector row_offset(d_num_procs + 1); row_offset[d_num_procs] = numDistributedRows(); @@ -1125,7 +1178,7 @@ Matrix::qr_factorize() const row_offset.data(), 1, MPI_INT, - MPI_COMM_WORLD) == MPI_SUCCESS); + d_comm) == MPI_SUCCESS); for (int i = d_num_procs - 1; i >= 0; i--) { row_offset[i] = row_offset[i + 1] - row_offset[i]; } @@ -1165,7 +1218,7 @@ Matrix::qr_factorize() const lqcompute(&QRmgr); Matrix* qr_factorized_matrix = new Matrix(row_offset[myid + 1] - row_offset[myid], - numColumns(), distributed()); + numColumns(), distributed(), d_comm); for (int rank = 0; rank < d_num_procs; ++rank) { gather_block(&qr_factorized_matrix->item(0, 0), QRmgr.A, 1, row_offset[rank] + 1, @@ -1204,6 +1257,8 @@ Matrix::qrcp_pivots_transpose_serial(int* row_pivot, int* row_pivot_owner, int pivots_requested) const { + CAROM_VERIFY(d_comm == MPI_COMM_WORLD); + // This method assumes this matrix is serial CAROM_VERIFY(!distributed()); @@ -1265,7 +1320,7 @@ Matrix::qrcp_pivots_transpose_serial(int* row_pivot, CAROM_VERIFY(MPI_Finalized(&is_mpi_finalized) == MPI_SUCCESS); int my_rank = 0; if(is_mpi_initialized && !is_mpi_finalized) { - const MPI_Comm my_comm = MPI_COMM_WORLD; + const MPI_Comm my_comm = d_comm; CAROM_VERIFY(MPI_Comm_rank(my_comm, &my_rank) == MPI_SUCCESS); } @@ -1312,6 +1367,7 @@ Matrix::qrcp_pivots_transpose_distributed_scalapack { // Check if distributed; otherwise, use serial implementation CAROM_VERIFY(distributed()); + CAROM_VERIFY(d_comm == MPI_COMM_WORLD); int num_total_rows = d_num_rows; CAROM_VERIFY(MPI_Allreduce(MPI_IN_PLACE, @@ -1319,13 +1375,13 @@ Matrix::qrcp_pivots_transpose_distributed_scalapack 1, MPI_INT, MPI_SUM, - MPI_COMM_WORLD) == MPI_SUCCESS); + d_comm) == MPI_SUCCESS); int *row_offset = new int[d_num_procs + 1]; row_offset[d_num_procs] = num_total_rows; int my_rank; - MPI_Comm_rank(MPI_COMM_WORLD, &my_rank); + MPI_Comm_rank(d_comm, &my_rank); row_offset[my_rank] = d_num_rows; @@ -1335,7 +1391,7 @@ Matrix::qrcp_pivots_transpose_distributed_scalapack row_offset, 1, MPI_INT, - MPI_COMM_WORLD) == MPI_SUCCESS); + d_comm) == MPI_SUCCESS); for (int i = d_num_procs - 1; i >= 0; i--) { row_offset[i] = row_offset[i + 1] - row_offset[i]; @@ -1378,7 +1434,7 @@ Matrix::qrcp_pivots_transpose_distributed_scalapack int *rcount = (my_rank == 0) ? new int[d_num_procs] : NULL; int *rdisp = (my_rank == 0) ? new int[d_num_procs] : NULL; - MPI_Gather(&scount, 1, MPI_INT, rcount, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Gather(&scount, 1, MPI_INT, rcount, 1, MPI_INT, 0, d_comm); if (my_rank == 0) { @@ -1390,7 +1446,7 @@ Matrix::qrcp_pivots_transpose_distributed_scalapack } MPI_Gatherv(mypivots, scount, MPI_INT, row_pivot, rcount, rdisp, MPI_INT, 0, - MPI_COMM_WORLD); + d_comm); delete [] mypivots; @@ -1490,7 +1546,8 @@ const CAROM_VERIFY(row_pivot_owner != NULL); // Compute total number of rows to set global sizes of matrix - const MPI_Comm comm = MPI_COMM_WORLD; + CAROM_VERIFY(d_comm == MPI_COMM_WORLD); + const MPI_Comm comm = d_comm; const int master_rank = 0; int num_total_rows = d_num_rows; @@ -1644,7 +1701,8 @@ const CAROM_VERIFY(row_pivot_owner != NULL); // Compute total number of rows to set global sizes of matrix - const MPI_Comm comm = MPI_COMM_WORLD; + CAROM_VERIFY(d_comm == MPI_COMM_WORLD); + const MPI_Comm comm = d_comm; const int master_rank = 0; int num_total_rows = d_num_rows; @@ -1803,7 +1861,7 @@ Matrix::orthogonalize(double zero_tol) if (d_distributed && d_num_procs > 1) { CAROM_VERIFY( MPI_Allreduce(MPI_IN_PLACE, &factor, 1, - MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD) == MPI_SUCCESS ); + MPI_DOUBLE, MPI_SUM, d_comm) == MPI_SUCCESS ); } for (int i = 0; i < d_num_rows; ++i) item(i, work) -= factor * item(i, col); @@ -1818,7 +1876,7 @@ Matrix::orthogonalize(double zero_tol) if (d_distributed && d_num_procs > 1) { CAROM_VERIFY( MPI_Allreduce(MPI_IN_PLACE, &norm, 1, MPI_DOUBLE, - MPI_SUM, MPI_COMM_WORLD) == MPI_SUCCESS ); + MPI_SUM, d_comm) == MPI_SUCCESS ); } if (norm > zero_tol) { @@ -1848,7 +1906,7 @@ Matrix::orthogonalize_last(int ncols, double zero_tol) if (d_distributed && d_num_procs > 1) { CAROM_VERIFY( MPI_Allreduce(MPI_IN_PLACE, &factor, 1, MPI_DOUBLE, - MPI_SUM, MPI_COMM_WORLD) == MPI_SUCCESS ); + MPI_SUM, d_comm) == MPI_SUCCESS ); } for (int i = 0; i < d_num_rows; ++i) item(i, last_col) -= factor * item(i, col); @@ -1863,7 +1921,7 @@ Matrix::orthogonalize_last(int ncols, double zero_tol) if (d_distributed && d_num_procs > 1) { CAROM_VERIFY( MPI_Allreduce(MPI_IN_PLACE, &norm, 1, MPI_DOUBLE, - MPI_SUM, MPI_COMM_WORLD) == MPI_SUCCESS ); + MPI_SUM, d_comm) == MPI_SUCCESS ); } if (norm > zero_tol) { @@ -1924,7 +1982,7 @@ Matrix::rescale_cols_max() if (d_distributed && d_num_procs > 1) { MPI_Allreduce(&local_max, &global_max, d_num_cols, MPI_DOUBLE, MPI_MAX, - MPI_COMM_WORLD); + d_comm); } else { @@ -1946,6 +2004,9 @@ Matrix::rescale_cols_max() Matrix outerProduct(const Vector &v, const Vector &w) { + CAROM_VERIFY(v.getComm() == w.getComm()); + MPI_Comm comm = v.getComm(); + /* * There are two cases of concern: * @@ -1957,7 +2018,7 @@ Matrix outerProduct(const Vector &v, const Vector &w) int result_num_rows = v.dim(); int result_num_cols; bool is_distributed = v.distributed(); - Vector gathered_w; + Vector gathered_w(comm); /* * Gather all of the entries in w on each process into a Vector stored @@ -1977,7 +2038,6 @@ Matrix outerProduct(const Vector &v, const Vector &w) // but the storage for std::vector containers must be contiguous // as defined by the C++ standard. int process_local_num_cols = w.dim(); - MPI_Comm comm = MPI_COMM_WORLD; MPI_Datatype num_cols_datatype = MPI_INT; int num_procs; MPI_Comm_size(comm, &num_procs); @@ -2030,7 +2090,7 @@ Matrix outerProduct(const Vector &v, const Vector &w) } /* Create the matrix */ - Matrix result(result_num_rows, result_num_cols, is_distributed); + Matrix result(result_num_rows, result_num_cols, is_distributed, comm); /* Compute the outer product using the gathered copy of w. */ for (int i = 0; i < result_num_rows; i++) @@ -2046,6 +2106,8 @@ Matrix outerProduct(const Vector &v, const Vector &w) Matrix DiagonalMatrixFactory(const Vector &v) { + const MPI_Comm comm = v.getComm(); + const int resultNumRows = v.dim(); int resultNumColumns, processRowStartIndex, processNumColumns; const bool isDistributed = v.distributed(); @@ -2065,7 +2127,6 @@ Matrix DiagonalMatrixFactory(const Vector &v) * counts must be summed to get the number of columns on each * process. */ - const MPI_Comm comm = MPI_COMM_WORLD; int numProcesses; MPI_Comm_size(comm, &numProcesses); std::vector @@ -2103,7 +2164,7 @@ Matrix DiagonalMatrixFactory(const Vector &v) * Create the diagonal matrix and assign process local entries in v * to process local entries in the diagonal matrix. */ - Matrix diagonalMatrix(resultNumRows, resultNumColumns, isDistributed); + Matrix diagonalMatrix(resultNumRows, resultNumColumns, isDistributed, comm); for (int i = 0; i < resultNumRows; i++) { for (int j = 0; j < resultNumColumns; j++) @@ -2129,6 +2190,8 @@ Matrix IdentityMatrixFactory(const Vector &v) struct EigenPair SymmetricRightEigenSolve(Matrix* A) { + CAROM_VERIFY(A->getComm() == MPI_COMM_WORLD); + char jobz = 'V', uplo = 'U'; int info; @@ -2176,6 +2239,9 @@ struct EigenPair SymmetricRightEigenSolve(Matrix* A) struct ComplexEigenPair NonSymmetricRightEigenSolve(Matrix* A) { + CAROM_VERIFY(A->getComm() == MPI_COMM_WORLD); + MPI_Comm comm = A->getComm(); + char jobvl = 'N', jobrl = 'V'; int info; @@ -2185,7 +2251,7 @@ struct ComplexEigenPair NonSymmetricRightEigenSolve(Matrix* A) double* e_real = new double [k]; double* e_imaginary = new double [k]; double* ev_l = NULL; - Matrix* ev_r = new Matrix(k, k, false); + Matrix* ev_r = new Matrix(k, k, false, comm); Matrix* A_copy = new Matrix(*A); // A now in a row major representation. Put it @@ -2213,8 +2279,8 @@ struct ComplexEigenPair NonSymmetricRightEigenSolve(Matrix* A) } ComplexEigenPair eigenpair; - eigenpair.ev_real = new Matrix(k, k, false); - eigenpair.ev_imaginary = new Matrix(k, k, false); + eigenpair.ev_real = new Matrix(k, k, false, comm); + eigenpair.ev_imaginary = new Matrix(k, k, false, comm); // Separate lapack eigenvector into real and imaginary parts for (int i = 0; i < k; ++i) @@ -2255,7 +2321,13 @@ void SerialSVD(Matrix* A, Vector* S, Matrix* V) { + const MPI_Comm comm = A->getComm(); + CAROM_VERIFY(comm == MPI_COMM_WORLD); + CAROM_VERIFY(comm == U->getComm()); + CAROM_VERIFY(comm == S->getComm()); + CAROM_VERIFY(comm == V->getComm()); CAROM_VERIFY(!A->distributed()); + int m = A->numRows(); int n = A->numColumns(); @@ -2335,6 +2407,10 @@ Matrix* SpaceTimeProduct(const CAROM::Matrix* As, const CAROM::Matrix* At, { // TODO: implement reduction in parallel for the spatial matrices CAROM_VERIFY(As->distributed() && Bs->distributed()); + const MPI_Comm comm = As->getComm(); + CAROM_VERIFY(comm == Bs->getComm()); + CAROM_VERIFY(comm == At->getComm()); + CAROM_VERIFY(comm == Bt->getComm()); const int AtOS = At0at0 ? 1 : 0; const int BtOS0 = Bt0at0 ? 1 : 0; @@ -2358,7 +2434,7 @@ Matrix* SpaceTimeProduct(const CAROM::Matrix* As, const CAROM::Matrix* At, CAROM_VERIFY(tscale == NULL || (tscale->size() == ntime-1 && k0 > 0)); - Matrix* p = new CAROM::Matrix(nrows, ncols, false); + Matrix* p = new CAROM::Matrix(nrows, ncols, false, comm); for (int i=0; i