From 8557797b64cadcde1ee7c2e00e19af6fd0d39b79 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Mon, 4 Nov 2024 13:33:50 +0100 Subject: [PATCH] Bugfix in MUMPSSolver: do numerical factorization on latest value of the matrix, in case it was regularized + added unit test on singularity detection for MA57 and MUMPS solvers --- uno/solvers/MUMPS/MUMPSSolver.cpp | 4 +-- unotest/MA57SolverTests.cpp | 42 +++++++++++++++++++----------- unotest/MUMPSSolverTests.cpp | 43 ++++++++++++++++++++----------- 3 files changed, 57 insertions(+), 32 deletions(-) diff --git a/uno/solvers/MUMPS/MUMPSSolver.cpp b/uno/solvers/MUMPS/MUMPSSolver.cpp index 7c3b901c..8a754c3d 100644 --- a/uno/solvers/MUMPS/MUMPSSolver.cpp +++ b/uno/solvers/MUMPS/MUMPSSolver.cpp @@ -54,9 +54,9 @@ namespace uno { dmumps_c(&this->mumps_structure); } - void MUMPSSolver::do_numerical_factorization(const SymmetricMatrix& /*matrix*/) { + void MUMPSSolver::do_numerical_factorization(const SymmetricMatrix& matrix) { this->mumps_structure.job = MUMPSSolver::JOB_FACTORIZATION; - this->mumps_structure.a = this->COO_matrix.data_pointer(); + this->mumps_structure.a = const_cast(matrix.data_pointer()); dmumps_c(&this->mumps_structure); } diff --git a/unotest/MA57SolverTests.cpp b/unotest/MA57SolverTests.cpp index d7664444..839dbd23 100644 --- a/unotest/MA57SolverTests.cpp +++ b/unotest/MA57SolverTests.cpp @@ -7,16 +7,9 @@ using namespace uno; -const size_t n = 5; -const size_t nnz = 7; -const std::array reference{1., 2., 3., 4., 5.}; - -Vector create_rhs() { - Vector rhs{8., 45., 31., 15., 17.}; - return rhs; -} - TEST(MA57Solver, SystemSize5) { + const size_t n = 5; + const size_t nnz = 7; SymmetricMatrix matrix(n, nnz, false, "COO"); matrix.insert(2., 0, 0); matrix.insert(3., 0, 1); @@ -25,9 +18,10 @@ TEST(MA57Solver, SystemSize5) { matrix.insert(1., 2, 2); matrix.insert(5., 2, 3); matrix.insert(1., 4, 4); - const Vector rhs = create_rhs(); + const Vector rhs{8., 45., 31., 15., 17.}; Vector result(n); result.fill(0.); + const std::array reference{1., 2., 3., 4., 5.}; MA57Solver solver(n, nnz); solver.do_symbolic_factorization(matrix); @@ -40,6 +34,8 @@ TEST(MA57Solver, SystemSize5) { } TEST(MA57Solver, Inertia) { + const size_t n = 5; + const size_t nnz = 7; SymmetricMatrix matrix(n, nnz, false, "COO"); matrix.insert(2., 0, 0); matrix.insert(3., 0, 1); @@ -48,17 +44,33 @@ TEST(MA57Solver, Inertia) { matrix.insert(1., 2, 2); matrix.insert(5., 2, 3); matrix.insert(1., 4, 4); - const Vector rhs = create_rhs(); - Vector result(n); - result.fill(0.); MA57Solver solver(n, nnz); solver.do_symbolic_factorization(matrix); solver.do_numerical_factorization(matrix); - solver.solve_indefinite_system(matrix, rhs, result); const auto [number_positive, number_negative, number_zero] = solver.get_inertia(); ASSERT_EQ(number_positive, 3); ASSERT_EQ(number_negative, 2); ASSERT_EQ(number_zero, 0); -} \ No newline at end of file +} + +TEST(MA57Solver, SingularMatrix) { + const size_t n = 4; + const size_t nnz = 7; + // comes from hs015 solved with byrd preset + SymmetricMatrix matrix(n, nnz, false, "COO"); + matrix.insert( -0.0198, 0, 0); + matrix.insert(0.625075, 0, 0); + matrix.insert(-0.277512, 0, 1); + matrix.insert(-0.624975, 1, 1); + matrix.insert(0.625075, 1, 1); + matrix.insert(0., 2, 2); + matrix.insert(0., 3, 3); + MA57Solver solver(n, nnz); + solver.do_symbolic_factorization(matrix); + solver.do_numerical_factorization(matrix); + + // expected inertia (1, 1, 2) + ASSERT_TRUE(solver.matrix_is_singular()); +} diff --git a/unotest/MUMPSSolverTests.cpp b/unotest/MUMPSSolverTests.cpp index 5586f6c0..44b3c297 100644 --- a/unotest/MUMPSSolverTests.cpp +++ b/unotest/MUMPSSolverTests.cpp @@ -7,17 +7,11 @@ using namespace uno; -const size_t n = 5; -const size_t nnz = 7; -const std::array reference{1., 2., 3., 4., 5.}; -const double tolerance = 1e-8; - -Vector create_my_rhs() { - Vector rhs{8., 45., 31., 15., 17.}; - return rhs; -} - TEST(MUMPSSolver, SystemSize5) { + const double tolerance = 1e-8; + + const size_t n = 5; + const size_t nnz = 7; SymmetricMatrix matrix(n, nnz, false, "COO"); matrix.insert(2., 0, 0); matrix.insert(3., 0, 1); @@ -26,9 +20,10 @@ TEST(MUMPSSolver, SystemSize5) { matrix.insert(1., 2, 2); matrix.insert(5., 2, 3); matrix.insert(1., 4, 4); - const Vector rhs = create_my_rhs(); + const Vector rhs{8., 45., 31., 15., 17.}; Vector result(n); result.fill(0.); + const std::array reference{1., 2., 3., 4., 5.}; MUMPSSolver solver(n, nnz); solver.do_symbolic_factorization(matrix); @@ -58,6 +53,8 @@ array([-7.83039207, 8.94059148, -3.50815575, 1.7888887 , 4.60906763]) */ TEST(MUMPSSolver, Inertia) { + const size_t n = 5; + const size_t nnz = 7; SymmetricMatrix matrix(n, nnz, false, "COO"); matrix.insert(2., 0, 0); matrix.insert(3., 0, 1); @@ -66,17 +63,33 @@ TEST(MUMPSSolver, Inertia) { matrix.insert(1., 2, 2); matrix.insert(5., 2, 3); matrix.insert(1., 4, 4); - const Vector rhs = create_my_rhs(); - Vector result(n); - result.fill(0.); MUMPSSolver solver(n, nnz); solver.do_symbolic_factorization(matrix); solver.do_numerical_factorization(matrix); - solver.solve_indefinite_system(matrix, rhs, result); const auto [number_positive, number_negative, number_zero] = solver.get_inertia(); ASSERT_EQ(number_positive, 3); ASSERT_EQ(number_negative, 2); ASSERT_EQ(number_zero, 0); +} + +TEST(MUMPSSolver, SingularMatrix) { + const size_t n = 4; + const size_t nnz = 7; + // comes from hs015 solved with byrd preset + SymmetricMatrix matrix(n, nnz, false, "COO"); + matrix.insert( -0.0198, 0, 0); + matrix.insert(0.625075, 0, 0); + matrix.insert(-0.277512, 0, 1); + matrix.insert(-0.624975, 1, 1); + matrix.insert(0.625075, 1, 1); + matrix.insert(0., 2, 2); + matrix.insert(0., 3, 3); + MUMPSSolver solver(n, nnz); + solver.do_symbolic_factorization(matrix); + solver.do_numerical_factorization(matrix); + + // expected inertia (1, 1, 2) + ASSERT_TRUE(solver.matrix_is_singular()); } \ No newline at end of file