Skip to content

Encapsulate MUMPS solver class #182

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
julianlitz wants to merge 15 commits into
base: main
Choose a base branch
from
Open

Encapsulate MUMPS solver class #182

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
ff125ae
extract
julianlitz Feb 18, 2026
8bb865c
...
julianlitz Feb 18, 2026
83f0815
Reorder MUMPS solver structure declaration
julianlitz Feb 18, 2026
489375c
Reorganize MUMPS solver structure declaration
julianlitz Feb 18, 2026
ff96285
Implement CooMumpsSolver for sparse matrix solving
julianlitz Feb 19, 2026
f2b8165
Refactor CooMumpsSolver class interface
julianlitz Feb 19, 2026
a47f84b
Add linear algebra sources to CMakeLists
julianlitz Feb 19, 2026
88eb792
Expand COO_MUMPS_Solver
julianlitz Feb 21, 2026
697740a
Formatting
julianlitz Feb 21, 2026
7421a1e
Fix mumps solver test
julianlitz Feb 21, 2026
d93ccaf
Update coo_mumps_solver.cpp
julianlitz Feb 22, 2026
8f88e53
Update smootherGive.cpp
julianlitz Feb 23, 2026
6dcf5d0
Update smootherTake.cpp
julianlitz Feb 23, 2026
7f62372
Update extrapolatedSmootherGive.cpp
julianlitz Feb 23, 2026
fdc909f
Update extrapolatedSmootherTake.cpp
julianlitz Feb 23, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
19 changes: 4 additions & 15 deletions include/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,14 +12,13 @@ class DirectSolver_COO_MUMPS_Give : public DirectSolver
const DensityProfileCoefficients& density_profile_coefficients,
bool DirBC_Interior, int num_omp_threads);

~DirectSolver_COO_MUMPS_Give() override;
// Note: The rhs (right-hand side) vector gets overwritten during the solution process.
void solveInPlace(Vector<double> solution) override;

private:
// Solver matrix and MUMPS solver structure
SparseMatrixCOO<double> solver_matrix_;
DMUMPS_STRUC_C mumps_solver_;
// MUMPS solver structure with the solver matrix initialized in the constructor.
// std::optional is used because CooMumpsSolver cannot be default-constructed.
std::optional<CooMumpsSolver> mumps_solver_;

// clang-format off
const Stencil stencil_interior_ = {
Expand Down Expand Up @@ -54,10 +53,6 @@ class DirectSolver_COO_MUMPS_Give : public DirectSolver
void buildSolverMatrixCircleSection(const int i_r, SparseMatrixCOO<double>& solver_matrix);
void buildSolverMatrixRadialSection(const int i_theta, SparseMatrixCOO<double>& solver_matrix);

// Initializes the MUMPS solver with the specified matrix.
// Converts to 1-based indexing.
void initializeMumpsSolver(DMUMPS_STRUC_C& mumps_solver, SparseMatrixCOO<double>& solver_matrix);

// Adjusts the right-hand side vector for symmetry corrections.
// This modifies the system from
// A * solution = rhs
Expand All @@ -69,12 +64,6 @@ class DirectSolver_COO_MUMPS_Give : public DirectSolver
void applySymmetryShiftInnerBoundary(Vector<double> x) const;
void applySymmetryShiftOuterBoundary(Vector<double> x) const;

// Solves the adjusted system symmetric(matrixA) * solution = rhs using the MUMPS solver.
void solveWithMumps(Vector<double> solution);

// Finalizes the MUMPS solver, releasing any allocated resources.
void finalizeMumpsSolver(DMUMPS_STRUC_C& mumps_solver);

// Returns the total number of non-zero elements in the solver matrix.
int getNonZeroCountSolverMatrix() const;

Expand All @@ -89,4 +78,4 @@ class DirectSolver_COO_MUMPS_Give : public DirectSolver
double detDF, double coeff_beta);
};

#endif
#endif
19 changes: 4 additions & 15 deletions include/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,14 +12,13 @@ class DirectSolver_COO_MUMPS_Take : public DirectSolver
const DensityProfileCoefficients& density_profile_coefficients,
bool DirBC_Interior, int num_omp_threads);

~DirectSolver_COO_MUMPS_Take() override;
// Note: The rhs (right-hand side) vector gets overwritten during the solution process.
void solveInPlace(Vector<double> solution) override;

private:
// Solver matrix and MUMPS solver structure
SparseMatrixCOO<double> solver_matrix_;
DMUMPS_STRUC_C mumps_solver_;
// MUMPS solver structure with the solver matrix initialized in the constructor.
// std::optional is used because CooMumpsSolver cannot be default-constructed.
std::optional<CooMumpsSolver> mumps_solver_;

// clang-format off
const Stencil stencil_interior_ = {
Expand Down Expand Up @@ -54,10 +53,6 @@ class DirectSolver_COO_MUMPS_Take : public DirectSolver
void buildSolverMatrixCircleSection(const int i_r, SparseMatrixCOO<double>& solver_matrix);
void buildSolverMatrixRadialSection(const int i_theta, SparseMatrixCOO<double>& solver_matrix);

// Initializes the MUMPS solver with the specified matrix.
// Converts to 1-based indexing.
void initializeMumpsSolver(DMUMPS_STRUC_C& mumps_solver, SparseMatrixCOO<double>& solver_matrix);

// Adjusts the right-hand side vector for symmetry corrections.
// This modifies the system from
// A * solution = rhs
Expand All @@ -69,12 +64,6 @@ class DirectSolver_COO_MUMPS_Take : public DirectSolver
void applySymmetryShiftInnerBoundary(Vector<double> x) const;
void applySymmetryShiftOuterBoundary(Vector<double> x) const;

// Solves the adjusted system symmetric(matrixA) * solution = rhs using the MUMPS solver.
void solveWithMumps(Vector<double> solution);

// Finalizes the MUMPS solver, releasing any allocated resources.
void finalizeMumpsSolver(DMUMPS_STRUC_C& mumps_solver);

// Returns the total number of non-zero elements in the solver matrix.
int getNonZeroCountSolverMatrix() const;

Expand All @@ -90,4 +79,4 @@ class DirectSolver_COO_MUMPS_Take : public DirectSolver
ConstVector<double>& coeff_beta);
};

#endif
#endif
1 change: 1 addition & 0 deletions include/DirectSolver/directSolver.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@ class Level;
#include "../LinearAlgebra/Matrix/coo_matrix.h"
#include "../LinearAlgebra/Matrix/csr_matrix.h"
#include "../LinearAlgebra/Solvers/csr_lu_solver.h"
#include "../LinearAlgebra/Solvers/coo_mumps_solver.h"
#include "../Stencil/stencil.h"

#ifdef GMGPOLAR_USE_MUMPS
Expand Down
17 changes: 3 additions & 14 deletions include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,9 +61,6 @@ class ExtrapolatedSmootherGive : public ExtrapolatedSmoother
const DensityProfileCoefficients& density_profile_coefficients,
bool DirBC_Interior, int num_omp_threads);

// If MUMPS is enabled, this cleans up the inner boundary solver.
~ExtrapolatedSmootherGive() override;

// Performs one full coupled extrapolated smoothing sweep:
// BC -> WC -> BR -> WR
// Parallel implementation using OpenMP:
Expand All @@ -89,7 +86,9 @@ class ExtrapolatedSmootherGive : public ExtrapolatedSmoother
// When using the in-house solver, the matrix is stored in CSR format.
#ifdef GMGPOLAR_USE_MUMPS
using MatrixType = SparseMatrixCOO<double>;
DMUMPS_STRUC_C inner_boundary_mumps_solver_;
// MUMPS solver structure with the solver matrix initialized in the constructor.
// std::optional is used because CooMumpsSolver cannot be default-constructed.
std::optional<CooMumpsSolver> inner_boundary_mumps_solver_;
#else
using MatrixType = SparseMatrixCSR<double>;
SparseLUSolver<double> inner_boundary_lu_solver_;
Expand Down Expand Up @@ -179,14 +178,4 @@ class ExtrapolatedSmootherGive : public ExtrapolatedSmoother
void solveWhiteCircleSection(Vector<double> x, Vector<double> temp);
void solveBlackRadialSection(Vector<double> x, Vector<double> temp);
void solveWhiteRadialSection(Vector<double> x, Vector<double> temp);

/* ----------------------------------- */
/* Initialize and destroy MUMPS solver */
/* ----------------------------------- */
#ifdef GMGPOLAR_USE_MUMPS
// Initialize sparse MUMPS solver with assembled COO matrix.
void initializeMumpsSolver(DMUMPS_STRUC_C& mumps_solver, SparseMatrixCOO<double>& solver_matrix);
// Release MUMPS internal memory and MPI structures.
void finalizeMumpsSolver(DMUMPS_STRUC_C& mumps_solver);
#endif
};
17 changes: 3 additions & 14 deletions include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/extrapolatedSmootherTake.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,9 +61,6 @@ class ExtrapolatedSmootherTake : public ExtrapolatedSmoother
const DensityProfileCoefficients& density_profile_coefficients,
bool DirBC_Interior, int num_omp_threads);

// If MUMPS is enabled, this cleans up the inner boundary solver.
~ExtrapolatedSmootherTake() override;

// Performs one full coupled extrapolated smoothing sweep:
// BC -> WC -> BR -> WR
// using temp as RHS workspace.
Expand All @@ -87,7 +84,9 @@ class ExtrapolatedSmootherTake : public ExtrapolatedSmoother
// When using the in-house solver, the matrix is stored in CSR format.
#ifdef GMGPOLAR_USE_MUMPS
using MatrixType = SparseMatrixCOO<double>;
DMUMPS_STRUC_C inner_boundary_mumps_solver_;
// MUMPS solver structure with the solver matrix initialized in the constructor.
// std::optional is used because CooMumpsSolver cannot be default-constructed.
std::optional<CooMumpsSolver> inner_boundary_mumps_solver_;
#else
using MatrixType = SparseMatrixCSR<double>;
SparseLUSolver<double> inner_boundary_lu_solver_;
Expand Down Expand Up @@ -176,14 +175,4 @@ class ExtrapolatedSmootherTake : public ExtrapolatedSmoother
void solveWhiteCircleSection(Vector<double> x, Vector<double> temp);
void solveBlackRadialSection(Vector<double> x, Vector<double> temp);
void solveWhiteRadialSection(Vector<double> x, Vector<double> temp);

/* ----------------------------------- */
/* Initialize and destroy MUMPS solver */
/* ----------------------------------- */
#ifdef GMGPOLAR_USE_MUMPS
// Initialize sparse MUMPS solver with assembled COO matrix.
void initializeMumpsSolver(DMUMPS_STRUC_C& mumps_solver, SparseMatrixCOO<double>& solver_matrix);
// Release MUMPS internal memory and MPI structures.
void finalizeMumpsSolver(DMUMPS_STRUC_C& mumps_solver);
#endif
};
1 change: 1 addition & 0 deletions include/ExtrapolatedSmoother/extrapolatedSmoother.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@ class Level;
#include "../LinearAlgebra/Matrix/coo_matrix.h"
#include "../LinearAlgebra/Matrix/csr_matrix.h"
#include "../LinearAlgebra/Solvers/csr_lu_solver.h"
#include "../LinearAlgebra/Solvers/coo_mumps_solver.h"
#include "../Stencil/stencil.h"

#ifdef GMGPOLAR_USE_MUMPS
Expand Down
78 changes: 78 additions & 0 deletions include/LinearAlgebra/Solvers/coo_mumps_solver.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,78 @@
#pragma once

#ifdef GMGPOLAR_USE_MUMPS

#include "dmumps_c.h"
#include "../../LinearAlgebra/Matrix/coo_matrix.h"
#include "../../LinearAlgebra/Vector/vector.h"

/*
* Wraps MUMPS for solving sparse linear systems given in COO format.
* For general matrices, all non-zero entries must be provided.
* For symmetric matrices (is_symmetric = true), only the upper or lower
* triangular entries should be provided, and the matrix is assumed to be
* positive definite. In GMGPolar this holds true since the domain mapping is invertible.
*/
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];
}

/* ----------------------------------- */
/* 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
17 changes: 3 additions & 14 deletions include/Smoother/SmootherGive/smootherGive.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -53,9 +53,6 @@ class SmootherGive : public Smoother
const DensityProfileCoefficients& density_profile_coefficients, bool DirBC_Interior,
int num_omp_threads);

// If MUMPS is enabled, this cleans up the inner boundary solver.
~SmootherGive() override;

// Performs one full coupled smoothing sweep:
// BC -> WC -> BR -> WR
// Parallel implementation using OpenMP:
Expand All @@ -81,7 +78,9 @@ class SmootherGive : public Smoother
// When using the in-house solver, the matrix is stored in CSR format.
#ifdef GMGPOLAR_USE_MUMPS
using MatrixType = SparseMatrixCOO<double>;
DMUMPS_STRUC_C inner_boundary_mumps_solver_;
// MUMPS solver structure with the solver matrix initialized in the constructor.
// std::optional is used because CooMumpsSolver cannot be default-constructed.
std::optional<CooMumpsSolver> inner_boundary_mumps_solver_;
#else
using MatrixType = SparseMatrixCSR<double>;
SparseLUSolver<double> inner_boundary_lu_solver_;
Expand Down Expand Up @@ -172,14 +171,4 @@ class SmootherGive : public Smoother
void solveWhiteCircleSection(Vector<double> x, Vector<double> temp);
void solveBlackRadialSection(Vector<double> x, Vector<double> temp);
void solveWhiteRadialSection(Vector<double> x, Vector<double> temp);

/* ----------------------------------- */
/* Initialize and destroy MUMPS solver */
/* ----------------------------------- */
#ifdef GMGPOLAR_USE_MUMPS
// Initialize sparse MUMPS solver with assembled COO matrix.
void initializeMumpsSolver(DMUMPS_STRUC_C& mumps_solver, SparseMatrixCOO<double>& solver_matrix);
// Release MUMPS internal memory and MPI structures.
void finalizeMumpsSolver(DMUMPS_STRUC_C& mumps_solver);
#endif
};
17 changes: 3 additions & 14 deletions include/Smoother/SmootherTake/smootherTake.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -53,9 +53,6 @@ class SmootherTake : public Smoother
const DensityProfileCoefficients& density_profile_coefficients, bool DirBC_Interior,
int num_omp_threads);

// If MUMPS is enabled, this cleans up the inner boundary solver.
~SmootherTake() override;

// Performs one full coupled smoothing sweep:
// BC -> WC -> BR -> WR
// using temp as RHS workspace.
Expand All @@ -79,7 +76,9 @@ class SmootherTake : public Smoother
// When using the in-house solver, the matrix is stored in CSR format.
#ifdef GMGPOLAR_USE_MUMPS
using MatrixType = SparseMatrixCOO<double>;
DMUMPS_STRUC_C inner_boundary_mumps_solver_;
// MUMPS solver structure with the solver matrix initialized in the constructor.
// std::optional is used because CooMumpsSolver cannot be default-constructed.
std::optional<CooMumpsSolver> inner_boundary_mumps_solver_;
#else
using MatrixType = SparseMatrixCSR<double>;
SparseLUSolver<double> inner_boundary_lu_solver_;
Expand Down Expand Up @@ -168,14 +167,4 @@ class SmootherTake : public Smoother
void solveWhiteCircleSection(Vector<double> x, Vector<double> temp);
void solveBlackRadialSection(Vector<double> x, Vector<double> temp);
void solveWhiteRadialSection(Vector<double> x, Vector<double> temp);

/* ----------------------------------- */
/* Initialize and destroy MUMPS solver */
/* ----------------------------------- */
#ifdef GMGPOLAR_USE_MUMPS
// Initialize sparse MUMPS solver with assembled COO matrix.
void initializeMumpsSolver(DMUMPS_STRUC_C& mumps_solver, SparseMatrixCOO<double>& solver_matrix);
// Release MUMPS internal memory and MPI structures.
void finalizeMumpsSolver(DMUMPS_STRUC_C& mumps_solver);
#endif
};
1 change: 1 addition & 0 deletions include/Smoother/smoother.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@ class Level;
#include "../LinearAlgebra/Matrix/coo_matrix.h"
#include "../LinearAlgebra/Matrix/csr_matrix.h"
#include "../LinearAlgebra/Solvers/csr_lu_solver.h"
#include "../LinearAlgebra/Solvers/coo_mumps_solver.h"
#include "../Stencil/stencil.h"

#ifdef GMGPOLAR_USE_MUMPS
Expand Down
Loading