Merge pull request #121 from cvanaret/fix_mumps2

Improve the MUMPSSolver interface
cvanaret · Dec 2, 2024 · d5defd6 · d5defd6
2 parents 6a7537a + 0c1fe24
commit d5defd6
Show file tree

Hide file tree

Showing 3 changed files with 21 additions and 19 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -184,8 +184,6 @@ else()
       list(APPEND LIBRARIES ${ScaLAPACK_LIBRARY})
       list(APPEND DIRECTORIES ${ScaLAPACK_INCLUDE_DIRS})
    else()
-      # sequential
-      add_definitions("-D MUMPS_SEQUENTIAL")
       # link dummy parallel library (provided by MUMPS distribution)
       link_to_uno(MUMPS_MPISEQ_LIBRARY ${MUMPS_MPISEQ_LIBRARY})
    endif()

diff --git a/uno/solvers/MUMPS/MUMPSSolver.cpp b/uno/solvers/MUMPS/MUMPSSolver.cpp
@@ -7,17 +7,18 @@
 #include "mpi.h"
 #endif
 
-#define USE_COMM_WORLD -987654
+#define USE_COMM_WORLD (-987654)
 
 namespace uno {
-   MUMPSSolver::MUMPSSolver(size_t dimension, size_t number_nonzeros) : DirectSymmetricIndefiniteLinearSolver<size_t, double>(dimension),
-         COO_matrix(dimension, number_nonzeros, false) {
+   MUMPSSolver::MUMPSSolver(size_t dimension, size_t number_nonzeros) : DirectSymmetricIndefiniteLinearSolver<size_t, double>(dimension) {
+      this->row_indices.reserve(number_nonzeros);
+      this->column_indices.reserve(number_nonzeros);
+
       this->mumps_structure.sym = MUMPSSolver::GENERAL_SYMMETRIC;
 #if defined(HAS_MPI) && defined(MUMPS_PARALLEL)
       // TODO load number of processes from option file
       this->mumps_structure.par = 1;
-#endif
-#ifdef MUMPS_SEQUENTIAL
+#else
       this->mumps_structure.par = 1;
 #endif
       this->mumps_structure.job = MUMPSSolver::JOB_INIT;
@@ -38,14 +39,14 @@ namespace uno {
    }
 
    void MUMPSSolver::do_symbolic_analysis(const SymmetricMatrix<size_t, double>& matrix) {
-      this->save_matrix_to_local_format(matrix);
+      this->save_sparsity_to_local_format(matrix);
       this->mumps_structure.n = static_cast<int>(matrix.dimension());
       this->mumps_structure.nnz = static_cast<int>(matrix.number_nonzeros());
       this->mumps_structure.job = MUMPSSolver::JOB_ANALYSIS;
       this->mumps_structure.a = nullptr;
       // connect the local COO matrix with the pointers in the structure
-      this->mumps_structure.irn = this->COO_matrix.row_indices_pointer();
-      this->mumps_structure.jcn = this->COO_matrix.column_indices_pointer();
+      this->mumps_structure.irn = this->row_indices.data();
+      this->mumps_structure.jcn = this->column_indices.data();
       dmumps_c(&this->mumps_structure);
    }
 
@@ -87,11 +88,13 @@ namespace uno {
       return this->dimension - this->number_zero_eigenvalues();
    }
 
-   void MUMPSSolver::save_matrix_to_local_format(const SymmetricMatrix<size_t, double>& matrix) {
+   void MUMPSSolver::save_sparsity_to_local_format(const SymmetricMatrix<size_t, double>& matrix) {
       // build the internal matrix representation
-      this->COO_matrix.reset();
-      for (const auto [row_index, column_index, element]: matrix) {
-         this->COO_matrix.insert(element, static_cast<int>(row_index + this->fortran_shift), static_cast<int>(column_index + this->fortran_shift));
+      this->row_indices.clear();
+      this->column_indices.clear();
+      for (const auto [row_index, column_index, _]: matrix) {
+         this->row_indices.emplace_back(row_index + this->fortran_shift);
+         this->column_indices.emplace_back(column_index + this->fortran_shift);
       }
    }
 } // namespace
diff --git a/uno/solvers/MUMPS/MUMPSSolver.hpp b/uno/solvers/MUMPS/MUMPSSolver.hpp
@@ -4,8 +4,8 @@
 #ifndef UNO_MUMPSSOLVER_H
 #define UNO_MUMPSSOLVER_H
 
+#include <vector>
 #include "solvers/DirectSymmetricIndefiniteLinearSolver.hpp"
-#include "linear_algebra/COOSparseStorage.hpp"
 #include "dmumps_c.h"
 
 namespace uno {
@@ -16,8 +16,7 @@ namespace uno {
 
       void do_symbolic_analysis(const SymmetricMatrix<size_t, double>& matrix) override;
       void do_numerical_factorization(const SymmetricMatrix<size_t, double>& matrix) override;
-      void solve_indefinite_system(const SymmetricMatrix<size_t, double>& matrix, const Vector<double>& rhs,
-            Vector<double>& result) override;
+      void solve_indefinite_system(const SymmetricMatrix<size_t, double>& matrix, const Vector<double>& rhs, Vector<double>& result) override;
 
       [[nodiscard]] std::tuple<size_t, size_t, size_t> get_inertia() const override;
       [[nodiscard]] size_t number_negative_eigenvalues() const override;
@@ -28,7 +27,9 @@ namespace uno {
 
    protected:
       DMUMPS_STRUC_C mumps_structure{};
-      COOSparseStorage<int, double> COO_matrix;
+      // matrix sparsity
+      std::vector<int> row_indices{};
+      std::vector<int> column_indices{};
 
       static const int JOB_INIT = -1;
       static const int JOB_END = -2;
@@ -39,7 +40,7 @@ namespace uno {
       static const int GENERAL_SYMMETRIC = 2;
 
       const size_t fortran_shift{1};
-      void save_matrix_to_local_format(const SymmetricMatrix<size_t, double>& row_index);
+      void save_sparsity_to_local_format(const SymmetricMatrix<size_t, double>& matrix);
    };
 } // namespace