Merge pull request #69 from cvanaret/mumps_bugfix

Bug fix in MUMPSSolver
cvanaret · Nov 4, 2024 · b6159bb · b6159bb
2 parents baae5c1 + 8557797
commit b6159bb
Show file tree

Hide file tree

Showing 3 changed files with 57 additions and 32 deletions.
diff --git 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<size_t, double>& /*matrix*/) {
+   void MUMPSSolver::do_numerical_factorization(const SymmetricMatrix<size_t, double>& matrix) {
       this->mumps_structure.job = MUMPSSolver::JOB_FACTORIZATION;
-      this->mumps_structure.a = this->COO_matrix.data_pointer();
+      this->mumps_structure.a = const_cast<double*>(matrix.data_pointer());
       dmumps_c(&this->mumps_structure);
    }
 

diff --git 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<double, n> reference{1., 2., 3., 4., 5.};
-
-Vector<double> create_rhs() {
-   Vector<double> rhs{8., 45., 31., 15., 17.};
-   return rhs;
-}
-
 TEST(MA57Solver, SystemSize5) {
+   const size_t n = 5;
+   const size_t nnz = 7;
    SymmetricMatrix<size_t, double> 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<double> rhs = create_rhs();
+   const Vector<double> rhs{8., 45., 31., 15., 17.};
    Vector<double> result(n);
    result.fill(0.);
+   const std::array<double, n> 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<size_t, double> 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<double> rhs = create_rhs();
-   Vector<double> 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);
-}
+}
+
+TEST(MA57Solver, SingularMatrix) {
+   const size_t n = 4;
+   const size_t nnz = 7;
+   // comes from hs015 solved with byrd preset
+   SymmetricMatrix<size_t, double> 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
@@ -7,17 +7,11 @@
 
 using namespace uno;
 
-const size_t n = 5;
-const size_t nnz = 7;
-const std::array<double, n> reference{1., 2., 3., 4., 5.};
-const double tolerance = 1e-8;
-
-Vector<double> create_my_rhs() {
-   Vector<double> 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<size_t, double> 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<double> rhs = create_my_rhs();
+   const Vector<double> rhs{8., 45., 31., 15., 17.};
    Vector<double> result(n);
    result.fill(0.);
+   const std::array<double, n> 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<size_t, double> 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<double> rhs = create_my_rhs();
-   Vector<double> 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<size_t, double> 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());
 }