FESADev/src/fesa/solvers/linear/mkl_pardiso_tangent_solver.cpp

#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"
#include "math/dense_blas_internal.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) {
    const dense_blas_internal::ScopedSingleEquationThreads thread_scope{
        static_cast<std::size_t>(equation_count_)};
    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