KOKO\Mimi
128 lines
5.1 KiB
C++
128 lines
5.1 KiB
C++
#include "fesa/solvers/linear/linear_solver.h"
|
|
|
|
#include <gtest/gtest.h>
|
|
|
|
#include <type_traits>
|
|
#include <utility>
|
|
#include <vector>
|
|
|
|
#include "fesa/fem/dof_manager.h"
|
|
#include "fesa/math/sparse_matrix.h"
|
|
#include "fesa/solvers/linear/mkl_pardiso_solver.h"
|
|
|
|
namespace {
|
|
|
|
fesa::SparseMatrix MakeDenseCsr(const std::size_t rows,
|
|
const std::size_t columns,
|
|
const std::vector<double>& values) {
|
|
EXPECT_EQ(values.size(), rows * columns);
|
|
|
|
fesa::SparsePattern pattern;
|
|
std::vector<fesa::CooContribution> contributions;
|
|
pattern.row_offsets.reserve(rows + 1U);
|
|
pattern.row_offsets.push_back(0U);
|
|
for (std::size_t row = 0U; row < rows; ++row) {
|
|
for (std::size_t column = 0U; column < columns; ++column) {
|
|
pattern.column_indices.push_back(column);
|
|
contributions.push_back(
|
|
{row, column, values[row * columns + column], row, column});
|
|
}
|
|
pattern.row_offsets.push_back(pattern.column_indices.size());
|
|
}
|
|
|
|
auto matrix = fesa::SparseMatrix::FromCoo(rows, columns,
|
|
std::move(contributions), pattern);
|
|
EXPECT_TRUE(matrix.HasValue());
|
|
return std::move(matrix.Value());
|
|
}
|
|
|
|
void ExpectSolverFailure(const fesa::Status& status) {
|
|
EXPECT_FALSE(status.IsOk());
|
|
EXPECT_EQ(status.Category(), fesa::FailureCategory::kSolver);
|
|
ASSERT_FALSE(status.Diagnostics().empty());
|
|
EXPECT_EQ(status.Diagnostics().front().severity, fesa::Severity::kError);
|
|
}
|
|
|
|
} // namespace
|
|
|
|
TEST(MklPardisoSolver, RejectsInvalidCsrStateAndDimensions) {
|
|
static_assert(std::is_base_of_v<fesa::LinearSolver, fesa::MklPardisoSolver>);
|
|
static_assert(std::has_virtual_destructor_v<fesa::LinearSolver>);
|
|
|
|
fesa::MklPardisoSolver solver;
|
|
fesa::Vector untouched{2U};
|
|
untouched[0U] = 17.0;
|
|
untouched[1U] = -4.0;
|
|
const auto before_factorize = solver.Solve(fesa::Vector{2U, 1.0}, untouched);
|
|
ExpectSolverFailure(before_factorize);
|
|
EXPECT_EQ(before_factorize.Diagnostics().front().code,
|
|
"solver-not-factorized");
|
|
EXPECT_DOUBLE_EQ(untouched[0U], 17.0);
|
|
EXPECT_DOUBLE_EQ(untouched[1U], -4.0);
|
|
|
|
// A fully constrained model has a valid 0x0 Kff. It still observes the
|
|
// factorize-then-solve lifecycle without invoking a numerical backend.
|
|
const auto empty = MakeDenseCsr(0U, 0U, {});
|
|
ASSERT_TRUE(solver.Factorize(empty).IsOk());
|
|
fesa::Vector empty_solution{0U};
|
|
EXPECT_TRUE(solver.Solve(fesa::Vector{0U}, empty_solution).IsOk());
|
|
EXPECT_EQ(empty_solution.Size(), 0U);
|
|
|
|
// Refactorization from the trivial state must establish ordinary PARDISO
|
|
// state rather than retaining a zero-equation shortcut.
|
|
const auto spd = MakeDenseCsr(2U, 2U, {4.0, 1.0, 1.0, 3.0});
|
|
ASSERT_TRUE(solver.Factorize(spd).IsOk());
|
|
fesa::Vector solution{2U};
|
|
ASSERT_TRUE(solver.Solve(fesa::Vector{2U, 1.0}, solution).IsOk());
|
|
EXPECT_NEAR(solution[0U], 2.0 / 11.0, 1.0e-14);
|
|
EXPECT_NEAR(solution[1U], 3.0 / 11.0, 1.0e-14);
|
|
|
|
const double solved_first = solution[0U];
|
|
const double solved_second = solution[1U];
|
|
ExpectSolverFailure(solver.Solve(fesa::Vector{1U, 1.0}, solution));
|
|
EXPECT_DOUBLE_EQ(solution[0U], solved_first);
|
|
EXPECT_DOUBLE_EQ(solution[1U], solved_second);
|
|
|
|
fesa::Vector wrong_solution{1U};
|
|
wrong_solution[0U] = 41.0;
|
|
ExpectSolverFailure(solver.Solve(fesa::Vector{2U, 1.0}, wrong_solution));
|
|
EXPECT_DOUBLE_EQ(wrong_solution[0U], 41.0);
|
|
|
|
const auto rectangular = MakeDenseCsr(2U, 3U, {2.0, 0.0, 0.0, 0.0, 3.0, 0.0});
|
|
const auto rectangular_status = solver.Factorize(rectangular);
|
|
ExpectSolverFailure(rectangular_status);
|
|
EXPECT_EQ(rectangular_status.Diagnostics().front().code,
|
|
"solver-matrix-not-square");
|
|
|
|
fesa::SparsePattern invalid_pattern{{0U, 2U}, {0U}};
|
|
auto invalid_csr = fesa::SparseMatrix::FromCoo(
|
|
1U, 1U, {{0U, 0U, 1.0, 0U, 0U}}, invalid_pattern);
|
|
EXPECT_FALSE(invalid_csr.HasValue());
|
|
|
|
const auto nonsymmetric = MakeDenseCsr(2U, 2U, {2.0, 1.0, 0.0, 3.0});
|
|
const auto nonsymmetric_status = solver.Factorize(nonsymmetric);
|
|
ExpectSolverFailure(nonsymmetric_status);
|
|
EXPECT_EQ(nonsymmetric_status.Diagnostics().front().code,
|
|
"solver-matrix-not-symmetric");
|
|
|
|
const auto scaled_nonsymmetric =
|
|
MakeDenseCsr(2U, 2U, {2.0e-20, 1.0e-20, 1.1e-20, 3.0e-20});
|
|
const auto scaled_nonsymmetric_status = solver.Factorize(scaled_nonsymmetric);
|
|
// Stop this case before inspecting diagnostics when the production code
|
|
// incorrectly accepts the matrix; this keeps the RED failure deterministic.
|
|
ASSERT_FALSE(scaled_nonsymmetric_status.IsOk());
|
|
ExpectSolverFailure(scaled_nonsymmetric_status);
|
|
EXPECT_EQ(scaled_nonsymmetric_status.Diagnostics().front().code,
|
|
"solver-matrix-not-symmetric");
|
|
|
|
fesa::SparsePattern no_diagonal_pattern{{0U, 1U, 2U}, {1U, 0U}};
|
|
auto no_diagonal = fesa::SparseMatrix::FromCoo(
|
|
2U, 2U, {{0U, 1U, 1.0, 0U, 0U}, {1U, 0U, 1.0, 1U, 0U}},
|
|
no_diagonal_pattern);
|
|
ASSERT_TRUE(no_diagonal.HasValue());
|
|
const auto no_diagonal_status = solver.Factorize(no_diagonal.Value());
|
|
ExpectSolverFailure(no_diagonal_status);
|
|
EXPECT_EQ(no_diagonal_status.Diagnostics().front().code,
|
|
"solver-missing-diagonal");
|
|
}
|