diff --git a/CMakeLists.txt b/CMakeLists.txt index 423765fd..78bb4408 100644 --- 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 index e1883aee..6a26f257 100644 --- 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(dimension), - COO_matrix(dimension, number_nonzeros, false) { + MUMPSSolver::MUMPSSolver(size_t dimension, size_t number_nonzeros) : DirectSymmetricIndefiniteLinearSolver(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& matrix) { - this->save_matrix_to_local_format(matrix); + this->save_sparsity_to_local_format(matrix); this->mumps_structure.n = static_cast(matrix.dimension()); this->mumps_structure.nnz = static_cast(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& matrix) { + void MUMPSSolver::save_sparsity_to_local_format(const SymmetricMatrix& 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(row_index + this->fortran_shift), static_cast(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 index 0c7e6ca4..da7ebd98 100644 --- 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 #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& matrix) override; void do_numerical_factorization(const SymmetricMatrix& matrix) override; - void solve_indefinite_system(const SymmetricMatrix& matrix, const Vector& rhs, - Vector& result) override; + void solve_indefinite_system(const SymmetricMatrix& matrix, const Vector& rhs, Vector& result) override; [[nodiscard]] std::tuple 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 COO_matrix; + // matrix sparsity + std::vector row_indices{}; + std::vector 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& row_index); + void save_sparsity_to_local_format(const SymmetricMatrix& matrix); }; } // namespace