From ff125ae3dfeecf71683f9d8ad7fb6f960a1df349 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Wed, 18 Feb 2026 21:35:49 +0100 Subject: [PATCH 01/15] extract --- .../directSolverGive.h | 16 +- .../directSolverTake.h | 16 +- include/DirectSolver/directSolver.h | 1 + .../LinearAlgebra/Solvers/coo_mumps_solver.h | 224 ++++++++++++++++++ src/CMakeLists.txt | 10 +- .../applySymmetryShift.cpp | 17 +- .../buildSolverMatrix.cpp | 2 +- .../directSolverGive.cpp | 10 +- .../initializeMumps.cpp | 125 ---------- .../applySymmetryShift.cpp | 17 +- .../directSolverTake.cpp | 10 +- .../initializeMumps.cpp | 125 ---------- 12 files changed, 246 insertions(+), 327 deletions(-) create mode 100644 include/LinearAlgebra/Solvers/coo_mumps_solver.h delete mode 100644 src/DirectSolver/DirectSolver-COO-MUMPS-Give/initializeMumps.cpp delete mode 100644 src/DirectSolver/DirectSolver-COO-MUMPS-Take/initializeMumps.cpp diff --git a/include/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.h b/include/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.h index 4417e578..8bd3e86f 100644 --- a/include/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.h +++ b/include/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.h @@ -12,14 +12,12 @@ 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 solution) override; private: - // Solver matrix and MUMPS solver structure - SparseMatrixCOO solver_matrix_; - DMUMPS_STRUC_C mumps_solver_; + // MUMPS solver structure with the solver matrix initialized in the constructor. + CooMumpsSolver mumps_solver_; // clang-format off const Stencil stencil_interior_ = { @@ -54,10 +52,6 @@ class DirectSolver_COO_MUMPS_Give : public DirectSolver void buildSolverMatrixCircleSection(const int i_r, SparseMatrixCOO& solver_matrix); void buildSolverMatrixRadialSection(const int i_theta, SparseMatrixCOO& solver_matrix); - // Initializes the MUMPS solver with the specified matrix. - // Converts to 1-based indexing. - void initializeMumpsSolver(DMUMPS_STRUC_C& mumps_solver, SparseMatrixCOO& solver_matrix); - // Adjusts the right-hand side vector for symmetry corrections. // This modifies the system from // A * solution = rhs @@ -69,12 +63,6 @@ class DirectSolver_COO_MUMPS_Give : public DirectSolver void applySymmetryShiftInnerBoundary(Vector x) const; void applySymmetryShiftOuterBoundary(Vector x) const; - // Solves the adjusted system symmetric(matrixA) * solution = rhs using the MUMPS solver. - void solveWithMumps(Vector 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; diff --git a/include/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.h b/include/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.h index 55ae6067..1d494bdf 100644 --- a/include/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.h +++ b/include/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.h @@ -12,14 +12,12 @@ 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 solution) override; private: - // Solver matrix and MUMPS solver structure - SparseMatrixCOO solver_matrix_; - DMUMPS_STRUC_C mumps_solver_; + // MUMPS solver structure with the solver matrix initialized in the constructor. + CooMumpsSolver mumps_solver_; // clang-format off const Stencil stencil_interior_ = { @@ -54,10 +52,6 @@ class DirectSolver_COO_MUMPS_Take : public DirectSolver void buildSolverMatrixCircleSection(const int i_r, SparseMatrixCOO& solver_matrix); void buildSolverMatrixRadialSection(const int i_theta, SparseMatrixCOO& solver_matrix); - // Initializes the MUMPS solver with the specified matrix. - // Converts to 1-based indexing. - void initializeMumpsSolver(DMUMPS_STRUC_C& mumps_solver, SparseMatrixCOO& solver_matrix); - // Adjusts the right-hand side vector for symmetry corrections. // This modifies the system from // A * solution = rhs @@ -69,12 +63,6 @@ class DirectSolver_COO_MUMPS_Take : public DirectSolver void applySymmetryShiftInnerBoundary(Vector x) const; void applySymmetryShiftOuterBoundary(Vector x) const; - // Solves the adjusted system symmetric(matrixA) * solution = rhs using the MUMPS solver. - void solveWithMumps(Vector 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; diff --git a/include/DirectSolver/directSolver.h b/include/DirectSolver/directSolver.h index 252e03d3..f0edd0cf 100644 --- a/include/DirectSolver/directSolver.h +++ b/include/DirectSolver/directSolver.h @@ -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 diff --git a/include/LinearAlgebra/Solvers/coo_mumps_solver.h b/include/LinearAlgebra/Solvers/coo_mumps_solver.h new file mode 100644 index 00000000..12bd356e --- /dev/null +++ b/include/LinearAlgebra/Solvers/coo_mumps_solver.h @@ -0,0 +1,224 @@ +#pragma once + +#ifdef GMGPOLAR_USE_MUMPS + + #include + #include + + #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 matrix) + : matrix_(std::move(matrix)) + { + initialize(); + } + + ~CooMumpsSolver() + { + finalize(); + } + + // rhs is overwritten in-place with the solution on return. + void solve(Vector& rhs) + { + assert(std::ssize(rhs) == mumps_solver_.n); + + mumps_solver_.job = JOB_COMPUTE_SOLUTION; + mumps_solver_.nrhs = 1; + mumps_solver_.lrhs = mumps_solver_.n; // leading dimension: must equal n for dense centralized RHS + mumps_solver_.rhs = rhs.data(); // in: RHS, out: solution (overwritten in-place) + + dmumps_c(&mumps_solver_); + + if (INFOG(1) != 0) { + throw std::runtime_error("MUMPS reported an error during solution phase " + "(INFOG(1) = " + + std::to_string(INFOG(1)) + ")."); + } + } + +private: + void initialize() + { + assert(matrix_.rows() == matrix_.columns()); + + /* + * MUMPS uses 1-based indexing; our COO matrix uses 0-based indexing. + * Adjust row and column indices to match MUMPS' requirements. + */ + 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) { + throw std::runtime_error("MUMPS reported an error during analysis/factorization " + "(INFOG(1) = " + + std::to_string(INFOG(1)) + ")."); + } + + if (mumps_solver_.sym == SYM_POSITIVE_DEFINITE && INFOG(12) != 0) { + throw std::runtime_error("Matrix declared positive definite, " + "but negative pivots were encountered during factorization " + "(INFOG(12) = " + + std::to_string(INFOG(12)) + ")."); + } + } + + void finalize() + { + mumps_solver_.job = JOB_END; + dmumps_c(&mumps_solver_); + } + + void 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 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 + } + +private: + SparseMatrixCOO matrix_; + DMUMPS_STRUC_C mumps_solver_ = {}; + + /* ------------------------------------------------ */ + /* 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; +}; + +#endif // GMGPOLAR_USE_MUMPS \ No newline at end of file diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 0b00e56f..c1b9967c 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -66,26 +66,24 @@ set(DIRECT_SOLVER_SOURCES # Main DirectSolver files ${CMAKE_CURRENT_SOURCE_DIR}/DirectSolver/directSolver.cpp - # DirectSolverGive + # DirectSolver-COO-MUMPS-Give ${CMAKE_CURRENT_SOURCE_DIR}/DirectSolver/DirectSolver-COO-MUMPS-Give/applySymmetryShift.cpp ${CMAKE_CURRENT_SOURCE_DIR}/DirectSolver/DirectSolver-COO-MUMPS-Give/buildSolverMatrix.cpp ${CMAKE_CURRENT_SOURCE_DIR}/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/DirectSolver/DirectSolver-COO-MUMPS-Give/initializeMumps.cpp ${CMAKE_CURRENT_SOURCE_DIR}/DirectSolver/DirectSolver-COO-MUMPS-Give/matrixStencil.cpp - # DirectSolverGiveCustomLU + # DirectSolver-CSR-LU-Give ${CMAKE_CURRENT_SOURCE_DIR}/DirectSolver/DirectSolver-CSR-LU-Give/buildSolverMatrix.cpp ${CMAKE_CURRENT_SOURCE_DIR}/DirectSolver/DirectSolver-CSR-LU-Give/directSolverGive.cpp ${CMAKE_CURRENT_SOURCE_DIR}/DirectSolver/DirectSolver-CSR-LU-Give/matrixStencil.cpp - # DirectSolverTake + # DirectSolver-COO-MUMPS-Take ${CMAKE_CURRENT_SOURCE_DIR}/DirectSolver/DirectSolver-COO-MUMPS-Take/applySymmetryShift.cpp ${CMAKE_CURRENT_SOURCE_DIR}/DirectSolver/DirectSolver-COO-MUMPS-Take/buildSolverMatrix.cpp ${CMAKE_CURRENT_SOURCE_DIR}/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/DirectSolver/DirectSolver-COO-MUMPS-Take/initializeMumps.cpp ${CMAKE_CURRENT_SOURCE_DIR}/DirectSolver/DirectSolver-COO-MUMPS-Take/matrixStencil.cpp - # DirectSolverTakeCustomLU + # DirectSolver-CSR-LU-Take ${CMAKE_CURRENT_SOURCE_DIR}/DirectSolver/DirectSolver-CSR-LU-Take/buildSolverMatrix.cpp ${CMAKE_CURRENT_SOURCE_DIR}/DirectSolver/DirectSolver-CSR-LU-Take/directSolverTake.cpp ${CMAKE_CURRENT_SOURCE_DIR}/DirectSolver/DirectSolver-CSR-LU-Take/matrixStencil.cpp diff --git a/src/DirectSolver/DirectSolver-COO-MUMPS-Give/applySymmetryShift.cpp b/src/DirectSolver/DirectSolver-COO-MUMPS-Give/applySymmetryShift.cpp index b1dc034f..c91a8129 100644 --- a/src/DirectSolver/DirectSolver-COO-MUMPS-Give/applySymmetryShift.cpp +++ b/src/DirectSolver/DirectSolver-COO-MUMPS-Give/applySymmetryShift.cpp @@ -125,20 +125,11 @@ void DirectSolver_COO_MUMPS_Give::applySymmetryShift(Vector x) const assert(std::ssize(x) == grid_.numberOfNodes()); assert(grid_.nr() >= 4); - #pragma omp parallel sections num_threads(num_omp_threads_) - { - #pragma omp section - { - if (DirBC_Interior_) { - applySymmetryShiftInnerBoundary(x); - } - } - - #pragma omp section - { - applySymmetryShiftOuterBoundary(x); - } + if (DirBC_Interior_) { + applySymmetryShiftInnerBoundary(x); } + + applySymmetryShiftOuterBoundary(x); } #endif diff --git a/src/DirectSolver/DirectSolver-COO-MUMPS-Give/buildSolverMatrix.cpp b/src/DirectSolver/DirectSolver-COO-MUMPS-Give/buildSolverMatrix.cpp index 7cd0d75a..240804ad 100644 --- a/src/DirectSolver/DirectSolver-COO-MUMPS-Give/buildSolverMatrix.cpp +++ b/src/DirectSolver/DirectSolver-COO-MUMPS-Give/buildSolverMatrix.cpp @@ -920,7 +920,7 @@ SparseMatrixCOO DirectSolver_COO_MUMPS_Give::buildSolverMatrix() if (row <= col) { symmetric_solver_matrix.row_index(current_nz) = row; symmetric_solver_matrix.col_index(current_nz) = col; - symmetric_solver_matrix.value(current_nz) = std::move(solver_matrix.value(nz_index)); + symmetric_solver_matrix.value(current_nz) = solver_matrix.value(nz_index); current_nz++; } } diff --git a/src/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.cpp b/src/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.cpp index 0a3e025c..509ba004 100644 --- a/src/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.cpp +++ b/src/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.cpp @@ -7,9 +7,8 @@ DirectSolver_COO_MUMPS_Give::DirectSolver_COO_MUMPS_Give(const PolarGrid& grid, const DensityProfileCoefficients& density_profile_coefficients, bool DirBC_Interior, int num_omp_threads) : DirectSolver(grid, level_cache, domain_geometry, density_profile_coefficients, DirBC_Interior, num_omp_threads) + , mumps_solver_(buildSolverMatrix()) { - solver_matrix_ = buildSolverMatrix(); - initializeMumpsSolver(mumps_solver_, solver_matrix_); } void DirectSolver_COO_MUMPS_Give::solveInPlace(Vector solution) @@ -21,12 +20,7 @@ void DirectSolver_COO_MUMPS_Give::solveInPlace(Vector solution) // ensuring that the solution at the boundary is correctly adjusted and maintains the required symmetry. applySymmetryShift(solution); // Solves the adjusted system symmetric(matrixA) * solution = rhs using the MUMPS solver. - solveWithMumps(solution); -} - -DirectSolver_COO_MUMPS_Give::~DirectSolver_COO_MUMPS_Give() -{ - finalizeMumpsSolver(mumps_solver_); + mumps_solver_.solve(solution); } #endif diff --git a/src/DirectSolver/DirectSolver-COO-MUMPS-Give/initializeMumps.cpp b/src/DirectSolver/DirectSolver-COO-MUMPS-Give/initializeMumps.cpp deleted file mode 100644 index 2ebf5530..00000000 --- a/src/DirectSolver/DirectSolver-COO-MUMPS-Give/initializeMumps.cpp +++ /dev/null @@ -1,125 +0,0 @@ -#include "../../../include/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.h" - -#ifdef GMGPOLAR_USE_MUMPS - -void DirectSolver_COO_MUMPS_Give::initializeMumpsSolver(DMUMPS_STRUC_C& mumps_solver, - SparseMatrixCOO& solver_matrix) -{ - /* - * MUMPS (a parallel direct solver) uses 1-based indexing, - * whereas the input matrix follows 0-based indexing. - * Adjust row and column indices to match MUMPS' requirements. - */ - for (int i = 0; i < solver_matrix.non_zero_size(); i++) { - solver_matrix.row_index(i) += 1; - solver_matrix.col_index(i) += 1; - } - - mumps_solver.job = JOB_INIT; - mumps_solver.par = PAR_PARALLEL; - /* The matrix is positive definite for invertible mappings. */ - /* Therefore we use SYM_POSITIVE_DEFINITE instead of SYM_GENERAL_SYMMETRIC. */ - mumps_solver.sym = (solver_matrix.is_symmetric() ? SYM_POSITIVE_DEFINITE : SYM_UNSYMMETRIC); - mumps_solver.comm_fortran = USE_COMM_WORLD; - dmumps_c(&mumps_solver); - - ICNTL(mumps_solver, 1) = 0; // Output stream for error messages. - ICNTL(mumps_solver, 2) = 0; // Output stream for diagnostic printing and statistics local to each MPI process. - ICNTL(mumps_solver, 3) = 0; // Output stream for global information, collected on the host - ICNTL(mumps_solver, 4) = 0; // Level of printing for error, warning, and diagnostic messages. - ICNTL(mumps_solver, 5) = 0; // Controls the matrix input format - ICNTL(mumps_solver, 6) = 7; // Permutes the matrix to a zero-free diagonal and/or scale the matrix - ICNTL(mumps_solver, 7) = - 5; // Computes a symmetric permutation (ordering) to determine the pivot order to be used for the factorization in case of sequential analysis - ICNTL(mumps_solver, 8) = 77; // Describes the scaling strategy - ICNTL(mumps_solver, 9) = 1; // Computes the solution using A or A^T - ICNTL(mumps_solver, 10) = 0; // Applies the iterative refinement to the computed solution - ICNTL(mumps_solver, 11) = 0; // Computes statistics related to an error analysis of the linear system solved - ICNTL(mumps_solver, 12) = 0; // Defines an ordering strategy for symmetric matrices and is used - ICNTL(mumps_solver, 13) = 0; // Controls the parallelism of the root node - ICNTL(mumps_solver, 14) = // Controls the percentage increase in the estimated working space - (solver_matrix.is_symmetric() ? 5 : 20); - ICNTL(mumps_solver, 15) = 0; // Exploits compression of the input matrix resulting from a block format - ICNTL(mumps_solver, 16) = 0; // Controls the setting of the number of OpenMP threads - // ICNTL(17) Doesn't exist - ICNTL(mumps_solver, 18) = 0; // Defines the strategy for the distributed input matrix - ICNTL(mumps_solver, 19) = 0; // Computes the Schur complement matrix - ICNTL(mumps_solver, 20) = 0; // Determines the format (dense, sparse, or distributed) of the right-hand sides - ICNTL(mumps_solver, 21) = 0; // Determines the distribution (centralized or distributed) of the solution vectors. - ICNTL(mumps_solver, 22) = 0; // Controls the in-core/out-of-core (OOC) factorization and solve. - ICNTL(mumps_solver, 23) = 0; // Corresponds to the maximum size of the working memory in MegaBytes that MUMPS can - // allocate per working process - ICNTL(mumps_solver, 24) = 0; // Controls the detection of “null pivot rows”. - ICNTL(mumps_solver, 25) = - 0; // Allows the computation of a solution of a deficient matrix and also of a null space basis - ICNTL(mumps_solver, 26) = 0; // Drives the solution phase if a Schur complement matrix has been computed - ICNTL(mumps_solver, 27) = -32; // Controls the blocking size for multiple right-hand sides. - ICNTL(mumps_solver, 28) = 0; // Determines whether a sequential or parallel computation of the ordering is performed - ICNTL(mumps_solver, 29) = - 0; // Defines the parallel ordering tool (when ICNTL(28)=1) to be used to compute the fill-in reducing permutation. - ICNTL(mumps_solver, 30) = 0; // Computes a user-specified set of entries in the inverse A^−1 of the original matrix - ICNTL(mumps_solver, 31) = 0; // Indicates which factors may be discarded during the factorization. - ICNTL(mumps_solver, 32) = 0; // Performs the forward elimination of the right-hand sides during the factorization - ICNTL(mumps_solver, 33) = 0; // Computes the determinant of the input matrix. - ICNTL(mumps_solver, 34) = 0; // Controls the conservation of the OOC files during JOB= –3 - ICNTL(mumps_solver, 35) = 0; // Controls the activation of the BLR feature - ICNTL(mumps_solver, 36) = 0; // Controls the choice of BLR factorization variant - ICNTL(mumps_solver, 37) = 0; // Controls the BLR compression of the contribution blocks - ICNTL(mumps_solver, 38) = 600; // Estimates compression rate of LU factors - ICNTL(mumps_solver, 39) = 500; // Estimates compression rate of contribution blocks - // ICNTL(40-47) Don't exist - ICNTL(mumps_solver, 48) = 0; // Multithreading with tree parallelism - ICNTL(mumps_solver, 49) = 0; // Compact workarray id%S at the end of factorization phase - // ICNTL(50-55) Don't exist - ICNTL(mumps_solver, 56) = - 0; // Detects pseudo-singularities during factorization and factorizes the root node with a rankrevealing method - // ICNTL(57) Doesn't exist - ICNTL(mumps_solver, 58) = 2; // Defines options for symbolic factorization - // ICNTL(59-60) Don't exist - - CNTL(mumps_solver, 1) = -1.0; // Relative threshold for numerical pivoting - CNTL(mumps_solver, 2) = -1.0; // Stopping criterion for iterative refinement - CNTL(mumps_solver, 3) = 0.0; // Determine null pivot rows - CNTL(mumps_solver, 4) = -1.0; // Determines the threshold for static pivoting - CNTL(mumps_solver, 5) = - 0.0; // Defines the fixation for null pivots and is effective only when null pivot row detection is active - // CNTL(6) Doesn't exist - CNTL(mumps_solver, 7) = 0.0; // Defines the precision of the dropping parameter used during BLR compression - // CNTL(8-15) Don't exist - - mumps_solver.job = JOB_ANALYSIS_AND_FACTORIZATION; - assert(solver_matrix.rows() == solver_matrix.columns()); - mumps_solver.n = solver_matrix.rows(); - mumps_solver.nz = solver_matrix.non_zero_size(); - mumps_solver.irn = solver_matrix.row_indices_data(); - mumps_solver.jcn = solver_matrix.column_indices_data(); - mumps_solver.a = solver_matrix.values_data(); - dmumps_c(&mumps_solver); - - if (mumps_solver.sym == SYM_POSITIVE_DEFINITE && INFOG(mumps_solver, 12) != 0) { - std::cout - << "Warning: DirectSolver matrix is not positive definite: Negative pivots in the factorization phase." - << std::endl; - } -} - -void DirectSolver_COO_MUMPS_Give::solveWithMumps(Vector result_rhs) -{ - mumps_solver_.job = JOB_COMPUTE_SOLUTION; - mumps_solver_.nrhs = 1; - mumps_solver_.nz_rhs = result_rhs.size(); - mumps_solver_.rhs = result_rhs.data(); - mumps_solver_.lrhs = result_rhs.size(); - dmumps_c(&mumps_solver_); - if (mumps_solver_.info[0] != 0) { - std::cerr << "Error solving the direct system: " << mumps_solver_.info[0] << std::endl; - } -} - -void DirectSolver_COO_MUMPS_Give::finalizeMumpsSolver(DMUMPS_STRUC_C& mumps_solver) -{ - mumps_solver.job = JOB_END; - dmumps_c(&mumps_solver); -} - -#endif diff --git a/src/DirectSolver/DirectSolver-COO-MUMPS-Take/applySymmetryShift.cpp b/src/DirectSolver/DirectSolver-COO-MUMPS-Take/applySymmetryShift.cpp index 2b7d8e6e..c81352e0 100644 --- a/src/DirectSolver/DirectSolver-COO-MUMPS-Take/applySymmetryShift.cpp +++ b/src/DirectSolver/DirectSolver-COO-MUMPS-Take/applySymmetryShift.cpp @@ -83,20 +83,11 @@ void DirectSolver_COO_MUMPS_Take::applySymmetryShift(Vector x) const assert(std::ssize(x) == grid_.numberOfNodes()); assert(grid_.nr() >= 4); - #pragma omp parallel sections num_threads(num_omp_threads_) - { - #pragma omp section - { - if (DirBC_Interior_) { - applySymmetryShiftInnerBoundary(x); - } - } - - #pragma omp section - { - applySymmetryShiftOuterBoundary(x); - } + if (DirBC_Interior_) { + applySymmetryShiftInnerBoundary(x); } + + applySymmetryShiftOuterBoundary(x); } #endif diff --git a/src/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.cpp b/src/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.cpp index acf8a942..3408fbb7 100644 --- a/src/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.cpp +++ b/src/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.cpp @@ -7,9 +7,8 @@ DirectSolver_COO_MUMPS_Take::DirectSolver_COO_MUMPS_Take(const PolarGrid& grid, const DensityProfileCoefficients& density_profile_coefficients, bool DirBC_Interior, int num_omp_threads) : DirectSolver(grid, level_cache, domain_geometry, density_profile_coefficients, DirBC_Interior, num_omp_threads) + , mumps_solver_(buildSolverMatrix()) { - solver_matrix_ = buildSolverMatrix(); - initializeMumpsSolver(mumps_solver_, solver_matrix_); } void DirectSolver_COO_MUMPS_Take::solveInPlace(Vector solution) @@ -21,12 +20,7 @@ void DirectSolver_COO_MUMPS_Take::solveInPlace(Vector solution) // ensuring that the solution at the boundary is correctly adjusted and maintains the required symmetry. applySymmetryShift(solution); // Solves the adjusted system symmetric(matrixA) * solution = rhs using the MUMPS solver. - solveWithMumps(solution); -} - -DirectSolver_COO_MUMPS_Take::~DirectSolver_COO_MUMPS_Take() -{ - finalizeMumpsSolver(mumps_solver_); + mumps_solver_.solve(solution); } #endif diff --git a/src/DirectSolver/DirectSolver-COO-MUMPS-Take/initializeMumps.cpp b/src/DirectSolver/DirectSolver-COO-MUMPS-Take/initializeMumps.cpp deleted file mode 100644 index c19743d7..00000000 --- a/src/DirectSolver/DirectSolver-COO-MUMPS-Take/initializeMumps.cpp +++ /dev/null @@ -1,125 +0,0 @@ -#include "../../../include/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.h" - -#ifdef GMGPOLAR_USE_MUMPS - -void DirectSolver_COO_MUMPS_Take::initializeMumpsSolver(DMUMPS_STRUC_C& mumps_solver, - SparseMatrixCOO& solver_matrix) -{ - /* - * MUMPS (a parallel direct solver) uses 1-based indexing, - * whereas the input matrix follows 0-based indexing. - * Adjust row and column indices to match MUMPS' requirements. - */ - for (int i = 0; i < solver_matrix.non_zero_size(); i++) { - solver_matrix.row_index(i) += 1; - solver_matrix.col_index(i) += 1; - } - - mumps_solver.job = JOB_INIT; - mumps_solver.par = PAR_PARALLEL; - /* The matrix is positive definite for invertible mappings. */ - /* Therefore we use SYM_POSITIVE_DEFINITE instead of SYM_GENERAL_SYMMETRIC. */ - mumps_solver.sym = (solver_matrix.is_symmetric() ? SYM_POSITIVE_DEFINITE : SYM_UNSYMMETRIC); - mumps_solver.comm_fortran = USE_COMM_WORLD; - dmumps_c(&mumps_solver); - - ICNTL(mumps_solver, 1) = 0; // Output stream for error messages. - ICNTL(mumps_solver, 2) = 0; // Output stream for diagnostic printing and statistics local to each MPI process. - ICNTL(mumps_solver, 3) = 0; // Output stream for global information, collected on the host - ICNTL(mumps_solver, 4) = 0; // Level of printing for error, warning, and diagnostic messages. - ICNTL(mumps_solver, 5) = 0; // Controls the matrix input format - ICNTL(mumps_solver, 6) = 7; // Permutes the matrix to a zero-free diagonal and/or scale the matrix - ICNTL(mumps_solver, 7) = - 5; // Computes a symmetric permutation (ordering) to determine the pivot order to be used for the factorization in case of sequential analysis - ICNTL(mumps_solver, 8) = 77; // Describes the scaling strategy - ICNTL(mumps_solver, 9) = 1; // Computes the solution using A or A^T - ICNTL(mumps_solver, 10) = 0; // Applies the iterative refinement to the computed solution - ICNTL(mumps_solver, 11) = 0; // Computes statistics related to an error analysis of the linear system solved - ICNTL(mumps_solver, 12) = 0; // Defines an ordering strategy for symmetric matrices and is used - ICNTL(mumps_solver, 13) = 0; // Controls the parallelism of the root node - ICNTL(mumps_solver, 14) = // Controls the percentage increase in the estimated working space - (solver_matrix.is_symmetric() ? 5 : 20); - ICNTL(mumps_solver, 15) = 0; // Exploits compression of the input matrix resulting from a block format - ICNTL(mumps_solver, 16) = 0; // Controls the setting of the number of OpenMP threads - // ICNTL(17) Doesn't exist - ICNTL(mumps_solver, 18) = 0; // Defines the strategy for the distributed input matrix - ICNTL(mumps_solver, 19) = 0; // Computes the Schur complement matrix - ICNTL(mumps_solver, 20) = 0; // Determines the format (dense, sparse, or distributed) of the right-hand sides - ICNTL(mumps_solver, 21) = 0; // Determines the distribution (centralized or distributed) of the solution vectors. - ICNTL(mumps_solver, 22) = 0; // Controls the in-core/out-of-core (OOC) factorization and solve. - ICNTL(mumps_solver, 23) = 0; // Corresponds to the maximum size of the working memory in MegaBytes that MUMPS can - // allocate per working process - ICNTL(mumps_solver, 24) = 0; // Controls the detection of “null pivot rows”. - ICNTL(mumps_solver, 25) = - 0; // Allows the computation of a solution of a deficient matrix and also of a null space basis - ICNTL(mumps_solver, 26) = 0; // Drives the solution phase if a Schur complement matrix has been computed - ICNTL(mumps_solver, 27) = -32; // Controls the blocking size for multiple right-hand sides. - ICNTL(mumps_solver, 28) = 0; // Determines whether a sequential or parallel computation of the ordering is performed - ICNTL(mumps_solver, 29) = - 0; // Defines the parallel ordering tool (when ICNTL(28)=1) to be used to compute the fill-in reducing permutation. - ICNTL(mumps_solver, 30) = 0; // Computes a user-specified set of entries in the inverse A^−1 of the original matrix - ICNTL(mumps_solver, 31) = 0; // Indicates which factors may be discarded during the factorization. - ICNTL(mumps_solver, 32) = 0; // Performs the forward elimination of the right-hand sides during the factorization - ICNTL(mumps_solver, 33) = 0; // Computes the determinant of the input matrix. - ICNTL(mumps_solver, 34) = 0; // Controls the conservation of the OOC files during JOB= –3 - ICNTL(mumps_solver, 35) = 0; // Controls the activation of the BLR feature - ICNTL(mumps_solver, 36) = 0; // Controls the choice of BLR factorization variant - ICNTL(mumps_solver, 37) = 0; // Controls the BLR compression of the contribution blocks - ICNTL(mumps_solver, 38) = 600; // Estimates compression rate of LU factors - ICNTL(mumps_solver, 39) = 500; // Estimates compression rate of contribution blocks - // ICNTL(40-47) Don't exist - ICNTL(mumps_solver, 48) = 0; // Multithreading with tree parallelism - ICNTL(mumps_solver, 49) = 0; // Compact workarray id%S at the end of factorization phase - // ICNTL(50-55) Don't exist - ICNTL(mumps_solver, 56) = - 0; // Detects pseudo-singularities during factorization and factorizes the root node with a rankrevealing method - // ICNTL(57) Doesn't exist - ICNTL(mumps_solver, 58) = 2; // Defines options for symbolic factorization - // ICNTL(59-60) Don't exist - - CNTL(mumps_solver, 1) = -1.0; // Relative threshold for numerical pivoting - CNTL(mumps_solver, 2) = -1.0; // Stopping criterion for iterative refinement - CNTL(mumps_solver, 3) = 0.0; // Determine null pivot rows - CNTL(mumps_solver, 4) = -1.0; // Determines the threshold for static pivoting - CNTL(mumps_solver, 5) = - 0.0; // Defines the fixation for null pivots and is effective only when null pivot row detection is active - // CNTL(6) Doesn't exist - CNTL(mumps_solver, 7) = 0.0; // Defines the precision of the dropping parameter used during BLR compression - // CNTL(8-15) Don't exist - - mumps_solver.job = JOB_ANALYSIS_AND_FACTORIZATION; - assert(solver_matrix.rows() == solver_matrix.columns()); - mumps_solver.n = solver_matrix.rows(); - mumps_solver.nz = solver_matrix.non_zero_size(); - mumps_solver.irn = solver_matrix.row_indices_data(); - mumps_solver.jcn = solver_matrix.column_indices_data(); - mumps_solver.a = solver_matrix.values_data(); - dmumps_c(&mumps_solver); - - if (mumps_solver.sym == SYM_POSITIVE_DEFINITE && INFOG(mumps_solver, 12) != 0) { - std::cout - << "Warning: DirectSolver matrix is not positive definite: Negative pivots in the factorization phase." - << std::endl; - } -} - -void DirectSolver_COO_MUMPS_Take::solveWithMumps(Vector result_rhs) -{ - mumps_solver_.job = JOB_COMPUTE_SOLUTION; - mumps_solver_.nrhs = 1; - mumps_solver_.nz_rhs = result_rhs.size(); - mumps_solver_.rhs = result_rhs.data(); - mumps_solver_.lrhs = result_rhs.size(); - dmumps_c(&mumps_solver_); - if (mumps_solver_.info[0] != 0) { - std::cerr << "Error solving the direct system: " << mumps_solver_.info[0] << std::endl; - } -} - -void DirectSolver_COO_MUMPS_Take::finalizeMumpsSolver(DMUMPS_STRUC_C& mumps_solver) -{ - mumps_solver.job = JOB_END; - dmumps_c(&mumps_solver); -} - -#endif From 8bb865c5834e4412c2c5c5f4451ba5ca9a70dc01 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Wed, 18 Feb 2026 21:39:21 +0100 Subject: [PATCH 02/15] ... --- .../DirectSolver-COO-MUMPS-Take/buildSolverMatrix.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/DirectSolver/DirectSolver-COO-MUMPS-Take/buildSolverMatrix.cpp b/src/DirectSolver/DirectSolver-COO-MUMPS-Take/buildSolverMatrix.cpp index 42566e19..f0886fc1 100644 --- a/src/DirectSolver/DirectSolver-COO-MUMPS-Take/buildSolverMatrix.cpp +++ b/src/DirectSolver/DirectSolver-COO-MUMPS-Take/buildSolverMatrix.cpp @@ -547,7 +547,7 @@ SparseMatrixCOO DirectSolver_COO_MUMPS_Take::buildSolverMatrix() if (row <= col) { symmetric_solver_matrix.row_index(current_nz) = row; symmetric_solver_matrix.col_index(current_nz) = col; - symmetric_solver_matrix.value(current_nz) = std::move(solver_matrix.value(nz_index)); + symmetric_solver_matrix.value(current_nz) = solver_matrix.value(nz_index); current_nz++; } } From 83f0815e293e05881f316871a384d7d5d12e69ef Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Thu, 19 Feb 2026 00:28:32 +0100 Subject: [PATCH 03/15] Reorder MUMPS solver structure declaration Moved MUMPS solver structure declaration below stencils for proper initialization. --- .../DirectSolver-COO-MUMPS-Give/directSolverGive.h | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/include/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.h b/include/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.h index 8bd3e86f..7d532bd5 100644 --- a/include/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.h +++ b/include/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.h @@ -16,9 +16,6 @@ class DirectSolver_COO_MUMPS_Give : public DirectSolver void solveInPlace(Vector solution) override; private: - // MUMPS solver structure with the solver matrix initialized in the constructor. - CooMumpsSolver mumps_solver_; - // clang-format off const Stencil stencil_interior_ = { 7, 4, 8, @@ -47,6 +44,10 @@ class DirectSolver_COO_MUMPS_Give : public DirectSolver }; // clang-format on + // MUMPS solver structure with the solver matrix initialized in the constructor. + // Defined below stencils to ensure that the solver matrix is built after the stencils are defined. + CooMumpsSolver mumps_solver_; + // Constructs a symmetric solver matrix. SparseMatrixCOO buildSolverMatrix(); void buildSolverMatrixCircleSection(const int i_r, SparseMatrixCOO& solver_matrix); @@ -77,4 +78,4 @@ class DirectSolver_COO_MUMPS_Give : public DirectSolver double detDF, double coeff_beta); }; -#endif \ No newline at end of file +#endif From 489375ce1fe5cefcfbac88592a25a8086e3cb461 Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Thu, 19 Feb 2026 00:28:49 +0100 Subject: [PATCH 04/15] Reorganize MUMPS solver structure declaration Reordered MUMPS solver structure declaration to follow stencils. --- .../DirectSolver-COO-MUMPS-Take/directSolverTake.h | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/include/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.h b/include/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.h index 1d494bdf..0151fe61 100644 --- a/include/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.h +++ b/include/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.h @@ -16,9 +16,6 @@ class DirectSolver_COO_MUMPS_Take : public DirectSolver void solveInPlace(Vector solution) override; private: - // MUMPS solver structure with the solver matrix initialized in the constructor. - CooMumpsSolver mumps_solver_; - // clang-format off const Stencil stencil_interior_ = { 7, 4, 8, @@ -47,6 +44,10 @@ class DirectSolver_COO_MUMPS_Take : public DirectSolver }; // clang-format on + // MUMPS solver structure with the solver matrix initialized in the constructor. + // Defined below stencils to ensure that the solver matrix is built after the stencils are defined. + CooMumpsSolver mumps_solver_; + // Constructs a symmetric solver matrix. SparseMatrixCOO buildSolverMatrix(); void buildSolverMatrixCircleSection(const int i_r, SparseMatrixCOO& solver_matrix); @@ -78,4 +79,4 @@ class DirectSolver_COO_MUMPS_Take : public DirectSolver ConstVector& coeff_beta); }; -#endif \ No newline at end of file +#endif From ff962852fb183c867dc99073d8fd1a925b2f8481 Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Fri, 20 Feb 2026 00:01:26 +0100 Subject: [PATCH 05/15] Implement CooMumpsSolver for sparse matrix solving --- .../Solvers/coo_mumps_solver.cpp | 157 ++++++++++++++++++ 1 file changed, 157 insertions(+) create mode 100644 src/LinearAlgebra/Solvers/coo_mumps_solver.cpp diff --git a/src/LinearAlgebra/Solvers/coo_mumps_solver.cpp b/src/LinearAlgebra/Solvers/coo_mumps_solver.cpp new file mode 100644 index 00000000..b3ac63a3 --- /dev/null +++ b/src/LinearAlgebra/Solvers/coo_mumps_solver.cpp @@ -0,0 +1,157 @@ +#include "../../../include/LinearAlgebra/Solvers/coo_mumps_solver.h" + +#ifdef GMGPOLAR_USE_MUMPS + + #include + #include + #include + +CooMumpsSolver::CooMumpsSolver(SparseMatrixCOO matrix) + : matrix_(std::move(matrix)) +{ + initialize(); +} + +CooMumpsSolver::~CooMumpsSolver() +{ + finalize(); +} + +void CooMumpsSolver::solve(Vector& rhs) +{ + assert(std::ssize(rhs) == mumps_solver_.n); + + mumps_solver_.job = JOB_COMPUTE_SOLUTION; + mumps_solver_.nrhs = 1; + mumps_solver_.lrhs = mumps_solver_.n; // leading dimension: must equal n for dense centralized RHS + mumps_solver_.rhs = rhs.data(); // in: RHS, out: solution (overwritten in-place) + 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; our COO matrix uses 0-based indexing. + * Adjust row and column indices to match MUMPS' requirements. + */ + 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 +} + +#endif // GMGPOLAR_USE_MUMPS From f2b8165b5396a4c5e9e2c658d25d64a4a2cf37bf Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Fri, 20 Feb 2026 00:02:01 +0100 Subject: [PATCH 06/15] Refactor CooMumpsSolver class interface --- .../LinearAlgebra/Solvers/coo_mumps_solver.h | 168 ++---------------- 1 file changed, 11 insertions(+), 157 deletions(-) diff --git a/include/LinearAlgebra/Solvers/coo_mumps_solver.h b/include/LinearAlgebra/Solvers/coo_mumps_solver.h index 12bd356e..2521f8ff 100644 --- a/include/LinearAlgebra/Solvers/coo_mumps_solver.h +++ b/include/LinearAlgebra/Solvers/coo_mumps_solver.h @@ -2,11 +2,7 @@ #ifdef GMGPOLAR_USE_MUMPS - #include - #include - #include "dmumps_c.h" - #include "../../LinearAlgebra/Matrix/coo_matrix.h" #include "../../LinearAlgebra/Vector/vector.h" @@ -20,162 +16,17 @@ class CooMumpsSolver { public: - explicit CooMumpsSolver(SparseMatrixCOO matrix) - : matrix_(std::move(matrix)) - { - initialize(); - } - - ~CooMumpsSolver() - { - finalize(); - } + explicit CooMumpsSolver(SparseMatrixCOO matrix); + ~CooMumpsSolver(); // rhs is overwritten in-place with the solution on return. - void solve(Vector& rhs) - { - assert(std::ssize(rhs) == mumps_solver_.n); - - mumps_solver_.job = JOB_COMPUTE_SOLUTION; - mumps_solver_.nrhs = 1; - mumps_solver_.lrhs = mumps_solver_.n; // leading dimension: must equal n for dense centralized RHS - mumps_solver_.rhs = rhs.data(); // in: RHS, out: solution (overwritten in-place) - - dmumps_c(&mumps_solver_); - - if (INFOG(1) != 0) { - throw std::runtime_error("MUMPS reported an error during solution phase " - "(INFOG(1) = " + - std::to_string(INFOG(1)) + ")."); - } - } - -private: - void initialize() - { - assert(matrix_.rows() == matrix_.columns()); - - /* - * MUMPS uses 1-based indexing; our COO matrix uses 0-based indexing. - * Adjust row and column indices to match MUMPS' requirements. - */ - 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) { - throw std::runtime_error("MUMPS reported an error during analysis/factorization " - "(INFOG(1) = " + - std::to_string(INFOG(1)) + ")."); - } - - if (mumps_solver_.sym == SYM_POSITIVE_DEFINITE && INFOG(12) != 0) { - throw std::runtime_error("Matrix declared positive definite, " - "but negative pivots were encountered during factorization " - "(INFOG(12) = " + - std::to_string(INFOG(12)) + ")."); - } - } - - void finalize() - { - mumps_solver_.job = JOB_END; - dmumps_c(&mumps_solver_); - } - - void 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 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 - } + void solve(Vector& rhs); private: - SparseMatrixCOO matrix_; - DMUMPS_STRUC_C mumps_solver_ = {}; + void initialize(); + void finalize(); + void configureICNTL(); + void configureCNTL(); /* ------------------------------------------------ */ /* MUMPS uses 1-based indexing in the documentation */ @@ -219,6 +70,9 @@ class CooMumpsSolver static constexpr int SYM_UNSYMMETRIC = 0; static constexpr int SYM_POSITIVE_DEFINITE = 1; static constexpr int SYM_GENERAL_SYMMETRIC = 2; + + SparseMatrixCOO matrix_; + DMUMPS_STRUC_C mumps_solver_ = {}; }; -#endif // GMGPOLAR_USE_MUMPS \ No newline at end of file +#endif // GMGPOLAR_USE_MUMPS From a47f84b0fdcbe188b071ae610e24777fab969db0 Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Fri, 20 Feb 2026 00:02:31 +0100 Subject: [PATCH 07/15] Add linear algebra sources to CMakeLists --- src/CMakeLists.txt | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index c1b9967c..87b9c7fb 100644 --- a/src/CMakeLists.txt +++ b/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 @@ -166,6 +173,7 @@ set(CONFIG_PARSER_SOURCES # Create the main library add_library(GMGPolarLib STATIC ${POLAR_GRID_SOURCES} + ${LINEAR_ALGEBRA_SOURCES} ${GMG_POLAR_SOURCES} ${MULTIGRID_METHODS_SOURCES} ${LEVEL_SOURCES} From 88eb7924e9b9db2c2b227a05014ac85d054a44cd Mon Sep 17 00:00:00 2001 From: julianlitz Date: Sun, 22 Feb 2026 00:05:08 +0100 Subject: [PATCH 08/15] Expand COO_MUMPS_Solver --- .../directSolverGive.h | 8 +- .../directSolverTake.h | 8 +- .../extrapolatedSmootherGive.h | 17 +-- .../extrapolatedSmootherTake.h | 17 +-- .../extrapolatedSmoother.h | 1 + include/Smoother/SmootherGive/smootherGive.h | 17 +-- include/Smoother/SmootherTake/smootherTake.h | 17 +-- include/Smoother/smoother.h | 1 + src/CMakeLists.txt | 4 - .../directSolverGive.cpp | 5 +- .../directSolverTake.cpp | 5 +- .../buildAscMatrices.cpp | 3 - .../extrapolatedSmootherGive.cpp | 13 +- .../initializeMumps.cpp | 112 ------------------ .../solveAscSystem.cpp | 28 ++--- .../buildAscMatrices.cpp | 3 - .../extrapolatedSmootherTake.cpp | 13 +- .../initializeMumps.cpp | 112 ------------------ .../solveAscSystem.cpp | 28 ++--- src/Smoother/SmootherGive/buildMatrix.cpp | 3 - src/Smoother/SmootherGive/initializeMumps.cpp | 111 ----------------- src/Smoother/SmootherGive/smootherGive.cpp | 13 +- src/Smoother/SmootherGive/solveAscSystem.cpp | 28 ++--- src/Smoother/SmootherTake/buildMatrix.cpp | 3 - src/Smoother/SmootherTake/initializeMumps.cpp | 111 ----------------- src/Smoother/SmootherTake/smootherTake.cpp | 13 +- src/Smoother/SmootherTake/solveAscSystem.cpp | 28 ++--- tests/CMakeLists.txt | 1 + .../Solvers/coo_mumps_solver.cpp | 111 +++++++++++++++++ 29 files changed, 192 insertions(+), 642 deletions(-) delete mode 100644 src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/initializeMumps.cpp delete mode 100644 src/ExtrapolatedSmoother/ExtrapolatedSmootherTake/initializeMumps.cpp delete mode 100644 src/Smoother/SmootherGive/initializeMumps.cpp delete mode 100644 src/Smoother/SmootherTake/initializeMumps.cpp create mode 100644 tests/LinearAlgebra/Solvers/coo_mumps_solver.cpp diff --git a/include/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.h b/include/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.h index 7d532bd5..4fbeda06 100644 --- a/include/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.h +++ b/include/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.h @@ -16,6 +16,10 @@ class DirectSolver_COO_MUMPS_Give : public DirectSolver void solveInPlace(Vector solution) override; private: + // MUMPS solver structure with the solver matrix initialized in the constructor. + // std::optional is used because CooMumpsSolver cannot be default-constructed. + std::optional mumps_solver_; + // clang-format off const Stencil stencil_interior_ = { 7, 4, 8, @@ -44,10 +48,6 @@ class DirectSolver_COO_MUMPS_Give : public DirectSolver }; // clang-format on - // MUMPS solver structure with the solver matrix initialized in the constructor. - // Defined below stencils to ensure that the solver matrix is built after the stencils are defined. - CooMumpsSolver mumps_solver_; - // Constructs a symmetric solver matrix. SparseMatrixCOO buildSolverMatrix(); void buildSolverMatrixCircleSection(const int i_r, SparseMatrixCOO& solver_matrix); diff --git a/include/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.h b/include/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.h index 0151fe61..fbfa6cd2 100644 --- a/include/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.h +++ b/include/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.h @@ -16,6 +16,10 @@ class DirectSolver_COO_MUMPS_Take : public DirectSolver void solveInPlace(Vector solution) override; private: + // MUMPS solver structure with the solver matrix initialized in the constructor. + // std::optional is used because CooMumpsSolver cannot be default-constructed. + std::optional mumps_solver_; + // clang-format off const Stencil stencil_interior_ = { 7, 4, 8, @@ -44,10 +48,6 @@ class DirectSolver_COO_MUMPS_Take : public DirectSolver }; // clang-format on - // MUMPS solver structure with the solver matrix initialized in the constructor. - // Defined below stencils to ensure that the solver matrix is built after the stencils are defined. - CooMumpsSolver mumps_solver_; - // Constructs a symmetric solver matrix. SparseMatrixCOO buildSolverMatrix(); void buildSolverMatrixCircleSection(const int i_r, SparseMatrixCOO& solver_matrix); diff --git a/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.h b/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.h index 9d664335..6d753945 100644 --- a/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.h +++ b/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.h @@ -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: @@ -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; - 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 inner_boundary_mumps_solver_; #else using MatrixType = SparseMatrixCSR; SparseLUSolver inner_boundary_lu_solver_; @@ -179,14 +178,4 @@ class ExtrapolatedSmootherGive : public ExtrapolatedSmoother void solveWhiteCircleSection(Vector x, Vector temp); void solveBlackRadialSection(Vector x, Vector temp); void solveWhiteRadialSection(Vector x, Vector 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& solver_matrix); - // Release MUMPS internal memory and MPI structures. - void finalizeMumpsSolver(DMUMPS_STRUC_C& mumps_solver); -#endif }; diff --git a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/extrapolatedSmootherTake.h b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/extrapolatedSmootherTake.h index b400e7d4..b94863ec 100644 --- a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/extrapolatedSmootherTake.h +++ b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/extrapolatedSmootherTake.h @@ -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. @@ -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; - 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 inner_boundary_mumps_solver_; #else using MatrixType = SparseMatrixCSR; SparseLUSolver inner_boundary_lu_solver_; @@ -176,14 +175,4 @@ class ExtrapolatedSmootherTake : public ExtrapolatedSmoother void solveWhiteCircleSection(Vector x, Vector temp); void solveBlackRadialSection(Vector x, Vector temp); void solveWhiteRadialSection(Vector x, Vector 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& solver_matrix); - // Release MUMPS internal memory and MPI structures. - void finalizeMumpsSolver(DMUMPS_STRUC_C& mumps_solver); -#endif }; diff --git a/include/ExtrapolatedSmoother/extrapolatedSmoother.h b/include/ExtrapolatedSmoother/extrapolatedSmoother.h index e6c4f5eb..544d1c4c 100644 --- a/include/ExtrapolatedSmoother/extrapolatedSmoother.h +++ b/include/ExtrapolatedSmoother/extrapolatedSmoother.h @@ -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 diff --git a/include/Smoother/SmootherGive/smootherGive.h b/include/Smoother/SmootherGive/smootherGive.h index 0d656c6d..c413ba52 100644 --- a/include/Smoother/SmootherGive/smootherGive.h +++ b/include/Smoother/SmootherGive/smootherGive.h @@ -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: @@ -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; - 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 inner_boundary_mumps_solver_; #else using MatrixType = SparseMatrixCSR; SparseLUSolver inner_boundary_lu_solver_; @@ -172,14 +171,4 @@ class SmootherGive : public Smoother void solveWhiteCircleSection(Vector x, Vector temp); void solveBlackRadialSection(Vector x, Vector temp); void solveWhiteRadialSection(Vector x, Vector 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& solver_matrix); - // Release MUMPS internal memory and MPI structures. - void finalizeMumpsSolver(DMUMPS_STRUC_C& mumps_solver); -#endif }; diff --git a/include/Smoother/SmootherTake/smootherTake.h b/include/Smoother/SmootherTake/smootherTake.h index 20107bdc..1f5c7179 100644 --- a/include/Smoother/SmootherTake/smootherTake.h +++ b/include/Smoother/SmootherTake/smootherTake.h @@ -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. @@ -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; - 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 inner_boundary_mumps_solver_; #else using MatrixType = SparseMatrixCSR; SparseLUSolver inner_boundary_lu_solver_; @@ -168,14 +167,4 @@ class SmootherTake : public Smoother void solveWhiteCircleSection(Vector x, Vector temp); void solveBlackRadialSection(Vector x, Vector temp); void solveWhiteRadialSection(Vector x, Vector 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& solver_matrix); - // Release MUMPS internal memory and MPI structures. - void finalizeMumpsSolver(DMUMPS_STRUC_C& mumps_solver); -#endif }; diff --git a/include/Smoother/smoother.h b/include/Smoother/smoother.h index ba2fb9ed..2d47de74 100644 --- a/include/Smoother/smoother.h +++ b/include/Smoother/smoother.h @@ -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 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 87b9c7fb..429f0a8d 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -124,7 +124,6 @@ set(SMOOTHER_SOURCES # SmootherGive ${CMAKE_CURRENT_SOURCE_DIR}/Smoother/SmootherGive/applyAscOrtho.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Smoother/SmootherGive/buildMatrix.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/Smoother/SmootherGive/initializeMumps.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Smoother/SmootherGive/matrixStencil.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Smoother/SmootherGive/smootherGive.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Smoother/SmootherGive/solveAscSystem.cpp @@ -132,7 +131,6 @@ set(SMOOTHER_SOURCES # SmootherTake ${CMAKE_CURRENT_SOURCE_DIR}/Smoother/SmootherTake/applyAscOrtho.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Smoother/SmootherTake/buildMatrix.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/Smoother/SmootherTake/initializeMumps.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Smoother/SmootherTake/matrixStencil.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Smoother/SmootherTake/smootherTake.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Smoother/SmootherTake/solveAscSystem.cpp @@ -150,7 +148,6 @@ set(EXTRAPOLATED_SMOOTHER_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/ExtrapolatedSmoother/ExtrapolatedSmootherGive/applyAscOrtho.cpp ${CMAKE_CURRENT_SOURCE_DIR}/ExtrapolatedSmoother/ExtrapolatedSmootherGive/buildAscMatrices.cpp ${CMAKE_CURRENT_SOURCE_DIR}/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/ExtrapolatedSmoother/ExtrapolatedSmootherGive/initializeMumps.cpp ${CMAKE_CURRENT_SOURCE_DIR}/ExtrapolatedSmoother/ExtrapolatedSmootherGive/smootherStencil.cpp ${CMAKE_CURRENT_SOURCE_DIR}/ExtrapolatedSmoother/ExtrapolatedSmootherGive/solveAscSystem.cpp @@ -158,7 +155,6 @@ set(EXTRAPOLATED_SMOOTHER_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/ExtrapolatedSmoother/ExtrapolatedSmootherTake/applyAscOrtho.cpp ${CMAKE_CURRENT_SOURCE_DIR}/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildAscMatrices.cpp ${CMAKE_CURRENT_SOURCE_DIR}/ExtrapolatedSmoother/ExtrapolatedSmootherTake/extrapolatedSmootherTake.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/ExtrapolatedSmoother/ExtrapolatedSmootherTake/initializeMumps.cpp ${CMAKE_CURRENT_SOURCE_DIR}/ExtrapolatedSmoother/ExtrapolatedSmootherTake/smootherStencil.cpp ${CMAKE_CURRENT_SOURCE_DIR}/ExtrapolatedSmoother/ExtrapolatedSmootherTake/solveAscSystem.cpp ) diff --git a/src/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.cpp b/src/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.cpp index 509ba004..d4f1c97b 100644 --- a/src/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.cpp +++ b/src/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.cpp @@ -7,8 +7,9 @@ DirectSolver_COO_MUMPS_Give::DirectSolver_COO_MUMPS_Give(const PolarGrid& grid, const DensityProfileCoefficients& density_profile_coefficients, bool DirBC_Interior, int num_omp_threads) : DirectSolver(grid, level_cache, domain_geometry, density_profile_coefficients, DirBC_Interior, num_omp_threads) - , mumps_solver_(buildSolverMatrix()) { + SparseMatrixCOO solver_matrix = buildSolverMatrix(); + mumps_solver_.emplace(std::move(solver_matrix)); } void DirectSolver_COO_MUMPS_Give::solveInPlace(Vector solution) @@ -20,7 +21,7 @@ void DirectSolver_COO_MUMPS_Give::solveInPlace(Vector solution) // ensuring that the solution at the boundary is correctly adjusted and maintains the required symmetry. applySymmetryShift(solution); // Solves the adjusted system symmetric(matrixA) * solution = rhs using the MUMPS solver. - mumps_solver_.solve(solution); + mumps_solver_->solve(solution); } #endif diff --git a/src/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.cpp b/src/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.cpp index 3408fbb7..71543f0a 100644 --- a/src/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.cpp +++ b/src/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.cpp @@ -7,8 +7,9 @@ DirectSolver_COO_MUMPS_Take::DirectSolver_COO_MUMPS_Take(const PolarGrid& grid, const DensityProfileCoefficients& density_profile_coefficients, bool DirBC_Interior, int num_omp_threads) : DirectSolver(grid, level_cache, domain_geometry, density_profile_coefficients, DirBC_Interior, num_omp_threads) - , mumps_solver_(buildSolverMatrix()) { + SparseMatrixCOO solver_matrix = buildSolverMatrix(); + mumps_solver_.emplace(std::move(solver_matrix)); } void DirectSolver_COO_MUMPS_Take::solveInPlace(Vector solution) @@ -20,7 +21,7 @@ void DirectSolver_COO_MUMPS_Take::solveInPlace(Vector solution) // ensuring that the solution at the boundary is correctly adjusted and maintains the required symmetry. applySymmetryShift(solution); // Solves the adjusted system symmetric(matrixA) * solution = rhs using the MUMPS solver. - mumps_solver_.solve(solution); + mumps_solver_->solve(solution); } #endif diff --git a/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/buildAscMatrices.cpp b/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/buildAscMatrices.cpp index ccd0cf46..983b5ec6 100644 --- a/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/buildAscMatrices.cpp +++ b/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/buildAscMatrices.cpp @@ -1332,9 +1332,6 @@ void ExtrapolatedSmootherGive::buildAscMatrices() } } - circle_tridiagonal_solver_.setup(); - radial_tridiagonal_solver_.setup(); - #ifdef GMGPOLAR_USE_MUMPS /* ------------------------------------------------------------------- */ /* Part 3: Convert inner_boundary_circle_matrix_ to a symmetric matrix */ diff --git a/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.cpp b/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.cpp index c3c50019..1a34f678 100644 --- a/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.cpp +++ b/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.cpp @@ -10,20 +10,17 @@ ExtrapolatedSmootherGive::ExtrapolatedSmootherGive(const PolarGrid& grid, const , radial_tridiagonal_solver_(grid.lengthSmootherRadial(), grid.ntheta(), false) { buildAscMatrices(); + + circle_tridiagonal_solver_.setup(); + radial_tridiagonal_solver_.setup(); + #ifdef GMGPOLAR_USE_MUMPS - initializeMumpsSolver(inner_boundary_mumps_solver_, inner_boundary_circle_matrix_); + inner_boundary_mumps_solver_.emplace(std::move(inner_boundary_circle_matrix_)); #else inner_boundary_lu_solver_ = SparseLUSolver(inner_boundary_circle_matrix_); #endif } -ExtrapolatedSmootherGive::~ExtrapolatedSmootherGive() -{ -#ifdef GMGPOLAR_USE_MUMPS - finalizeMumpsSolver(inner_boundary_mumps_solver_); -#endif -} - // The smoothing solves linear systems of the form: // A_sc * u_sc = f_sc − A_sc^ortho * u_sc^ortho // where: diff --git a/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/initializeMumps.cpp b/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/initializeMumps.cpp deleted file mode 100644 index 69889eaa..00000000 --- a/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/initializeMumps.cpp +++ /dev/null @@ -1,112 +0,0 @@ -#include "../../../include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.h" - -#ifdef GMGPOLAR_USE_MUMPS - -void ExtrapolatedSmootherGive::initializeMumpsSolver(DMUMPS_STRUC_C& mumps_solver, - SparseMatrixCOO& solver_matrix) -{ - /* - * MUMPS (a parallel direct solver) uses 1-based indexing, - * whereas the input matrix follows 0-based indexing. - * Adjust row and column indices to match MUMPS' requirements. - */ - for (int i = 0; i < solver_matrix.non_zero_size(); i++) { - solver_matrix.row_index(i) += 1; - solver_matrix.col_index(i) += 1; - } - - mumps_solver.job = JOB_INIT; - mumps_solver.par = PAR_PARALLEL; - /* The matrix is positive definite for invertible mappings. */ - /* Therefore we use SYM_POSITIVE_DEFINITE instead of SYM_GENERAL_SYMMETRIC. */ - mumps_solver.sym = (solver_matrix.is_symmetric() ? SYM_POSITIVE_DEFINITE : SYM_UNSYMMETRIC); - mumps_solver.comm_fortran = USE_COMM_WORLD; - dmumps_c(&mumps_solver); - - ICNTL(mumps_solver, 1) = 0; // Output stream for error messages. - ICNTL(mumps_solver, 2) = 0; // Output stream for diagnostic printing and statistics local to each MPI process. - ICNTL(mumps_solver, 3) = 0; // Output stream for global information, collected on the host - ICNTL(mumps_solver, 4) = 0; // Level of printing for error, warning, and diagnostic messages. - ICNTL(mumps_solver, 5) = 0; // Controls the matrix input format - ICNTL(mumps_solver, 6) = 7; // Permutes the matrix to a zero-free diagonal and/or scale the matrix - ICNTL(mumps_solver, 7) = - 5; // Computes a symmetric permutation (ordering) to determine the pivot order to be used for the factorization in case of sequential analysis - ICNTL(mumps_solver, 8) = 77; // Describes the scaling strategy - ICNTL(mumps_solver, 9) = 1; // Computes the solution using A or A^T - ICNTL(mumps_solver, 10) = 0; // Applies the iterative refinement to the computed solution - ICNTL(mumps_solver, 11) = 0; // Computes statistics related to an error analysis of the linear system solved - ICNTL(mumps_solver, 12) = 0; // Defines an ordering strategy for symmetric matrices and is used - ICNTL(mumps_solver, 13) = 0; // Controls the parallelism of the root node - ICNTL(mumps_solver, 14) = // Controls the percentage increase in the estimated working space - (solver_matrix.is_symmetric() ? 5 : 20); - ICNTL(mumps_solver, 15) = 0; // Exploits compression of the input matrix resulting from a block format - ICNTL(mumps_solver, 16) = 0; // Controls the setting of the number of OpenMP threads - // ICNTL(17) Doesn't exist - ICNTL(mumps_solver, 18) = 0; // Defines the strategy for the distributed input matrix - ICNTL(mumps_solver, 19) = 0; // Computes the Schur complement matrix - ICNTL(mumps_solver, 20) = 0; // Determines the format (dense, sparse, or distributed) of the right-hand sides - ICNTL(mumps_solver, 21) = 0; // Determines the distribution (centralized or distributed) of the solution vectors. - ICNTL(mumps_solver, 22) = 0; // Controls the in-core/out-of-core (OOC) factorization and solve. - ICNTL(mumps_solver, 23) = 0; // Corresponds to the maximum size of the working memory in MegaBytes that MUMPS can - // allocate per working process - ICNTL(mumps_solver, 24) = 0; // Controls the detection of “null pivot rows”. - ICNTL(mumps_solver, 25) = - 0; // Allows the computation of a solution of a deficient matrix and also of a null space basis - ICNTL(mumps_solver, 26) = 0; // Drives the solution phase if a Schur complement matrix has been computed - ICNTL(mumps_solver, 27) = -32; // Controls the blocking size for multiple right-hand sides. - ICNTL(mumps_solver, 28) = 0; // Determines whether a sequential or parallel computation of the ordering is performed - ICNTL(mumps_solver, 29) = - 0; // Defines the parallel ordering tool (when ICNTL(28)=1) to be used to compute the fill-in reducing permutation. - ICNTL(mumps_solver, 30) = 0; // Computes a user-specified set of entries in the inverse A^−1 of the original matrix - ICNTL(mumps_solver, 31) = 0; // Indicates which factors may be discarded during the factorization. - ICNTL(mumps_solver, 32) = 0; // Performs the forward elimination of the right-hand sides during the factorization - ICNTL(mumps_solver, 33) = 0; // Computes the determinant of the input matrix. - ICNTL(mumps_solver, 34) = 0; // Controls the conservation of the OOC files during JOB= –3 - ICNTL(mumps_solver, 35) = 0; // Controls the activation of the BLR feature - ICNTL(mumps_solver, 36) = 0; // Controls the choice of BLR factorization variant - ICNTL(mumps_solver, 37) = 0; // Controls the BLR compression of the contribution blocks - ICNTL(mumps_solver, 38) = 600; // Estimates compression rate of LU factors - ICNTL(mumps_solver, 39) = 500; // Estimates compression rate of contribution blocks - // ICNTL(40-47) Don't exist - ICNTL(mumps_solver, 48) = 0; // Multithreading with tree parallelism - ICNTL(mumps_solver, 49) = 0; // Compact workarray id%S at the end of factorization phase - // ICNTL(50-55) Don't exist - ICNTL(mumps_solver, 56) = - 0; // Detects pseudo-singularities during factorization and factorizes the root node with a rankrevealing method - // ICNTL(57) Doesn't exist - ICNTL(mumps_solver, 58) = 2; // Defines options for symbolic factorization - // ICNTL(59-60) Don't exist - - CNTL(mumps_solver, 1) = -1.0; // Relative threshold for numerical pivoting - CNTL(mumps_solver, 2) = -1.0; // Stopping criterion for iterative refinement - CNTL(mumps_solver, 3) = 0.0; // Determine null pivot rows - CNTL(mumps_solver, 4) = -1.0; // Determines the threshold for static pivoting - CNTL(mumps_solver, 5) = - 0.0; // Defines the fixation for null pivots and is effective only when null pivot row detection is active - // CNTL(6) Doesn't exist - CNTL(mumps_solver, 7) = 0.0; // Defines the precision of the dropping parameter used during BLR compression - // CNTL(8-15) Don't exist - - mumps_solver.job = JOB_ANALYSIS_AND_FACTORIZATION; - assert(solver_matrix.rows() == solver_matrix.columns()); - mumps_solver.n = solver_matrix.rows(); - mumps_solver.nz = solver_matrix.non_zero_size(); - mumps_solver.irn = solver_matrix.row_indices_data(); - mumps_solver.jcn = solver_matrix.column_indices_data(); - mumps_solver.a = solver_matrix.values_data(); - dmumps_c(&mumps_solver); - - if (mumps_solver.sym == SYM_POSITIVE_DEFINITE && INFOG(mumps_solver, 12) != 0) { - std::cout << "Warning: ExtrapolatedSmoother inner boundary matrix is not positive definite: Negative pivots in " - "the factorization phase." - << std::endl; - } -} - -void ExtrapolatedSmootherGive::finalizeMumpsSolver(DMUMPS_STRUC_C& mumps_solver) -{ - mumps_solver.job = JOB_END; - dmumps_c(&mumps_solver); -} - -#endif diff --git a/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/solveAscSystem.cpp b/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/solveAscSystem.cpp index 68df770e..a438d610 100644 --- a/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/solveAscSystem.cpp +++ b/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/solveAscSystem.cpp @@ -18,18 +18,12 @@ void ExtrapolatedSmootherGive::solveBlackCircleSection(Vector x, Vector< int batch_stride = 2; circle_tridiagonal_solver_.solve_diagonal(circle_section, batch_offset, batch_stride); + Vector inner_boundary = Kokkos::subview(temp, Kokkos::make_pair(0, grid_.ntheta())); + #ifdef GMGPOLAR_USE_MUMPS - inner_boundary_mumps_solver_.job = JOB_COMPUTE_SOLUTION; - inner_boundary_mumps_solver_.nrhs = 1; // single rhs vector - inner_boundary_mumps_solver_.nz_rhs = grid_.ntheta(); // non-zeros in rhs - inner_boundary_mumps_solver_.rhs = circle_section.data(); - inner_boundary_mumps_solver_.lrhs = grid_.ntheta(); // leading dimension of rhs - dmumps_c(&inner_boundary_mumps_solver_); - if (inner_boundary_mumps_solver_.info[0] != 0) { - std::cerr << "Error solving the system: " << inner_boundary_mumps_solver_.info[0] << std::endl; - } + inner_boundary_mumps_solver_->solve(inner_boundary); #else - inner_boundary_lu_solver_.solveInPlace(circle_section.data()); + inner_boundary_lu_solver_.solveInPlace(inner_boundary); #endif } @@ -61,18 +55,12 @@ void ExtrapolatedSmootherGive::solveWhiteCircleSection(Vector x, Vector< int batch_stride = 2; circle_tridiagonal_solver_.solve_diagonal(circle_section, batch_offset, batch_stride); + Vector inner_boundary = Kokkos::subview(temp, Kokkos::make_pair(0, grid_.ntheta())); + #ifdef GMGPOLAR_USE_MUMPS - inner_boundary_mumps_solver_.job = JOB_COMPUTE_SOLUTION; - inner_boundary_mumps_solver_.nrhs = 1; // single rhs vector - inner_boundary_mumps_solver_.nz_rhs = grid_.ntheta(); // non-zeros in rhs - inner_boundary_mumps_solver_.rhs = circle_section.data(); - inner_boundary_mumps_solver_.lrhs = grid_.ntheta(); // leading dimension of rhs - dmumps_c(&inner_boundary_mumps_solver_); - if (inner_boundary_mumps_solver_.info[0] != 0) { - std::cerr << "Error solving the system: " << inner_boundary_mumps_solver_.info[0] << std::endl; - } + inner_boundary_mumps_solver_->solve(inner_boundary); #else - inner_boundary_lu_solver_.solveInPlace(circle_section.data()); + inner_boundary_lu_solver_.solveInPlace(inner_boundary); #endif } diff --git a/src/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildAscMatrices.cpp b/src/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildAscMatrices.cpp index cfb1fa24..7636df8a 100644 --- a/src/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildAscMatrices.cpp +++ b/src/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildAscMatrices.cpp @@ -690,9 +690,6 @@ void ExtrapolatedSmootherTake::buildAscMatrices() } } - circle_tridiagonal_solver_.setup(); - radial_tridiagonal_solver_.setup(); - #ifdef GMGPOLAR_USE_MUMPS /* ------------------------------------------------------------------- */ /* Part 3: Convert inner_boundary_circle_matrix_ to a symmetric matrix */ diff --git a/src/ExtrapolatedSmoother/ExtrapolatedSmootherTake/extrapolatedSmootherTake.cpp b/src/ExtrapolatedSmoother/ExtrapolatedSmootherTake/extrapolatedSmootherTake.cpp index 0c822f5c..525d7007 100644 --- a/src/ExtrapolatedSmoother/ExtrapolatedSmootherTake/extrapolatedSmootherTake.cpp +++ b/src/ExtrapolatedSmoother/ExtrapolatedSmootherTake/extrapolatedSmootherTake.cpp @@ -10,20 +10,17 @@ ExtrapolatedSmootherTake::ExtrapolatedSmootherTake(const PolarGrid& grid, const , radial_tridiagonal_solver_(grid.lengthSmootherRadial(), grid.ntheta(), false) { buildAscMatrices(); + + circle_tridiagonal_solver_.setup(); + radial_tridiagonal_solver_.setup(); + #ifdef GMGPOLAR_USE_MUMPS - initializeMumpsSolver(inner_boundary_mumps_solver_, inner_boundary_circle_matrix_); + inner_boundary_mumps_solver_.emplace(std::move(inner_boundary_circle_matrix_)); #else inner_boundary_lu_solver_ = SparseLUSolver(inner_boundary_circle_matrix_); #endif } -ExtrapolatedSmootherTake::~ExtrapolatedSmootherTake() -{ -#ifdef GMGPOLAR_USE_MUMPS - finalizeMumpsSolver(inner_boundary_mumps_solver_); -#endif -} - // The smoothing solves linear systems of the form: // A_sc * u_sc = f_sc − A_sc^ortho * u_sc^ortho // where: diff --git a/src/ExtrapolatedSmoother/ExtrapolatedSmootherTake/initializeMumps.cpp b/src/ExtrapolatedSmoother/ExtrapolatedSmootherTake/initializeMumps.cpp deleted file mode 100644 index 4ce13919..00000000 --- a/src/ExtrapolatedSmoother/ExtrapolatedSmootherTake/initializeMumps.cpp +++ /dev/null @@ -1,112 +0,0 @@ -#include "../../../include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/extrapolatedSmootherTake.h" - -#ifdef GMGPOLAR_USE_MUMPS - -void ExtrapolatedSmootherTake::initializeMumpsSolver(DMUMPS_STRUC_C& mumps_solver, - SparseMatrixCOO& solver_matrix) -{ - /* - * MUMPS (a parallel direct solver) uses 1-based indexing, - * whereas the input matrix follows 0-based indexing. - * Adjust row and column indices to match MUMPS' requirements. - */ - for (int i = 0; i < solver_matrix.non_zero_size(); i++) { - solver_matrix.row_index(i) += 1; - solver_matrix.col_index(i) += 1; - } - - mumps_solver.job = JOB_INIT; - mumps_solver.par = PAR_PARALLEL; - /* The matrix is positive definite for invertible mappings. */ - /* Therefore we use SYM_POSITIVE_DEFINITE instead of SYM_GENERAL_SYMMETRIC. */ - mumps_solver.sym = (solver_matrix.is_symmetric() ? SYM_POSITIVE_DEFINITE : SYM_UNSYMMETRIC); - mumps_solver.comm_fortran = USE_COMM_WORLD; - dmumps_c(&mumps_solver); - - ICNTL(mumps_solver, 1) = 0; // Output stream for error messages. - ICNTL(mumps_solver, 2) = 0; // Output stream for diagnostic printing and statistics local to each MPI process. - ICNTL(mumps_solver, 3) = 0; // Output stream for global information, collected on the host - ICNTL(mumps_solver, 4) = 0; // Level of printing for error, warning, and diagnostic messages. - ICNTL(mumps_solver, 5) = 0; // Controls the matrix input format - ICNTL(mumps_solver, 6) = 7; // Permutes the matrix to a zero-free diagonal and/or scale the matrix - ICNTL(mumps_solver, 7) = - 5; // Computes a symmetric permutation (ordering) to determine the pivot order to be used for the factorization in case of sequential analysis - ICNTL(mumps_solver, 8) = 77; // Describes the scaling strategy - ICNTL(mumps_solver, 9) = 1; // Computes the solution using A or A^T - ICNTL(mumps_solver, 10) = 0; // Applies the iterative refinement to the computed solution - ICNTL(mumps_solver, 11) = 0; // Computes statistics related to an error analysis of the linear system solved - ICNTL(mumps_solver, 12) = 0; // Defines an ordering strategy for symmetric matrices and is used - ICNTL(mumps_solver, 13) = 0; // Controls the parallelism of the root node - ICNTL(mumps_solver, 14) = // Controls the percentage increase in the estimated working space - (solver_matrix.is_symmetric() ? 5 : 20); - ICNTL(mumps_solver, 15) = 0; // Exploits compression of the input matrix resulting from a block format - ICNTL(mumps_solver, 16) = 0; // Controls the setting of the number of OpenMP threads - // ICNTL(17) Doesn't exist - ICNTL(mumps_solver, 18) = 0; // Defines the strategy for the distributed input matrix - ICNTL(mumps_solver, 19) = 0; // Computes the Schur complement matrix - ICNTL(mumps_solver, 20) = 0; // Determines the format (dense, sparse, or distributed) of the right-hand sides - ICNTL(mumps_solver, 21) = 0; // Determines the distribution (centralized or distributed) of the solution vectors. - ICNTL(mumps_solver, 22) = 0; // Controls the in-core/out-of-core (OOC) factorization and solve. - ICNTL(mumps_solver, 23) = 0; // Corresponds to the maximum size of the working memory in MegaBytes that MUMPS can - // allocate per working process - ICNTL(mumps_solver, 24) = 0; // Controls the detection of “null pivot rows”. - ICNTL(mumps_solver, 25) = - 0; // Allows the computation of a solution of a deficient matrix and also of a null space basis - ICNTL(mumps_solver, 26) = 0; // Drives the solution phase if a Schur complement matrix has been computed - ICNTL(mumps_solver, 27) = -32; // Controls the blocking size for multiple right-hand sides. - ICNTL(mumps_solver, 28) = 0; // Determines whether a sequential or parallel computation of the ordering is performed - ICNTL(mumps_solver, 29) = - 0; // Defines the parallel ordering tool (when ICNTL(28)=1) to be used to compute the fill-in reducing permutation. - ICNTL(mumps_solver, 30) = 0; // Computes a user-specified set of entries in the inverse A^−1 of the original matrix - ICNTL(mumps_solver, 31) = 0; // Indicates which factors may be discarded during the factorization. - ICNTL(mumps_solver, 32) = 0; // Performs the forward elimination of the right-hand sides during the factorization - ICNTL(mumps_solver, 33) = 0; // Computes the determinant of the input matrix. - ICNTL(mumps_solver, 34) = 0; // Controls the conservation of the OOC files during JOB= –3 - ICNTL(mumps_solver, 35) = 0; // Controls the activation of the BLR feature - ICNTL(mumps_solver, 36) = 0; // Controls the choice of BLR factorization variant - ICNTL(mumps_solver, 37) = 0; // Controls the BLR compression of the contribution blocks - ICNTL(mumps_solver, 38) = 600; // Estimates compression rate of LU factors - ICNTL(mumps_solver, 39) = 500; // Estimates compression rate of contribution blocks - // ICNTL(40-47) Don't exist - ICNTL(mumps_solver, 48) = 0; // Multithreading with tree parallelism - ICNTL(mumps_solver, 49) = 0; // Compact workarray id%S at the end of factorization phase - // ICNTL(50-55) Don't exist - ICNTL(mumps_solver, 56) = - 0; // Detects pseudo-singularities during factorization and factorizes the root node with a rankrevealing method - // ICNTL(57) Doesn't exist - ICNTL(mumps_solver, 58) = 2; // Defines options for symbolic factorization - // ICNTL(59-60) Don't exist - - CNTL(mumps_solver, 1) = -1.0; // Relative threshold for numerical pivoting - CNTL(mumps_solver, 2) = -1.0; // Stopping criterion for iterative refinement - CNTL(mumps_solver, 3) = 0.0; // Determine null pivot rows - CNTL(mumps_solver, 4) = -1.0; // Determines the threshold for static pivoting - CNTL(mumps_solver, 5) = - 0.0; // Defines the fixation for null pivots and is effective only when null pivot row detection is active - // CNTL(6) Doesn't exist - CNTL(mumps_solver, 7) = 0.0; // Defines the precision of the dropping parameter used during BLR compression - // CNTL(8-15) Don't exist - - mumps_solver.job = JOB_ANALYSIS_AND_FACTORIZATION; - assert(solver_matrix.rows() == solver_matrix.columns()); - mumps_solver.n = solver_matrix.rows(); - mumps_solver.nz = solver_matrix.non_zero_size(); - mumps_solver.irn = solver_matrix.row_indices_data(); - mumps_solver.jcn = solver_matrix.column_indices_data(); - mumps_solver.a = solver_matrix.values_data(); - dmumps_c(&mumps_solver); - - if (mumps_solver.sym == SYM_POSITIVE_DEFINITE && INFOG(mumps_solver, 12) != 0) { - std::cout << "Warning: ExtrapolatedSmoother inner boundary matrix is not positive definite: Negative pivots in " - "the factorization phase." - << std::endl; - } -} - -void ExtrapolatedSmootherTake::finalizeMumpsSolver(DMUMPS_STRUC_C& mumps_solver) -{ - mumps_solver.job = JOB_END; - dmumps_c(&mumps_solver); -} - -#endif diff --git a/src/ExtrapolatedSmoother/ExtrapolatedSmootherTake/solveAscSystem.cpp b/src/ExtrapolatedSmoother/ExtrapolatedSmootherTake/solveAscSystem.cpp index cd0ed6bd..d8b49742 100644 --- a/src/ExtrapolatedSmoother/ExtrapolatedSmootherTake/solveAscSystem.cpp +++ b/src/ExtrapolatedSmoother/ExtrapolatedSmootherTake/solveAscSystem.cpp @@ -18,18 +18,12 @@ void ExtrapolatedSmootherTake::solveBlackCircleSection(Vector x, Vector< int batch_stride = 2; circle_tridiagonal_solver_.solve_diagonal(circle_section, batch_offset, batch_stride); + Vector inner_boundary = Kokkos::subview(temp, Kokkos::make_pair(0, grid_.ntheta())); + #ifdef GMGPOLAR_USE_MUMPS - inner_boundary_mumps_solver_.job = JOB_COMPUTE_SOLUTION; - inner_boundary_mumps_solver_.nrhs = 1; // single rhs vector - inner_boundary_mumps_solver_.nz_rhs = grid_.ntheta(); // non-zeros in rhs - inner_boundary_mumps_solver_.rhs = circle_section.data(); - inner_boundary_mumps_solver_.lrhs = grid_.ntheta(); // leading dimension of rhs - dmumps_c(&inner_boundary_mumps_solver_); - if (inner_boundary_mumps_solver_.info[0] != 0) { - std::cerr << "Error solving the system: " << inner_boundary_mumps_solver_.info[0] << std::endl; - } + inner_boundary_mumps_solver_->solve(inner_boundary); #else - inner_boundary_lu_solver_.solveInPlace(circle_section.data()); + inner_boundary_lu_solver_.solveInPlace(inner_boundary); #endif } @@ -61,18 +55,12 @@ void ExtrapolatedSmootherTake::solveWhiteCircleSection(Vector x, Vector< int batch_stride = 2; circle_tridiagonal_solver_.solve_diagonal(circle_section, batch_offset, batch_stride); + Vector inner_boundary = Kokkos::subview(temp, Kokkos::make_pair(0, grid_.ntheta())); + #ifdef GMGPOLAR_USE_MUMPS - inner_boundary_mumps_solver_.job = JOB_COMPUTE_SOLUTION; - inner_boundary_mumps_solver_.nrhs = 1; // single rhs vector - inner_boundary_mumps_solver_.nz_rhs = grid_.ntheta(); // non-zeros in rhs - inner_boundary_mumps_solver_.rhs = circle_section.data(); - inner_boundary_mumps_solver_.lrhs = grid_.ntheta(); // leading dimension of rhs - dmumps_c(&inner_boundary_mumps_solver_); - if (inner_boundary_mumps_solver_.info[0] != 0) { - std::cerr << "Error solving the system: " << inner_boundary_mumps_solver_.info[0] << std::endl; - } + inner_boundary_mumps_solver_->solve(inner_boundary); #else - inner_boundary_lu_solver_.solveInPlace(circle_section.data()); + inner_boundary_lu_solver_.solveInPlace(inner_boundary); #endif } diff --git a/src/Smoother/SmootherGive/buildMatrix.cpp b/src/Smoother/SmootherGive/buildMatrix.cpp index d563beaf..8ac7e859 100644 --- a/src/Smoother/SmootherGive/buildMatrix.cpp +++ b/src/Smoother/SmootherGive/buildMatrix.cpp @@ -750,9 +750,6 @@ void SmootherGive::buildAscMatrices() } } - circle_tridiagonal_solver_.setup(); - radial_tridiagonal_solver_.setup(); - #ifdef GMGPOLAR_USE_MUMPS /* ------------------------------------------------------------------ */ /* Part 3: Convert inner_boundary_circle_matrix to a symmetric matrix */ diff --git a/src/Smoother/SmootherGive/initializeMumps.cpp b/src/Smoother/SmootherGive/initializeMumps.cpp deleted file mode 100644 index 7b96b1c4..00000000 --- a/src/Smoother/SmootherGive/initializeMumps.cpp +++ /dev/null @@ -1,111 +0,0 @@ -#include "../../../include/Smoother/SmootherGive/smootherGive.h" - -#ifdef GMGPOLAR_USE_MUMPS - -void SmootherGive::initializeMumpsSolver(DMUMPS_STRUC_C& mumps_solver, SparseMatrixCOO& solver_matrix) -{ - /* - * MUMPS (a parallel direct solver) uses 1-based indexing, - * whereas the input matrix follows 0-based indexing. - * Adjust row and column indices to match MUMPS' requirements. - */ - for (int i = 0; i < solver_matrix.non_zero_size(); i++) { - solver_matrix.row_index(i) += 1; - solver_matrix.col_index(i) += 1; - } - - mumps_solver.job = JOB_INIT; - mumps_solver.par = PAR_PARALLEL; - /* The matrix is positive definite for invertible mappings. */ - /* Therefore we use SYM_POSITIVE_DEFINITE instead of SYM_GENERAL_SYMMETRIC. */ - mumps_solver.sym = (solver_matrix.is_symmetric() ? SYM_POSITIVE_DEFINITE : SYM_UNSYMMETRIC); - mumps_solver.comm_fortran = USE_COMM_WORLD; - dmumps_c(&mumps_solver); - - ICNTL(mumps_solver, 1) = 0; // Output stream for error messages. - ICNTL(mumps_solver, 2) = 0; // Output stream for diagnostic printing and statistics local to each MPI process. - ICNTL(mumps_solver, 3) = 0; // Output stream for global information, collected on the host - ICNTL(mumps_solver, 4) = 0; // Level of printing for error, warning, and diagnostic messages. - ICNTL(mumps_solver, 5) = 0; // Controls the matrix input format - ICNTL(mumps_solver, 6) = 7; // Permutes the matrix to a zero-free diagonal and/or scale the matrix - ICNTL(mumps_solver, 7) = - 5; // Computes a symmetric permutation (ordering) to determine the pivot order to be used for the factorization in case of sequential analysis - ICNTL(mumps_solver, 8) = 77; // Describes the scaling strategy - ICNTL(mumps_solver, 9) = 1; // Computes the solution using A or A^T - ICNTL(mumps_solver, 10) = 0; // Applies the iterative refinement to the computed solution - ICNTL(mumps_solver, 11) = 0; // Computes statistics related to an error analysis of the linear system solved - ICNTL(mumps_solver, 12) = 0; // Defines an ordering strategy for symmetric matrices and is used - ICNTL(mumps_solver, 13) = 0; // Controls the parallelism of the root node - ICNTL(mumps_solver, 14) = // Controls the percentage increase in the estimated working space - (solver_matrix.is_symmetric() ? 5 : 20); - ICNTL(mumps_solver, 15) = 0; // Exploits compression of the input matrix resulting from a block format - ICNTL(mumps_solver, 16) = 0; // Controls the setting of the number of OpenMP threads - // ICNTL(17) Doesn't exist - ICNTL(mumps_solver, 18) = 0; // Defines the strategy for the distributed input matrix - ICNTL(mumps_solver, 19) = 0; // Computes the Schur complement matrix - ICNTL(mumps_solver, 20) = 0; // Determines the format (dense, sparse, or distributed) of the right-hand sides - ICNTL(mumps_solver, 21) = 0; // Determines the distribution (centralized or distributed) of the solution vectors. - ICNTL(mumps_solver, 22) = 0; // Controls the in-core/out-of-core (OOC) factorization and solve. - ICNTL(mumps_solver, 23) = 0; // Corresponds to the maximum size of the working memory in MegaBytes that MUMPS can - // allocate per working process - ICNTL(mumps_solver, 24) = 0; // Controls the detection of “null pivot rows”. - ICNTL(mumps_solver, 25) = - 0; // Allows the computation of a solution of a deficient matrix and also of a null space basis - ICNTL(mumps_solver, 26) = 0; // Drives the solution phase if a Schur complement matrix has been computed - ICNTL(mumps_solver, 27) = -32; // Controls the blocking size for multiple right-hand sides. - ICNTL(mumps_solver, 28) = 0; // Determines whether a sequential or parallel computation of the ordering is performed - ICNTL(mumps_solver, 29) = - 0; // Defines the parallel ordering tool (when ICNTL(28)=1) to be used to compute the fill-in reducing permutation. - ICNTL(mumps_solver, 30) = 0; // Computes a user-specified set of entries in the inverse A^−1 of the original matrix - ICNTL(mumps_solver, 31) = 0; // Indicates which factors may be discarded during the factorization. - ICNTL(mumps_solver, 32) = 0; // Performs the forward elimination of the right-hand sides during the factorization - ICNTL(mumps_solver, 33) = 0; // Computes the determinant of the input matrix. - ICNTL(mumps_solver, 34) = 0; // Controls the conservation of the OOC files during JOB= –3 - ICNTL(mumps_solver, 35) = 0; // Controls the activation of the BLR feature - ICNTL(mumps_solver, 36) = 0; // Controls the choice of BLR factorization variant - ICNTL(mumps_solver, 37) = 0; // Controls the BLR compression of the contribution blocks - ICNTL(mumps_solver, 38) = 600; // Estimates compression rate of LU factors - ICNTL(mumps_solver, 39) = 500; // Estimates compression rate of contribution blocks - // ICNTL(40-47) Don't exist - ICNTL(mumps_solver, 48) = 0; // Multithreading with tree parallelism - ICNTL(mumps_solver, 49) = 0; // Compact workarray id%S at the end of factorization phase - // ICNTL(50-55) Don't exist - ICNTL(mumps_solver, 56) = - 0; // Detects pseudo-singularities during factorization and factorizes the root node with a rankrevealing method - // ICNTL(57) Doesn't exist - ICNTL(mumps_solver, 58) = 2; // Defines options for symbolic factorization - // ICNTL(59-60) Don't exist - - CNTL(mumps_solver, 1) = -1.0; // Relative threshold for numerical pivoting - CNTL(mumps_solver, 2) = -1.0; // Stopping criterion for iterative refinement - CNTL(mumps_solver, 3) = 0.0; // Determine null pivot rows - CNTL(mumps_solver, 4) = -1.0; // Determines the threshold for static pivoting - CNTL(mumps_solver, 5) = - 0.0; // Defines the fixation for null pivots and is effective only when null pivot row detection is active - // CNTL(6) Doesn't exist - CNTL(mumps_solver, 7) = 0.0; // Defines the precision of the dropping parameter used during BLR compression - // CNTL(8-15) Don't exist - - mumps_solver.job = JOB_ANALYSIS_AND_FACTORIZATION; - assert(solver_matrix.rows() == solver_matrix.columns()); - mumps_solver.n = solver_matrix.rows(); - mumps_solver.nz = solver_matrix.non_zero_size(); - mumps_solver.irn = solver_matrix.row_indices_data(); - mumps_solver.jcn = solver_matrix.column_indices_data(); - mumps_solver.a = solver_matrix.values_data(); - dmumps_c(&mumps_solver); - - if (mumps_solver.sym == SYM_POSITIVE_DEFINITE && INFOG(mumps_solver, 12) != 0) { - std::cout << "Warning: Smoother inner boundary matrix is not positive definite: Negative pivots in the " - "factorization phase." - << std::endl; - } -} - -void SmootherGive::finalizeMumpsSolver(DMUMPS_STRUC_C& mumps_solver) -{ - mumps_solver.job = JOB_END; - dmumps_c(&mumps_solver); -} - -#endif diff --git a/src/Smoother/SmootherGive/smootherGive.cpp b/src/Smoother/SmootherGive/smootherGive.cpp index db458ad0..7b3eeb7a 100644 --- a/src/Smoother/SmootherGive/smootherGive.cpp +++ b/src/Smoother/SmootherGive/smootherGive.cpp @@ -8,20 +8,17 @@ SmootherGive::SmootherGive(const PolarGrid& grid, const LevelCache& level_cache, , radial_tridiagonal_solver_(grid.lengthSmootherRadial(), grid.ntheta(), false) { buildAscMatrices(); + + circle_tridiagonal_solver_.setup(); + radial_tridiagonal_solver_.setup(); + #ifdef GMGPOLAR_USE_MUMPS - initializeMumpsSolver(inner_boundary_mumps_solver_, inner_boundary_circle_matrix_); + inner_boundary_mumps_solver_.emplace(std::move(inner_boundary_circle_matrix_)); #else inner_boundary_lu_solver_ = SparseLUSolver(inner_boundary_circle_matrix_); #endif } -SmootherGive::~SmootherGive() -{ -#ifdef GMGPOLAR_USE_MUMPS - finalizeMumpsSolver(inner_boundary_mumps_solver_); -#endif -} - // The smoothing solves linear systems of the form: // A_sc * u_sc = f_sc − A_sc^ortho * u_sc^ortho // where: diff --git a/src/Smoother/SmootherGive/solveAscSystem.cpp b/src/Smoother/SmootherGive/solveAscSystem.cpp index b84b9e80..ef2a779e 100644 --- a/src/Smoother/SmootherGive/solveAscSystem.cpp +++ b/src/Smoother/SmootherGive/solveAscSystem.cpp @@ -13,18 +13,12 @@ void SmootherGive::solveBlackCircleSection(Vector x, Vector temp circle_tridiagonal_solver_.solve(circle_section, batch_offset, batch_stride); if (is_inner_circle_black) { + Vector inner_boundary = Kokkos::subview(temp, Kokkos::make_pair(0, grid_.ntheta())); + #ifdef GMGPOLAR_USE_MUMPS - inner_boundary_mumps_solver_.job = JOB_COMPUTE_SOLUTION; - inner_boundary_mumps_solver_.nrhs = 1; // single rhs vector - inner_boundary_mumps_solver_.nz_rhs = grid_.ntheta(); // non-zeros in rhs - inner_boundary_mumps_solver_.rhs = circle_section.data(); - inner_boundary_mumps_solver_.lrhs = grid_.ntheta(); // leading dimension of rhs - dmumps_c(&inner_boundary_mumps_solver_); - if (inner_boundary_mumps_solver_.info[0] != 0) { - std::cerr << "Error solving the system: " << inner_boundary_mumps_solver_.info[0] << std::endl; - } + inner_boundary_mumps_solver_->solve(inner_boundary); #else - inner_boundary_lu_solver_.solveInPlace(circle_section.data()); + inner_boundary_lu_solver_.solveInPlace(inner_boundary); #endif } @@ -51,18 +45,12 @@ void SmootherGive::solveWhiteCircleSection(Vector x, Vector temp circle_tridiagonal_solver_.solve(circle_section, batch_offset, batch_stride); if (is_inner_circle_white) { + Vector inner_boundary = Kokkos::subview(temp, Kokkos::make_pair(0, grid_.ntheta())); + #ifdef GMGPOLAR_USE_MUMPS - inner_boundary_mumps_solver_.job = JOB_COMPUTE_SOLUTION; - inner_boundary_mumps_solver_.nrhs = 1; // single rhs vector - inner_boundary_mumps_solver_.nz_rhs = grid_.ntheta(); // non-zeros in rhs - inner_boundary_mumps_solver_.rhs = circle_section.data(); - inner_boundary_mumps_solver_.lrhs = grid_.ntheta(); // leading dimension of rhs - dmumps_c(&inner_boundary_mumps_solver_); - if (inner_boundary_mumps_solver_.info[0] != 0) { - std::cerr << "Error solving the system: " << inner_boundary_mumps_solver_.info[0] << std::endl; - } + inner_boundary_mumps_solver_->solve(inner_boundary); #else - inner_boundary_lu_solver_.solveInPlace(circle_section.data()); + inner_boundary_lu_solver_.solveInPlace(inner_boundary); #endif } diff --git a/src/Smoother/SmootherTake/buildMatrix.cpp b/src/Smoother/SmootherTake/buildMatrix.cpp index 4a50bcac..703b939e 100644 --- a/src/Smoother/SmootherTake/buildMatrix.cpp +++ b/src/Smoother/SmootherTake/buildMatrix.cpp @@ -436,9 +436,6 @@ void SmootherTake::buildAscMatrices() } } - circle_tridiagonal_solver_.setup(); - radial_tridiagonal_solver_.setup(); - #ifdef GMGPOLAR_USE_MUMPS /* ------------------------------------------------------------------ */ /* Part 3: Convert inner_boundary_circle_matrix to a symmetric matrix */ diff --git a/src/Smoother/SmootherTake/initializeMumps.cpp b/src/Smoother/SmootherTake/initializeMumps.cpp deleted file mode 100644 index d9bbbdb2..00000000 --- a/src/Smoother/SmootherTake/initializeMumps.cpp +++ /dev/null @@ -1,111 +0,0 @@ -#include "../../../include/Smoother/SmootherTake/smootherTake.h" - -#ifdef GMGPOLAR_USE_MUMPS - -void SmootherTake::initializeMumpsSolver(DMUMPS_STRUC_C& mumps_solver, SparseMatrixCOO& solver_matrix) -{ - /* - * MUMPS (a parallel direct solver) uses 1-based indexing, - * whereas the input matrix follows 0-based indexing. - * Adjust row and column indices to match MUMPS' requirements. - */ - for (int i = 0; i < solver_matrix.non_zero_size(); i++) { - solver_matrix.row_index(i) += 1; - solver_matrix.col_index(i) += 1; - } - - mumps_solver.job = JOB_INIT; - mumps_solver.par = PAR_PARALLEL; - /* The matrix is positive definite for invertible mappings. */ - /* Therefore we use SYM_POSITIVE_DEFINITE instead of SYM_GENERAL_SYMMETRIC. */ - mumps_solver.sym = (solver_matrix.is_symmetric() ? SYM_POSITIVE_DEFINITE : SYM_UNSYMMETRIC); - mumps_solver.comm_fortran = USE_COMM_WORLD; - dmumps_c(&mumps_solver); - - ICNTL(mumps_solver, 1) = 0; // Output stream for error messages. - ICNTL(mumps_solver, 2) = 0; // Output stream for diagnostic printing and statistics local to each MPI process. - ICNTL(mumps_solver, 3) = 0; // Output stream for global information, collected on the host - ICNTL(mumps_solver, 4) = 0; // Level of printing for error, warning, and diagnostic messages. - ICNTL(mumps_solver, 5) = 0; // Controls the matrix input format - ICNTL(mumps_solver, 6) = 7; // Permutes the matrix to a zero-free diagonal and/or scale the matrix - ICNTL(mumps_solver, 7) = - 5; // Computes a symmetric permutation (ordering) to determine the pivot order to be used for the factorization in case of sequential analysis - ICNTL(mumps_solver, 8) = 77; // Describes the scaling strategy - ICNTL(mumps_solver, 9) = 1; // Computes the solution using A or A^T - ICNTL(mumps_solver, 10) = 0; // Applies the iterative refinement to the computed solution - ICNTL(mumps_solver, 11) = 0; // Computes statistics related to an error analysis of the linear system solved - ICNTL(mumps_solver, 12) = 0; // Defines an ordering strategy for symmetric matrices and is used - ICNTL(mumps_solver, 13) = 0; // Controls the parallelism of the root node - ICNTL(mumps_solver, 14) = // Controls the percentage increase in the estimated working space - (solver_matrix.is_symmetric() ? 5 : 20); - ICNTL(mumps_solver, 15) = 0; // Exploits compression of the input matrix resulting from a block format - ICNTL(mumps_solver, 16) = 0; // Controls the setting of the number of OpenMP threads - // ICNTL(17) Doesn't exist - ICNTL(mumps_solver, 18) = 0; // Defines the strategy for the distributed input matrix - ICNTL(mumps_solver, 19) = 0; // Computes the Schur complement matrix - ICNTL(mumps_solver, 20) = 0; // Determines the format (dense, sparse, or distributed) of the right-hand sides - ICNTL(mumps_solver, 21) = 0; // Determines the distribution (centralized or distributed) of the solution vectors. - ICNTL(mumps_solver, 22) = 0; // Controls the in-core/out-of-core (OOC) factorization and solve. - ICNTL(mumps_solver, 23) = 0; // Corresponds to the maximum size of the working memory in MegaBytes that MUMPS can - // allocate per working process - ICNTL(mumps_solver, 24) = 0; // Controls the detection of “null pivot rows”. - ICNTL(mumps_solver, 25) = - 0; // Allows the computation of a solution of a deficient matrix and also of a null space basis - ICNTL(mumps_solver, 26) = 0; // Drives the solution phase if a Schur complement matrix has been computed - ICNTL(mumps_solver, 27) = -32; // Controls the blocking size for multiple right-hand sides. - ICNTL(mumps_solver, 28) = 0; // Determines whether a sequential or parallel computation of the ordering is performed - ICNTL(mumps_solver, 29) = - 0; // Defines the parallel ordering tool (when ICNTL(28)=1) to be used to compute the fill-in reducing permutation. - ICNTL(mumps_solver, 30) = 0; // Computes a user-specified set of entries in the inverse A^−1 of the original matrix - ICNTL(mumps_solver, 31) = 0; // Indicates which factors may be discarded during the factorization. - ICNTL(mumps_solver, 32) = 0; // Performs the forward elimination of the right-hand sides during the factorization - ICNTL(mumps_solver, 33) = 0; // Computes the determinant of the input matrix. - ICNTL(mumps_solver, 34) = 0; // Controls the conservation of the OOC files during JOB= –3 - ICNTL(mumps_solver, 35) = 0; // Controls the activation of the BLR feature - ICNTL(mumps_solver, 36) = 0; // Controls the choice of BLR factorization variant - ICNTL(mumps_solver, 37) = 0; // Controls the BLR compression of the contribution blocks - ICNTL(mumps_solver, 38) = 600; // Estimates compression rate of LU factors - ICNTL(mumps_solver, 39) = 500; // Estimates compression rate of contribution blocks - // ICNTL(40-47) Don't exist - ICNTL(mumps_solver, 48) = 0; // Multithreading with tree parallelism - ICNTL(mumps_solver, 49) = 0; // Compact workarray id%S at the end of factorization phase - // ICNTL(50-55) Don't exist - ICNTL(mumps_solver, 56) = - 0; // Detects pseudo-singularities during factorization and factorizes the root node with a rankrevealing method - // ICNTL(57) Doesn't exist - ICNTL(mumps_solver, 58) = 2; // Defines options for symbolic factorization - // ICNTL(59-60) Don't exist - - CNTL(mumps_solver, 1) = -1.0; // Relative threshold for numerical pivoting - CNTL(mumps_solver, 2) = -1.0; // Stopping criterion for iterative refinement - CNTL(mumps_solver, 3) = 0.0; // Determine null pivot rows - CNTL(mumps_solver, 4) = -1.0; // Determines the threshold for static pivoting - CNTL(mumps_solver, 5) = - 0.0; // Defines the fixation for null pivots and is effective only when null pivot row detection is active - // CNTL(6) Doesn't exist - CNTL(mumps_solver, 7) = 0.0; // Defines the precision of the dropping parameter used during BLR compression - // CNTL(8-15) Don't exist - - mumps_solver.job = JOB_ANALYSIS_AND_FACTORIZATION; - assert(solver_matrix.rows() == solver_matrix.columns()); - mumps_solver.n = solver_matrix.rows(); - mumps_solver.nz = solver_matrix.non_zero_size(); - mumps_solver.irn = solver_matrix.row_indices_data(); - mumps_solver.jcn = solver_matrix.column_indices_data(); - mumps_solver.a = solver_matrix.values_data(); - dmumps_c(&mumps_solver); - - if (mumps_solver.sym == SYM_POSITIVE_DEFINITE && INFOG(mumps_solver, 12) != 0) { - std::cout << "Warning: Smoother inner boundary matrix is not positive definite: Negative pivots in the " - "factorization phase." - << std::endl; - } -} - -void SmootherTake::finalizeMumpsSolver(DMUMPS_STRUC_C& mumps_solver) -{ - mumps_solver.job = JOB_END; - dmumps_c(&mumps_solver); -} - -#endif diff --git a/src/Smoother/SmootherTake/smootherTake.cpp b/src/Smoother/SmootherTake/smootherTake.cpp index 9bdff112..8ba590e6 100644 --- a/src/Smoother/SmootherTake/smootherTake.cpp +++ b/src/Smoother/SmootherTake/smootherTake.cpp @@ -8,20 +8,17 @@ SmootherTake::SmootherTake(const PolarGrid& grid, const LevelCache& level_cache, , radial_tridiagonal_solver_(grid.lengthSmootherRadial(), grid.ntheta(), false) { buildAscMatrices(); + + circle_tridiagonal_solver_.setup(); + radial_tridiagonal_solver_.setup(); + #ifdef GMGPOLAR_USE_MUMPS - initializeMumpsSolver(inner_boundary_mumps_solver_, inner_boundary_circle_matrix_); + inner_boundary_mumps_solver_.emplace(std::move(inner_boundary_circle_matrix_)); #else inner_boundary_lu_solver_ = SparseLUSolver(inner_boundary_circle_matrix_); #endif } -SmootherTake::~SmootherTake() -{ -#ifdef GMGPOLAR_USE_MUMPS - finalizeMumpsSolver(inner_boundary_mumps_solver_); -#endif -} - // The smoothing solves linear systems of the form: // A_sc * u_sc = f_sc − A_sc^ortho * u_sc^ortho // where: diff --git a/src/Smoother/SmootherTake/solveAscSystem.cpp b/src/Smoother/SmootherTake/solveAscSystem.cpp index 9917ac83..64cd249c 100644 --- a/src/Smoother/SmootherTake/solveAscSystem.cpp +++ b/src/Smoother/SmootherTake/solveAscSystem.cpp @@ -13,18 +13,12 @@ void SmootherTake::solveBlackCircleSection(Vector x, Vector temp circle_tridiagonal_solver_.solve(circle_section, batch_offset, batch_stride); if (is_inner_circle_black) { + Vector inner_boundary = Kokkos::subview(temp, Kokkos::make_pair(0, grid_.ntheta())); + #ifdef GMGPOLAR_USE_MUMPS - inner_boundary_mumps_solver_.job = JOB_COMPUTE_SOLUTION; - inner_boundary_mumps_solver_.nrhs = 1; // single rhs vector - inner_boundary_mumps_solver_.nz_rhs = grid_.ntheta(); // non-zeros in rhs - inner_boundary_mumps_solver_.rhs = circle_section.data(); - inner_boundary_mumps_solver_.lrhs = grid_.ntheta(); // leading dimension of rhs - dmumps_c(&inner_boundary_mumps_solver_); - if (inner_boundary_mumps_solver_.info[0] != 0) { - std::cerr << "Error solving the system: " << inner_boundary_mumps_solver_.info[0] << std::endl; - } + inner_boundary_mumps_solver_->solve(inner_boundary); #else - inner_boundary_lu_solver_.solveInPlace(circle_section.data()); + inner_boundary_lu_solver_.solveInPlace(inner_boundary); #endif } @@ -51,18 +45,12 @@ void SmootherTake::solveWhiteCircleSection(Vector x, Vector temp circle_tridiagonal_solver_.solve(circle_section, batch_offset, batch_stride); if (is_inner_circle_white) { + Vector inner_boundary = Kokkos::subview(temp, Kokkos::make_pair(0, grid_.ntheta())); + #ifdef GMGPOLAR_USE_MUMPS - inner_boundary_mumps_solver_.job = JOB_COMPUTE_SOLUTION; - inner_boundary_mumps_solver_.nrhs = 1; // single rhs vector - inner_boundary_mumps_solver_.nz_rhs = grid_.ntheta(); // non-zeros in rhs - inner_boundary_mumps_solver_.rhs = circle_section.data(); - inner_boundary_mumps_solver_.lrhs = grid_.ntheta(); // leading dimension of rhs - dmumps_c(&inner_boundary_mumps_solver_); - if (inner_boundary_mumps_solver_.info[0] != 0) { - std::cerr << "Error solving the system: " << inner_boundary_mumps_solver_.info[0] << std::endl; - } + inner_boundary_mumps_solver_->solve(inner_boundary); #else - inner_boundary_lu_solver_.solveInPlace(circle_section.data()); + inner_boundary_lu_solver_.solveInPlace(inner_boundary); #endif } diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index f5699a4e..f8192ee2 100644 --- a/tests/CMakeLists.txt +++ b/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 diff --git a/tests/LinearAlgebra/Solvers/coo_mumps_solver.cpp b/tests/LinearAlgebra/Solvers/coo_mumps_solver.cpp new file mode 100644 index 00000000..b320fd43 --- /dev/null +++ b/tests/LinearAlgebra/Solvers/coo_mumps_solver.cpp @@ -0,0 +1,111 @@ +#pragma once +#ifdef GMGPOLAR_USE_MUMPS + + #include + #include "coo_mumps_solver.h" + +// ----------------------------------------------------------------------- +// Test 1: General (non-symmetric) 4x4 system +// +// Matrix A (non-symmetric): +// [ 4 1 0 0 ] +// [ 2 5 1 0 ] +// [ 0 1 6 2 ] +// [ 0 0 3 7 ] +// +// RHS b = [5, 8, 9, 10]^T +// +// Expected solution computed from A*x = b. +// ----------------------------------------------------------------------- +TEST(CooMumpsSolverTest, GeneralNonSymmetric4x4) +{ + using triplet = SparseMatrixCOO::triplet_type; + + // All non-zero entries (0-based row/col indices) + std::vector entries = {{0, 0, 4.0}, {0, 1, 1.0}, {1, 0, 2.0}, {1, 1, 5.0}, {1, 2, 1.0}, + {2, 1, 1.0}, {2, 2, 6.0}, {2, 3, 2.0}, {3, 2, 3.0}, {3, 3, 7.0}}; + + SparseMatrixCOO mat(4, 4, entries); + mat.is_symmetric(false); + + CooMumpsSolver solver(std::move(mat)); + + Vector rhs("rhs", 4); + rhs(0) = 5.0; + rhs(1) = 8.0; + rhs(2) = 9.0; + rhs(3) = 10.0; + + solver.solve(rhs); + + // Verify A*x = b by back-substitution check + // Reference solution (computed analytically / via numpy): + // x ~ [0.9526, 0.2105, 0.9298, 0.8319] + const double tol = 1e-10; + EXPECT_NEAR(rhs(0), 4.0 * 0.9526 + 1.0 * 0.2105, 1e-3); + + // More robust: re-multiply and check residual + std::vector x = {rhs(0), rhs(1), rhs(2), rhs(3)}; + + // A*x + double Ax0 = 4 * x[0] + 1 * x[1]; + double Ax1 = 2 * x[0] + 5 * x[1] + 1 * x[2]; + double Ax2 = 1 * x[1] + 6 * x[2] + 2 * x[3]; + double Ax3 = 3 * x[2] + 7 * x[3]; + + EXPECT_NEAR(Ax0, 5.0, tol); + EXPECT_NEAR(Ax1, 8.0, tol); + EXPECT_NEAR(Ax2, 9.0, tol); + EXPECT_NEAR(Ax3, 10.0, tol); +} + +// ----------------------------------------------------------------------- +// Test 2: Symmetric positive-definite 4x4 system (lower triangle only) +// +// Matrix A (SPD): +// [ 4 2 0 0 ] +// [ 2 5 1 0 ] +// [ 0 1 6 2 ] +// [ 0 0 2 7 ] +// +// Only lower triangular entries provided (including diagonal). +// RHS b = [6, 8, 9, 11]^T +// ----------------------------------------------------------------------- +TEST(CooMumpsSolverTest, SymmetricPositiveDefinite4x4) +{ + using triplet = SparseMatrixCOO::triplet_type; + + // Lower triangular entries only (0-based row/col indices) + std::vector entries = {{0, 0, 4.0}, {1, 0, 2.0}, {1, 1, 5.0}, {2, 1, 1.0}, + {2, 2, 6.0}, {3, 2, 2.0}, {3, 3, 7.0}}; + + SparseMatrixCOO mat(4, 4, entries); + mat.is_symmetric(true); // SPD, half-entries only + + CooMumpsSolver solver(std::move(mat)); + + Vector rhs("rhs", 4); + rhs(0) = 6.0; + rhs(1) = 8.0; + rhs(2) = 9.0; + rhs(3) = 11.0; + + solver.solve(rhs); + + // Verify residual: A*x = b using the full (symmetric) matrix + std::vector x = {rhs(0), rhs(1), rhs(2), rhs(3)}; + + const double tol = 1e-10; + + double Ax0 = 4 * x[0] + 2 * x[1]; + double Ax1 = 2 * x[0] + 5 * x[1] + 1 * x[2]; + double Ax2 = 1 * x[1] + 6 * x[2] + 2 * x[3]; + double Ax3 = 2 * x[2] + 7 * x[3]; + + EXPECT_NEAR(Ax0, 6.0, tol); + EXPECT_NEAR(Ax1, 8.0, tol); + EXPECT_NEAR(Ax2, 9.0, tol); + EXPECT_NEAR(Ax3, 11.0, tol); +} + +#endif // GMGPOLAR_USE_MUMPS \ No newline at end of file From 697740a9e63916467ea8065fca658020fe51ff68 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Sun, 22 Feb 2026 00:09:05 +0100 Subject: [PATCH 09/15] Formatting --- src/LinearAlgebra/Solvers/coo_mumps_solver.cpp | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/LinearAlgebra/Solvers/coo_mumps_solver.cpp b/src/LinearAlgebra/Solvers/coo_mumps_solver.cpp index b3ac63a3..75aba2eb 100644 --- a/src/LinearAlgebra/Solvers/coo_mumps_solver.cpp +++ b/src/LinearAlgebra/Solvers/coo_mumps_solver.cpp @@ -28,8 +28,7 @@ void CooMumpsSolver::solve(Vector& rhs) dmumps_c(&mumps_solver_); if (INFOG(1) != 0) { - std::cerr << "MUMPS reported an error during solution phase " - << "(INFOG(1) = " << INFOG(1) << ").\n"; + std::cerr << "MUMPS reported an error during solution phase " << "(INFOG(1) = " << INFOG(1) << ").\n"; } } @@ -64,15 +63,14 @@ void CooMumpsSolver::initialize() dmumps_c(&mumps_solver_); if (INFOG(1) != 0) { - std::cerr << "MUMPS reported an error during analysis/factorization " - << "(INFOG(1) = " << INFOG(1) << ").\n"; + 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"; + << "but negative pivots were encountered during factorization " << "(INFOG(12) = " << INFOG(12) + << ").\n"; } } From 7421a1e88c6aadfb19c90946c571e3b0455debfd Mon Sep 17 00:00:00 2001 From: julianlitz Date: Sun, 22 Feb 2026 00:38:11 +0100 Subject: [PATCH 10/15] Fix mumps solver test --- .../Solvers/coo_mumps_solver.cpp | 98 +++++++++---------- 1 file changed, 46 insertions(+), 52 deletions(-) diff --git a/tests/LinearAlgebra/Solvers/coo_mumps_solver.cpp b/tests/LinearAlgebra/Solvers/coo_mumps_solver.cpp index b320fd43..02cb0a52 100644 --- a/tests/LinearAlgebra/Solvers/coo_mumps_solver.cpp +++ b/tests/LinearAlgebra/Solvers/coo_mumps_solver.cpp @@ -2,28 +2,33 @@ #ifdef GMGPOLAR_USE_MUMPS #include - #include "coo_mumps_solver.h" + + #include "../../../include/LinearAlgebra/Vector/vector.h" + #include "../../../include/LinearAlgebra/Matrix/coo_matrix.h" + #include "../../../include/LinearAlgebra/Solvers/coo_mumps_solver.h" // ----------------------------------------------------------------------- // Test 1: General (non-symmetric) 4x4 system // // Matrix A (non-symmetric): -// [ 4 1 0 0 ] -// [ 2 5 1 0 ] -// [ 0 1 6 2 ] -// [ 0 0 3 7 ] +// [ 1 0 2 0 ] +// [ 3 0 4 5 ] +// [ 0 6 7 0 ] +// [ 0 8 0 9 ] // -// RHS b = [5, 8, 9, 10]^T +// RHS b = [2, 4, 6, 8] // // Expected solution computed from A*x = b. +// +// x = [140.0 / 43.0, 149.0 / 86.0, -27.0/43.0, -28.0/43.0] // ----------------------------------------------------------------------- TEST(CooMumpsSolverTest, GeneralNonSymmetric4x4) { using triplet = SparseMatrixCOO::triplet_type; // All non-zero entries (0-based row/col indices) - std::vector entries = {{0, 0, 4.0}, {0, 1, 1.0}, {1, 0, 2.0}, {1, 1, 5.0}, {1, 2, 1.0}, - {2, 1, 1.0}, {2, 2, 6.0}, {2, 3, 2.0}, {3, 2, 3.0}, {3, 3, 7.0}}; + std::vector entries = {{0, 0, 1.0}, {0, 2, 2.0}, {1, 0, 3.0}, {1, 2, 4.0}, {1, 3, 5.0}, + {2, 1, 6.0}, {2, 2, 7.0}, {3, 1, 8.0}, {3, 3, 9.0}}; SparseMatrixCOO mat(4, 4, entries); mat.is_symmetric(false); @@ -31,42 +36,34 @@ TEST(CooMumpsSolverTest, GeneralNonSymmetric4x4) CooMumpsSolver solver(std::move(mat)); Vector rhs("rhs", 4); - rhs(0) = 5.0; - rhs(1) = 8.0; - rhs(2) = 9.0; - rhs(3) = 10.0; + rhs(0) = 2.0; + rhs(1) = 4.0; + rhs(2) = 6.0; + rhs(3) = 8.0; solver.solve(rhs); - // Verify A*x = b by back-substitution check - // Reference solution (computed analytically / via numpy): - // x ~ [0.9526, 0.2105, 0.9298, 0.8319] - const double tol = 1e-10; - EXPECT_NEAR(rhs(0), 4.0 * 0.9526 + 1.0 * 0.2105, 1e-3); - - // More robust: re-multiply and check residual - std::vector x = {rhs(0), rhs(1), rhs(2), rhs(3)}; + Vector solution("solution", 4); + solution(0) = 140.0 / 43.0; + solution(1) = 149.0 / 86.0; + solution(2) = -27.0 / 43.0; + solution(3) = -28.0 / 43.0; - // A*x - double Ax0 = 4 * x[0] + 1 * x[1]; - double Ax1 = 2 * x[0] + 5 * x[1] + 1 * x[2]; - double Ax2 = 1 * x[1] + 6 * x[2] + 2 * x[3]; - double Ax3 = 3 * x[2] + 7 * x[3]; - - EXPECT_NEAR(Ax0, 5.0, tol); - EXPECT_NEAR(Ax1, 8.0, tol); - EXPECT_NEAR(Ax2, 9.0, tol); - EXPECT_NEAR(Ax3, 10.0, tol); + const double tol = 1e-10; + EXPECT_NEAR(rhs(0), solution(0), tol); + EXPECT_NEAR(rhs(1), solution(1), tol); + EXPECT_NEAR(rhs(2), solution(2), tol); + EXPECT_NEAR(rhs(3), solution(3), tol); } // ----------------------------------------------------------------------- // Test 2: Symmetric positive-definite 4x4 system (lower triangle only) // // Matrix A (SPD): -// [ 4 2 0 0 ] -// [ 2 5 1 0 ] -// [ 0 1 6 2 ] -// [ 0 0 2 7 ] +// [ 4 0 2 0 ] +// [ 0 5 1 3 ] +// [ 2 1 6 2 ] +// [ 0 3 2 7 ] // // Only lower triangular entries provided (including diagonal). // RHS b = [6, 8, 9, 11]^T @@ -76,8 +73,8 @@ TEST(CooMumpsSolverTest, SymmetricPositiveDefinite4x4) using triplet = SparseMatrixCOO::triplet_type; // Lower triangular entries only (0-based row/col indices) - std::vector entries = {{0, 0, 4.0}, {1, 0, 2.0}, {1, 1, 5.0}, {2, 1, 1.0}, - {2, 2, 6.0}, {3, 2, 2.0}, {3, 3, 7.0}}; + std::vector entries = {{0, 0, 4.0}, {1, 1, 5.0}, {2, 0, 2.0}, {2, 1, 1.0}, + {2, 2, 6.0}, {3, 1, 3.0}, {3, 2, 2.0}, {3, 3, 7.0}}; SparseMatrixCOO mat(4, 4, entries); mat.is_symmetric(true); // SPD, half-entries only @@ -85,27 +82,24 @@ TEST(CooMumpsSolverTest, SymmetricPositiveDefinite4x4) CooMumpsSolver solver(std::move(mat)); Vector rhs("rhs", 4); - rhs(0) = 6.0; - rhs(1) = 8.0; - rhs(2) = 9.0; - rhs(3) = 11.0; + rhs(0) = 2.0; + rhs(1) = 4.0; + rhs(2) = 6.0; + rhs(3) = 8.0; solver.solve(rhs); - // Verify residual: A*x = b using the full (symmetric) matrix - std::vector x = {rhs(0), rhs(1), rhs(2), rhs(3)}; + Vector solution("solution", 4); + solution(0) = 9.0 / 46.0; + solution(1) = 3.0 / 23.0; + solution(2) = 14.0 / 23.0; + solution(3) = 21.0 / 23.0; const double tol = 1e-10; - - double Ax0 = 4 * x[0] + 2 * x[1]; - double Ax1 = 2 * x[0] + 5 * x[1] + 1 * x[2]; - double Ax2 = 1 * x[1] + 6 * x[2] + 2 * x[3]; - double Ax3 = 2 * x[2] + 7 * x[3]; - - EXPECT_NEAR(Ax0, 6.0, tol); - EXPECT_NEAR(Ax1, 8.0, tol); - EXPECT_NEAR(Ax2, 9.0, tol); - EXPECT_NEAR(Ax3, 11.0, tol); + EXPECT_NEAR(rhs(0), solution(0), tol); + EXPECT_NEAR(rhs(1), solution(1), tol); + EXPECT_NEAR(rhs(2), solution(2), tol); + EXPECT_NEAR(rhs(3), solution(3), tol); } #endif // GMGPOLAR_USE_MUMPS \ No newline at end of file From d93ccaf733440a906065e55e1ddcdb032dedf544 Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Sun, 22 Feb 2026 17:01:42 +0100 Subject: [PATCH 11/15] Update coo_mumps_solver.cpp --- tests/LinearAlgebra/Solvers/coo_mumps_solver.cpp | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/tests/LinearAlgebra/Solvers/coo_mumps_solver.cpp b/tests/LinearAlgebra/Solvers/coo_mumps_solver.cpp index 02cb0a52..9ce60ab8 100644 --- a/tests/LinearAlgebra/Solvers/coo_mumps_solver.cpp +++ b/tests/LinearAlgebra/Solvers/coo_mumps_solver.cpp @@ -16,11 +16,7 @@ // [ 0 6 7 0 ] // [ 0 8 0 9 ] // -// RHS b = [2, 4, 6, 8] -// -// Expected solution computed from A*x = b. -// -// x = [140.0 / 43.0, 149.0 / 86.0, -27.0/43.0, -28.0/43.0] +// RHS b = [2, 4, 6, 8]^T // ----------------------------------------------------------------------- TEST(CooMumpsSolverTest, GeneralNonSymmetric4x4) { @@ -66,7 +62,7 @@ TEST(CooMumpsSolverTest, GeneralNonSymmetric4x4) // [ 0 3 2 7 ] // // Only lower triangular entries provided (including diagonal). -// RHS b = [6, 8, 9, 11]^T +// RHS b = [2, 4, 6, 8]^T // ----------------------------------------------------------------------- TEST(CooMumpsSolverTest, SymmetricPositiveDefinite4x4) { From 8f88e53a389f08ec5772ba87c8a45ea04f852310 Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Mon, 23 Feb 2026 11:42:58 +0100 Subject: [PATCH 12/15] Update smootherGive.cpp --- src/Smoother/SmootherGive/smootherGive.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Smoother/SmootherGive/smootherGive.cpp b/src/Smoother/SmootherGive/smootherGive.cpp index 7b3eeb7a..595aeee3 100644 --- a/src/Smoother/SmootherGive/smootherGive.cpp +++ b/src/Smoother/SmootherGive/smootherGive.cpp @@ -13,7 +13,7 @@ SmootherGive::SmootherGive(const PolarGrid& grid, const LevelCache& level_cache, radial_tridiagonal_solver_.setup(); #ifdef GMGPOLAR_USE_MUMPS - inner_boundary_mumps_solver_.emplace(std::move(inner_boundary_circle_matrix_)); + inner_boundary_mumps_solver_.emplace(inner_boundary_circle_matrix_); #else inner_boundary_lu_solver_ = SparseLUSolver(inner_boundary_circle_matrix_); #endif From 6dcf5d09b8f02dc56d90598b9face1ec38b0bd40 Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Mon, 23 Feb 2026 11:43:24 +0100 Subject: [PATCH 13/15] Update smootherTake.cpp --- src/Smoother/SmootherTake/smootherTake.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Smoother/SmootherTake/smootherTake.cpp b/src/Smoother/SmootherTake/smootherTake.cpp index 8ba590e6..76dc8b89 100644 --- a/src/Smoother/SmootherTake/smootherTake.cpp +++ b/src/Smoother/SmootherTake/smootherTake.cpp @@ -13,7 +13,7 @@ SmootherTake::SmootherTake(const PolarGrid& grid, const LevelCache& level_cache, radial_tridiagonal_solver_.setup(); #ifdef GMGPOLAR_USE_MUMPS - inner_boundary_mumps_solver_.emplace(std::move(inner_boundary_circle_matrix_)); + inner_boundary_mumps_solver_.emplace(inner_boundary_circle_matrix_); #else inner_boundary_lu_solver_ = SparseLUSolver(inner_boundary_circle_matrix_); #endif From 7f623725f59a67bb5e4348632893f96d8401e4df Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Mon, 23 Feb 2026 11:43:48 +0100 Subject: [PATCH 14/15] Update extrapolatedSmootherGive.cpp --- .../ExtrapolatedSmootherGive/extrapolatedSmootherGive.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.cpp b/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.cpp index 1a34f678..f7a27781 100644 --- a/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.cpp +++ b/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.cpp @@ -15,7 +15,7 @@ ExtrapolatedSmootherGive::ExtrapolatedSmootherGive(const PolarGrid& grid, const radial_tridiagonal_solver_.setup(); #ifdef GMGPOLAR_USE_MUMPS - inner_boundary_mumps_solver_.emplace(std::move(inner_boundary_circle_matrix_)); + inner_boundary_mumps_solver_.emplace(inner_boundary_circle_matrix_); #else inner_boundary_lu_solver_ = SparseLUSolver(inner_boundary_circle_matrix_); #endif From fdc909f3dd9aa61f010358eb929371a68d7c231e Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Mon, 23 Feb 2026 11:44:11 +0100 Subject: [PATCH 15/15] Update extrapolatedSmootherTake.cpp --- .../ExtrapolatedSmootherTake/extrapolatedSmootherTake.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ExtrapolatedSmoother/ExtrapolatedSmootherTake/extrapolatedSmootherTake.cpp b/src/ExtrapolatedSmoother/ExtrapolatedSmootherTake/extrapolatedSmootherTake.cpp index 525d7007..dfa4b963 100644 --- a/src/ExtrapolatedSmoother/ExtrapolatedSmootherTake/extrapolatedSmootherTake.cpp +++ b/src/ExtrapolatedSmoother/ExtrapolatedSmootherTake/extrapolatedSmootherTake.cpp @@ -15,7 +15,7 @@ ExtrapolatedSmootherTake::ExtrapolatedSmootherTake(const PolarGrid& grid, const radial_tridiagonal_solver_.setup(); #ifdef GMGPOLAR_USE_MUMPS - inner_boundary_mumps_solver_.emplace(std::move(inner_boundary_circle_matrix_)); + inner_boundary_mumps_solver_.emplace(inner_boundary_circle_matrix_); #else inner_boundary_lu_solver_ = SparseLUSolver(inner_boundary_circle_matrix_); #endif