Skip to content

Commit

Permalink
Bugfix in MUMPSSolver: do numerical factorization on latest value of …
Browse files Browse the repository at this point in the history 
…the matrix, in case it was regularized + added unit test on singularity detection for MA57 and MUMPS solvers
  • Loading branch information
@cvanaret
cvanaret committed Nov 4, 2024
1 parent baae5c1 commit 8557797
Show file tree
Hide file tree
Showing 3 changed files with 57 additions and 32 deletions.
4 changes: 2 additions & 2 deletions uno/solvers/MUMPS/MUMPSSolver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}

Expand Down
42 changes: 27 additions & 15 deletions unotest/MA57SolverTests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -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);
Expand All @@ -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);
Expand All @@ -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());
}
43 changes: 28 additions & 15 deletions unotest/MUMPSSolverTests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -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);
Expand Down Expand Up @@ -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);
Expand All @@ -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());
}

0 comments on commit 8557797

Please sign in to comment.