From e3878959511c0acf8328354d1914623f0ffa0a46 Mon Sep 17 00:00:00 2001 From: Marco Feder Date: Sat, 24 May 2025 18:07:41 +0200 Subject: [PATCH] Add MPI communicator to SparseDirectMUMPS constructor. Use integers from MUMPS --- include/deal.II/lac/sparse_direct.h | 29 ++++++++++++++++++++++++----- source/lac/sparse_direct.cc | 22 ++++++++++------------ 2 files changed, 34 insertions(+), 17 deletions(-) diff --git a/include/deal.II/lac/sparse_direct.h b/include/deal.II/lac/sparse_direct.h index 236230f1ae..c1ee01d1c7 100644 --- a/include/deal.II/lac/sparse_direct.h +++ b/include/deal.II/lac/sparse_direct.h @@ -48,6 +48,16 @@ namespace types #else using suitesparse_index = long int; #endif + +#ifdef DEAL_II_WITH_MUMPS + using mumps_index = MUMPS_INT; + using mumps_nnz = MUMPS_INT8; +#else + using mumps_index = int; + using mumps_nnz = std::size_t; +#endif + + } // namespace types /** @@ -475,12 +485,12 @@ public: : blr_ucfs(blr_ucfs) , lowrank_threshold(lowrank_threshold) {} - /** * If true, the Block Low-Rank approximation is used with the UCFS * algorithm, an algorithm with higher compression than the standard one. */ bool blr_ucfs; + /** * Threshold for the low-rank truncation of the blocks. */ @@ -510,11 +520,13 @@ public: * If true, the MUMPS solver will print out error statistics. */ bool error_statistics; + /* * If true, the MUMPS solver will use the symmetric factorization. This is * only possible if the matrix is symmetric. */ bool symmetric; + /* * If true, the MUMPS solver will use the positive definite factorization. * This is only possible if the matrix is symmetric and positive definite. @@ -536,7 +548,8 @@ public: /** * Constructor, takes AdditionalData to control MUMPS execution. */ - SparseDirectMUMPS(const AdditionalData &additional_data = AdditionalData()); + SparseDirectMUMPS(const AdditionalData &additional_data = AdditionalData(), + const MPI_Comm &communicator = MPI_COMM_WORLD); /** * Destructor. @@ -590,6 +603,7 @@ public: private: #ifdef DEAL_II_WITH_MUMPS mutable DMUMPS_STRUC_C id; + #endif // DEAL_II_WITH_MUMPS /** @@ -606,12 +620,12 @@ private: /** * irn contains the row indices of the non-zero entries of the matrix. */ - std::unique_ptr irn; + std::unique_ptr irn; /** * jcn contains the column indices of the non-zero entries of the matrix. */ - std::unique_ptr jcn; + std::unique_ptr jcn; /** * The number of rows of the matrix. The matrix is square. @@ -621,7 +635,7 @@ private: /** * The number of non-zero entries in the matrix. */ - std::size_t nnz; + types::mumps_nnz nnz; /** * This function hands over to MUMPS the system's matrix. @@ -646,6 +660,11 @@ private: * Struct that holds the additional data for the MUMPS solver. */ AdditionalData additional_data; + + /** + * MPI_Comm object for the MUMPS solver. + */ + const MPI_Comm mpi_communicator; }; DEAL_II_NAMESPACE_CLOSE diff --git a/source/lac/sparse_direct.cc b/source/lac/sparse_direct.cc index f4d6e6377f..b4b695b75c 100644 --- a/source/lac/sparse_direct.cc +++ b/source/lac/sparse_direct.cc @@ -861,9 +861,11 @@ SparseDirectUMFPACK::n() const #ifdef DEAL_II_WITH_MUMPS -SparseDirectMUMPS::SparseDirectMUMPS(const AdditionalData &data) +SparseDirectMUMPS::SparseDirectMUMPS(const AdditionalData &data, + const MPI_Comm &communicator) + : additional_data(data) + , mpi_communicator(communicator) { - this->additional_data = data; // Initialize MUMPS instance: id.job = -1; id.par = 1; @@ -883,8 +885,7 @@ SparseDirectMUMPS::SparseDirectMUMPS(const AdditionalData &data) else id.sym = 0; - // Use MPI_COMM_WORLD as communicator - id.comm_fortran = (MUMPS_INT)MPI_Comm_c2f(MPI_COMM_WORLD); + id.comm_fortran = (MUMPS_INT)MPI_Comm_c2f(mpi_communicator); dmumps_c(&id); if (additional_data.output_details == false) @@ -926,12 +927,9 @@ void SparseDirectMUMPS::initialize_matrix(const Matrix &matrix) { Assert(matrix.n() == matrix.m(), ExcMessage("Matrix needs to be square.")); - // Here we should be checking if the matrix is respecting the symmetry given - //(I.E. sym = 0 for non-symmetric matrix, sym = 1 for posdef matrix, sym = 2 - // for general symmetric). // Hand over matrix to MUMPS as centralized assembled matrix - if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0) { // Set number of unknowns n = matrix.n(); @@ -945,8 +943,8 @@ SparseDirectMUMPS::initialize_matrix(const Matrix &matrix) // matrix indices pointing to the row and column dimensions // respectively of the matrix representation above (a): ie. a[k] is // the matrix element (irn[k], jcn[k]) - irn = std::make_unique(nnz); - jcn = std::make_unique(nnz); + irn = std::make_unique(nnz); + jcn = std::make_unique(nnz); size_type n_non_zero_elements = 0; @@ -1001,7 +999,7 @@ SparseDirectMUMPS::copy_rhs_to_mumps(const Vector &new_rhs) const Assert(n == new_rhs.size(), ExcMessage("Matrix size and rhs length must be equal.")); - if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0) { rhs.resize(n); for (size_type i = 0; i < n; ++i) @@ -1022,7 +1020,7 @@ SparseDirectMUMPS::copy_solution(Vector &vector) const ExcMessage("Class not initialized with a rhs vector.")); // Copy solution into the given vector - if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0) { for (size_type i = 0; i < n; ++i) vector(i) = rhs[i]; -- 2.39.5