Integrate Refactor solver sources and dual MSVC CMake builds

src/fesa/math/compensated_sum.h

@@ -0,0 +1,30 @@
+#ifndef FESA_MATH_COMPENSATED_SUM_H_
+#define FESA_MATH_COMPENSATED_SUM_H_
+
+#include <cmath>
+
+namespace fesa::math_internal {
+
+/// @brief Preserves the compensated scalar accumulation graph and caller order.
+inline void AddCompensated(const double term, double* sum, double* correction) {
+  const double next = *sum + term;
+  *correction += std::abs(*sum) >= std::abs(term)
+                     ? (*sum - next) + term
+                     : (term - next) + *sum;
+  *sum = next;
+}
+
+/// @brief Owns one sum and its correction without changing reduction order.
+class CompensatedSum {
+ public:
+  void Add(const double term) { AddCompensated(term, &sum_, &correction_); }
+  double Value() const { return sum_ + correction_; }
+
+ private:
+  double sum_ = 0.0;
+  double correction_ = 0.0;
+};
+
+}  // namespace fesa::math_internal
+
+#endif  // FESA_MATH_COMPENSATED_SUM_H_

src/fesa/math/dense_blas_internal.cpp

@@ -7,8 +7,271 @@
 #include <utility>
 #include <vector>
 
+#include "math/finite_values.h"
+
 namespace fesa::dense_blas_internal {
 
+namespace {
+
+using math_internal::IsFinite;
+
+/// @brief Keeps dense failures backend-neutral for caller-specific diagnostics.
+Status FactorFailure(const char* code, const char* message) {
+  return Status::Failure({{Severity::kError, code, {}, "", "", message}});
+}
+
+/// @brief Preserves each value owner's size exception at the private boundary.
+MKL_INT MklSizeOrThrow(const std::size_t size, const char* message) {
+  const auto converted = ToMklSize(size);
+  if (!converted.HasValue()) throw std::length_error{message};
+  return converted.Value();
+}
+
+/// @brief Converts a contiguous BLAS vector length with its existing contract.
+MKL_INT VectorSizeOrThrow(const std::size_t size) {
+  return MklSizeOrThrow(size,
+                       "Dense vector size exceeds the MKL integer range.");
+}
+
+/// @brief Converts a matrix dimension with its existing exception contract.
+MKL_INT MatrixSizeOrThrow(const std::size_t size) {
+  return MklSizeOrThrow(size,
+                       "Dense matrix dimension exceeds the MKL integer range.");
+}
+
+/// @brief Rejects Matrix product inputs before allocation or empty returns.
+void RequireFinite(const double* values, const std::size_t size) {
+  for (std::size_t i = 0; i < size; ++i)
+    if (!std::isfinite(values[i]))
+      throw std::invalid_argument{"Matrix product inputs must be finite."};
+}
+
+}  // namespace
+
+DenseFactor::DenseFactor(const Matrix& coefficients,
+                         const bool positive_definite)
+    : factors_{coefficients},
+      pivots_(positive_definite ? 0U : coefficients.Rows()),
+      positive_definite_{positive_definite} {
+  for (std::size_t row = 0; row < coefficients.Rows(); ++row) {
+    double row_sum = 0.0;
+    for (std::size_t column = 0; column < coefficients.Columns(); ++column)
+      row_sum += std::abs(coefficients(row, column));
+    coefficient_norm_inf_ = (std::max)(coefficient_norm_inf_, row_sum);
+  }
+}
+
+Result<DenseFactor> DenseFactor::Lu(const Matrix& coefficients) {
+  return Create(coefficients, false);
+}
+
+Result<DenseFactor> DenseFactor::Cholesky(const Matrix& coefficients) {
+  return Create(coefficients, true);
+}
+
+Result<DenseFactor> DenseFactor::Create(const Matrix& coefficients,
+                                        const bool positive_definite) {
+  if (coefficients.Rows() != coefficients.Columns())
+    return Result<DenseFactor>::Failure(FactorFailure(
+        "dense-factor-shape", "Dense factor requires a square matrix."));
+  const auto size = ToMklSize(coefficients.Rows());
+  if (!size.HasValue()) return Result<DenseFactor>::Failure(size.GetStatus());
+  if (!IsFinite(coefficients))
+    return Result<DenseFactor>::Failure(FactorFailure(
+        "dense-factor-nonfinite", "Dense factor inputs must be finite."));
+  DenseFactor candidate{coefficients, positive_definite};
+  if (size.Value() != 0) {
+    const auto info =
+        positive_definite
+            ? LAPACKE_dpotrf(LAPACK_ROW_MAJOR, 'L', size.Value(),
+                             candidate.factors_.Data(), size.Value())
+            : LAPACKE_dgetrf(LAPACK_ROW_MAJOR, size.Value(), size.Value(),
+                             candidate.factors_.Data(), size.Value(),
+                             candidate.pivots_.data());
+    if (info != 0 || !IsFinite(candidate.factors_))
+      return Result<DenseFactor>::Failure(FactorFailure(
+          "dense-factor-pivot", "Dense numerical factorization failed."));
+  }
+  return Result<DenseFactor>::Success(std::move(candidate));
+}
+
+Result<Matrix> DenseFactor::Solve(const Matrix& rhs,
+                                  const bool transpose) const {
+  return SolveImpl(rhs, transpose, false, false);
+}
+
+Result<Matrix> DenseFactor::SolveRight(const Matrix& rhs,
+                                       const bool transpose) const {
+  return SolveImpl(rhs, transpose, true, false);
+}
+
+Result<Matrix> DenseFactor::SolveLower(const Matrix& rhs,
+                                       const bool transpose) const {
+  return SolveImpl(rhs, transpose, false, true);
+}
+
+Result<Matrix> DenseFactor::SolveLowerRight(const Matrix& rhs,
+                                            const bool transpose) const {
+  return SolveImpl(rhs, transpose, true, true);
+}
+
+Result<Matrix> DenseFactor::SolveImpl(const Matrix& rhs, const bool transpose,
+                                      const bool right,
+                                      const bool lower_only) const {
+  const auto order = factors_.Rows();
+  if ((right ? rhs.Columns() : rhs.Rows()) != order ||
+      (lower_only && !positive_definite_))
+    return Result<Matrix>::Failure(FactorFailure(
+        "dense-solve-shape", "Dense solve shape or factor property differs."));
+  const auto count = ToMklSize(right ? rhs.Rows() : rhs.Columns());
+  if (!count.HasValue()) return Result<Matrix>::Failure(count.GetStatus());
+  if (!IsFinite(rhs))
+    return Result<Matrix>::Failure(FactorFailure(
+        "dense-solve-nonfinite", "Dense solve RHS must be finite."));
+  if (order == 0U || count.Value() == 0)
+    return Result<Matrix>::Success(Matrix{rhs});
+  Matrix candidate = right && !lower_only ? Matrix{order, rhs.Rows()}
+                                         : Matrix{rhs};
+  const auto size = ToMklSize(order).Value();
+  if (lower_only) {
+    cblas_dtrsm(CblasRowMajor, right ? CblasRight : CblasLeft, CblasLower,
+                transpose ? CblasTrans : CblasNoTrans, CblasNonUnit,
+                right ? count.Value() : size, right ? size : count.Value(), 1.0,
+                factors_.Data(), size, candidate.Data(),
+                right ? size : count.Value());
+  } else {
+    // X*op(A)=B becomes op(A)^T*X^T=B^T. GETRS uses a compact n-by-nrhs
+    // row-major RHS with ldb=nrhs, including nonsquare multi-RHS storage.
+    if (right) {
+      for (std::size_t i = 0; i < rhs.Rows(); ++i)
+        for (std::size_t j = 0; j < order; ++j) candidate(j, i) = rhs(i, j);
+    }
+    const auto info =
+        positive_definite_
+            ? LAPACKE_dpotrs(LAPACK_ROW_MAJOR, 'L', size, count.Value(),
+                             factors_.Data(), size, candidate.Data(),
+                             count.Value())
+            : LAPACKE_dgetrs(LAPACK_ROW_MAJOR, transpose != right ? 'T' : 'N',
+                             size, count.Value(), factors_.Data(), size,
+                             pivots_.data(), candidate.Data(), count.Value());
+    if (info != 0)
+      return Result<Matrix>::Failure(
+          FactorFailure("dense-solve-backend", "Dense backend solve failed."));
+    if (right) {
+      Matrix transposed{rhs.Rows(), order};
+      for (std::size_t i = 0; i < rhs.Rows(); ++i)
+        for (std::size_t j = 0; j < order; ++j)
+          transposed(i, j) = candidate(j, i);
+      candidate = std::move(transposed);
+    }
+  }
+  if (!IsFinite(candidate))
+    return Result<Matrix>::Failure(FactorFailure(
+        "dense-solve-nonfinite", "Dense solve result is nonfinite."));
+  return Result<Matrix>::Success(std::move(candidate));
+}
+
+Result<Matrix> DenseFactor::Lower() const {
+  if (!positive_definite_)
+    return Result<Matrix>::Failure(FactorFailure(
+        "dense-factor-property", "An LU factor has no Cholesky lower factor."));
+  Matrix lower{factors_.Rows(), factors_.Columns()};
+  for (std::size_t i = 0; i < factors_.Rows(); ++i)
+    for (std::size_t j = 0; j <= i; ++j) lower(i, j) = factors_(i, j);
+  return Result<Matrix>::Success(std::move(lower));
+}
+
+double DenseFactor::Determinant() const {
+  double determinant = 1.0;
+  for (std::size_t i = 0; i < factors_.Rows(); ++i) {
+    determinant *= factors_(i, i);
+    if (positive_definite_)
+      determinant *= factors_(i, i);
+    else if (pivots_[i] != static_cast<MKL_INT>(i + 1U))
+      determinant = -determinant;
+  }
+  return determinant;
+}
+
+Result<double> DenseFactor::ReciprocalConditionInfinity() const {
+  if (!(coefficient_norm_inf_ > 0.0) || !std::isfinite(coefficient_norm_inf_))
+    return Result<double>::Failure(FactorFailure(
+        "dense-condition-norm", "Dense coefficient infinity norm is invalid."));
+  const auto order = factors_.Rows();
+  std::vector<double> row_sums(order, 0.0);
+  Matrix rhs{order, 1U};
+  for (std::size_t column = 0; column < order; ++column) {
+    if (column != 0U) rhs(column - 1U, 0U) = 0.0;
+    rhs(column, 0U) = 1.0;
+    const auto solved = Solve(rhs);
+    if (!solved.HasValue()) return Result<double>::Failure(solved.GetStatus());
+    for (std::size_t row = 0; row < order; ++row)
+      row_sums[row] += std::abs(solved.Value()(row, 0U));
+  }
+  double inverse_norm = 0.0;
+  for (const auto row_sum : row_sums)
+    inverse_norm = (std::max)(inverse_norm, row_sum);
+  if (!(inverse_norm > 0.0) || !std::isfinite(inverse_norm))
+    return Result<double>::Failure(FactorFailure(
+        "dense-condition-norm", "Dense inverse infinity norm is invalid."));
+  return Result<double>::Success((1.0 / coefficient_norm_inf_) / inverse_norm);
+}
+
+double DotValues(const double* left, const std::size_t size,
+                  const double* right) {
+  if (size == 0U) return 0.0;
+  return cblas_ddot(VectorSizeOrThrow(size), left, 1, right, 1);
+}
+
+double NormValues(const double* values, const std::size_t size) {
+  if (size == 0U) return 0.0;
+  return cblas_dnrm2(VectorSizeOrThrow(size), values, 1);
+}
+
+void ScaleValues(const double alpha, double* values, const std::size_t size) {
+  if (size == 0U) return;
+  cblas_dscal(VectorSizeOrThrow(size), alpha, values, 1);
+}
+
+void AxpyValues(const double alpha, const double* x, const std::size_t size,
+                double* y) {
+  if (size == 0U) return;
+  cblas_daxpy(VectorSizeOrThrow(size), alpha, x, 1, y, 1);
+}
+
+Vector Multiply(const Matrix& left, const Vector& right) {
+  if (left.Columns() != right.Size())
+    throw std::invalid_argument{
+        "Matrix-vector multiplication has incompatible dimensions."};
+  RequireFinite(left.Data(), left.Rows() * left.Columns());
+  RequireFinite(right.Data(), right.Size());
+  Vector result{left.Rows()};
+  if (left.Rows() == 0U || left.Columns() == 0U) return result;
+  const auto rows = MatrixSizeOrThrow(left.Rows());
+  const auto columns = MatrixSizeOrThrow(left.Columns());
+  cblas_dgemv(CblasRowMajor, CblasNoTrans, rows, columns, 1.0, left.Data(),
+              columns, right.Data(), 1, 0.0, result.Data(), 1);
+  return result;
+}
+
+Matrix Multiply(const Matrix& left, const Matrix& right) {
+  if (left.Columns() != right.Rows())
+    throw std::invalid_argument{
+        "Matrix multiplication has incompatible dimensions."};
+  RequireFinite(left.Data(), left.Rows() * left.Columns());
+  RequireFinite(right.Data(), right.Rows() * right.Columns());
+  Matrix result{left.Rows(), right.Columns()};
+  if (left.Rows() == 0U || left.Columns() == 0U || right.Columns() == 0U)
+    return result;
+  const auto rows = MatrixSizeOrThrow(left.Rows());
+  const auto columns = MatrixSizeOrThrow(right.Columns());
+  const auto inner = MatrixSizeOrThrow(left.Columns());
+  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, rows, columns, inner,
+              1.0, left.Data(), inner, right.Data(), columns, 0.0,
+              result.Data(), columns);
+  return result;
+}
+
 ScopedSingleEquationThreads::ScopedSingleEquationThreads(
     const std::size_t equation_count) noexcept {
   // The measured 1x1 PARDISO case has backend/thread overhead and an
@@ -45,60 +308,6 @@ LocalDensePolicy ReadScopedLocalDensePolicy() noexcept {
   return {mkl_domain_get_max_threads(MKL_DOMAIN_BLAS), mkl_get_dynamic()};
 }
 
-bool CholeskyLowerAndInverse(const Matrix& matrix, Matrix& lower,
-                             Matrix& inverse_lower, CholeskyFailure* failure) {
-  if (failure != nullptr) *failure = CholeskyFailure::kNone;
-  if (matrix.Rows() != matrix.Columns())
-    throw std::invalid_argument{"Cholesky requires a square matrix."};
-  const auto size = ToMklSize(matrix.Rows());
-  if (!size.HasValue())
-    throw std::length_error{"Cholesky size exceeds the MKL integer range."};
-  Matrix candidate_lower{matrix.Rows(), matrix.Columns()};
-  Matrix candidate_inverse{matrix.Rows(), matrix.Columns()};
-  if (matrix.Rows() == 0U) {
-    lower = std::move(candidate_lower);
-    inverse_lower = std::move(candidate_inverse);
-    return true;
-  }
-  std::vector<double> values(matrix.Rows() * matrix.Columns());
-  for (std::size_t row = 0; row < matrix.Rows(); ++row)
-    for (std::size_t column = 0; column < matrix.Columns(); ++column) {
-      if (!std::isfinite(matrix(row, column)))
-        throw std::invalid_argument{"Cholesky inputs must be finite."};
-      values[row * matrix.Columns() + column] = matrix(row, column);
-    }
-  if (LAPACKE_dpotrf(LAPACK_ROW_MAJOR, 'L', size.Value(), values.data(),
-                     size.Value()) != 0) {
-    if (failure != nullptr) *failure = CholeskyFailure::kPivot;
-    return false;
-  }
-  for (std::size_t row = 0; row < matrix.Rows(); ++row)
-    for (std::size_t column = 0; column <= row; ++column) {
-      candidate_lower(row, column) = values[row * matrix.Columns() + column];
-      if (!std::isfinite(candidate_lower(row, column))) {
-        if (failure != nullptr) *failure = CholeskyFailure::kPivot;
-        return false;
-      }
-    }
-  if (LAPACKE_dtrtri(LAPACK_ROW_MAJOR, 'L', 'N', size.Value(), values.data(),
-                     size.Value()) != 0) {
-    if (failure != nullptr) *failure = CholeskyFailure::kInverse;
-    return false;
-  }
-  for (std::size_t row = 0; row < matrix.Rows(); ++row)
-    for (std::size_t column = 0; column <= row; ++column) {
-      const auto value = values[row * matrix.Columns() + column];
-      if (!std::isfinite(value)) {
-        if (failure != nullptr) *failure = CholeskyFailure::kInverse;
-        return false;
-      }
-      candidate_inverse(row, column) = value;
-    }
-  lower = std::move(candidate_lower);
-  inverse_lower = std::move(candidate_inverse);
-  return true;
-}
-
 Vector MultiplyDiagonal(const Vector& diagonal, const Vector& rhs) {
   if (diagonal.Size() != rhs.Size())
     throw std::invalid_argument{"Diagonal product dimensions differ."};
@@ -122,11 +331,19 @@ Matrix MultiplyTransposeLeft(const Matrix& left, const Matrix& right) {
   if (left.Rows() != right.Rows()) {
     throw std::invalid_argument{"Transpose product dimensions differ."};
   }
-  Matrix transposed{left.Columns(), left.Rows()};
-  for (std::size_t row = 0; row < left.Rows(); ++row)
-    for (std::size_t column = 0; column < left.Columns(); ++column)
-      transposed(column, row) = left(row, column);
-  return transposed.Multiply(right);
+  RequireFinite(left.Data(), left.Rows() * left.Columns());
+  RequireFinite(right.Data(), right.Rows() * right.Columns());
+  Matrix result{left.Columns(), right.Columns()};
+  if (left.Columns() == 0U || left.Rows() == 0U || right.Columns() == 0U)
+    return result;
+  const auto rows = MatrixSizeOrThrow(left.Columns());
+  const auto columns = MatrixSizeOrThrow(right.Columns());
+  const auto inner = MatrixSizeOrThrow(left.Rows());
+  // Row-major A is stored inner-by-rows even when GEMM consumes A-transpose.
+  cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans, rows, columns, inner,
+              1.0, left.Data(), rows, right.Data(), columns, 0.0,
+              result.Data(), columns);
+  return result;
 }
 
 bool OuterProduct(const double* left, const std::size_t rows,

src/fesa/math/dense_blas_internal.h

@@ -2,6 +2,7 @@
 #define FESA_SRC_FESA_MATH_DENSE_BLAS_INTERNAL_H_
 
 #include <cstddef>
+#include <vector>
 
 #include <mkl.h>
 
@@ -11,6 +12,73 @@
 
 namespace fesa::dense_blas_internal {
 
+/// @brief Owns a double-precision numerical factor of a coefficient snapshot.
+/// @note Repeated RHS may reuse this immutable factor concurrently. Changing
+/// coefficients, dimension, matrix property or precision requires a new factor;
+/// matching sparsity alone never permits numerical reuse. No symbolic cache or
+/// borrowed coefficient storage exists. Backend thread policy stays
+/// caller-owned.
+class DenseFactor {
+ public:
+  /// @brief Factors a finite square general matrix with partial-pivoting LU.
+  /// @return An independent factor or a shape, finite, or singular failure.
+  static Result<DenseFactor> Lu(const Matrix& coefficients);
+
+  /// @brief Factors the lower triangle of a positive definite matrix.
+  /// @pre The caller establishes symmetry, including its own tolerance policy.
+  /// @return An independent factor or a shape, finite, or pivot failure.
+  /// @note All entries are checked for finiteness; upper values are not
+  /// averaged.
+  static Result<DenseFactor> Cholesky(const Matrix& coefficients);
+
+  /// @brief Solves op(A)*X=rhs with one or more row-major RHS columns.
+  /// @return A finite owned solution or failure, leaving factor and rhs intact.
+  Result<Matrix> Solve(const Matrix& rhs, bool transpose = false) const;
+
+  /// @brief Solves X*op(A)=rhs with one or more row-major RHS rows.
+  /// @return A finite owned solution or failure, leaving factor and rhs intact.
+  Result<Matrix> SolveRight(const Matrix& rhs, bool transpose = false) const;
+
+  /// @brief Solves op(L)*X=rhs using this Cholesky lower factor.
+  /// @return Failure for an LU factor, incompatible shape or nonfinite result.
+  Result<Matrix> SolveLower(const Matrix& rhs, bool transpose = false) const;
+
+  /// @brief Solves X*op(L)=rhs using this Cholesky lower factor.
+  /// @return Failure for an LU factor, incompatible shape or nonfinite result.
+  Result<Matrix> SolveLowerRight(const Matrix& rhs,
+                                 bool transpose = false) const;
+
+  /// @brief Returns an owned lower Cholesky factor with an exact zero upper
+  /// part.
+  /// @return Failure for an LU factor.
+  Result<Matrix> Lower() const;
+
+  /// @brief Computes the signed determinant from pivots and diagonal factors.
+  /// @note Empty determinant is one; double overflow/underflow propagates.
+  double Determinant() const;
+
+  /// @brief Computes (1/||A||inf)/||A^-1||inf without retaining an inverse.
+  /// @note Solves canonical RHS in ascending order and streams row absolute
+  /// sums. No estimator, acceptance threshold or regularization is applied.
+  /// @return Failure for an empty matrix or nonfinite norm/solution arithmetic.
+  Result<double> ReciprocalConditionInfinity() const;
+
+ private:
+  /// @brief Constructs a candidate before any backend factorization succeeds.
+  DenseFactor(const Matrix& coefficients, bool positive_definite);
+  /// @brief Implements factor creation while keeping failed candidates local.
+  static Result<DenseFactor> Create(const Matrix& coefficients,
+                                    bool positive_definite);
+  /// @brief Solves into a candidate with a shared finite/shape failure gate.
+  Result<Matrix> SolveImpl(const Matrix& rhs, bool transpose, bool right,
+                           bool lower_only) const;
+
+  Matrix factors_;
+  std::vector<MKL_INT> pivots_;
+  double coefficient_norm_inf_{0.0};
+  bool positive_definite_;
+};
+
 /// @brief Limits only default-policy one-equation backend calls to one thread.
 /// @note Existing caller TLS overrides are retained. Destruction restores the
 /// default request; process-wide MKL/OpenMP settings are never changed.
@@ -54,6 +122,37 @@ LocalDensePolicy ReadScopedLocalDensePolicy() noexcept;
 /// @return The converted length or an uncategorized overflow failure.
 Result<MKL_INT> ToMklSize(std::size_t size);
 
+/// @brief Computes a contiguous dot product without rejecting nonfinite input.
+/// @pre Both arrays contain size readable entries; null is valid for size zero.
+/// @throws std::length_error if size exceeds the backend integer range.
+double DotValues(const double* left, std::size_t size, const double* right);
+
+/// @brief Computes a robust contiguous Euclidean norm with BLAS propagation.
+/// @pre values contains size readable entries; null is valid for size zero.
+/// @throws std::length_error if size exceeds the backend integer range.
+double NormValues(const double* values, std::size_t size);
+
+/// @brief Scales contiguous values with BLAS nonfinite propagation.
+/// @pre values contains size writable entries; null is valid for size zero.
+/// @throws std::length_error if size exceeds the backend integer range.
+void ScaleValues(double alpha, double* values, std::size_t size);
+
+/// @brief Accumulates alpha*x into y, allowing x and y to be identical.
+/// @pre Arrays contain size entries; other overlap is unsupported.
+/// @note Null storage is valid for size zero; nonfinite values follow BLAS.
+/// @throws std::length_error if size exceeds the backend integer range.
+void AxpyValues(double alpha, const double* x, std::size_t size, double* y);
+
+/// @brief Multiplies owned row-major values with the Matrix failure contract.
+/// @throws std::invalid_argument for incompatible shapes or nonfinite input.
+/// @throws std::length_error if dimensions exceed the backend integer range.
+Vector Multiply(const Matrix& left, const Vector& right);
+
+/// @brief Multiplies owned row-major matrices with the Matrix failure contract.
+/// @throws std::invalid_argument for incompatible shapes or nonfinite input.
+/// @throws std::length_error if dimensions exceed the backend integer range.
+Matrix Multiply(const Matrix& left, const Matrix& right);
+
 /// @brief Computes left-transpose times right through ordinary dense MKL.
 /// @throws std::invalid_argument for mismatched rows or nonfinite operands.
 Matrix MultiplyTransposeLeft(const Matrix& left, const Matrix& right);
@@ -83,19 +182,6 @@ bool MultiplyTransposeRight(const double* left, std::size_t rows,
                             const double* right, std::size_t columns,
                             std::size_t inner, double* destination);
 
-/// @brief Distinguishes positive-factor failure from inverse failure.
-enum class CholeskyFailure { kNone, kPivot, kInverse };
-
-/// @brief Computes a lower Cholesky factor and its inverse with MKL LAPACK.
-/// @return False for a nonpositive pivot or nonfinite factor/inverse.
-/// @note Outputs are replaced only after both finite candidates succeed.
-/// The optional failure output is initialized even when validation throws.
-/// @throws std::invalid_argument for a nonsquare or nonfinite input.
-/// @throws std::length_error for a backend size overflow.
-bool CholeskyLowerAndInverse(const Matrix& matrix, Matrix& lower,
-                             Matrix& inverse_lower,
-                             CholeskyFailure* failure = nullptr);
-
 /// @brief Copies contiguous double values through the private BLAS backend.
 /// @param source Readable storage for size values when size is nonzero.
 /// @param size Number of contiguous values to copy.

src/fesa/math/finite_values.h

@@ -0,0 +1,40 @@
+#ifndef FESA_MATH_FINITE_VALUES_H_
+#define FESA_MATH_FINITE_VALUES_H_
+
+#include <array>
+#include <cmath>
+#include <cstddef>
+
+#include "fesa/math/matrix.h"
+#include "fesa/math/vector.h"
+
+namespace fesa::math_internal {
+
+/// @brief Checks ordinary doubles in ascending order without dereferencing empty storage.
+inline bool IsFiniteValues(const double* values, const std::size_t size) {
+  for (std::size_t index = 0U; index < size; ++index) {
+    if (!std::isfinite(values[index])) {
+      return false;
+    }
+  }
+  return true;
+}
+
+template <std::size_t Size>
+bool IsFinite(const std::array<double, Size>& values) {
+  return IsFiniteValues(values.data(), Size);
+}
+
+inline bool IsFinite(const Vector& values) {
+  return IsFiniteValues(values.Data(), values.Size());
+}
+
+inline bool IsFinite(const Matrix& values) {
+  const auto size = values.Rows() * values.Columns();
+  const double* const data = values.Data();
+  return IsFiniteValues(data, size);
+}
+
+}  // namespace fesa::math_internal
+
+#endif  // FESA_MATH_FINITE_VALUES_H_

src/fesa/math/high_low_arithmetic.h

@@ -88,6 +88,12 @@ inline double Rounded(HighLow value) {
 /// @brief Adds normalized expansions; rounding of discarded tail is bounded.
 inline HighLow Add(HighLow left, HighLow right) {
   if (!IsNormalized(left) || !IsNormalized(right)) return InvalidHighLow();
+  // Exact-zero additions only normalize signed zero in the retained pair.
+  // The nonzero graph and its discarded-tail association stay unchanged.
+  if (left.high == 0 && left.low == 0)
+    return {right.high, right.low == 0 ? 0.0 : right.low};
+  if (right.high == 0 && right.low == 0)
+    return {left.high, left.low == 0 ? 0.0 : left.low};
   const auto head = TwoSum(left.high, right.high);
   const auto tail = TwoSum(left.low, right.low);
   const auto middle = TwoSum(head.low, tail.high);
@@ -164,9 +170,16 @@ inline HighLow TwoProduct(double left, double right) {
 inline HighLow Multiply(HighLow left, HighLow right) {
   if (!IsNormalized(left) || !IsNormalized(right)) return InvalidHighLow();
   auto product = TwoProduct(left.high, right.high);
-  product = Add(product, TwoProduct(left.high, right.low));
-  product = Add(product, TwoProduct(left.low, right.high));
-  return Add(product, TwoProduct(left.low, right.low));
+  // Omitted terms are exact zero; retain the zero normalization that preceded
+  // a nonzero low-high term and the final canonicalization/invalid check.
+  if (right.low != 0)
+    product = Add(product, TwoProduct(left.high, right.low));
+  else if (left.low != 0)
+    product = Add(product, HighLow{});
+  if (left.low != 0)
+    product = Add(product, TwoProduct(left.low, right.high));
+  return Add(product, left.low != 0 && right.low != 0
+                          ? TwoProduct(left.low, right.low) : HighLow{});
 }
 /// @brief Orders canonical values; 2 means unordered/unsupported input.
 inline int Compare(HighLow left, HighLow right) {

src/fesa/math/jet_arithmetic.h

@@ -0,0 +1,38 @@
+#ifndef FESA_SRC_FESA_MATH_JET_ARITHMETIC_H_
+#define FESA_SRC_FESA_MATH_JET_ARITHMETIC_H_
+
+#include <array>
+#include <cstddef>
+
+namespace fesa::math_internal {
+/// @brief Multiplies fixed Jet matrices in row/column/index order.
+/// @note The caller retains its Jet precision, derivative storage, product rule
+/// and constant provenance. Each result starts at its actual constant zero.
+template <class Jet, class Constant>
+std::array<std::array<Jet, 3>, 3> MultiplyJetMatrix(
+    const std::array<std::array<Jet, 3>, 3>& left,
+    const std::array<std::array<Jet, 3>, 3>& right, Constant constant) {
+  std::array<std::array<Jet, 3>, 3> result{};
+  for (std::size_t row = 0; row < 3; ++row)
+    for (std::size_t column = 0; column < 3; ++column) {
+      result[row][column] = constant(0.0);
+      for (std::size_t index = 0; index < 3; ++index)
+        result[row][column] = result[row][column] +
+                              left[row][index] * right[index][column];
+    }
+  return result;
+}
+/// @brief Multiplies a fixed Jet matrix and actual vector RHS in row order.
+/// @note The zero-start component accumulation matches both Jet consumers.
+template <class Jet, class Constant>
+std::array<Jet, 3> MultiplyJetMatrixVector(
+    const std::array<std::array<Jet, 3>, 3>& matrix,
+    const std::array<Jet, 3>& vector, Constant constant) {
+  std::array<Jet, 3> result{constant(0.0),constant(0.0),constant(0.0)};
+  for (std::size_t row = 0; row < 3; ++row)
+    for (std::size_t column = 0; column < 3; ++column)
+      result[row] = result[row] + matrix[row][column] * vector[column];
+  return result;
+}
+}  // namespace fesa::math_internal
+#endif  // FESA_SRC_FESA_MATH_JET_ARITHMETIC_H_

src/fesa/math/jet_linear_solve.h

@@ -0,0 +1,85 @@
+#ifndef FESA_SRC_FESA_MATH_JET_LINEAR_SOLVE_H_
+#define FESA_SRC_FESA_MATH_JET_LINEAR_SOLVE_H_
+
+#include <array>
+#include <cstddef>
+#include <utility>
+
+#include "dense_blas_internal.h"
+
+namespace fesa::math_internal {
+/// @brief Solves one three-component Jet action by increasing derivative degree.
+/// @pre factor is LU of exactly the value coefficients of a. Channel zero of
+/// degree zero is the value; each later degree contains every retained channel.
+/// The reader/writer preserve the caller's derivative or Taylor convention.
+/// @note Each degree solves A0*x_degree=(b-A*x_lower)_degree using the same
+/// immutable factor. The caller's product supplies every product-rule term;
+/// the not-yet-written degree is zero. Distinct Jet storage/order stays owned
+/// by the caller. No inverse or primal-only approximation is formed.
+/// @return False on a nonfinite solve, leaving output and inputs unchanged.
+template <class Jet, class Constant, class Read, class Write>
+bool SolveJetAction(const std::array<std::array<Jet, 3>, 3>& a,
+                    const std::array<Jet, 3>& b,
+                    const dense_blas_internal::DenseFactor& factor,
+                    const std::array<std::size_t, 4>& channels,
+                    Constant constant, Read read, Write write,
+                    std::array<Jet, 3>& output) {
+  std::array<Jet, 3> candidate{constant(0.0), constant(0.0), constant(0.0)};
+  for (std::size_t degree = 0; degree < channels.size(); ++degree) {
+    Matrix rhs{3, channels[degree]};
+    for (std::size_t row = 0; row < 3; ++row) {
+      Jet residual = b[row];
+      // At degree zero the entire candidate is zero. Avoid 0*nonfinite higher
+      // derivatives affecting a value-only RHS before their own finite gate.
+      if (degree != 0)
+        for (std::size_t column = 0; column < 3; ++column)
+          residual = residual - a[row][column] * candidate[column];
+      for (std::size_t channel = 0; channel < channels[degree]; ++channel)
+        rhs(row, channel) = read(residual, degree, channel);
+    }
+    const auto solved = factor.Solve(rhs);
+    if (!solved.HasValue()) return false;
+    for (std::size_t row = 0; row < 3; ++row)
+      for (std::size_t channel = 0; channel < channels[degree]; ++channel)
+        write(candidate[row], degree, channel, solved.Value()(row, channel));
+  }
+  output = std::move(candidate);
+  return true;
+}
+/// @brief Applies a fixed double factor to every channel of one Jet RHS.
+/// @note Fixed coefficients have no variations, so the differentiated equation
+/// is A0*x_channel=b_channel with no lower-order product terms. Packing never
+/// rounds, contracts or omits a derivative channel. Candidate publication is
+/// atomic and the factor remains immutable.
+template <class Jet, class Constant, class Read, class Write>
+bool SolveConstantJetAction(const std::array<Jet, 3>& b,
+                            const dense_blas_internal::DenseFactor& factor,
+                            bool transpose,
+                            const std::array<std::size_t, 4>& channels,
+                            Constant constant, Read read, Write write,
+                            std::array<Jet, 3>& output) {
+  std::size_t count = 0;
+  for (const auto size : channels) count += size;
+  Matrix rhs{3, count};
+  std::size_t offset = 0;
+  for (std::size_t degree = 0; degree < channels.size(); ++degree) {
+    for (std::size_t row = 0; row < 3; ++row)
+      for (std::size_t channel = 0; channel < channels[degree]; ++channel)
+        rhs(row, offset + channel) = read(b[row], degree, channel);
+    offset += channels[degree];
+  }
+  const auto solved = factor.Solve(rhs, transpose);
+  if (!solved.HasValue()) return false;
+  std::array<Jet, 3> candidate{constant(0), constant(0), constant(0)};
+  offset = 0;
+  for (std::size_t degree = 0; degree < channels.size(); ++degree) {
+    for (std::size_t row = 0; row < 3; ++row)
+      for (std::size_t channel = 0; channel < channels[degree]; ++channel)
+        write(candidate[row], degree, channel, solved.Value()(row, offset+channel));
+    offset += channels[degree];
+  }
+  output = std::move(candidate);
+  return true;
+}
+}  // namespace fesa::math_internal
+#endif  // FESA_SRC_FESA_MATH_JET_LINEAR_SOLVE_H_

src/fesa/math/matrix.cpp

@@ -1,8 +1,5 @@
 #include "fesa/math/matrix.h"
 
-#include <cmath>
-#include <mkl.h>
-
 #include <limits>
 #include <stdexcept>
 #include <utility>
@@ -24,17 +21,6 @@ std::size_t CheckedStorageSize(const std::size_t rows,
   return rows * columns;
 }
 
-/// @brief Preserves the Matrix exception contract around private size
-/// conversion.
-MKL_INT MklSizeOrThrow(const std::size_t size) {
-  const auto converted = dense_blas_internal::ToMklSize(size);
-  if (!converted.HasValue()) {
-    throw std::length_error{
-        "Dense matrix dimension exceeds the MKL integer range."};
-  }
-  return converted.Value();
-}
-
 }  // namespace
 
 Matrix::Matrix(const std::size_t rows, const std::size_t columns,
@@ -88,6 +74,18 @@ std::size_t Matrix::Rows() const noexcept { return rows_; }
 
 std::size_t Matrix::Columns() const noexcept { return columns_; }
 
+double* Matrix::Data() noexcept { return values_.data(); }
+
+const double* Matrix::Data() const noexcept { return values_.data(); }
+
+double Matrix::Norm() const {
+  return dense_blas_internal::NormValues(Data(), values_.size());
+}
+
+void Matrix::Scale(const double alpha) {
+  dense_blas_internal::ScaleValues(alpha, Data(), values_.size());
+}
+
 double& Matrix::operator()(const std::size_t row, const std::size_t column) {
   if (row >= rows_ || column >= columns_) {
     throw std::out_of_range{"Matrix index is outside its dimensions."};
@@ -104,53 +102,11 @@ const double& Matrix::operator()(const std::size_t row,
 }
 
 Vector Matrix::Multiply(const Vector& rhs) const {
-  if (columns_ != rhs.Size()) {
-    throw std::invalid_argument{
-        "Matrix-vector multiplication has incompatible dimensions."};
-  }
-  for (const double value : values_)
-    if (!std::isfinite(value))
-      throw std::invalid_argument{"Matrix product inputs must be finite."};
-  for (std::size_t i = 0; i < rhs.Size(); ++i)
-    if (!std::isfinite(rhs.Data()[i]))
-      throw std::invalid_argument{"Matrix product inputs must be finite."};
-
-  Vector result{rows_};
-  if (rows_ == 0 || columns_ == 0) {
-    return result;
-  }
-
-  // The owned layout is row-major, so the leading dimension is the column
-  // count for the adapter call and remains invisible to public consumers.
-  cblas_dgemv(CblasRowMajor, CblasNoTrans, MklSizeOrThrow(rows_),
-              MklSizeOrThrow(columns_), 1.0, values_.data(),
-              MklSizeOrThrow(columns_), rhs.Data(), 1, 0.0, result.Data(), 1);
-  return result;
+  return dense_blas_internal::Multiply(*this, rhs);
 }
 
 Matrix Matrix::Multiply(const Matrix& rhs) const {
-  if (columns_ != rhs.rows_) {
-    throw std::invalid_argument{
-        "Matrix multiplication has incompatible dimensions."};
-  }
-  for (const double value : values_)
-    if (!std::isfinite(value))
-      throw std::invalid_argument{"Matrix product inputs must be finite."};
-  for (const double value : rhs.values_)
-    if (!std::isfinite(value))
-      throw std::invalid_argument{"Matrix product inputs must be finite."};
-
-  Matrix result{rows_, rhs.columns_};
-  if (rows_ == 0 || columns_ == 0 || rhs.columns_ == 0) {
-    return result;
-  }
-
-  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, MklSizeOrThrow(rows_),
-              MklSizeOrThrow(rhs.columns_), MklSizeOrThrow(columns_), 1.0,
-              values_.data(), MklSizeOrThrow(columns_), rhs.values_.data(),
-              MklSizeOrThrow(rhs.columns_), 0.0, result.values_.data(),
-              MklSizeOrThrow(rhs.columns_));
-  return result;
+  return dense_blas_internal::Multiply(*this, rhs);
 }
 
 }  // namespace fesa

src/fesa/math/sparse_matrix.cpp

@@ -10,7 +10,6 @@
 
 #include <mkl.h>
 
-#include "fesa/fem/dof_manager.h"
 #include "math/dense_blas_internal.h"
 
 namespace fesa {

src/fesa/math/vector.cpp

@@ -1,26 +1,11 @@
 #include "fesa/math/vector.h"
 
-#include <mkl.h>
-
 #include <stdexcept>
 #include <utility>
 
 #include "math/dense_blas_internal.h"
 
 namespace fesa {
-namespace {
-
-/// @brief Preserves the Vector exception contract around private size
-/// conversion.
-MKL_INT MklSizeOrThrow(const std::size_t size) {
-  const auto converted = dense_blas_internal::ToMklSize(size);
-  if (!converted.HasValue()) {
-    throw std::length_error{"Dense vector size exceeds the MKL integer range."};
-  }
-  return converted.Value();
-}
-
-}  // namespace
 
 Vector::Vector(const std::size_t size, const double value)
     : values_(size, value) {}
@@ -71,38 +56,22 @@ double Vector::Dot(const Vector& rhs) const {
     throw std::invalid_argument{
         "Vector dot product requires equal dimensions."};
   }
-  if (values_.empty()) {
-    return 0.0;
-  }
-
-  return cblas_ddot(MklSizeOrThrow(Size()), Data(), 1, rhs.Data(), 1);
+  return dense_blas_internal::DotValues(Data(), Size(), rhs.Data());
 }
 
 double Vector::Norm() const {
-  if (values_.empty()) {
-    return 0.0;
-  }
-
-  return cblas_dnrm2(MklSizeOrThrow(Size()), Data(), 1);
+  return dense_blas_internal::NormValues(Data(), Size());
 }
 
 void Vector::Scale(const double alpha) {
-  if (values_.empty()) {
-    return;
-  }
-
-  cblas_dscal(MklSizeOrThrow(Size()), alpha, Data(), 1);
+  dense_blas_internal::ScaleValues(alpha, Data(), Size());
 }
 
 void Vector::Axpy(const double alpha, const Vector& x) {
   if (Size() != x.Size()) {
     throw std::invalid_argument{"Vector axpy requires equal dimensions."};
   }
-  if (values_.empty()) {
-    return;
-  }
-
-  cblas_daxpy(MklSizeOrThrow(Size()), alpha, x.Data(), 1, Data(), 1);
+  dense_blas_internal::AxpyValues(alpha, x.Data(), Size(), Data());
 }
 
 }  // namespace fesa

src/fesa/math/vector3_operations.h

@@ -0,0 +1,57 @@
+#ifndef FESA_MATH_VECTOR3_OPERATIONS_H_
+#define FESA_MATH_VECTOR3_OPERATIONS_H_
+
+#include <array>
+#include <cstddef>
+
+#include "fesa/math/matrix.h"
+#include "fesa/math/vector3.h"
+
+namespace fesa::math_internal {
+
+/// @brief Reduces four source-ordered geometric vectors with scalar-left products.
+inline Vector3 WeightedSum(const std::array<double, 4>& weights,
+                           const std::array<Vector3, 4>& values) {
+  Vector3 result{};
+  for (std::size_t node = 0U; node < 4U; ++node) {
+    result = result + weights[node] * values[node];
+  }
+  return result;
+}
+
+/// @brief Cancels the first-node anchor before the same ordered reduction.
+inline Vector3 DerivativeSum(const std::array<double, 4>& derivatives,
+                             const std::array<Vector3, 4>& values) {
+  std::array<Vector3, 4> relative{};
+  for (std::size_t node = 0U; node < 4U; ++node) {
+    relative[node] = values[node] - values[0U];
+  }
+  return WeightedSum(derivatives, relative);
+}
+
+/// @brief Forms the positive cross-product matrix in fixed Cartesian coordinates.
+inline Matrix Skew(const Vector3& vector) {
+  Matrix result{3U, 3U};
+  result(0U, 1U) = -vector.Z();
+  result(0U, 2U) = vector.Y();
+  result(1U, 0U) = vector.Z();
+  result(1U, 2U) = -vector.X();
+  result(2U, 0U) = -vector.Y();
+  result(2U, 1U) = vector.X();
+  return result;
+}
+
+/// @brief Applies a fixed 3-by-3 geometric frame with the original component graph.
+inline Vector3 Multiply(const Matrix& matrix, const Vector3& vector) {
+  return Vector3{
+      matrix(0U, 0U) * vector.X() + matrix(0U, 1U) * vector.Y() +
+          matrix(0U, 2U) * vector.Z(),
+      matrix(1U, 0U) * vector.X() + matrix(1U, 1U) * vector.Y() +
+          matrix(1U, 2U) * vector.Z(),
+      matrix(2U, 0U) * vector.X() + matrix(2U, 1U) * vector.Y() +
+          matrix(2U, 2U) * vector.Z()};
+}
+
+}  // namespace fesa::math_internal
+
+#endif  // FESA_MATH_VECTOR3_OPERATIONS_H_