feat(nonlinear-static-newton-raphson): step 6 - backend-neutral-tangent-solver

include/fesa/solvers/linear/mkl_pardiso_tangent_solver.h

@@ -0,0 +1,32 @@
+#ifndef FESA_SOLVERS_LINEAR_MKL_PARDISO_TANGENT_SOLVER_H_
+#define FESA_SOLVERS_LINEAR_MKL_PARDISO_TANGENT_SOLVER_H_
+
+#include <memory>
+
+#include "fesa/solvers/linear/tangent_solver.h"
+
+namespace fesa {
+
+/// @brief Adapts full-CSR unsymmetric oneMKL PARDISO tangent solves.
+class MklPardisoTangentSolver final : public TangentSolver {
+ public:
+  /// @brief Constructs an empty, unfactorized PARDISO tangent adapter.
+  MklPardisoTangentSolver();
+
+  /// @brief Releases retained PARDISO tangent backend state.
+  ~MklPardisoTangentSolver() override;
+
+  /// @brief Validates and factorizes one fresh nonlinear tangent matrix.
+  Status FactorizeFresh(const SparseMatrix& matrix) override;
+
+  /// @brief Solves one correction without refactorizing the tangent.
+  Status SolveCorrection(const Vector& rhs, Vector& correction) const override;
+
+ private:
+  class Impl;
+  std::unique_ptr<Impl> impl_;
+};
+
+}  // namespace fesa
+
+#endif  // FESA_SOLVERS_LINEAR_MKL_PARDISO_TANGENT_SOLVER_H_

include/fesa/solvers/linear/tangent_solver.h

@@ -0,0 +1,31 @@
+#ifndef FESA_SOLVERS_LINEAR_TANGENT_SOLVER_H_
+#define FESA_SOLVERS_LINEAR_TANGENT_SOLVER_H_
+
+#include "fesa/core/status.h"
+
+namespace fesa {
+
+class SparseMatrix;
+class Vector;
+
+/// @brief Factorizes one nonlinear tangent and solves its correction.
+class TangentSolver {
+ public:
+  /// @brief Destroys a backend-neutral nonlinear tangent solver.
+  virtual ~TangentSolver() = default;
+
+  /// @brief Replaces any retained state with a fresh tangent factorization.
+  /// @return Success or a structured solver failure.
+  virtual Status FactorizeFresh(const SparseMatrix& matrix) = 0;
+
+  /// @brief Solves one correction using the current tangent factorization.
+  /// @param rhs Immutable free-DOF residual right-hand side.
+  /// @param correction Updated only after successful finite substitution.
+  /// @return Success or a structured solver failure.
+  virtual Status SolveCorrection(const Vector& rhs, Vector& correction) const =
+      0;
+};
+
+}  // namespace fesa
+
+#endif  // FESA_SOLVERS_LINEAR_TANGENT_SOLVER_H_

src/fesa/CMakeLists.txt

@@ -53,6 +53,7 @@ add_library(
     results/result_recovery.cpp
     results/shell_result_recovery.cpp
     solvers/linear/mkl_pardiso_solver.cpp
+    solvers/linear/mkl_pardiso_tangent_solver.cpp
 )
 
 target_include_directories(

src/fesa/solvers/linear/mkl_pardiso_tangent_solver.cpp

@@ -0,0 +1,294 @@
+#include "fesa/solvers/linear/mkl_pardiso_tangent_solver.h"
+
+#include <mkl.h>
+
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <limits>
+#include <memory>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include "fesa/math/sparse_matrix.h"
+#include "fesa/math/vector.h"
+
+namespace fesa {
+namespace {
+
+/// @brief Builds a structured nonlinear tangent-solver failure.
+Status SolverFailure(const std::string& code, const std::string& identity,
+                     const std::string& message) {
+  return Status::Failure(
+      FailureCategory::kSolver,
+      {{Severity::kError, code, {{}, 0U}, "PARDISO", identity, message}});
+}
+
+/// @brief Translates a PARDISO phase error into the stable solver taxonomy.
+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 phase_text = std::to_string(phase);
+  const std::string error_text = std::to_string(error);
+  return SolverFailure(code, "phase=" + phase_text + ",error=" + error_text,
+                       "oneMKL PARDISO phase " + phase_text +
+                           " failed with error " + error_text + " (" + reason +
+                           ").");
+}
+
+/// @brief Reports whether a public size fits the private MKL integer type.
+bool ConvertsToMklInt(const std::size_t value) {
+  return value <=
+         static_cast<std::size_t>((std::numeric_limits<MKL_INT>::max)());
+}
+
+}  // namespace
+
+class MklPardisoTangentSolver::Impl {
+ public:
+  Impl() = default;
+
+  ~Impl() { static_cast<void>(Release()); }
+
+  /// @brief Replaces any retained state with one validated full-CSR tangent.
+  Status FactorizeFresh(const SparseMatrix& matrix) {
+    const MKL_INT release_error = Release();
+    if (release_error != 0) {
+      return PardisoFailure(-1, release_error);
+    }
+
+    const Status csr_status = matrix.Validate();
+    if (!csr_status.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) {
+      // PARDISO does not accept n=0, but a fully constrained nonlinear
+      // correction remains a valid factorize/solve lifecycle.
+      factorized_ = true;
+      return Status::Ok();
+    }
+    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 copy_status = CopyValidatedFullCsr(matrix);
+    if (!copy_status.IsOk()) {
+      ClearOwnedArrays();
+      return copy_status;
+    }
+
+    // The nonlinear tangent is generally nonsymmetric, so PARDISO consumes
+    // the complete public CSR with mtype=11 instead of an SPD triangle.
+    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<std::size_t>(equation_count_), 0);
+    owns_pardiso_state_ = 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<void>(Release());
+      return failure;
+    }
+
+    phase = 22;
+    error = 0;
+    CallPardiso(phase, nullptr, nullptr, error);
+    if (error != 0) {
+      const Status failure = PardisoFailure(phase, error);
+      static_cast<void>(Release());
+      return failure;
+    }
+    if (iparm_[13] != 0) {
+      const Status failure = SolverFailure(
+          "pardiso-perturbed-pivot", "phase=22",
+          "oneMKL PARDISO perturbed one or more tangent pivots.");
+      static_cast<void>(Release());
+      return failure;
+    }
+
+    factorized_ = true;
+    return Status::Ok();
+  }
+
+  /// @brief Solves one correction while preserving caller output on failure.
+  Status SolveCorrection(const Vector& rhs, Vector& correction) {
+    if (!factorized_) {
+      return SolverFailure(
+          "solver-not-factorized", "factorization-state",
+          "Substitution requires a successful retained factorization.");
+    }
+    const std::size_t size = static_cast<std::size_t>(equation_count_);
+    if (rhs.Size() != size || correction.Size() != size) {
+      return SolverFailure(
+          "solver-vector-dimension-mismatch", "rhs-or-correction",
+          "RHS and correction 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.");
+      }
+    }
+    if (size == 0U) {
+      return Status::Ok();
+    }
+
+    std::vector<double> rhs_copy(rhs.Data(), rhs.Data() + rhs.Size());
+    Vector candidate{size};
+    MKL_INT phase = 33;
+    MKL_INT error = 0;
+    CallPardiso(phase, rhs_copy.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.");
+      }
+    }
+
+    correction = std::move(candidate);
+    return Status::Ok();
+  }
+
+ private:
+  /// @brief Copies validated public CSR without symmetry or diagonal filtering.
+  Status CopyValidatedFullCsr(const SparseMatrix& matrix) {
+    equation_count_ = static_cast<MKL_INT>(matrix.Rows());
+    row_offsets_.clear();
+    column_indices_.clear();
+    values_.clear();
+    row_offsets_.reserve(matrix.RowOffsets().size());
+    column_indices_.reserve(matrix.ColumnIndices().size());
+    values_.reserve(matrix.Values().size());
+
+    for (const std::size_t offset : matrix.RowOffsets()) {
+      if (!ConvertsToMklInt(offset)) {
+        return SolverFailure("solver-dimension-overflow", "row-offsets",
+                             "CSR row offsets exceed the oneMKL integer "
+                             "range.");
+      }
+      row_offsets_.push_back(static_cast<MKL_INT>(offset));
+    }
+    for (const std::size_t column : matrix.ColumnIndices()) {
+      if (!ConvertsToMklInt(column)) {
+        return SolverFailure("solver-dimension-overflow", "column-indices",
+                             "CSR indices exceed the oneMKL integer range.");
+      }
+      column_indices_.push_back(static_cast<MKL_INT>(column));
+    }
+    values_ = matrix.Values();
+    return Status::Ok();
+  }
+
+  /// @brief Calls PARDISO with retained private arrays and phase state.
+  void CallPardiso(const MKL_INT phase, double* rhs, double* correction,
+                   MKL_INT& error) {
+    pardiso(pt_.data(), &max_factorizations_, &matrix_number_, &mtype_, &phase,
+            &equation_count_, values_.data(), row_offsets_.data(),
+            column_indices_.data(), permutation_.data(), &rhs_count_,
+            iparm_.data(), &message_level_, rhs, correction, &error);
+  }
+
+  /// @brief Releases backend state and resets every retained array.
+  MKL_INT Release() noexcept {
+    MKL_INT error = 0;
+    if (owns_pardiso_state_) {
+      const MKL_INT phase = -1;
+      double placeholder = 0.0;
+      CallPardiso(phase, &placeholder, &placeholder, error);
+    }
+    owns_pardiso_state_ = false;
+    factorized_ = false;
+    pt_.fill(nullptr);
+    iparm_.fill(0);
+    permutation_.clear();
+    ClearOwnedArrays();
+    return error;
+  }
+
+  /// @brief Clears FESA-owned CSR arrays without touching backend state.
+  void ClearOwnedArrays() noexcept {
+    equation_count_ = 0;
+    row_offsets_.clear();
+    column_indices_.clear();
+    values_.clear();
+  }
+
+  std::array<void*, 64U> pt_{};
+  std::array<MKL_INT, 64U> iparm_{};
+  std::vector<MKL_INT> row_offsets_;
+  std::vector<MKL_INT> column_indices_;
+  std::vector<MKL_INT> permutation_;
+  std::vector<double> values_;
+  MKL_INT equation_count_{0};
+  MKL_INT max_factorizations_{1};
+  MKL_INT matrix_number_{1};
+  MKL_INT mtype_{11};
+  MKL_INT rhs_count_{1};
+  MKL_INT message_level_{0};
+  bool owns_pardiso_state_{false};
+  bool factorized_{false};
+};
+
+MklPardisoTangentSolver::MklPardisoTangentSolver()
+    : impl_{std::make_unique<Impl>()} {}
+
+MklPardisoTangentSolver::~MklPardisoTangentSolver() = default;
+
+Status MklPardisoTangentSolver::FactorizeFresh(const SparseMatrix& matrix) {
+  return impl_->FactorizeFresh(matrix);
+}
+
+Status MklPardisoTangentSolver::SolveCorrection(const Vector& rhs,
+                                                 Vector& correction) const {
+  return impl_->SolveCorrection(rhs, correction);
+}
+
+}  // namespace fesa

tests/CMakeLists.txt

@@ -48,6 +48,7 @@ add_executable(
     unit/results/results_writer_test.cpp
     unit/solvers/linear/linear_solver_test.cpp
     unit/solvers/linear/mkl_pardiso_solver_test.cpp
+    unit/solvers/linear/tangent_solver_test.cpp
 )
 
 target_include_directories(

tests/unit/solvers/linear/mkl_pardiso_solver_test.cpp

@@ -123,6 +123,15 @@ TEST(MklPardisoSolver, ReusesOneFactorizationForRepeatedRhs) {
   EXPECT_LE(NormalizedResidual(matrix, second, rhs_second), 1.0e-10);
 }
 
+TEST(MklPardisoSolver, StillRejectsNonsymmetricPublicCsr) {
+  const auto nonsymmetric = MakeDenseCsr(2U, {2.0, 1.0, 0.0, 3.0});
+  fesa::MklPardisoSolver solver;
+
+  const auto status = solver.Factorize(nonsymmetric);
+  ExpectStructuredSolverFailure(status);
+  EXPECT_EQ(status.Diagnostics()[0U].code, "solver-matrix-not-symmetric");
+}
+
 TEST(MklPardisoSolver, RefactorizesWithoutLeakingState) {
   const auto first_matrix = MakeDenseCsr(2U, {4.0, 1.0, 1.0, 3.0});
   const auto second_matrix = MakeDenseCsr(2U, {2.0, 0.0, 0.0, 5.0});

tests/unit/solvers/linear/tangent_solver_test.cpp

@@ -0,0 +1,132 @@
+#include "fesa/solvers/linear/mkl_pardiso_tangent_solver.h"
+
+#include <gtest/gtest.h>
+
+#include <cstddef>
+#include <initializer_list>
+#include <limits>
+#include <type_traits>
+#include <utility>
+#include <vector>
+
+#include "fesa/fem/dof_manager.h"
+#include "fesa/math/sparse_matrix.h"
+#include "fesa/math/vector.h"
+#include "fesa/solvers/linear/tangent_solver.h"
+
+namespace {
+
+fesa::SparseMatrix MakeDenseCsr(const std::size_t size,
+                                const std::vector<double>& values) {
+  EXPECT_EQ(values.size(), size * size);
+
+  fesa::SparsePattern pattern;
+  std::vector<fesa::CooContribution> contributions;
+  pattern.row_offsets.reserve(size + 1U);
+  pattern.row_offsets.push_back(0U);
+  for (std::size_t row = 0U; row < size; ++row) {
+    for (std::size_t column = 0U; column < size; ++column) {
+      pattern.column_indices.push_back(column);
+      contributions.push_back(
+          {row, column, values[row * size + column], row, column});
+    }
+    pattern.row_offsets.push_back(pattern.column_indices.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<double> values) {
+  fesa::Vector result{values.size()};
+  std::size_t index = 0U;
+  for (const double value : values) {
+    result[index++] = value;
+  }
+  return result;
+}
+
+void ExpectSolverFailure(const fesa::Status& status) {
+  EXPECT_FALSE(status.IsOk());
+  EXPECT_EQ(status.Category(), fesa::FailureCategory::kSolver);
+  ASSERT_EQ(status.Diagnostics().size(), 1U);
+  EXPECT_EQ(status.Diagnostics().front().severity, fesa::Severity::kError);
+}
+
+TEST(TangentSolver, SolvesNonsymmetricKnownCorrection) {
+  static_assert(
+      std::is_base_of_v<fesa::TangentSolver, fesa::MklPardisoTangentSolver>);
+  static_assert(std::has_virtual_destructor_v<fesa::TangentSolver>);
+
+  const auto matrix = MakeDenseCsr(2U, {2.0, 1.0, 0.0, 3.0});
+  const auto rhs = MakeVector({4.0, 6.0});
+  fesa::MklPardisoTangentSolver concrete_solver;
+  fesa::TangentSolver& solver = concrete_solver;
+
+  ASSERT_TRUE(solver.FactorizeFresh(matrix).IsOk());
+  fesa::Vector correction{2U};
+  ASSERT_TRUE(solver.SolveCorrection(rhs, correction).IsOk());
+  EXPECT_NEAR(correction[0U], 1.0, 1.0e-14);
+  EXPECT_NEAR(correction[1U], 2.0, 1.0e-14);
+}
+
+TEST(TangentSolver, SolvesTrivialZeroEquationSystem) {
+  const auto matrix = MakeDenseCsr(0U, {});
+  fesa::MklPardisoTangentSolver solver;
+
+  ASSERT_TRUE(solver.FactorizeFresh(matrix).IsOk());
+  fesa::Vector correction{0U};
+  EXPECT_TRUE(solver.SolveCorrection(fesa::Vector{0U}, correction).IsOk());
+  EXPECT_EQ(correction.Size(), 0U);
+}
+
+TEST(TangentSolver, RejectsNonfiniteInputsWithoutModifyingCorrection) {
+  const fesa::SparsePattern pattern{{0U, 1U}, {0U}};
+  const auto nonfinite_matrix = fesa::SparseMatrix::FromCoo(
+      1U, 1U,
+      {{0U, 0U, (std::numeric_limits<double>::quiet_NaN)(), 0U, 0U}},
+      pattern);
+  EXPECT_FALSE(nonfinite_matrix.HasValue());
+  EXPECT_EQ(nonfinite_matrix.GetStatus().Diagnostics().front().code,
+            "nonfinite-sparse-value");
+
+  fesa::MklPardisoTangentSolver solver;
+  ASSERT_TRUE(solver.FactorizeFresh(MakeDenseCsr(2U, {2.0, 1.0, 0.0, 3.0}))
+                  .IsOk());
+  auto rhs = MakeVector({4.0, 6.0});
+  rhs[1U] = (std::numeric_limits<double>::infinity)();
+  fesa::Vector correction = MakeVector({17.0, -4.0});
+
+  const auto status = solver.SolveCorrection(rhs, correction);
+  ExpectSolverFailure(status);
+  EXPECT_EQ(status.Diagnostics().front().code, "nonfinite-solver-rhs");
+  EXPECT_DOUBLE_EQ(correction[0U], 17.0);
+  EXPECT_DOUBLE_EQ(correction[1U], -4.0);
+}
+
+TEST(TangentSolver, ReportsSingularFactorizationDeterministically) {
+  fesa::MklPardisoTangentSolver solver;
+
+  const auto status =
+      solver.FactorizeFresh(MakeDenseCsr(2U, {0.0, 0.0, 0.0, 0.0}));
+  ASSERT_FALSE(status.IsOk());
+  ExpectSolverFailure(status);
+  EXPECT_EQ(status.Diagnostics().front().code, "pardiso-perturbed-pivot");
+  EXPECT_EQ(status.Diagnostics().front().entity_identity, "phase=22");
+}
+
+TEST(TangentSolver, RejectsSolveBeforeSuccessfulFactorization) {
+  fesa::MklPardisoTangentSolver solver;
+  fesa::Vector correction = MakeVector({23.0, -9.0});
+
+  const auto status =
+      solver.SolveCorrection(MakeVector({1.0, 2.0}), correction);
+  ExpectSolverFailure(status);
+  EXPECT_EQ(status.Diagnostics().front().code, "solver-not-factorized");
+  EXPECT_DOUBLE_EQ(correction[0U], 23.0);
+  EXPECT_DOUBLE_EQ(correction[1U], -9.0);
+}
+
+}  // namespace