From 0016812ac4eda5bb97c8aa98c6004e389e10f556 Mon Sep 17 00:00:00 2001 From: Marco Feder Date: Fri, 6 Jun 2025 14:21:44 +0200 Subject: [PATCH] Address some comments --- include/deal.II/lac/sparse_direct.h | 6 ++-- source/lac/sparse_direct.cc | 45 +++++++++++++++++++---------- tests/mumps/mumps_06.cc | 4 +-- 3 files changed, 34 insertions(+), 21 deletions(-) diff --git a/include/deal.II/lac/sparse_direct.h b/include/deal.II/lac/sparse_direct.h index 6f0f46232e..2a75a0493c 100644 --- a/include/deal.II/lac/sparse_direct.h +++ b/include/deal.II/lac/sparse_direct.h @@ -526,12 +526,12 @@ public: {} /** - * Print details of the execution + * If true, the MUMPS solver will print out details of the execution. */ bool output_details; /** - * Print error statistics + * If true, the MUMPS solver will print out error statistics. */ bool error_statistics; @@ -672,7 +672,7 @@ private: /** * IndexSet storing the locally owned rows of the matrix. */ - IndexSet row_range; + IndexSet locally_owned_rows; /** * This function initializes a MUMPS instance and hands over the system's diff --git a/source/lac/sparse_direct.cc b/source/lac/sparse_direct.cc index 2d02cd554b..3bbf2ceb92 100644 --- a/source/lac/sparse_direct.cc +++ b/source/lac/sparse_direct.cc @@ -931,15 +931,14 @@ SparseDirectMUMPS::initialize_matrix(const Matrix &matrix) n = matrix.n(); id.n = n; - if constexpr (std::is_same_v, Matrix>) + if constexpr (std::is_same_v>) { // Serial matrix: hand over matrix to MUMPS as centralized assembled // matrix if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0) { // number of nonzero elements in matrix - if constexpr (std::is_same_v>) - nnz = matrix.n_actually_nonzero_elements(); + nnz = matrix.n_actually_nonzero_elements(); // representation of the matrix a = std::make_unique(nnz); @@ -994,9 +993,17 @@ SparseDirectMUMPS::initialize_matrix(const Matrix &matrix) id.a = a.get(); } } - else if constexpr (std::is_same_v || - std::is_same_v) + else if constexpr (std::is_same_v || + std::is_same_v) { + int result; + MPI_Comm_compare(mpi_communicator, + matrix.get_mpi_communicator(), + &result); + AssertThrow(result == MPI_IDENT, + ExcMessage("The matrix communicator must match the MUMPS " + "communicator.")); + // Distributed matrix case id.icntl[17] = 3; // distributed matrix assembly id.nnz = matrix.n_nonzero_elements(); @@ -1004,34 +1011,36 @@ SparseDirectMUMPS::initialize_matrix(const Matrix &matrix) size_type n_non_zero_local = 0; // Get the range of rows owned by this process - row_range = matrix.locally_owned_range_indices(); + locally_owned_rows = matrix.locally_owned_range_indices(); size_type local_non_zeros = 0; - if constexpr (std::is_same_v) + if constexpr (std::is_same_v) { const auto &trilinos_matrix = matrix.trilinos_matrix(); local_non_zeros = trilinos_matrix.NumMyNonzeros(); } - else if constexpr (std::is_same_v) + else if constexpr (std::is_same_v) { Mat &petsc_matrix = const_cast(matrix) .petsc_matrix(); MatInfo info; MatGetInfo(petsc_matrix, MAT_LOCAL, &info); - local_non_zeros = (PetscInt)info.nz_used; + local_non_zeros = (size_type)info.nz_used; } + // We allocate enough entries for the general, nonsymmetric, case. If case + // of a symmetric matrix, we will end up with fewer entries. irn = std::make_unique(local_non_zeros); jcn = std::make_unique(local_non_zeros); a = std::make_unique(local_non_zeros); - irhs_loc.resize(row_range.n_elements()); + irhs_loc.resize(locally_owned_rows.n_elements()); if (additional_data.symmetric == true) { - for (const auto &row : row_range) + for (const auto &row : locally_owned_rows) { for (auto it = matrix.begin(row); it != matrix.end(row); ++it) if (std::abs(it->value()) > 0.0 && it->column() >= row) @@ -1049,13 +1058,13 @@ SparseDirectMUMPS::initialize_matrix(const Matrix &matrix) } const types::global_cell_index local_index = - row_range.index_within_set(row); + locally_owned_rows.index_within_set(row); irhs_loc[local_index] = row + 1; } } else { - for (const auto &row : row_range) + for (const auto &row : locally_owned_rows) { // Loop over columns for (auto it = matrix.begin(row); it != matrix.end(row); ++it) @@ -1075,7 +1084,7 @@ SparseDirectMUMPS::initialize_matrix(const Matrix &matrix) } const types::global_cell_index local_index = - row_range.index_within_set(row); + locally_owned_rows.index_within_set(row); irhs_loc[local_index] = row + 1; } } @@ -1092,7 +1101,7 @@ SparseDirectMUMPS::initialize_matrix(const Matrix &matrix) id.icntl[20] = 0; // centralized solution, stored on rank 0 by MUMPS id.nrhs = 1; id.lrhs_loc = n; - id.nloc_rhs = row_range.n_elements(); + id.nloc_rhs = locally_owned_rows.n_elements(); } else { @@ -1220,6 +1229,10 @@ SparseDirectMUMPS::vmult(VectorType &dst, const VectorType &src) const rhs.resize(0); // remove rhs again } + else + { + DEAL_II_NOT_IMPLEMENTED(); + } } diff --git a/tests/mumps/mumps_06.cc b/tests/mumps/mumps_06.cc index b07422dee4..d899174489 100644 --- a/tests/mumps/mumps_06.cc +++ b/tests/mumps/mumps_06.cc @@ -56,7 +56,7 @@ solve_and_check(const MatrixType &M, data.output_details = false; data.symmetric = true; data.posdef = true; - SparseDirectMUMPS solver(data); + SparseDirectMUMPS solver(data, M.get_mpi_communicator()); solver.initialize(M); VectorType dst(rhs); solver.vmult(dst, rhs); @@ -110,7 +110,7 @@ test() SparseDirectMUMPS::AdditionalData data; data.output_details = true; data.error_statistics = true; - SparseDirectMUMPS Binv(data); + SparseDirectMUMPS Binv(data, B.get_mpi_communicator()); Binv.initialize(B); // for a number of different solution -- 2.39.5