From be5f4eb86df2b883c11912b471a69c0b2c2c79e7 Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Sun, 9 Aug 2026 20:20:37 +0900 Subject: [PATCH] feat(linear-static-3d-euler-beam): step 20 - mkl-pardiso-solver-review-fix --- ...tic-3d-euler-beam-implementation-report.md | 46 +++++++++++++ .../solvers/linear/mkl_pardiso_solver.cpp | 26 ++++++-- .../solvers/linear/linear_solver_test.cpp | 65 ++++++++++++++----- .../linear/mkl_pardiso_solver_test.cpp | 27 +++++++- 4 files changed, 139 insertions(+), 25 deletions(-) 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 88f459b..1a8202e 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 @@ -923,3 +923,49 @@ depending on MKL types. - concerns: none; no critical implementation, environment, backend, or upstream-contract conflict was found. + +### Step 20 Review Fix Round 1 — fully constrained lifecycle and normalized validation + +- review_classification: Important integration/correctness fixes. Step 19 + intentionally produces a valid `0x0 Kff` for an all-constrained model, so + rejecting that matrix contradicted the accumulated solver workflow. +- test_scope: the existing six `MklPardisoSolver` test names were retained. + The invalid-state test now covers trivial factorize/solve followed by a + nonempty refactor, unchanged caller output on failed solves, and a + small-scale asymmetric matrix. Numerical helpers now divide by the actual + nonzero RHS/expected norm and fail closed for zero or nonfinite norms. + +| stage | exact command | exit_code | observed_result | +| --- | --- | ---: | --- | +| RED-zero-equation | `cmake --build .harness/build --config Debug; ctest --test-dir .harness/build -C Debug --output-on-failure -R MklPardisoSolver` | 1 | Build passed; 5/6 tests passed and `RejectsInvalidCsrStateAndDimensions` failed because valid `0x0 Kff` factorization returned failure | +| TEST-hardening | same targeted command after replacing unit-floor norm helpers | 1 | The same zero-equation production defect remained the only failure; no artificial RED was claimed for already-correct numerical results | +| GREEN-zero-equation | same targeted command after the minimum production change | 0 | 6/6 passed; empty factorize/solve avoids PARDISO and a following nonempty refactor solves correctly | +| RED-scaled-symmetry | `cmake --build .harness/build --config Debug; ctest --test-dir .harness/build -C Debug --output-on-failure -R MklPardisoSolver` | 1 | Build passed; 5/6 tests passed and the existing invalid-state test showed that a scaled asymmetric matrix was incorrectly accepted | +| GREEN-scaled-symmetry | same targeted command after the minimum production change | 0 | 6/6 passed with symmetry normalized by the matrix's actual maximum absolute entry and no unit-size floor | +| 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 | Explicit MSVC x64 dependencies configured and generated successfully | +| VERIFY-build/targeted | `cmake --build .harness/build --config Debug; ctest --test-dir .harness/build -C Debug --output-on-failure -R MklPardisoSolver` | 0 | Debug build completed without a new warning; exact Step 20 suite passed 6/6 | +| VERIFY-discovery/full | `ctest --test-dir .harness/build -C Debug --show-only=json-v1; ctest --test-dir .harness/build -C Debug --output-on-failure` | 0 | 57 tests discovered and 57/57 passed | +| VERIFY-contract | backend API leak scan, exact-test count, forbidden fallback/regularization scan, `git diff --check`, and reference diff/status | 0 | backend API leaks 0; exact tests 6; forbidden implementation matches 0; whitespace clean; `reference/` unchanged | + +- production_fix: a successful zero-equation factorization is represented as + trivial private state without allocating or calling PARDISO; empty solve + succeeds only after that observable factorization. The ordinary `release()` + path safely clears this state before later refactorization or destruction. +- production_fix: symmetry requires every reverse structural pair and compares + `abs(aij-aji)` to `1e-12 * max(abs(A))`. An exactly zero matrix uses an exact + zero difference check; SparseMatrix validation remains responsible for + rejecting nonfinite values. +- failure_atomicity: solve-before-factorize, dimension mismatch, and nonfinite + RHS failures are asserted to leave the caller's solution values unchanged; + successful backend output is still copied from a private candidate only + after phase 33 and finite-result validation. +- evidence_hygiene: two intermediate symmetry RED attempts were discarded + because a concurrent clean build locked MSBuild files. The uncontended RED + above reproduced the production defect deterministically and is the recorded + evidence. +- phase_index: intentionally unchanged during this review fix; hash + `2cfb6ee5e0c49fbf6b06d27cb33ec9ac2c34dbc4`. +- supersession: this review section supersedes the original Step 20 statements + that the factorized CSR must be nonempty or that empty input is rejected. +- concerns: none; the transient build contention was resolved before recorded + RED/GREEN/VERIFY runs and no critical or upstream blocker remains. diff --git a/src/fesa/solvers/linear/mkl_pardiso_solver.cpp b/src/fesa/solvers/linear/mkl_pardiso_solver.cpp index 8f5e357..8b2fafa 100644 --- a/src/fesa/solvers/linear/mkl_pardiso_solver.cpp +++ b/src/fesa/solvers/linear/mkl_pardiso_solver.cpp @@ -110,10 +110,11 @@ public: "PARDISO factorization requires a square matrix."); } if (matrix.rows() == 0U) { - return solverFailure( - "solver-empty-matrix", - "matrix-shape", - "PARDISO factorization requires at least one equation."); + // A fully constrained model has no free equations. Preserve the + // observable factorize/solve lifecycle without creating backend + // state or calling PARDISO with its invalid n=0 input. + factorized_ = true; + return Status::ok(); } if (!convertsToMklInt(matrix.rows()) || !convertsToMklInt(matrix.values().size())) { @@ -183,6 +184,9 @@ public: "PARDISO RHS values must be finite."); } } + if (size == 0U) { + return Status::ok(); + } std::vector rhsCopy(rhs.data(), rhs.data() + rhs.size()); Vector candidate{size}; @@ -211,6 +215,11 @@ private: const auto& publicColumns = matrix.columnIndices(); const auto& publicValues = matrix.values(); + double matrixScale = 0.0; + for (const double value : publicValues) { + matrixScale = (std::max)(matrixScale, std::abs(value)); + } + for (std::size_t row = 0U; row < matrix.rows(); ++row) { for (std::size_t position = publicOffsets[row]; position < publicOffsets[row + 1U]; @@ -233,8 +242,13 @@ private: 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) { + const double difference = std::abs(left - right); + // The approved symmetry test is normalized by the matrix's + // actual nonzero scale, without an absolute unit-size floor. + const bool isSymmetric = matrixScale == 0.0 ? + difference == 0.0 : + difference <= 1.0e-12 * matrixScale; + if (!isSymmetric) { return solverFailure( "solver-matrix-not-symmetric", std::to_string(row) + ":" + std::to_string(column), diff --git a/tests/unit/solvers/linear/linear_solver_test.cpp b/tests/unit/solvers/linear/linear_solver_test.cpp index 4c07636..cb4b976 100644 --- a/tests/unit/solvers/linear/linear_solver_test.cpp +++ b/tests/unit/solvers/linear/linear_solver_test.cpp @@ -55,14 +55,45 @@ TEST(MklPardisoSolver, RejectsInvalidCsrStateAndDimensions) { static_assert(std::has_virtual_destructor_v); fesa::MklPardisoSolver solver; - fesa::Vector solution{2U}; - expectSolverFailure(solver.solve(fesa::Vector{2U, 1.0}, solution)); + fesa::Vector untouched{2U}; + untouched[0U] = 17.0; + untouched[1U] = -4.0; + const auto beforeFactorize = + solver.solve(fesa::Vector{2U, 1.0}, untouched); + expectSolverFailure(beforeFactorize); EXPECT_EQ( - solver.solve(fesa::Vector{2U, 1.0}, solution) - .diagnostics() - .front() - .code, + beforeFactorize.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 emptySolution{0U}; + EXPECT_TRUE(solver.solve(fesa::Vector{0U}, emptySolution).isOk()); + EXPECT_EQ(emptySolution.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 solvedFirst = solution[0U]; + const double solvedSecond = solution[1U]; + expectSolverFailure(solver.solve(fesa::Vector{1U, 1.0}, solution)); + EXPECT_DOUBLE_EQ(solution[0U], solvedFirst); + EXPECT_DOUBLE_EQ(solution[1U], solvedSecond); + + fesa::Vector wrongSolution{1U}; + wrongSolution[0U] = 41.0; + expectSolverFailure(solver.solve(fesa::Vector{2U, 1.0}, wrongSolution)); + EXPECT_DOUBLE_EQ(wrongSolution[0U], 41.0); const auto rectangular = makeDenseCsr( 2U, 3U, {2.0, 0.0, 0.0, 0.0, 3.0, 0.0}); @@ -72,11 +103,6 @@ TEST(MklPardisoSolver, RejectsInvalidCsrStateAndDimensions) { 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, @@ -92,6 +118,18 @@ TEST(MklPardisoSolver, RejectsInvalidCsrStateAndDimensions) { nonsymmetricStatus.diagnostics().front().code, "solver-matrix-not-symmetric"); + const auto scaledNonsymmetric = makeDenseCsr( + 2U, 2U, {2.0e-20, 1.0e-20, 1.1e-20, 3.0e-20}); + const auto scaledNonsymmetricStatus = + solver.factorize(scaledNonsymmetric); + // Stop this case before inspecting diagnostics when the production code + // incorrectly accepts the matrix; this keeps the RED failure deterministic. + ASSERT_FALSE(scaledNonsymmetricStatus.isOk()); + expectSolverFailure(scaledNonsymmetricStatus); + EXPECT_EQ( + scaledNonsymmetricStatus.diagnostics().front().code, + "solver-matrix-not-symmetric"); + fesa::SparsePattern noDiagonalPattern{{0U, 1U, 2U}, {1U, 0U}}; auto noDiagonal = fesa::SparseMatrix::fromCoo( 2U, @@ -105,9 +143,4 @@ TEST(MklPardisoSolver, RejectsInvalidCsrStateAndDimensions) { 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 index 4041d30..28c84c0 100644 --- a/tests/unit/solvers/linear/mkl_pardiso_solver_test.cpp +++ b/tests/unit/solvers/linear/mkl_pardiso_solver_test.cpp @@ -5,7 +5,6 @@ #include -#include #include #include #include @@ -58,7 +57,16 @@ double normalizedResidual( const fesa::Vector& rhs) { auto residual = matrix.multiply(solution); residual.axpy(-1.0, rhs); - return residual.norm() / (std::max)(1.0, rhs.norm()); + const double numerator = residual.norm(); + const double denominator = rhs.norm(); + if (!std::isfinite(numerator) || !std::isfinite(denominator)) { + return (std::numeric_limits::infinity)(); + } + if (denominator == 0.0) { + return numerator == 0.0 ? 0.0 : + (std::numeric_limits::infinity)(); + } + return numerator / denominator; } double relativeError( @@ -66,7 +74,16 @@ double relativeError( const fesa::Vector& expected) { auto difference = actual; difference.axpy(-1.0, expected); - return difference.norm() / (std::max)(1.0, expected.norm()); + const double numerator = difference.norm(); + const double denominator = expected.norm(); + if (!std::isfinite(numerator) || !std::isfinite(denominator)) { + return (std::numeric_limits::infinity)(); + } + if (denominator == 0.0) { + return numerator == 0.0 ? 0.0 : + (std::numeric_limits::infinity)(); + } + return numerator / denominator; } void expectStructuredSolverFailure(const fesa::Status& status) { @@ -174,9 +191,13 @@ TEST(MklPardisoSolver, ClassifiesSingularIndefiniteAndNonfiniteFailures) { auto rhs = makeVector({1.0, 2.0}); rhs[1U] = (std::numeric_limits::infinity)(); fesa::Vector solution{2U}; + solution[0U] = 23.0; + solution[1U] = -9.0; const auto rhsStatus = solver.solve(rhs, solution); expectStructuredSolverFailure(rhsStatus); EXPECT_EQ(rhsStatus.diagnostics()[0U].code, "nonfinite-solver-rhs"); + EXPECT_DOUBLE_EQ(solution[0U], 23.0); + EXPECT_DOUBLE_EQ(solution[1U], -9.0); fesa::SparsePattern pattern{{0U, 1U}, {0U}}; auto nonfiniteMatrix = fesa::SparseMatrix::fromCoo(