From 80f25e569a886f6a35659991fb6cbc980110f784 Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Sun, 9 Aug 2026 20:05:20 +0900 Subject: [PATCH] feat(linear-static-3d-euler-beam): step 20 - mkl-pardiso-solver --- ...tic-3d-euler-beam-implementation-report.md | 63 +++ include/fesa/solvers/linear/linear_solver.hpp | 18 + .../solvers/linear/mkl_pardiso_solver.hpp | 23 ++ src/fesa/CMakeLists.txt | 1 + .../solvers/linear/mkl_pardiso_solver.cpp | 369 ++++++++++++++++++ tests/CMakeLists.txt | 2 + .../solvers/linear/linear_solver_test.cpp | 113 ++++++ .../linear/mkl_pardiso_solver_test.cpp | 251 ++++++++++++ 8 files changed, 840 insertions(+) create mode 100644 include/fesa/solvers/linear/linear_solver.hpp create mode 100644 include/fesa/solvers/linear/mkl_pardiso_solver.hpp create mode 100644 src/fesa/solvers/linear/mkl_pardiso_solver.cpp create mode 100644 tests/unit/solvers/linear/linear_solver_test.cpp create mode 100644 tests/unit/solvers/linear/mkl_pardiso_solver_test.cpp diff --git a/docs/implementation-plans/linear-static-3d-euler-beam-implementation-report.md b/docs/implementation-plans/linear-static-3d-euler-beam-implementation-report.md index ca81a5f..88f459b 100644 --- a/docs/implementation-plans/linear-static-3d-euler-beam-implementation-report.md +++ b/docs/implementation-plans/linear-static-3d-euler-beam-implementation-report.md @@ -860,3 +860,66 @@ effective RHS without changing this module's mapping responsibility. - concerns: none; no critical implementation or upstream contract conflict was found. + +## Step 20 — mkl-pardiso-solver + +- task_id: `TASK-20` +- status: `completed` +- changed_files: `include/fesa/solvers/linear/linear_solver.hpp`, + `include/fesa/solvers/linear/mkl_pardiso_solver.hpp`, + `src/fesa/solvers/linear/mkl_pardiso_solver.cpp`, + `tests/unit/solvers/linear/linear_solver_test.cpp`, + `tests/unit/solvers/linear/mkl_pardiso_solver_test.cpp`, + `src/fesa/CMakeLists.txt`, `tests/CMakeLists.txt`, + `docs/implementation-plans/linear-static-3d-euler-beam-implementation-report.md`, + `phases/linear-static-3d-euler-beam/index.json` +- requirement_ids: `FESA-REQ-LS3DEB-025`, `FESA-REQ-LS3DEB-026`, + `FESA-REQ-LS3DEB-034`, `FESA-REQ-LS3DEB-035` +- test_ids: `T20-PARDISO-001`, `T20-PARDISO-002`, + `T20-PARDISO-003`, `T20-PARDISO-004`, `T20-PARDISO-005`, + `T20-PARDISO-006` + +| stage | exact command | exit_code | expected_or_observed_result | evidence_tail | +| --- | --- | ---: | --- | --- | +| RED-build | `cmake --build .harness/build --config Debug --target fesa_tests` | 1 | Exactly six planned tests were registered before production and both solver public APIs were absent | CMake regenerated with all dependencies resolved; MSVC C1083 reported missing `fesa/solvers/linear/linear_solver.hpp` and `fesa/solvers/linear/mkl_pardiso_solver.hpp` from the two new test files | +| GREEN-build | `cmake --build .harness/build --config Debug --target fesa_tests` | 0 | Minimum backend boundary, private PARDISO state, six tests, solver library, and unit executable compile and link | `mkl_pardiso_solver.cpp`, both tests, `fesa_solver.lib`, and `fesa_unit_tests.exe` built without a FESA warning under `/W4 /WX` | +| GREEN-test | `ctest --test-dir .harness/build -C Debug -R MklPardisoSolver --output-on-failure` | 0 | SPD solution, repeated RHS, refactorization, validation/failures, and conditioning behavior pass | 6/6 exact `MklPardisoSolver` tests passed | +| VERIFY-configure | `cmake -S . -B .harness/build -A x64 -DFESA_GTEST_SOURCE_DIR=C:/git/googletest "-DMKL_DIR=C:/Program Files (x86)/Intel/oneAPI/mkl/2026.1/lib/cmake/mkl" "-DTBB_DIR=C:/Program Files (x86)/Intel/oneAPI/tbb/2023.1/lib/cmake/tbb" "-DHDF5_DIR=C:/Program Files/HDF_Group/HDF5/2.1.1/cmake"` | 0 | Approved explicit-dependency MSVC x64 build tree generates | Windows SDK 10.0.26100.0, oneMKL 2026.1 ILP64/dynamic, oneTBB, and HDF5 resolved; configure/generate completed | +| VERIFY-build | `cmake --build .harness/build --config Debug` | 0 | Full Debug build passes without a new FESA warning | `fesa_solver.lib` and `fesa_unit_tests.exe` built under `/W4 /WX` | +| VERIFY-targeted | `ctest --test-dir .harness/build -C Debug -R MklPardisoSolver --output-on-failure` | 0 | Focused Step 20 suite remains green | 6/6 exact `MklPardisoSolver` tests passed | +| VERIFY-discovery | `ctest --test-dir .harness/build -C Debug --show-only=json-v1` | 0 | CTest discovers the accumulated suite and all six exact names | 57 tests discovered, including 6 `MklPardisoSolver` tests | +| VERIFY-full | `ctest --test-dir .harness/build -C Debug --output-on-failure` | 0 | Full accumulated C++ suite has zero failures | 57/57 tests passed | +| VERIFY-contract-scans | Fail-on-match public-header scan for `mkl.h`, `MKL_INT`, `pardiso`, oneTBB, and HDF5 symbols; exact-test-count, `git diff --check`, and reference diff/status checks | 0 | The solver backend remains private and Step scope/artifact invariants hold | public PARDISO leaks 0; solver-header backend leaks 0; exact tests 6; whitespace clean; reference unchanged | + +- contract_checks: `LinearSolver` exposes only the exact separate + `factorize(SparseMatrix)` and const `solve(rhs,solution)` boundary. + `MklPardisoSolver` exposes only its constructor, virtual destructor, + overrides, and `unique_ptr`; every MKL header/type/handle is confined + to the `.cpp` and private `Impl`. +- contract_checks: factorization first validates the full square, nonempty, + finite, sorted-unique public CSR, full symmetric structure/value tolerance, + every diagonal slot, and all `size_t` to `MKL_INT` conversions. It then + copies only the sorted upper triangle while retaining exact structural-zero + diagonals. +- contract_checks: the private adapter uses real SPD `mtype=2`, explicit + phases `11`, `22`, `33`, and `-1`, `iparm[26]=1` matrix checking, and + `iparm[34]=1` zero-based indexing. A successful factorization is retained + for repeated RHS; refactorization and destruction release the previous + PARDISO state through the same RAII path. No regularization or fallback is + present. +- contract_checks: solve-before-factorize, RHS/solution dimension mismatch, + nonfinite RHS/solution, nonsquare/empty/asymmetric/missing-diagonal input, + singular or indefinite SPD failure, checker/integer/backend errors, and + phase-specific `-4`/`-7` failures return one deterministic + `FailureCategory::solver` diagnostic. The known SPD system meets normalized + residual `1e-10` and relative analytical error `1e-9`; common-scale and + stiffness-ratio sweeps either meet both limits or return the structured + solver failure without adding a conditioning threshold. +- generated_evidence: `.harness/build/src/fesa/Debug/fesa_solver.lib`, + `.harness/build/tests/Debug/fesa_unit_tests.exe` +- reference_diff: unchanged; `git diff --exit-code -- reference/` exit 0 +- handoff: Step 21 can form the effective free RHS after this adapter has + factorized `Kff`, then call `solve` repeatedly without refactorization or + depending on MKL types. +- concerns: none; no critical implementation, environment, backend, or + upstream-contract conflict was found. diff --git a/include/fesa/solvers/linear/linear_solver.hpp b/include/fesa/solvers/linear/linear_solver.hpp new file mode 100644 index 0000000..cbf7e50 --- /dev/null +++ b/include/fesa/solvers/linear/linear_solver.hpp @@ -0,0 +1,18 @@ +#pragma once + +#include "fesa/core/status.hpp" + +namespace fesa { + +class SparseMatrix; +class Vector; + +// Separates reusable matrix factorization from right-hand-side substitution. +class LinearSolver { +public: + virtual ~LinearSolver() = default; + virtual Status factorize(const SparseMatrix& matrix) = 0; + virtual Status solve(const Vector& rhs, Vector& solution) const = 0; +}; + +} // namespace fesa diff --git a/include/fesa/solvers/linear/mkl_pardiso_solver.hpp b/include/fesa/solvers/linear/mkl_pardiso_solver.hpp new file mode 100644 index 0000000..f0d1eee --- /dev/null +++ b/include/fesa/solvers/linear/mkl_pardiso_solver.hpp @@ -0,0 +1,23 @@ +#pragma once + +#include "fesa/solvers/linear/linear_solver.hpp" + +#include + +namespace fesa { + +// Keeps every oneMKL type and the retained factorization in the private Impl. +class MklPardisoSolver final : public LinearSolver { +public: + MklPardisoSolver(); + ~MklPardisoSolver() override; + + Status factorize(const SparseMatrix& matrix) override; + Status solve(const Vector& rhs, Vector& solution) const override; + +private: + class Impl; + std::unique_ptr impl_; +}; + +} // namespace fesa diff --git a/src/fesa/CMakeLists.txt b/src/fesa/CMakeLists.txt index 29a765c..679f27a 100644 --- a/src/fesa/CMakeLists.txt +++ b/src/fesa/CMakeLists.txt @@ -17,6 +17,7 @@ add_library( math/sparse_matrix.cpp math/vector.cpp model/domain.cpp + solvers/linear/mkl_pardiso_solver.cpp ) target_include_directories( diff --git a/src/fesa/solvers/linear/mkl_pardiso_solver.cpp b/src/fesa/solvers/linear/mkl_pardiso_solver.cpp new file mode 100644 index 0000000..8f5e357 --- /dev/null +++ b/src/fesa/solvers/linear/mkl_pardiso_solver.cpp @@ -0,0 +1,369 @@ +#include "fesa/solvers/linear/mkl_pardiso_solver.hpp" + +#include "fesa/math/sparse_matrix.hpp" +#include "fesa/math/vector.hpp" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace fesa { +namespace { + +Status solverFailure( + const std::string& code, + const std::string& identity, + const std::string& message) { + return Status::failure( + FailureCategory::solver, + {{Severity::error, + code, + {{}, 0U}, + "PARDISO", + identity, + message}}); +} + +Status pardisoFailure( + const MKL_INT phase, + const MKL_INT error) { + std::string code; + std::string reason; + switch (error) { + case -4: + code = "pardiso-zero-or-negative-pivot"; + reason = "zero or negative pivot"; + break; + case -7: + code = "pardiso-singular-diagonal"; + reason = "singular diagonal"; + break; + case -8: + code = "pardiso-integer-overflow"; + reason = "32-bit backend integer overflow"; + break; + case 21: + case 22: + case 23: + case 24: + code = "pardiso-invalid-csr"; + reason = "matrix checker rejected the CSR indices"; + break; + default: + code = phase == 11 ? "pardiso-analysis-failed" : + phase == 22 ? "pardiso-factorization-failed" : + phase == 33 ? "pardiso-solve-failed" : + "pardiso-release-failed"; + reason = "backend error"; + break; + } + + const std::string phaseText = std::to_string(phase); + const std::string errorText = std::to_string(error); + return solverFailure( + code, + "phase=" + phaseText + ",error=" + errorText, + "oneMKL PARDISO phase " + phaseText + " failed with error " + + errorText + " (" + reason + ")."); +} + +bool convertsToMklInt(const std::size_t value) { + return value <= + static_cast((std::numeric_limits::max)()); +} + +} // namespace + +class MklPardisoSolver::Impl { +public: + Impl() = default; + + ~Impl() { + static_cast(release()); + } + + Status factorize(const SparseMatrix& matrix) { + const MKL_INT releaseError = release(); + if (releaseError != 0) { + return pardisoFailure(-1, releaseError); + } + + const Status csrStatus = matrix.validate(); + if (!csrStatus.isOk()) { + return solverFailure( + "solver-invalid-csr", + "public-csr", + "The public sparse matrix failed CSR validation."); + } + if (matrix.rows() != matrix.columns()) { + return solverFailure( + "solver-matrix-not-square", + "matrix-shape", + "PARDISO factorization requires a square matrix."); + } + if (matrix.rows() == 0U) { + return solverFailure( + "solver-empty-matrix", + "matrix-shape", + "PARDISO factorization requires at least one equation."); + } + if (!convertsToMklInt(matrix.rows()) || + !convertsToMklInt(matrix.values().size())) { + return solverFailure( + "solver-dimension-overflow", + "matrix-shape", + "Sparse matrix dimensions exceed the oneMKL integer range."); + } + + const Status copyStatus = copyValidatedUpperTriangle(matrix); + if (!copyStatus.isOk()) { + clearOwnedArrays(); + return copyStatus; + } + + // PARDISO owns internal memory behind pt after phase 11. Initialize + // once per factorization and retain it until refactorization/destruction. + pt_.fill(nullptr); + iparm_.fill(0); + pardisoinit(pt_.data(), &mtype_, iparm_.data()); + iparm_[26] = 1; // Validate sorted CSR integer arrays. + iparm_[34] = 1; // Consume the project's native zero-based CSR. + permutation_.assign(static_cast(equationCount_), 0); + ownsPardisoState_ = true; + + MKL_INT phase = 11; + MKL_INT error = 0; + callPardiso(phase, nullptr, nullptr, error); + if (error != 0) { + const Status failure = pardisoFailure(phase, error); + static_cast(release()); + return failure; + } + + phase = 22; + error = 0; + callPardiso(phase, nullptr, nullptr, error); + if (error != 0) { + const Status failure = pardisoFailure(phase, error); + static_cast(release()); + return failure; + } + + factorized_ = true; + return Status::ok(); + } + + Status solve(const Vector& rhs, Vector& solution) { + if (!factorized_) { + return solverFailure( + "solver-not-factorized", + "factorization-state", + "Substitution requires a successful retained factorization."); + } + const std::size_t size = static_cast(equationCount_); + if (rhs.size() != size || solution.size() != size) { + return solverFailure( + "solver-vector-dimension-mismatch", + "rhs-or-solution", + "RHS and solution dimensions must match the factorized matrix."); + } + for (std::size_t index = 0U; index < rhs.size(); ++index) { + if (!std::isfinite(rhs[index])) { + return solverFailure( + "nonfinite-solver-rhs", + std::to_string(index), + "PARDISO RHS values must be finite."); + } + } + + std::vector rhsCopy(rhs.data(), rhs.data() + rhs.size()); + Vector candidate{size}; + MKL_INT phase = 33; + MKL_INT error = 0; + callPardiso(phase, rhsCopy.data(), candidate.data(), error); + if (error != 0) { + return pardisoFailure(phase, error); + } + for (std::size_t index = 0U; index < candidate.size(); ++index) { + if (!std::isfinite(candidate[index])) { + return solverFailure( + "nonfinite-solver-solution", + std::to_string(index), + "PARDISO substitution produced a nonfinite solution."); + } + } + + solution = std::move(candidate); + return Status::ok(); + } + +private: + Status copyValidatedUpperTriangle(const SparseMatrix& matrix) { + const auto& publicOffsets = matrix.rowOffsets(); + const auto& publicColumns = matrix.columnIndices(); + const auto& publicValues = matrix.values(); + + for (std::size_t row = 0U; row < matrix.rows(); ++row) { + for (std::size_t position = publicOffsets[row]; + position < publicOffsets[row + 1U]; + ++position) { + const std::size_t column = publicColumns[position]; + const auto reverseBegin = publicColumns.begin() + + static_cast(publicOffsets[column]); + const auto reverseEnd = publicColumns.begin() + + static_cast(publicOffsets[column + 1U]); + const auto reverse = + std::lower_bound(reverseBegin, reverseEnd, row); + if (reverse == reverseEnd || *reverse != row) { + return solverFailure( + "solver-matrix-not-symmetric", + std::to_string(row) + ":" + std::to_string(column), + "The full public CSR must contain both symmetric entries."); + } + + const std::size_t reversePosition = static_cast( + std::distance(publicColumns.begin(), reverse)); + const double left = publicValues[position]; + const double right = publicValues[reversePosition]; + const double scale = (std::max)({1.0, std::abs(left), std::abs(right)}); + if (std::abs(left - right) > 1.0e-12 * scale) { + return solverFailure( + "solver-matrix-not-symmetric", + std::to_string(row) + ":" + std::to_string(column), + "The full public CSR values violate the approved symmetry tolerance."); + } + } + } + + equationCount_ = static_cast(matrix.rows()); + rowOffsets_.clear(); + columnIndices_.clear(); + values_.clear(); + rowOffsets_.reserve(matrix.rows() + 1U); + rowOffsets_.push_back(0); + + for (std::size_t row = 0U; row < matrix.rows(); ++row) { + bool hasDiagonal = false; + for (std::size_t position = publicOffsets[row]; + position < publicOffsets[row + 1U]; + ++position) { + const std::size_t column = publicColumns[position]; + if (column < row) { + continue; + } + if (!convertsToMklInt(column) || + !convertsToMklInt(columnIndices_.size())) { + return solverFailure( + "solver-dimension-overflow", + std::to_string(row) + ":" + std::to_string(column), + "CSR indices exceed the oneMKL integer range."); + } + hasDiagonal = hasDiagonal || column == row; + columnIndices_.push_back(static_cast(column)); + values_.push_back(publicValues[position]); + } + if (!hasDiagonal) { + return solverFailure( + "solver-missing-diagonal", + std::to_string(row), + "Every PARDISO SPD row must retain its diagonal slot."); + } + if (!convertsToMklInt(columnIndices_.size())) { + return solverFailure( + "solver-dimension-overflow", + std::to_string(row), + "CSR row offsets exceed the oneMKL integer range."); + } + rowOffsets_.push_back( + static_cast(columnIndices_.size())); + } + return Status::ok(); + } + + void callPardiso( + const MKL_INT phase, + double* rhs, + double* solution, + MKL_INT& error) { + pardiso( + pt_.data(), + &maxFactorizations_, + &matrixNumber_, + &mtype_, + &phase, + &equationCount_, + values_.data(), + rowOffsets_.data(), + columnIndices_.data(), + permutation_.data(), + &rhsCount_, + iparm_.data(), + &messageLevel_, + rhs, + solution, + &error); + } + + MKL_INT release() noexcept { + MKL_INT error = 0; + if (ownsPardisoState_) { + const MKL_INT phase = -1; + double placeholder = 0.0; + callPardiso(phase, &placeholder, &placeholder, error); + } + ownsPardisoState_ = false; + factorized_ = false; + pt_.fill(nullptr); + iparm_.fill(0); + permutation_.clear(); + clearOwnedArrays(); + return error; + } + + void clearOwnedArrays() noexcept { + equationCount_ = 0; + rowOffsets_.clear(); + columnIndices_.clear(); + values_.clear(); + } + + std::array pt_{}; + std::array iparm_{}; + std::vector rowOffsets_; + std::vector columnIndices_; + std::vector permutation_; + std::vector values_; + MKL_INT equationCount_{0}; + MKL_INT maxFactorizations_{1}; + MKL_INT matrixNumber_{1}; + MKL_INT mtype_{2}; + MKL_INT rhsCount_{1}; + MKL_INT messageLevel_{0}; + bool ownsPardisoState_{false}; + bool factorized_{false}; +}; + +MklPardisoSolver::MklPardisoSolver() + : impl_{std::make_unique()} {} + +MklPardisoSolver::~MklPardisoSolver() = default; + +Status MklPardisoSolver::factorize(const SparseMatrix& matrix) { + return impl_->factorize(matrix); +} + +Status MklPardisoSolver::solve( + const Vector& rhs, + Vector& solution) const { + return impl_->solve(rhs, solution); +} + +} // namespace fesa diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 256ab6b..36b02a9 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -22,6 +22,8 @@ add_executable( unit/model/domain_test.cpp unit/model/model_types_test.cpp unit/results/result_records_test.cpp + unit/solvers/linear/linear_solver_test.cpp + unit/solvers/linear/mkl_pardiso_solver_test.cpp ) target_link_libraries( diff --git a/tests/unit/solvers/linear/linear_solver_test.cpp b/tests/unit/solvers/linear/linear_solver_test.cpp new file mode 100644 index 0000000..4c07636 --- /dev/null +++ b/tests/unit/solvers/linear/linear_solver_test.cpp @@ -0,0 +1,113 @@ +#include "fesa/solvers/linear/linear_solver.hpp" +#include "fesa/solvers/linear/mkl_pardiso_solver.hpp" + +#include "fesa/fem/dof_manager.hpp" +#include "fesa/math/sparse_matrix.hpp" + +#include + +#include +#include +#include + +namespace { + +fesa::SparseMatrix makeDenseCsr( + const std::size_t rows, + const std::size_t columns, + const std::vector& values) { + EXPECT_EQ(values.size(), rows * columns); + + fesa::SparsePattern pattern; + std::vector contributions; + pattern.rowOffsets.reserve(rows + 1U); + pattern.rowOffsets.push_back(0U); + for (std::size_t row = 0U; row < rows; ++row) { + for (std::size_t column = 0U; column < columns; ++column) { + pattern.columnIndices.push_back(column); + contributions.push_back({ + row, + column, + values[row * columns + column], + row, + column}); + } + pattern.rowOffsets.push_back(pattern.columnIndices.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.failureCategory(), fesa::FailureCategory::solver); + ASSERT_FALSE(status.diagnostics().empty()); + EXPECT_EQ(status.diagnostics().front().severity, fesa::Severity::error); +} + +} // namespace + +TEST(MklPardisoSolver, RejectsInvalidCsrStateAndDimensions) { + static_assert(std::is_base_of_v); + static_assert(std::has_virtual_destructor_v); + + fesa::MklPardisoSolver solver; + fesa::Vector solution{2U}; + expectSolverFailure(solver.solve(fesa::Vector{2U, 1.0}, solution)); + EXPECT_EQ( + solver.solve(fesa::Vector{2U, 1.0}, solution) + .diagnostics() + .front() + .code, + "solver-not-factorized"); + + const auto rectangular = makeDenseCsr( + 2U, 3U, {2.0, 0.0, 0.0, 0.0, 3.0, 0.0}); + const auto rectangularStatus = solver.factorize(rectangular); + expectSolverFailure(rectangularStatus); + EXPECT_EQ( + rectangularStatus.diagnostics().front().code, + "solver-matrix-not-square"); + + const auto empty = makeDenseCsr(0U, 0U, {}); + const auto emptyStatus = solver.factorize(empty); + expectSolverFailure(emptyStatus); + EXPECT_EQ(emptyStatus.diagnostics().front().code, "solver-empty-matrix"); + + fesa::SparsePattern invalidPattern{{0U, 2U}, {0U}}; + auto invalidCsr = fesa::SparseMatrix::fromCoo( + 1U, + 1U, + {{0U, 0U, 1.0, 0U, 0U}}, + invalidPattern); + EXPECT_FALSE(invalidCsr.hasValue()); + + const auto nonsymmetric = makeDenseCsr(2U, 2U, {2.0, 1.0, 0.0, 3.0}); + const auto nonsymmetricStatus = solver.factorize(nonsymmetric); + expectSolverFailure(nonsymmetricStatus); + EXPECT_EQ( + nonsymmetricStatus.diagnostics().front().code, + "solver-matrix-not-symmetric"); + + fesa::SparsePattern noDiagonalPattern{{0U, 1U, 2U}, {1U, 0U}}; + auto noDiagonal = fesa::SparseMatrix::fromCoo( + 2U, + 2U, + {{0U, 1U, 1.0, 0U, 0U}, {1U, 0U, 1.0, 1U, 0U}}, + noDiagonalPattern); + ASSERT_TRUE(noDiagonal.hasValue()); + const auto noDiagonalStatus = solver.factorize(noDiagonal.value()); + expectSolverFailure(noDiagonalStatus); + EXPECT_EQ( + noDiagonalStatus.diagnostics().front().code, + "solver-missing-diagonal"); + + const auto spd = makeDenseCsr(2U, 2U, {4.0, 1.0, 1.0, 3.0}); + ASSERT_TRUE(solver.factorize(spd).isOk()); + expectSolverFailure(solver.solve(fesa::Vector{1U, 1.0}, solution)); + fesa::Vector wrongSolution{1U}; + expectSolverFailure(solver.solve(fesa::Vector{2U, 1.0}, wrongSolution)); +} diff --git a/tests/unit/solvers/linear/mkl_pardiso_solver_test.cpp b/tests/unit/solvers/linear/mkl_pardiso_solver_test.cpp new file mode 100644 index 0000000..4041d30 --- /dev/null +++ b/tests/unit/solvers/linear/mkl_pardiso_solver_test.cpp @@ -0,0 +1,251 @@ +#include "fesa/solvers/linear/mkl_pardiso_solver.hpp" + +#include "fesa/fem/dof_manager.hpp" +#include "fesa/math/sparse_matrix.hpp" + +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace { + +fesa::SparseMatrix makeDenseCsr( + const std::size_t size, + const std::vector& values) { + EXPECT_EQ(values.size(), size * size); + + fesa::SparsePattern pattern; + std::vector contributions; + pattern.rowOffsets.reserve(size + 1U); + pattern.rowOffsets.push_back(0U); + for (std::size_t row = 0U; row < size; ++row) { + for (std::size_t column = 0U; column < size; ++column) { + pattern.columnIndices.push_back(column); + contributions.push_back({ + row, + column, + values[row * size + column], + row, + column}); + } + pattern.rowOffsets.push_back(pattern.columnIndices.size()); + } + + auto matrix = fesa::SparseMatrix::fromCoo( + size, size, std::move(contributions), pattern); + EXPECT_TRUE(matrix.hasValue()); + return std::move(matrix.value()); +} + +fesa::Vector makeVector(const std::initializer_list values) { + fesa::Vector result{values.size()}; + std::size_t index = 0U; + for (const double value : values) { + result[index++] = value; + } + return result; +} + +double normalizedResidual( + const fesa::SparseMatrix& matrix, + const fesa::Vector& solution, + const fesa::Vector& rhs) { + auto residual = matrix.multiply(solution); + residual.axpy(-1.0, rhs); + return residual.norm() / (std::max)(1.0, rhs.norm()); +} + +double relativeError( + const fesa::Vector& actual, + const fesa::Vector& expected) { + auto difference = actual; + difference.axpy(-1.0, expected); + return difference.norm() / (std::max)(1.0, expected.norm()); +} + +void expectStructuredSolverFailure(const fesa::Status& status) { + EXPECT_FALSE(status.isOk()); + EXPECT_EQ(status.failureCategory(), fesa::FailureCategory::solver); + ASSERT_EQ(status.diagnostics().size(), 1U); + EXPECT_EQ(status.diagnostics()[0U].severity, fesa::Severity::error); + EXPECT_FALSE(status.diagnostics()[0U].code.empty()); + EXPECT_FALSE(status.diagnostics()[0U].message.empty()); +} + +} // namespace + +TEST(MklPardisoSolver, SolvesKnownSpdWithNormalizedResidual) { + const auto matrix = makeDenseCsr(3U, { + 6.0, 2.0, 1.0, + 2.0, 5.0, 2.0, + 1.0, 2.0, 4.0}); + const auto expected = makeVector({1.0, -2.0, 3.0}); + const auto rhs = matrix.multiply(expected); + + fesa::MklPardisoSolver concreteSolver; + fesa::LinearSolver& solver = concreteSolver; + ASSERT_TRUE(solver.factorize(matrix).isOk()); + + fesa::Vector solution{3U}; + ASSERT_TRUE(solver.solve(rhs, solution).isOk()); + EXPECT_LE(normalizedResidual(matrix, solution, rhs), 1.0e-10); + EXPECT_LE(relativeError(solution, expected), 1.0e-9); +} + +TEST(MklPardisoSolver, ReusesOneFactorizationForRepeatedRhs) { + const auto matrix = makeDenseCsr(2U, {4.0, 1.0, 1.0, 3.0}); + const auto expectedFirst = makeVector({1.0, 2.0}); + const auto expectedSecond = makeVector({-2.0, 0.5}); + const auto rhsFirst = matrix.multiply(expectedFirst); + const auto rhsSecond = matrix.multiply(expectedSecond); + + fesa::MklPardisoSolver solver; + ASSERT_TRUE(solver.factorize(matrix).isOk()); + fesa::Vector first{2U}; + fesa::Vector second{2U}; + ASSERT_TRUE(solver.solve(rhsFirst, first).isOk()); + ASSERT_TRUE(solver.solve(rhsSecond, second).isOk()); + + EXPECT_LE(relativeError(first, expectedFirst), 1.0e-9); + EXPECT_LE(relativeError(second, expectedSecond), 1.0e-9); + EXPECT_LE(normalizedResidual(matrix, first, rhsFirst), 1.0e-10); + EXPECT_LE(normalizedResidual(matrix, second, rhsSecond), 1.0e-10); +} + +TEST(MklPardisoSolver, RefactorizesWithoutLeakingState) { + const auto firstMatrix = makeDenseCsr(2U, {4.0, 1.0, 1.0, 3.0}); + const auto secondMatrix = makeDenseCsr(2U, {2.0, 0.0, 0.0, 5.0}); + const auto firstExpected = makeVector({1.0, 2.0}); + const auto secondExpected = makeVector({-3.0, 4.0}); + + fesa::MklPardisoSolver solver; + ASSERT_TRUE(solver.factorize(firstMatrix).isOk()); + fesa::Vector firstSolution{2U}; + ASSERT_TRUE( + solver.solve(firstMatrix.multiply(firstExpected), firstSolution).isOk()); + EXPECT_LE(relativeError(firstSolution, firstExpected), 1.0e-9); + + ASSERT_TRUE(solver.factorize(secondMatrix).isOk()); + fesa::Vector secondSolution{2U}; + const auto secondRhs = secondMatrix.multiply(secondExpected); + ASSERT_TRUE(solver.solve(secondRhs, secondSolution).isOk()); + EXPECT_LE(relativeError(secondSolution, secondExpected), 1.0e-9); + EXPECT_LE( + normalizedResidual(secondMatrix, secondSolution, secondRhs), 1.0e-10); +} + +TEST(MklPardisoSolver, ClassifiesSingularIndefiniteAndNonfiniteFailures) { + fesa::MklPardisoSolver solver; + + const auto singular = makeDenseCsr(2U, {1.0, 1.0, 1.0, 1.0}); + const auto singularStatus = solver.factorize(singular); + expectStructuredSolverFailure(singularStatus); + EXPECT_TRUE( + singularStatus.diagnostics()[0U].code == + "pardiso-zero-or-negative-pivot" || + singularStatus.diagnostics()[0U].code == + "pardiso-singular-diagonal"); + EXPECT_NE( + singularStatus.diagnostics()[0U].entityIdentity.find( + "phase=22,error="), + std::string::npos); + + const auto indefinite = makeDenseCsr(2U, {1.0, 2.0, 2.0, 1.0}); + const auto indefiniteStatus = solver.factorize(indefinite); + expectStructuredSolverFailure(indefiniteStatus); + EXPECT_TRUE( + indefiniteStatus.diagnostics()[0U].code == + "pardiso-zero-or-negative-pivot" || + indefiniteStatus.diagnostics()[0U].code == + "pardiso-singular-diagonal"); + EXPECT_NE( + indefiniteStatus.diagnostics()[0U].entityIdentity.find( + "phase=22,error="), + std::string::npos); + + const auto spd = makeDenseCsr(2U, {3.0, 1.0, 1.0, 2.0}); + ASSERT_TRUE(solver.factorize(spd).isOk()); + auto rhs = makeVector({1.0, 2.0}); + rhs[1U] = (std::numeric_limits::infinity)(); + fesa::Vector solution{2U}; + const auto rhsStatus = solver.solve(rhs, solution); + expectStructuredSolverFailure(rhsStatus); + EXPECT_EQ(rhsStatus.diagnostics()[0U].code, "nonfinite-solver-rhs"); + + fesa::SparsePattern pattern{{0U, 1U}, {0U}}; + auto nonfiniteMatrix = fesa::SparseMatrix::fromCoo( + 1U, + 1U, + {{0U, + 0U, + (std::numeric_limits::quiet_NaN)(), + 0U, + 0U}}, + pattern); + EXPECT_FALSE(nonfiniteMatrix.hasValue()); + EXPECT_EQ( + nonfiniteMatrix.status().diagnostics()[0U].code, + "nonfinite-sparse-value"); +} + +TEST(MklPardisoSolver, ConditioningSweepPassesResolvedCasesAndFailsUnresolvedCasesExplicitly) { + const std::vector commonScales{1.0e-12, 1.0, 1.0e12}; + for (const double scale : commonScales) { + const auto matrix = makeDenseCsr(2U, { + 4.0 * scale, 1.0 * scale, + 1.0 * scale, 3.0 * scale}); + const auto expected = makeVector({1.25, -0.75}); + const auto rhs = matrix.multiply(expected); + + fesa::MklPardisoSolver solver; + ASSERT_TRUE(solver.factorize(matrix).isOk()) << "scale=" << scale; + fesa::Vector solution{2U}; + ASSERT_TRUE(solver.solve(rhs, solution).isOk()) << "scale=" << scale; + EXPECT_LE(normalizedResidual(matrix, solution, rhs), 1.0e-10); + EXPECT_LE(relativeError(solution, expected), 1.0e-9); + } + + const double resolvedRatio = 1.0e-8; + const auto resolvedMatrix = makeDenseCsr( + 2U, {1.0, 0.0, 0.0, resolvedRatio}); + const auto resolvedExpected = makeVector({0.5, -2.0}); + const auto resolvedRhs = resolvedMatrix.multiply(resolvedExpected); + fesa::MklPardisoSolver resolvedSolver; + ASSERT_TRUE(resolvedSolver.factorize(resolvedMatrix).isOk()); + fesa::Vector resolvedSolution{2U}; + ASSERT_TRUE(resolvedSolver.solve(resolvedRhs, resolvedSolution).isOk()); + EXPECT_LE( + normalizedResidual(resolvedMatrix, resolvedSolution, resolvedRhs), + 1.0e-10); + EXPECT_LE(relativeError(resolvedSolution, resolvedExpected), 1.0e-9); + + const std::vector unresolvedCandidates{1.0e-16, 1.0e-300}; + for (const double ratio : unresolvedCandidates) { + const auto matrix = makeDenseCsr(2U, {1.0, 0.0, 0.0, ratio}); + const auto expected = makeVector({0.5, -2.0}); + const auto rhs = matrix.multiply(expected); + fesa::MklPardisoSolver solver; + + const auto factorStatus = solver.factorize(matrix); + if (!factorStatus.isOk()) { + expectStructuredSolverFailure(factorStatus); + continue; + } + + fesa::Vector solution{2U}; + const auto solveStatus = solver.solve(rhs, solution); + if (!solveStatus.isOk()) { + expectStructuredSolverFailure(solveStatus); + continue; + } + + EXPECT_LE(normalizedResidual(matrix, solution, rhs), 1.0e-10); + EXPECT_LE(relativeError(solution, expected), 1.0e-9); + } +}