Encapsulate MUMPS solver class

include/LinearAlgebra/Solvers/coo_mumps_solver.h

@@ -0,0 +1,97 @@
+#pragma once
+
+#ifdef GMGPOLAR_USE_MUMPS
+
+    #include "dmumps_c.h"
+    #include "../../LinearAlgebra/Matrix/coo_matrix.h"
+    #include "../../LinearAlgebra/Vector/vector.h"
+
+/*
+ * CooMumpsSolver - sparse direct solver backed by MUMPS (MUltifrontal Massively Parallel Solver).
+ *
+ * Accepts a matrix in COO (Coordinate) format and solves systems of the form Ax = b.
+ * The solver takes ownership of the matrix, as MUMPS may reorder and modify it in place
+ * during the analysis and factorization phases.
+ *
+ * Usage:
+ *   1. Construct with a fully assembled COO matrix (all non-zeros must be present,
+ *      including both triangles if the matrix is symmetric).
+ *   2. Call solve(rhs) one or more times. The right-hand side vector is overwritten
+ *      in place with the solution x on return.
+ *
+ * Symmetry:
+ *   If the matrix is marked symmetric (is_symmetric = true), it is assumed to be
+ *   symmetric positive definite and MUMPS is configured accordingly (SYM = 1).
+ *   In that case, only the upper triangle is passed to MUMPS internally.
+ *   This assumption is valid within GMGPolar because all system matrices arise from
+ *   an invertible domain mapping, guaranteeing positive definiteness.
+ *
+ * Thread safety:
+ *   Not thread-safe. Each instance manages its own MUMPS communicator and
+ *   must not be shared across threads.
+ */
+class CooMumpsSolver
+{
+public:
+    explicit CooMumpsSolver(SparseMatrixCOO<double> matrix);
+    ~CooMumpsSolver();
+
+    // rhs is overwritten in-place with the solution on return.
+    void solve(Vector<double>& rhs);
+
+private:
+    void initialize();
+    void finalize();
+    void configureICNTL();
+    void configureCNTL();
+
+    /* ------------------------------------------------ */
+    /* MUMPS uses 1-based indexing in the documentation */
+    /* ------------------------------------------------ */
+    int& ICNTL(int i)
+    {
+        return mumps_solver_.icntl[i - 1];
+    }
+    double& CNTL(int i)
+    {
+        return mumps_solver_.cntl[i - 1];
+    }
+    int& INFOG(int i)
+    {
+        return mumps_solver_.infog[i - 1];
+    }
+
+    SparseMatrixCOO<double> extractUpperTriangle(const SparseMatrixCOO<double>& matrix) const;
+
+    /* ----------------------------------- */
+    /* MUMPS jobs and constant definitions */
+    /* ----------------------------------- */
+    static constexpr int USE_COMM_WORLD   = -987654;
+    static constexpr int PAR_NOT_PARALLEL = 0;
+    static constexpr int PAR_PARALLEL     = 1;
+
+    static constexpr int JOB_INIT               = -1;
+    static constexpr int JOB_END                = -2;
+    static constexpr int JOB_REMOVE_SAVED_DATA  = -3;
+    static constexpr int JOB_FREE_INTERNAL_DATA = -4;
+    static constexpr int JOB_SUPPRESS_OOC_FILES = -200;
+
+    static constexpr int JOB_ANALYSIS_PHASE                  = 1;
+    static constexpr int JOB_FACTORIZATION_PHASE             = 2;
+    static constexpr int JOB_COMPUTE_SOLUTION                = 3;
+    static constexpr int JOB_ANALYSIS_AND_FACTORIZATION      = 4;
+    static constexpr int JOB_FACTORIZATION_AND_SOLUTION      = 5;
+    static constexpr int JOB_ANALYSIS_FACTORIZATION_SOLUTION = 6;
+    static constexpr int JOB_SAVE_INTERNAL_DATA              = 7;
+    static constexpr int JOB_RESTORE_INTERNAL_DATA           = 8;
+    static constexpr int JOB_DISTRIBUTE_RHS                  = 9;
+
+    static constexpr int SYM_UNSYMMETRIC       = 0;
+    static constexpr int SYM_POSITIVE_DEFINITE = 1;
+    static constexpr int SYM_GENERAL_SYMMETRIC = 2;
+
+    SparseMatrixCOO<double> matrix_;
+    DMUMPS_STRUC_C mumps_solver_ = {};
+};
+
+#endif // GMGPOLAR_USE_MUMPS

src/CMakeLists.txt

@@ -16,6 +16,13 @@ set(POLAR_GRID_SOURCES
     ${CMAKE_CURRENT_SOURCE_DIR}/PolarGrid/anisotropic_division.cpp
 )
 
+# Gather all source files
+# file(GLOB_RECURSE LINEAR_ALGEBRA_SOURCES 
+# ${CMAKE_CURRENT_SOURCE_DIR}/LinearAlgebra/Solvers/*.cpp)
+set(LINEAR_ALGEBRA_SOURCES
+    ${CMAKE_CURRENT_SOURCE_DIR}/LinearAlgebra/Solvers/coo_mumps_solver.cpp
+)
+
 # file(GLOB_RECURSE GMG_POLAR_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/*.cpp)
 set(GMG_POLAR_SOURCES
     ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/gmgpolar.cpp
@@ -47,6 +54,7 @@ set(CONFIG_PARSER_SOURCES
 # Create the main library
 add_library(GMGPolarLib STATIC
     ${POLAR_GRID_SOURCES}
+    ${LINEAR_ALGEBRA_SOURCES}
     ${GMG_POLAR_SOURCES}
     ${STENCIL_SOURCES}
     ${INTERPOLATION_SOURCES}

src/LinearAlgebra/Solvers/coo_mumps_solver.cpp

@@ -0,0 +1,196 @@
+#include "../../../include/LinearAlgebra/Solvers/coo_mumps_solver.h"
+
+#ifdef GMGPOLAR_USE_MUMPS
+
+    #include <iostream>
+    #include <cassert>
+    #include <stdexcept>
+
+CooMumpsSolver::CooMumpsSolver(SparseMatrixCOO<double> matrix)
+{
+    if (matrix.is_symmetric()) {
+        matrix_ = extractUpperTriangle(matrix);
+    }
+    else {
+        matrix_ = std::move(matrix);
+    }
+
+    initialize();
+}
+
+CooMumpsSolver::~CooMumpsSolver()
+{
+    finalize();
+}
+
+void CooMumpsSolver::solve(Vector<double>& rhs)
+{
+    assert(std::ssize(rhs) == mumps_solver_.n);
+
+    mumps_solver_.job  = JOB_COMPUTE_SOLUTION;
+    mumps_solver_.nrhs = 1;
+    mumps_solver_.lrhs = mumps_solver_.n;
+    mumps_solver_.rhs  = rhs.data();
+
+    dmumps_c(&mumps_solver_);
+
+    if (INFOG(1) != 0) {
+        std::cerr << "MUMPS reported an error during solution phase " << "(INFOG(1) = " << INFOG(1) << ").\n";
+    }
+}
+
+void CooMumpsSolver::initialize()
+{
+    assert(matrix_.rows() == matrix_.columns());
+
+    /*
+     * MUMPS uses 1-based indexing; convert from 0-based.
+     */
+    for (int i = 0; i < matrix_.non_zero_size(); i++) {
+        matrix_.row_index(i) += 1;
+        matrix_.col_index(i) += 1;
+    }
+
+    mumps_solver_.job          = JOB_INIT;
+    mumps_solver_.par          = PAR_PARALLEL;
+    mumps_solver_.sym          = matrix_.is_symmetric() ? SYM_POSITIVE_DEFINITE : SYM_UNSYMMETRIC;
+    mumps_solver_.comm_fortran = USE_COMM_WORLD;
+
+    dmumps_c(&mumps_solver_);
+
+    configureICNTL();
+    configureCNTL();
+
+    mumps_solver_.job = JOB_ANALYSIS_AND_FACTORIZATION;
+    mumps_solver_.n   = matrix_.rows();
+    mumps_solver_.nz  = matrix_.non_zero_size();
+    mumps_solver_.irn = matrix_.row_indices_data();
+    mumps_solver_.jcn = matrix_.column_indices_data();
+    mumps_solver_.a   = matrix_.values_data();
+
+    dmumps_c(&mumps_solver_);
+
+    if (INFOG(1) != 0) {
+        std::cerr << "MUMPS reported an error during analysis/factorization " << "(INFOG(1) = " << INFOG(1) << ").\n";
+        return;
+    }
+
+    if (mumps_solver_.sym == SYM_POSITIVE_DEFINITE && INFOG(12) != 0) {
+        std::cerr << "Matrix declared positive definite, "
+                  << "but negative pivots were encountered during factorization " << "(INFOG(12) = " << INFOG(12)
+                  << ").\n";
+    }
+}
+
+void CooMumpsSolver::finalize()
+{
+    mumps_solver_.job = JOB_END;
+    dmumps_c(&mumps_solver_);
+}
+
+void CooMumpsSolver::configureICNTL()
+{
+    // All ICNTL values are left at their defaults unless noted below.
+    // ICNTL(1) = 0: suppress error message output
+    // ICNTL(3) = 0: suppress global information output
+    // ICNTL(6) = 7: automatically choose permutation/scaling strategy
+    // ICNTL(7) = 5: use METIS for fill-reducing ordering
+    // ICNTL(48) = 0: disable tree parallelism (conflicts with OpenMP in newer MUMPS releases)
+
+    ICNTL(1)  = 0; // Output stream for error messages
+    ICNTL(2)  = 0; // Output stream for diagnostic printing and statistics local to each MPI process
+    ICNTL(3)  = 0; // Output stream for global information, collected on the host
+    ICNTL(4)  = 0; // Level of printing for error, warning, and diagnostic messages
+    ICNTL(5)  = 0; // Controls the matrix input format
+    ICNTL(6)  = 7; // Permutes the matrix to a zero-free diagonal and/or scales the matrix
+    ICNTL(7)  = 5; // Symmetric permutation (ordering) to determine pivot order for sequential analysis
+    ICNTL(8)  = 77; // Scaling strategy
+    ICNTL(9)  = 1; // Computes the solution using A or A^T
+    ICNTL(10) = 0; // Iterative refinement steps applied to the computed solution
+    ICNTL(11) = 0; // Error analysis statistics
+    ICNTL(12) = 0; // Ordering strategy for symmetric matrices
+    ICNTL(13) = 0; // Controls the parallelism of the root node
+    ICNTL(14) = matrix_.is_symmetric() ? 5 : 20; // Percentage increase in estimated working space
+    ICNTL(15) = 0; // Exploits compression of the input matrix resulting from a block format
+    ICNTL(16) = 0; // Controls the setting of the number of OpenMP threads
+    // ICNTL(17) does not exist
+    ICNTL(18) = 0; // Strategy for the distributed input matrix
+    ICNTL(19) = 0; // Computes the Schur complement matrix
+    ICNTL(20) = 0; // Format of the right-hand sides (dense, sparse, or distributed)
+    ICNTL(21) = 0; // Distribution of the solution vectors (centralized or distributed)
+    ICNTL(22) = 0; // In-core/out-of-core (OOC) factorization and solve
+    ICNTL(23) = 0; // Maximum working memory in MegaBytes per working process
+    ICNTL(24) = 0; // Detection of null pivot rows
+    ICNTL(25) = 0; // Solution of a deficient matrix and null space basis computation
+    ICNTL(26) = 0; // Solution phase when Schur complement has been computed
+    ICNTL(27) = -32; // Blocking size for multiple right-hand sides
+    ICNTL(28) = 0; // Sequential or parallel computation of the ordering
+    ICNTL(29) = 0; // Parallel ordering tool when ICNTL(28)=1
+    ICNTL(30) = 0; // User-specified entries in the inverse A^-1
+    ICNTL(31) = 0; // Which factors may be discarded during factorization
+    ICNTL(32) = 0; // Forward elimination of the right-hand sides during factorization
+    ICNTL(33) = 0; // Computes the determinant of the input matrix
+    ICNTL(34) = 0; // Conservation of OOC files during JOB=-3
+    ICNTL(35) = 0; // Activation of the BLR feature
+    ICNTL(36) = 0; // Choice of BLR factorization variant
+    ICNTL(37) = 0; // BLR compression of the contribution blocks
+    ICNTL(38) = 600; // Estimated compression rate of LU factors
+    ICNTL(39) = 500; // Estimated compression rate of contribution blocks
+    // ICNTL(40-47) do not exist
+    ICNTL(48) = 0; // Multithreading with tree parallelism
+    ICNTL(49) = 0; // Compact workarray id%S at end of factorization phase
+    // ICNTL(50-55) do not exist
+    ICNTL(56) = 0; // Detects pseudo-singularities; rank-revealing factorization of root node
+    // ICNTL(57) does not exist
+    ICNTL(58) = 2; // Options for symbolic factorization
+    // ICNTL(59-60) do not exist
+}
+
+void CooMumpsSolver::configureCNTL()
+{
+    // All CNTL values are left at their defaults unless noted below.
+
+    CNTL(1) = -1.0; // Relative threshold for numerical pivoting
+    CNTL(2) = -1.0; // Stopping criterion for iterative refinement
+    CNTL(3) = 0.0; // Threshold for null pivot row detection
+    CNTL(4) = -1.0; // Threshold for static pivoting
+    CNTL(5) = 0.0; // Fixation for null pivots (effective only when null pivot detection is active)
+    // CNTL(6) does not exist
+    CNTL(7) = 0.0; // Dropping parameter precision for BLR compression
+    // CNTL(8-15) do not exist
+}
+
+SparseMatrixCOO<double> CooMumpsSolver::extractUpperTriangle(const SparseMatrixCOO<double>& matrix) const
+{
+    const int full_nnz = matrix.non_zero_size();
+    const int rows     = matrix.rows();
+    const int cols     = matrix.columns();
+
+    int symmetric_nnz = 0;
+
+    for (int i = 0; i < full_nnz; ++i) {
+        if (matrix.row_index(i) <= matrix.col_index(i))
+            symmetric_nnz++;
+    }
+
+    SparseMatrixCOO<double> sym(rows, cols, symmetric_nnz);
+    sym.is_symmetric(true);
+
+    int current = 0;
+
+    for (int i = 0; i < full_nnz; ++i) {
+        const int row    = matrix.row_index(i);
+        const int column = matrix.col_index(i);
+
+        if (row <= column) {
+            sym.row_index(current) = row;
+            sym.col_index(current) = column;
+            sym.value(current)     = matrix.value(i);
+            current++;
+        }
+    }
+
+    return sym;
+}
+
+#endif

tests/CMakeLists.txt

@@ -12,6 +12,7 @@ add_executable(gmgpolar_tests
     LinearAlgebra/Matrix/coo_matrix.cpp
     LinearAlgebra/Matrix/csr_matrix.cpp
     LinearAlgebra/Solvers/csr_lu_solver.cpp
+    LinearAlgebra/Solvers/coo_mumps_solver.cpp
     LinearAlgebra/Solvers/tridiagonal_solver.cpp
     PolarGrid/polargrid.cpp
     Interpolation/prolongation.cpp