From e216660d518eff0a10cfe549c4af8d1ca6e0eb7c Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Fri, 31 Jul 2026 16:16:28 +0900 Subject: [PATCH] =?UTF-8?q?feat(equation-and-linear-solve):=20step=202=20?= =?UTF-8?q?=E2=80=94=20pardiso-linear-solver?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- CMakeLists.txt | 2 + include/fesa/solvers/linear/linear_solver.hpp | 26 ++ .../solvers/linear/pardiso_linear_solver.hpp | 19 ++ .../solvers/linear/pardiso_linear_solver.cpp | 316 ++++++++++++++++++ tests/CMakeLists.txt | 30 ++ .../linear/pardiso_linear_solver_test.cpp | 148 ++++++++ 6 files changed, 541 insertions(+) create mode 100644 include/fesa/solvers/linear/linear_solver.hpp create mode 100644 include/fesa/solvers/linear/pardiso_linear_solver.hpp create mode 100644 src/fesa/solvers/linear/pardiso_linear_solver.cpp create mode 100644 tests/unit/solvers/linear/pardiso_linear_solver_test.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index a33195d..eaba369 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -38,6 +38,7 @@ add_library(fesa_core STATIC src/fesa/io/abaqus/semantic_mapper.cpp src/fesa/model/domain.cpp src/fesa/model/domain_builder.cpp + src/fesa/solvers/linear/pardiso_linear_solver.cpp ) target_include_directories(fesa_core @@ -47,6 +48,7 @@ target_include_directories(fesa_core target_compile_features(fesa_core PUBLIC cxx_std_20) target_compile_options(fesa_core PRIVATE /W4 /permissive- /EHsc) +target_link_libraries(fesa_core PRIVATE MKL::MKL) add_executable(fesa src/fesa/cli/main.cpp diff --git a/include/fesa/solvers/linear/linear_solver.hpp b/include/fesa/solvers/linear/linear_solver.hpp new file mode 100644 index 0000000..e258861 --- /dev/null +++ b/include/fesa/solvers/linear/linear_solver.hpp @@ -0,0 +1,26 @@ +#pragma once + +#include +#include + +#include +#include + +namespace fesa { + +struct LinearSolveResult final { + std::vector solution; + double relative_residual{}; + std::vector diagnostics; +}; + +class LinearSolver { +public: + virtual ~LinearSolver() = default; + + [[nodiscard]] virtual LinearSolveResult solve( + const SymmetricCsr& matrix, + std::span rhs) = 0; +}; + +} // namespace fesa diff --git a/include/fesa/solvers/linear/pardiso_linear_solver.hpp b/include/fesa/solvers/linear/pardiso_linear_solver.hpp new file mode 100644 index 0000000..d7de39e --- /dev/null +++ b/include/fesa/solvers/linear/pardiso_linear_solver.hpp @@ -0,0 +1,19 @@ +#pragma once + +#include + +#include + +namespace fesa { + +class PardisoLinearSolver final : public LinearSolver { +public: + PardisoLinearSolver(); + ~PardisoLinearSolver() override; + + [[nodiscard]] LinearSolveResult solve( + const SymmetricCsr& matrix, + std::span rhs) override; +}; + +} // namespace fesa diff --git a/src/fesa/solvers/linear/pardiso_linear_solver.cpp b/src/fesa/solvers/linear/pardiso_linear_solver.cpp new file mode 100644 index 0000000..172a2c5 --- /dev/null +++ b/src/fesa/solvers/linear/pardiso_linear_solver.cpp @@ -0,0 +1,316 @@ +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace fesa { +namespace { + +static_assert( + std::is_same_v, + "FESA requires the oneMKL LP64 interface."); + +Diagnostic solver_error(std::string code, std::string message) { + return { + DiagnosticStage::solver, + Severity::error, + std::move(code), + std::move(message), + std::nullopt, + }; +} + +std::optional validate_matrix( + const SymmetricCsr& matrix) { + if (matrix.order == 0) { + return "PARDISO requires a nonempty reduced system."; + } + if (matrix.order > + static_cast( + std::numeric_limits::max())) { + return "Matrix order exceeds the oneMKL LP64 index range."; + } + if (matrix.row_offsets.size() != matrix.order + 1 || + matrix.row_offsets.front() != 0) { + return "CSR row offsets must contain order + 1 entries " + "starting at zero."; + } + if (matrix.column_indices.size() != matrix.values.size()) { + return "CSR column and value counts must match."; + } + + MKL_INT previous_offset = 0; + for (const MKL_INT offset : matrix.row_offsets) { + if (offset < previous_offset || offset < 0 || + static_cast(offset) > + matrix.column_indices.size()) { + return "CSR row offsets must be nondecreasing and in range."; + } + previous_offset = offset; + } + if (static_cast(matrix.row_offsets.back()) != + matrix.column_indices.size()) { + return "The final CSR row offset must equal the entry count."; + } + + for (std::size_t row = 0; row < matrix.order; ++row) { + MKL_INT previous_column = -1; + bool has_diagonal = false; + const std::size_t begin = + static_cast(matrix.row_offsets[row]); + const std::size_t end = + static_cast(matrix.row_offsets[row + 1]); + for (std::size_t entry = begin; entry < end; ++entry) { + const MKL_INT column = matrix.column_indices[entry]; + if (column < static_cast(row) || + column >= static_cast(matrix.order) || + column <= previous_column) { + return "CSR rows must contain sorted unique " + "upper-triangle columns."; + } + if (!std::isfinite(matrix.values[entry])) { + return "CSR values must be finite."; + } + has_diagonal = + has_diagonal || column == static_cast(row); + previous_column = column; + } + if (!has_diagonal) { + return "Every CSR row must contain its diagonal entry."; + } + } + return std::nullopt; +} + +class PardisoSession final { +public: + PardisoSession() { + pardisoinit(handles_.data(), &matrix_type_, parameters_.data()); + parameters_[26] = 1; + parameters_[34] = 1; + } + + PardisoSession(const PardisoSession&) = delete; + PardisoSession& operator=(const PardisoSession&) = delete; + + ~PardisoSession() noexcept { + if (!active_) { + return; + } + + constexpr MKL_INT release_all = -1; + MKL_INT error = 0; + pardiso( + handles_.data(), + &max_factorizations_, + &matrix_number_, + &matrix_type_, + &release_all, + &order_, + matrix_->values.data(), + matrix_->row_offsets.data(), + matrix_->column_indices.data(), + permutation_.data(), + &right_hand_side_count_, + parameters_.data(), + &message_level_, + right_hand_side_->data(), + solution_->data(), + &error); + } + + MKL_INT execute( + const MKL_INT phase, + const SymmetricCsr& matrix, + std::vector& right_hand_side, + std::vector& solution) { + order_ = static_cast(matrix.order); + matrix_ = &matrix; + right_hand_side_ = &right_hand_side; + solution_ = &solution; + permutation_.resize(matrix.order); + active_ = true; + + MKL_INT error = 0; + pardiso( + handles_.data(), + &max_factorizations_, + &matrix_number_, + &matrix_type_, + &phase, + &order_, + matrix.values.data(), + matrix.row_offsets.data(), + matrix.column_indices.data(), + permutation_.data(), + &right_hand_side_count_, + parameters_.data(), + &message_level_, + right_hand_side.data(), + solution.data(), + &error); + return error; + } + +private: + std::array handles_{}; + std::array parameters_{}; + std::vector permutation_; + const SymmetricCsr* matrix_{}; + std::vector* right_hand_side_{}; + std::vector* solution_{}; + MKL_INT order_{}; + MKL_INT max_factorizations_{1}; + MKL_INT matrix_number_{1}; + MKL_INT matrix_type_{2}; + MKL_INT right_hand_side_count_{1}; + MKL_INT message_level_{}; + bool active_{}; +}; + +double relative_residual( + const SymmetricCsr& matrix, + const std::span rhs, + const std::span solution) { + std::vector residual(rhs.begin(), rhs.end()); + for (double& value : residual) { + value = -value; + } + + for (std::size_t row = 0; row < matrix.order; ++row) { + const std::size_t begin = + static_cast(matrix.row_offsets[row]); + const std::size_t end = + static_cast(matrix.row_offsets[row + 1]); + for (std::size_t entry = begin; entry < end; ++entry) { + const std::size_t column = + static_cast(matrix.column_indices[entry]); + const double value = matrix.values[entry]; + residual[row] += value * solution[column]; + if (column != row) { + residual[column] += value * solution[row]; + } + } + } + + double residual_squared = 0.0; + double rhs_squared = 0.0; + for (std::size_t index = 0; index < rhs.size(); ++index) { + residual_squared += residual[index] * residual[index]; + rhs_squared += rhs[index] * rhs[index]; + } + const double residual_norm = std::sqrt(residual_squared); + const double rhs_norm = std::sqrt(rhs_squared); + return rhs_norm == 0.0 ? residual_norm : residual_norm / rhs_norm; +} + +std::string pardiso_failure( + const std::string_view phase, + const MKL_INT error) { + return "PARDISO " + std::string{phase} + + " failed with error " + std::to_string(error) + "."; +} + +} // namespace + +PardisoLinearSolver::PardisoLinearSolver() = default; + +PardisoLinearSolver::~PardisoLinearSolver() = default; + +LinearSolveResult PardisoLinearSolver::solve( + const SymmetricCsr& matrix, + const std::span rhs) { + LinearSolveResult result; + if (const auto error = validate_matrix(matrix); + error.has_value()) { + result.diagnostics.push_back( + solver_error("solver.invalid_csr", *error)); + return result; + } + if (rhs.size() != matrix.order) { + result.diagnostics.push_back(solver_error( + "solver.dimension_mismatch", + "Right-hand side size must equal the matrix order.")); + return result; + } + if (!std::ranges::all_of(rhs, [](const double value) { + return std::isfinite(value); + })) { + result.diagnostics.push_back(solver_error( + "solver.invalid_rhs", + "Right-hand side values must be finite.")); + return result; + } + + std::vector right_hand_side(rhs.begin(), rhs.end()); + std::vector solution(matrix.order, 0.0); + { + PardisoSession session; + + constexpr MKL_INT analyze = 11; + if (const MKL_INT error = + session.execute( + analyze, matrix, right_hand_side, solution); + error != 0) { + result.diagnostics.push_back(solver_error( + "solver.analysis_failed", + pardiso_failure("analysis", error))); + return result; + } + + constexpr MKL_INT factorize = 22; + if (const MKL_INT error = + session.execute( + factorize, matrix, right_hand_side, solution); + error != 0) { + result.diagnostics.push_back(solver_error( + "solver.factorization_failed", + pardiso_failure("factorization", error))); + return result; + } + + constexpr MKL_INT solve_system = 33; + if (const MKL_INT error = + session.execute( + solve_system, matrix, right_hand_side, solution); + error != 0) { + result.diagnostics.push_back(solver_error( + "solver.solve_failed", + pardiso_failure("solve", error))); + return result; + } + } + if (!std::ranges::all_of(solution, [](const double value) { + return std::isfinite(value); + })) { + result.diagnostics.push_back(solver_error( + "solver.nonfinite_solution", + "PARDISO produced a nonfinite solution.")); + return result; + } + + result.relative_residual = + relative_residual(matrix, rhs, solution); + if (!std::isfinite(result.relative_residual)) { + result.diagnostics.push_back(solver_error( + "solver.nonfinite_residual", + "The independently computed relative residual is nonfinite.")); + return result; + } + result.solution = std::move(solution); + return result; +} + +} // namespace fesa diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 113a566..fe4c5c9 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -345,3 +345,33 @@ add_test( COMMAND "$" --gtest_filter=Reaction.* ) + +add_executable(fesa_linear_solver_tests + unit/solvers/linear/pardiso_linear_solver_test.cpp +) + +target_compile_features(fesa_linear_solver_tests PRIVATE cxx_std_20) +target_compile_options( + fesa_linear_solver_tests + PRIVATE + /W4 + /permissive- + /EHsc +) + +target_link_libraries(fesa_linear_solver_tests + PRIVATE + fesa_core + GTest::gtest_main +) + +add_test( + NAME PardisoLinearSolver + COMMAND "$" +) + +set_property( + TEST PardisoLinearSolver + PROPERTY ENVIRONMENT_MODIFICATION + ${FESA_DEPENDENCY_RUNTIME_MODIFICATIONS} +) diff --git a/tests/unit/solvers/linear/pardiso_linear_solver_test.cpp b/tests/unit/solvers/linear/pardiso_linear_solver_test.cpp new file mode 100644 index 0000000..dc03630 --- /dev/null +++ b/tests/unit/solvers/linear/pardiso_linear_solver_test.cpp @@ -0,0 +1,148 @@ +#include + +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace { + +fesa::SymmetricCsr spd_matrix() { + return { + 3, + {0, 2, 4, 5}, + {0, 1, 1, 2, 2}, + {4.0, 1.0, 3.0, 1.0, 2.0}, + }; +} + +double independent_relative_residual( + const fesa::SymmetricCsr& matrix, + const std::span rhs, + const std::span solution) { + std::vector residual(rhs.begin(), rhs.end()); + for (double& value : residual) { + value = -value; + } + + for (std::size_t row = 0; row < matrix.order; ++row) { + const std::size_t begin = + static_cast(matrix.row_offsets[row]); + const std::size_t end = + static_cast(matrix.row_offsets[row + 1]); + for (std::size_t entry = begin; entry < end; ++entry) { + const std::size_t column = + static_cast(matrix.column_indices[entry]); + const double value = matrix.values[entry]; + residual[row] += value * solution[column]; + if (column != row) { + residual[column] += value * solution[row]; + } + } + } + + double residual_squared = 0.0; + double rhs_squared = 0.0; + for (std::size_t index = 0; index < rhs.size(); ++index) { + residual_squared += residual[index] * residual[index]; + rhs_squared += rhs[index] * rhs[index]; + } + const double residual_norm = std::sqrt(residual_squared); + const double rhs_norm = std::sqrt(rhs_squared); + return rhs_norm == 0.0 ? residual_norm : residual_norm / rhs_norm; +} + +bool has_diagnostic( + const fesa::LinearSolveResult& result, + const std::string_view code) { + return std::ranges::find( + result.diagnostics, code, &fesa::Diagnostic::code) != + result.diagnostics.end(); +} + +TEST(PardisoLinearSolver, SolvesThreeByThreeSpdSystem) { + fesa::PardisoLinearSolver pardiso; + fesa::LinearSolver& solver = pardiso; + const fesa::SymmetricCsr matrix = spd_matrix(); + const std::vector rhs{6.0, 10.0, 8.0}; + + const fesa::LinearSolveResult result = solver.solve(matrix, rhs); + + ASSERT_TRUE(result.diagnostics.empty()); + ASSERT_EQ(result.solution.size(), 3); + EXPECT_NEAR(result.solution[0], 1.0, 1.0e-12); + EXPECT_NEAR(result.solution[1], 2.0, 1.0e-12); + EXPECT_NEAR(result.solution[2], 3.0, 1.0e-12); + const double independent = + independent_relative_residual(matrix, rhs, result.solution); + EXPECT_NEAR(result.relative_residual, independent, 1.0e-15); + EXPECT_LT(result.relative_residual, 1.0e-12); +} + +TEST(PardisoLinearSolver, SupportsRepeatedSolve) { + fesa::PardisoLinearSolver solver; + const fesa::SymmetricCsr matrix = spd_matrix(); + + const fesa::LinearSolveResult first = + solver.solve(matrix, std::vector{6.0, 10.0, 8.0}); + const fesa::LinearSolveResult second = + solver.solve(matrix, std::vector{-3.5, 2.5, 4.5}); + + ASSERT_TRUE(first.diagnostics.empty()); + ASSERT_TRUE(second.diagnostics.empty()); + ASSERT_EQ(second.solution.size(), 3); + EXPECT_NEAR(second.solution[0], -1.0, 1.0e-12); + EXPECT_NEAR(second.solution[1], 0.5, 1.0e-12); + EXPECT_NEAR(second.solution[2], 2.0, 1.0e-12); +} + +TEST(LinearSolver, RejectsInvalidUpperTriangleCsr) { + fesa::PardisoLinearSolver solver; + fesa::SymmetricCsr invalid = spd_matrix(); + invalid.column_indices[2] = 0; + + const fesa::LinearSolveResult result = + solver.solve(invalid, std::vector{6.0, 10.0, 8.0}); + + EXPECT_TRUE(result.solution.empty()); + EXPECT_TRUE(has_diagnostic(result, "solver.invalid_csr")); +} + +TEST(LinearSolver, RejectsRhsDimensionMismatch) { + fesa::PardisoLinearSolver solver; + + const fesa::LinearSolveResult result = + solver.solve(spd_matrix(), std::vector{1.0, 2.0}); + + EXPECT_TRUE(result.solution.empty()); + EXPECT_TRUE(has_diagnostic(result, "solver.dimension_mismatch")); +} + +TEST(PardisoLinearSolver, ReportsSingularMatrix) { + fesa::PardisoLinearSolver solver; + const fesa::SymmetricCsr singular{ + 2, + {0, 2, 3}, + {0, 1, 1}, + {1.0, 1.0, 1.0}, + }; + + const fesa::LinearSolveResult result = + solver.solve(singular, std::vector{2.0, 2.0}); + + EXPECT_TRUE(result.solution.empty()); + EXPECT_FALSE(result.diagnostics.empty()); + EXPECT_TRUE(std::ranges::all_of( + result.diagnostics, + [](const fesa::Diagnostic& diagnostic) { + return diagnostic.stage == fesa::DiagnosticStage::solver && + diagnostic.severity == fesa::Severity::error; + })); +} + +} // namespace