Skip to content

feat: Add MUMPS direct sparse solver support#212

Open
afbrbza wants to merge 14 commits intolabmec:developfrom
afbrbza:develop
Open

feat: Add MUMPS direct sparse solver support#212
afbrbza wants to merge 14 commits intolabmec:developfrom
afbrbza:develop

Conversation

@afbrbza
Copy link

@afbrbza afbrbza commented Feb 24, 2026

Summary

This PR integrates the MUMPS (MUltifrontal Massively Parallel sparse direct Solver) library into NeoPZ as a new direct solver option for sparse linear systems.

What was added

New solver and matrix classes

  • TPZMumpsSolver (Solvers/): A new TPZMatrixSolver subclass that wraps the MUMPS C interface (dmumps_c). Supports symmetric positive definite, symmetric indefinite, and non-symmetric systems. Exposes ICNTL parameter control for fine-tuning
    MUMPS behavior.
  • TPZFYsmpMatrixMumps (Matrix/): Non-symmetric sparse matrix (Yale format) backed by MUMPS for factorization and solve.
  • TPZSYsmpMatrixMumps (Matrix/): Symmetric sparse matrix (Yale format) backed by MUMPS. Both classes store COO format arrays (1-based indexing) as required by MUMPS.

Structural matrix integration

  • TPZSSpStructMatrix and TPZSpStructMatrix were updated to optionally use TPZFYsmpMatrixMumps / TPZSYsmpMatrixMumps when MUMPS is enabled.

CMake build system

  • Added cmake/EnableMUMPS.cmake to detect and link MUMPS and its dependencies (BLAS, LAPACK, ScaLAPACK, OpenMP).
  • Cross-platform support: tested on Ubuntu (static and shared libs) and macOS (static and shared libs), including OpenMP path resolution for macOS.
  • MUMPS support is opt-in via the USING_MUMPS CMake option.

Testing

  • Verified working on Ubuntu with both static and shared MUMPS libraries.
  • Verified working on macOS (Apple Silicon) with both static and shared MUMPS libraries.

@orlandini
Copy link
Member

Hello,

This feature is quite useful, specially when using ARM macs.

I didn't have the time to look through everything, but I would like to share a few things that caught my attention:

  1. In cmake/enable_MUMPS.cmake, I think that lines 34-35 should be removed, otherwise OpenMP_ROOT will be set even if the library was not found.
  2. In TPZ[S]YSMPMumps.h (and non-symmetric variant as well), the description on the header file should be reviewed.
  3. If the C00 arrays will indeed be stored (do they really need to? are they modified by MUMPS calls?), I think that calling TPZ[S]YSMPMumps::UpdateCOOFormat (depending on the value of TPZ[S]YSMPMumps::fCOOValid) at the beginning of TPZ[S]YSMPMumps::Decompose and TPZ[S]YSMPMumps::SolveDirect calls would be better than delegating the responsibility to other classes.
  4. Regarding fCOOValid : if the user calls TPZ[S]YSMPMatrix::SetData, it should be set to false, right?
  5. MUMPS has both complex and real-valued versions, I don't think this is being checked at CMake level.
  6. If NeoPZ will only support real-valued MUMPS (as it seems to be the case, given the unit tests), this should be documented at the class level. And, maybe, aborting the execution before any MUMPS operation if a complex type is being used.

What do you think, @philippedevloo , @giavancini , @afbrbza ?

@afbrbza
Copy link
Author

afbrbza commented Feb 27, 2026

Hi, thank you for the thorough review! Here is a point-by-point response:

1. cmake/EnableMUMPS.cmakeOpenMP_ROOT set even if library was not found
Fixed. The target_compile_definitions call that propagated OpenMP_ROOT as a compile macro was unconditional. It is now guarded with if(OpenMP_ROOT), so the macro is only defined when OpenMP was actually located.

2. Header descriptions should be reviewed
The @brief and @note fields in both TPZSYSMPMumps.h and TPZFYsmpMatrixMumps.h (non-symmetric variant) were updated to accurately describe each class and its limitations.

3. UpdateCOOFormat should be called from Decompose/SolveDirect, not delegated to other classes
Done. Both symmetric and non-symmetric variants now call UpdateCOOFormat() lazily at the beginning of Decompose (guarded by if (!fCOOValid)), so the COO arrays are always up-to-date before being passed to MUMPS, regardless of how the matrix was populated. The COO arrays are freed after Decompose and fCOOValid is reset to false.

4. fCOOValid should be set to false when SetData is called
Done. Both SetData overloads in the symmetric and non-symmetric classes set fCOOValid = false after delegating to the parent class. Additionally, the raw-pointer overload of GetData(int64_t* &IA, int64_t* &JA, TVar* &A) was overridden in both (@giavancini's suggestion) TPZSYsmpMatrixMumps and TPZFYsmpMatrixMumps to also set fCOOValid = false. The reason is that this overload hands the caller direct pointers into the internal fIA/fJA/fA arrays, so values can be modified without the class being aware — making the COO arrays stale. The TPZVec<> overload was intentionally not overridden since it returns copies of the data.

5. MUMPS complex/real support not checked at CMake level
Added a check in EnableMUMPS.cmake after find_package(MUMPS). It reads the MUMPS_DOUBLE / MUMPS_d_FOUND variables that MUMPSConfig.cmake exposes and emits a WARNING if double-precision real MUMPS (dmumps) is not available, since that is what NeoPZ requires. A comment was also added noting that complex MUMPS support (zmumps/cmumps) is planned for a future release.

6. Document real-only support; abort at runtime for complex types
Done on both fronts:

  • Both headers now carry a @note stating that only real-valued types are supported and that complex support is planned for a future release.
  • TPZMumpsSolver::Decompose uses if constexpr (is_complex<TVar>) to print a descriptive error message and call DebugStop() before any MUMPS operation is attempted with a complex type.

Additional note — MUMPS and Pardiso can be used together
It is worth mentioning that MUMPS and Pardiso (MKL) are not mutually exclusive in this implementation. Both can be enabled simultaneously (-DUSING_MUMPS=ON -DUSING_MKL=ON), and the user can choose either solver at runtime through the analysis object, which makes it straightforward to compare or switch between them.

Unit test comparing both solvers
To validate the integration and as a regression baseline, a unit test (TestSolverComparison) was added under UnitTest_PZ/. It is conditionally compiled only when both USING_MKL and USING_MUMPS are enabled. The test solves the same 3D Darcy flow problem with both solvers — in symmetric (SPD, using TPZSYsmpMatrix / TPZSYsmpMatrixMumps) and non-symmetric (TPZFYsmpMatrix / TPZFYsmpMatrixMumps) configurations — and asserts that the relative L2 difference between the two solutions is below a tight tolerance (1e-8). We would appreciate a review of this test as part of this PR.

CI fix — MUMPS-specific StrMatrix sources compiled unconditionally
A build failure was also found and fixed: TPZSpStructMatrixMumps.cpp/.h and TPZSSpStructMatrixMumps.cpp/.h were being added to the pz target unconditionally in StrMatrix/CMakeLists.txt, causing a missing dmumps_c.h error on systems without MUMPS. These files are now added inside if(USING_MUMPS), consistent with how Matrix/CMakeLists.txt already handled TPZYSMPMumps and TPZSYSMPMumps.

@afbrbza
Copy link
Author

afbrbza commented Mar 5, 2026

The ZMUMPS (complex<double>) support has been fully enabled. The DebugStop() guard introduced in point 6 has been removed, as zmumps_c was already implemented via

MumpsStrucSelector<complex<double>> → ZMUMPS_STRUC_C.

References to complex support as a planned future release have been removed from headers and comments accordingly.

To validate this, the existing TestSolverComparison unit test was extended with a new case ("Complex symmetric solvers: ZMUMPS vs Pardiso – 3D L2 projection"), guarded by MUMPS_HAVE_COMPLEX16. The real-STATE test cases were additionally wrapped with #if defined(MUMPS_HAVE_DOUBLE) to keep the file compilable when only the complex flavour of MUMPS is built.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants