Parallel reading of serial basis files (#316)

dylan-copeland · web-flow · commit 80aa762b128b · 2025-09-16T13:27:39.000-07:00
diff --git a/lib/linalg/BasisReader.cpp b/lib/linalg/BasisReader.cpp
@@ -23,19 +23,22 @@ namespace CAROM {
 BasisReader::BasisReader(
     const std::string& base_file_name,
     Database::formats db_format,
-    const int dim) :
+    const int dim,
+    MPI_Comm comm) :
     d_dim(dim),
     full_file_name(""),
     base_file_name_(base_file_name),
     d_format(db_format)
 {
-    CAROM_ASSERT(!base_file_name.empty());
+    CAROM_VERIFY(!base_file_name.empty());
+    CAROM_VERIFY(comm == MPI_COMM_NULL || comm == MPI_COMM_WORLD);
+    d_distributed = comm != MPI_COMM_NULL;
 
     int mpi_init;
     MPI_Initialized(&mpi_init);
     int rank;
-    if (mpi_init) {
-        MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+    if (mpi_init && d_distributed) {
+        MPI_Comm_rank(comm, &rank);
     }
     else {
         rank = 0;
@@ -55,15 +58,15 @@ BasisReader::BasisReader(
             We allow 0 local dimension. (global dimension still needs to be positive)
         */
         std::vector<int> tmp;
-        d_global_dim = get_global_offsets(d_dim, tmp, MPI_COMM_WORLD);
+        d_global_dim = get_global_offsets(d_dim, tmp, comm);
         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", MPI_COMM_WORLD);
+    d_database->open(full_file_name, "r", comm);
 }
 
 BasisReader::~BasisReader()
@@ -78,12 +81,12 @@ BasisReader::getSpatialBasis()
     int num_rows = getDim("basis");
     int num_cols = getNumSamples("basis");
 
-    Matrix* spatial_basis_vectors = new Matrix(num_rows, num_cols, true);
+    Matrix* spatial_basis_vectors = new Matrix(num_rows, num_cols, d_distributed);
 
     d_database->getDoubleArray("spatial_basis",
                                &spatial_basis_vectors->item(0, 0),
                                num_rows*num_cols,
-                               true);
+                               d_distributed);
     return std::unique_ptr<Matrix>(spatial_basis_vectors);
 }
 
@@ -107,15 +110,16 @@ BasisReader::getSpatialBasis(
     CAROM_VERIFY(start_col <= end_col && end_col <= num_cols);
     int num_cols_to_read = end_col - start_col + 1;
 
-    Matrix* spatial_basis_vectors = new Matrix(num_rows, num_cols_to_read, true);
+    Matrix* spatial_basis_vectors = new Matrix(num_rows, num_cols_to_read,
+            d_distributed);
     sprintf(tmp, "spatial_basis");
     d_database->getDoubleArray(tmp,
                                &spatial_basis_vectors->item(0, 0),
                                num_rows*num_cols_to_read,
                                start_col - 1,
                                num_cols_to_read,
                                num_cols,
-                               true);
+                               d_distributed);
     return std::unique_ptr<Matrix>(spatial_basis_vectors);
 }
 
@@ -156,7 +160,7 @@ BasisReader::getTemporalBasis()
     sprintf(tmp, "temporal_basis");
     d_database->getDoubleArray(tmp,
                                &temporal_basis_vectors->item(0, 0),
-                               num_rows*num_cols);
+                               num_rows*num_cols, d_distributed);
     return std::unique_ptr<Matrix>(temporal_basis_vectors);
 }
 
@@ -187,7 +191,8 @@ BasisReader::getTemporalBasis(
                                num_rows*num_cols_to_read,
                                start_col - 1,
                                num_cols_to_read,
-                               num_cols);
+                               num_cols,
+                               d_distributed);
     return std::unique_ptr<Matrix>(temporal_basis_vectors);
 }
 
@@ -229,7 +234,8 @@ BasisReader::getSingularValues()
     sprintf(tmp, "singular_value");
     d_database->getDoubleArray(tmp,
                                &singular_values->item(0),
-                               size);
+                               size,
+                               d_distributed);
     return std::unique_ptr<Vector>(singular_values);
 }
 
@@ -324,12 +330,12 @@ BasisReader::getSnapshotMatrix()
     int num_cols = getNumSamples("snapshot");
 
     char tmp[100];
-    Matrix* snapshots = new Matrix(num_rows, num_cols, true);
+    Matrix* snapshots = new Matrix(num_rows, num_cols, d_distributed);
     sprintf(tmp, "snapshot_matrix");
     d_database->getDoubleArray(tmp,
                                &snapshots->item(0, 0),
                                num_rows*num_cols,
-                               true);
+                               d_distributed);
     return std::unique_ptr<Matrix>(snapshots);
 }
 
@@ -353,15 +359,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, true);
+    Matrix* snapshots = new Matrix(num_rows, num_cols_to_read, d_distributed);
     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,
-                               true);
+                               d_distributed);
     return std::unique_ptr<Matrix>(snapshots);
 }
 }
diff --git a/lib/linalg/BasisReader.h b/lib/linalg/BasisReader.h
@@ -51,7 +51,8 @@ class BasisReader {
     BasisReader(
         const std::string& base_file_name,
         Database::formats db_format = Database::formats::HDF5,
-        const int dim = -1);
+        const int dim = -1,
+        MPI_Comm comm = MPI_COMM_WORLD);
 
     /**
      * @brief Destructor.
@@ -314,6 +315,8 @@ class BasisReader {
      * If negative, use the dimension from the rank-specific local file.
      */
     int d_global_dim;
+
+    bool d_distributed;
 };
 
 }
diff --git a/unit_tests/CMakeLists.txt b/unit_tests/CMakeLists.txt
@@ -88,6 +88,8 @@ endforeach()
 # Copy testing inputs and scripts to build/tests/
 file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/baselines/basis_conversion/
      DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/baselines/basis_conversion/")
+file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/basis_data/
+     DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/basis_data/")
 file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/s_opt_data/
      DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/s_opt_data/")
 
diff --git a/unit_tests/basis_data/basis.000000 b/unit_tests/basis_data/basis.000000
diff --git a/unit_tests/test_HDFDatabase.cpp b/unit_tests/test_HDFDatabase.cpp
@@ -543,6 +543,35 @@ TEST(HDF5, Test_selective_parallel_reading)
 
 #endif
 
+TEST(HDF5, Test_serial_file_parallel_reading)
+{
+    // Read the matrix on all ranks.
+    CAROM::BasisReader reader("./basis_data/basis.000000",
+                              CAROM::Database::formats::HDF5, -1,
+                              MPI_COMM_NULL);
+    std::unique_ptr<CAROM::Matrix> B = reader.getSpatialBasis();
+
+    EXPECT_EQ(B->numRows(), 4);
+    EXPECT_EQ(B->numColumns(), 2);
+    EXPECT_FALSE(B->distributed());
+
+    double error = 0.0;
+    for (int i=0; i<B->numRows(); ++i)
+        for (int j=0; j<B->numColumns(); ++j)
+        {
+            const double B_ij = (i == 1 && j == 0) || (i == 2 && j == 1) ? -1.0 : 0.0;
+            error = std::max(error, std::abs(B_ij - (*B)(i, j)));
+        }
+
+    std::unique_ptr<CAROM::Vector> sv = reader.getSingularValues();
+    EXPECT_EQ(sv->dim(), 2);
+
+    for (int i=0; i<sv->dim(); ++i)
+        error = std::max(error, std::abs((*sv)(i) - 1.0));
+
+    EXPECT_NEAR(error, 0.0, threshold);
+}
+
 TEST(BasisGeneratorIO, HDFDatabase)
 {
     // Get the rank of this process, and the number of processors.
