// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

// solvers.hpp
//
// Solvers for systems of linear equations
//
// This file implements various solvers for systems of linear equations.
// It includes both direct and iterative solvers.
//
// Main features and selection guide:
// - Direct solvers
//   * SolverType::LU        — general dense matrix, O(n^3), most versatile
//   * SolverType::Cholesky  — for symmetric positive definite matrices, about half the cost of LU
//   * SolverType::QR        — least squares for overdetermined systems (m > n), stable even when rank-deficient
//   * SolverType::SVD       — handles rank-deficient cases, most stable but slowest
//   * SolverType::LDL       — symmetric (also non-positive-definite), more general than Cholesky
//   * SolverType::Auto      — inspects matrix properties and selects automatically
// - Iterative solvers
//   * conjugate_gradient()  — for symmetric positive definite matrices, suited to large sparse matrices
//   * Jacobi / Gauss-Seidel — convergence guaranteed for diagonally dominant matrices
// - Solvers for special matrices
//   * tridiagonal matrix (Thomas method), band matrix solver
//
// See: decomposition.hpp (LU/Cholesky/QR/SVD implementations),
//       eigenvalues.hpp (eigenvalue problems),
//       sparse_lu.hpp (sparse matrix LU)
//
// Note: for large-scale problems, CG/MINRES/GMRES in
//     sangi::linalg::iterative_solvers.hpp are recommended

#ifndef SANGI_SOLVERS_HPP
#define SANGI_SOLVERS_HPP

#include <math/core/matrix.hpp>
#include <math/core/vector.hpp>
#include <math/core/traits.hpp>
#include <math/linalg/decomposition.hpp>
#include <cassert>

#include <string>
#include <stdexcept>
#include <limits>
#include <cmath>
#include <vector>
#include <utility>

namespace sangi {
    namespace algorithms {

        // Forward declaration
        template<typename MA, typename VB>
            requires sangi::BaseMatrixLike<MA>
                  && sangi::BaseVectorLike<VB>
                  && std::same_as<sangi::element_t<MA>, sangi::element_t<VB>>
        Vector<sangi::element_t<MA>> conjugate_gradient(const MA& a, const VB& b,
            sangi::element_t<MA> tolerance = std::numeric_limits<sangi::element_t<MA>>::epsilon() * 10,
            std::size_t max_iterations = 1000);

        //-----------------------------------------------------------------------------
        // Utilities for solver selection
        //-----------------------------------------------------------------------------

        // Enumeration of solver types
        enum class SolverType {
            Auto,     // Automatically select the optimal solver
            LU,       // LU decomposition solver
            Cholesky, // Cholesky decomposition solver
            QR,       // QR decomposition solver
            SVD,      // Singular value decomposition solver
            LDL,      // LDL decomposition solver
            Iterative // Iterative solver
        };

        // Convert a string to a solver type
        inline SolverType parse_solver_type(const std::string& solver_name) {
            if (solver_name == "auto") return SolverType::Auto;
            if (solver_name == "lu") return SolverType::LU;
            if (solver_name == "cholesky") return SolverType::Cholesky;
            if (solver_name == "qr") return SolverType::QR;
            if (solver_name == "svd") return SolverType::SVD;
            if (solver_name == "ldl") return SolverType::LDL;
            if (solver_name == "iterative") return SolverType::Iterative;

            // The default is automatic selection
            return SolverType::Auto;
        }

        //-----------------------------------------------------------------------------
        // Direct solvers
        //-----------------------------------------------------------------------------

        // General-purpose solver interface (Ax = b)
        template<typename T>
        Vector<T> solve(const BaseMatrix<T>& a, const BaseVector<T>& b,
            SolverType solver_type = SolverType::Auto) {
            // Dimension check
            if (a.rows() != b.size()) {
                assert(false && "DimensionError: solve: matrix and vector dimensions mismatch");
                throw DimensionError("solve: matrix and vector dimensions mismatch");
            }

            // Branch the dispatch on whether T is a floating-point-like type.
            // For types closed within Q such as sangi::Rational, Cholesky/QR/SVD/CG
            // which use sqrt cannot be instantiated, so restrict to LU only.
            if constexpr (numeric_traits<T>::is_floating_point) {
                // Solver selection logic
                if (solver_type == SolverType::Auto) {
                    if (a.is_square()) {
                        // Symmetry check
                        bool is_symmetric = true;
                        for (typename Matrix<T>::size_type i = 0; i < a.rows() && is_symmetric; ++i) {
                            for (typename Matrix<T>::size_type j = i + 1; j < a.cols() && is_symmetric; ++j) {
                                if (std::abs(a(i, j) - a(j, i)) > std::numeric_limits<T>::epsilon() * 10) {
                                    is_symmetric = false;
                                }
                            }
                        }

                        if (is_symmetric) {
                            // Positive-definiteness check (as a simple method, check whether the diagonal elements are positive)
                            bool is_positive_definite = true;
                            for (typename Matrix<T>::size_type i = 0; i < a.rows() && is_positive_definite; ++i) {
                                if (a(i, i) <= 0) {
                                    is_positive_definite = false;
                                }
                            }

                            if (is_positive_definite) {
                                // Cholesky decomposition is optimal for symmetric positive definite matrices
                                return cholesky_solve(a, b);
                            }
                            else {
                                // LDL decomposition is suitable for symmetric matrices
                                return ldl_solve(a, b);
                            }
                        }
                        else {
                            // LU decomposition is the general choice for non-symmetric matrices
                            return lu_solve(a, b);
                        }
                    }
                    else {
                        // For non-square matrices, least-squares solution (QR or SVD)
                        return qr_solve(a, b);
                    }
                }

                // Explicit solver selection
                switch (solver_type) {
                case SolverType::LU:
                    return lu_solve(a, b);
                case SolverType::Cholesky:
                    return cholesky_solve(a, b);
                case SolverType::QR:
                    return qr_solve(a, b);
                case SolverType::SVD:
                    return svd_solve(a, b);
                case SolverType::LDL:
                    return ldl_solve(a, b);
                case SolverType::Iterative:
                    return conjugate_gradient(a, b);
                default:
                    // Unreachable
                    throw std::runtime_error("solve: unknown solver type");
                }
            } else {
                // Exact types such as Rational / Int: only LU is closed within Q.
                // Cholesky/QR/SVD/CG require real-number operations for sqrt and
                // iterative convergence tests, so they are unsuitable. Solve with LU
                // either when SolverType::LU is given explicitly or for Auto.
                if (solver_type != SolverType::Auto && solver_type != SolverType::LU) {
                    throw std::runtime_error(
                        "solve: non-floating-point types (Rational, etc.) support only the LU solver");
                }
                return lu_solve(a, b);
            }
        }

        // String-based solver selection interface
        template<typename T>
        Vector<T> solve(const BaseMatrix<T>& a, const BaseVector<T>& b, const std::string& solver_name) {
            return solve(a, b, parse_solver_type(solver_name));
        }

        //-----------------------------------------------------------------------------
        // solve_to_precision: Iterative Refinement (Wilkinson 1948)
        //
        // Purpose: when the user specifies the required number of digits (target_digits),
        //       automatically estimate the algorithm/matrix condition number and perform
        //       iterative refinement until the required precision is met.
        //
        // Behavior:
        //   1. Obtain an initial solution x via LU solve
        //   2. Measure the residual r = b - A·x
        //   3. Return if ||r||/||b|| < 10^(-target_digits)
        //   4. Otherwise double sangi::Float::defaultPrecision and solve again
        //   5. If it has not converged after max_iter iterations, return the final solution with a warning
        //
        // Constraints:
        //   - Only meaningful for T = sangi::Float (float/double have fixed precision)
        //   - Modifies the global defaultPrecision, so beware of multi-threading
        //     (saved & restored within the function, but interferes if other threads use it concurrently)
        //
        // See: Wilkinson 1948 "Iterative Refinement", Higham 2002 "Accuracy and
        //       Stability of Numerical Algorithms" §12.
        //-----------------------------------------------------------------------------

        struct SolveToPrecisionResult {
            int iterations = 0;
            int final_precision_digits = 0;
            double relative_residual = 0.0;
            bool converged = false;
        };

        template<typename T>
        Vector<T> solve_to_precision(const BaseMatrix<T>& A, const BaseVector<T>& b,
                                      int target_digits,
                                      SolveToPrecisionResult* report = nullptr,
                                      int max_iter = 20) {
            if (A.rows() != b.size()) {
                throw DimensionError("solve_to_precision: dim mismatch");
            }

            if constexpr (std::is_floating_point_v<T>) {
                // float/double: precision is fixed. Just solve once and measure the residual.
                Vector<T> x = lu_solve(A, b);
                if (report) {
                    Vector<T> r_vec(b.size());
                    for (std::size_t i = 0; i < b.size(); ++i) {
                        T s = T(0);
                        for (std::size_t j = 0; j < A.cols(); ++j) s += A(i, j) * x[j];
                        r_vec[i] = b[i] - s;
                    }
                    T r_norm = r_vec.norm();  // via member function (Vector<T>)
                    // b is a BaseVector — does it have a member norm()?
                    T b_norm_sq = T(0);
                    for (std::size_t i = 0; i < b.size(); ++i) b_norm_sq = b_norm_sq + b[i] * b[i];
                    T b_norm = std::sqrt(b_norm_sq);
                    report->iterations = 1;
                    report->final_precision_digits = std::numeric_limits<T>::digits10;
                    report->relative_residual = static_cast<double>((b_norm > T(0)) ? r_norm / b_norm : r_norm);
                    report->converged = (report->relative_residual < std::pow(10.0, -target_digits));
                }
                return x;
            } else {
                // Arbitrary-precision types such as sangi::Float: iterative refinement (precision doubling)
                const int initial_p = T::defaultPrecision();
                // Initial precision: target + safety margin (16 digits ≈ Wilkinson's log2(n) + log2(κ) upper bound)
                int p = std::max(initial_p, target_digits + 16);

                Vector<T> x(b.size());
                double rel_residual = 1.0;
                int iter = 0;
                bool converged = false;

                for (iter = 0; iter < max_iter; ++iter) {
                    // Solve at the current precision
                    int saved = T::defaultPrecision();
                    T::setDefaultPrecision(p);
                    try {
                        x = lu_solve(A, b);
                    } catch (...) {
                        T::setDefaultPrecision(saved);
                        throw;
                    }

                    // Residual r = b - A·x (measured at the current precision)
                    Vector<T> r_vec(b.size());
                    for (std::size_t i = 0; i < b.size(); ++i) {
                        T s = T(0);
                        for (std::size_t j = 0; j < A.cols(); ++j) s += A(i, j) * x[j];
                        r_vec[i] = b[i] - s;
                    }
                    T r_norm = r_vec.norm();  // via member function (Vector<T>)
                    // b is a BaseVector — does it have a member norm()?
                    T b_norm_sq = T(0);
                    for (std::size_t i = 0; i < b.size(); ++i) b_norm_sq = b_norm_sq + b[i] * b[i];
                    T b_norm = std::sqrt(b_norm_sq);
                    T rel = (b_norm > T(0)) ? r_norm / b_norm : r_norm;
                    rel_residual = rel.toDouble();

                    T::setDefaultPrecision(saved);

                    // Convergence test: ||r||/||b|| < 10^(-target_digits)
                    if (rel_residual < std::pow(10.0, -target_digits)) {
                        converged = true;
                        ++iter;
                        break;
                    }

                    // Double the precision
                    p *= 2;
                }

                if (report) {
                    report->iterations = iter;
                    report->final_precision_digits = p;
                    report->relative_residual = rel_residual;
                    report->converged = converged;
                }
                return x;
            }
        }

    } // namespace algorithms

    //-------------------------------------------------------------------------
    // Matrix<T>::solve — declared in math/core/matrix.hpp, defined here so that
    // matrix.hpp does not pick up a dependency on the algorithms layer.
    //
    // Wraps algorithms::solve(*this, b) in a try-block and returns
    // std::nullopt on dimension mismatch or singular system, so that
    // sangi::Matrix satisfies the
    //     { m.solve(v) } -> std::convertible_to<std::optional<V>>
    // concept used by sangi::newton_raphson_nd and other generic
    // solver-driven algorithms.
    //
    // Users must include <math/linalg/solvers.hpp> in addition to
    // <math/core/matrix.hpp> to make this member callable.
    //-------------------------------------------------------------------------
    template<typename T>
    [[nodiscard]] std::optional<Vector<T>> Matrix<T>::solve(const Vector<T>& b) const {
        try {
            return algorithms::solve(*this, b);
        }
        catch (...) {
            return std::nullopt;
        }
    }

    namespace algorithms {

        // Solver for multiple right-hand-side vectors (AX = B)
        template<typename T>
        Matrix<T> solve_multi(const BaseMatrix<T>& a, const BaseMatrix<T>& b,
            SolverType solver_type = SolverType::Auto) {
            // Dimension check
            if (a.rows() != b.rows()) {
                assert(false && "DimensionError: solve_multi: matrix dimensions mismatch");
                throw DimensionError("solve_multi: matrix dimensions mismatch");
            }

            const auto m = a.rows();
            const auto n = a.cols();
            const auto nrhs = b.cols();

#if SANGI_HAS_MKL
            if constexpr (std::is_same_v<T, float> || std::is_same_v<T, double>) {
                if (a.is_square() && (solver_type == SolverType::Auto || solver_type == SolverType::LU)) {
                    // LAPACKE_gesv: solve multiple right-hand sides at once
                    Matrix<T> Acopy = a;
                    Matrix<T> x = b;
                    std::vector<lapack_int> ipiv(n);
                    lapack_int info;
                    if constexpr (std::is_same_v<T, float>)
                        info = LAPACKE_sgesv(LAPACK_ROW_MAJOR, (lapack_int)n, (lapack_int)nrhs,
                                             Acopy.data(), (lapack_int)n, ipiv.data(),
                                             x.data(), (lapack_int)nrhs);
                    else
                        info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, (lapack_int)n, (lapack_int)nrhs,
                                             Acopy.data(), (lapack_int)n, ipiv.data(),
                                             x.data(), (lapack_int)nrhs);
                    if (info > 0)
                        throw MathError("solve_multi: singular matrix detected");
                    return x;
                }
            }
#endif

            // Initialize the result matrix
            Matrix<T> x(n, nrhs);

            // Apply the solver to each right-hand-side vector
            Vector<T> bi(m);
            Vector<T> xi(n);

            for (typename Matrix<T>::size_type j = 0; j < nrhs; ++j) {
                // Extract the right-hand-side vector
                for (typename Matrix<T>::size_type i = 0; i < m; ++i) {
                    bi[i] = b(i, j);
                }

                // Solve the equation
                xi = solve(a, bi, solver_type);

                // Store the result
                for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                    x(i, j) = xi[i];
                }
            }

            return x;
        }

        // String-based multiple right-hand-side solver interface
        template<typename T>
        Matrix<T> solve_multi(const BaseMatrix<T>& a, const BaseMatrix<T>& b, const std::string& solver_name) {
            return solve_multi(a, b, parse_solver_type(solver_name));
        }

        //-----------------------------------------------------------------------------
        // Iterative solvers
        //-----------------------------------------------------------------------------

        // Conjugate gradient method
        template<typename MA, typename VB>
            requires sangi::BaseMatrixLike<MA>
                  && sangi::BaseVectorLike<VB>
                  && std::same_as<sangi::element_t<MA>, sangi::element_t<VB>>
        Vector<sangi::element_t<MA>> conjugate_gradient(const MA& a, const VB& b,
            sangi::element_t<MA> tolerance,
            std::size_t max_iterations) {
            using T = sangi::element_t<MA>;
            SANGI_STATIC_ASSERT_SQUARE(MA);
            SANGI_STATIC_ASSERT_MATRIX_ROWS_EQ_VECTOR_SIZE(MA, VB);

            // Dimension check
            if (a.rows() != a.cols()) {
                assert(false && "DimensionError: conjugate_gradient: matrix must be square");
                throw DimensionError("conjugate_gradient: matrix must be square");
            }

            if (a.rows() != b.size()) {
                assert(false && "DimensionError: conjugate_gradient: matrix and vector dimensions mismatch");
                throw DimensionError("conjugate_gradient: matrix and vector dimensions mismatch");
            }

            const auto n = a.rows();
            Vector<T> x(n, T(0));     // Set the initial solution to zero
            Vector<T> r = b;          // Initial residual r = b - A*x = b
            Vector<T> p = r;          // Initial search direction
            Vector<T> ap(n);          // A*p buffer (allocated once outside the loop)

            T r_norm = dot(r, r);     // Residual norm (squared)
            const T b_norm = r_norm;  // Right-hand-side vector norm (squared) = r = b at the initial point

            // Zero-vector check
            if (b_norm < std::numeric_limits<T>::epsilon()) {
                return x;  // If the right-hand side is zero, the zero vector is the solution
            }

            // Relative tolerance (squared): ‖r‖² ≤ (tolerance * ‖b‖)²
            const T rel_tol = tolerance * tolerance * b_norm;

            // Iteration
            for (typename Matrix<T>::size_type iter = 0; iter < max_iterations; ++iter) {
                // Convergence check
                if (r_norm <= rel_tol) {
                    break;
                }

                // Compute A*p (ap already allocated outside the loop)
                for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                    T s = T(0);
                    for (typename Matrix<T>::size_type j = 0; j < n; ++j)
                        s += a(i, j) * p[j];
                    ap[i] = s;
                }

                // Step size α = r^T*r / (p^T*A*p)
                const T pAp = dot(p, ap);
                if (std::abs(pAp) < std::numeric_limits<T>::epsilon()) break;
                const T alpha = r_norm / pAp;

                // Update the solution and residual (fused loop)
                axpy(alpha, p, x);
                axpy(-alpha, ap, r);

                // New residual norm
                const T r_next_norm = dot(r, r);

                // Update the direction vector: p = r + beta * p
                if (r_norm < std::numeric_limits<T>::epsilon()) break;  // Avoid division by zero
                const T beta = r_next_norm / r_norm;
                axpby(T(1), r, beta, p);

                // Update the residual norm
                r_norm = r_next_norm;
            }

            return x;
        }

        // Jacobi method
        template<typename MA, typename VB>
            requires sangi::BaseMatrixLike<MA>
                  && sangi::BaseVectorLike<VB>
                  && std::same_as<sangi::element_t<MA>, sangi::element_t<VB>>
        Vector<sangi::element_t<MA>> jacobi_method(const MA& a, const VB& b,
            sangi::element_t<MA> tolerance = std::numeric_limits<sangi::element_t<MA>>::epsilon() * 10,
            typename Matrix<sangi::element_t<MA>>::size_type max_iterations = 1000) {
            using T = sangi::element_t<MA>;
            SANGI_STATIC_ASSERT_SQUARE(MA);
            SANGI_STATIC_ASSERT_MATRIX_ROWS_EQ_VECTOR_SIZE(MA, VB);

            // Dimension check
            if (a.rows() != a.cols()) {
                assert(false && "DimensionError: jacobi_method: matrix must be square");
                throw DimensionError("jacobi_method: matrix must be square");
            }

            if (a.rows() != b.size()) {
                assert(false && "DimensionError: jacobi_method: matrix and vector dimensions mismatch");
                throw DimensionError("jacobi_method: matrix and vector dimensions mismatch");
            }

            const auto n = a.rows();
            Vector<T> x(n, T(0));    // Set the initial solution to zero
            Vector<T> x_new(n);      // New solution

            // Zero check on the diagonal elements
            for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                if (std::abs(a(i, i)) < std::numeric_limits<T>::epsilon()) {
                    throw MathError("jacobi_method: zero diagonal element");
                }
            }

            // Iteration
            for (typename Matrix<T>::size_type iter = 0; iter < max_iterations; ++iter) {
                // Compute the new solution
                for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                    T sum = b[i];
                    for (typename Matrix<T>::size_type j = 0; j < n; ++j) {
                        if (i != j) {
                            sum -= a(i, j) * x[j];
                        }
                    }
                    x_new[i] = sum / a(i, i);
                }

                // Convergence check
                T diff_norm = 0;
                for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                    diff_norm += (x_new[i] - x[i]) * (x_new[i] - x[i]);
                }
                diff_norm = std::sqrt(diff_norm);

                if (diff_norm < tolerance) {
                    break;
                }

                // Update the solution
                std::swap(x, x_new);
            }

            return x;
        }

        // Gauss-Seidel method
        template<typename MA, typename VB>
            requires sangi::BaseMatrixLike<MA>
                  && sangi::BaseVectorLike<VB>
                  && std::same_as<sangi::element_t<MA>, sangi::element_t<VB>>
        Vector<sangi::element_t<MA>> gauss_seidel_method(const MA& a, const VB& b,
            sangi::element_t<MA> tolerance = std::numeric_limits<sangi::element_t<MA>>::epsilon() * 10,
            typename Matrix<sangi::element_t<MA>>::size_type max_iterations = 1000) {
            using T = sangi::element_t<MA>;
            SANGI_STATIC_ASSERT_SQUARE(MA);
            SANGI_STATIC_ASSERT_MATRIX_ROWS_EQ_VECTOR_SIZE(MA, VB);

            // Dimension check
            if (a.rows() != a.cols()) {
                assert(false && "DimensionError: gauss_seidel_method: matrix must be square");
                throw DimensionError("gauss_seidel_method: matrix must be square");
            }

            if (a.rows() != b.size()) {
                assert(false && "DimensionError: gauss_seidel_method: matrix and vector dimensions mismatch");
                throw DimensionError("gauss_seidel_method: matrix and vector dimensions mismatch");
            }

            const auto n = a.rows();
            Vector<T> x(n, T(0));   // Set the initial solution to zero
            Vector<T> x_old(n);     // Previous solution

            // Zero check on the diagonal elements
            for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                if (std::abs(a(i, i)) < std::numeric_limits<T>::epsilon()) {
                    throw MathError("gauss_seidel_method: zero diagonal element");
                }
            }

            // Iteration
            for (typename Matrix<T>::size_type iter = 0; iter < max_iterations; ++iter) {
                // Save the previous solution
                std::copy(x.begin(), x.end(), x_old.begin());

                // Compute the new solution
                for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                    T sum = b[i];

                    // Elements below i (already updated)
                    for (typename Matrix<T>::size_type j = 0; j < i; ++j) {
                        sum -= a(i, j) * x[j];
                    }

                    // Elements above i (not yet updated)
                    for (typename Matrix<T>::size_type j = i + 1; j < n; ++j) {
                        sum -= a(i, j) * x_old[j];
                    }

                    x[i] = sum / a(i, i);
                }

                // Convergence check
                T diff_norm = 0;
                for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                    diff_norm += (x[i] - x_old[i]) * (x[i] - x_old[i]);
                }
                diff_norm = std::sqrt(diff_norm);

                if (diff_norm < tolerance) {
                    break;
                }
            }

            return x;
        }

        // SOR method (successive over-relaxation)
        template<typename MA, typename VB>
            requires sangi::BaseMatrixLike<MA>
                  && sangi::BaseVectorLike<VB>
                  && std::same_as<sangi::element_t<MA>, sangi::element_t<VB>>
        Vector<sangi::element_t<MA>> sor_method(const MA& a, const VB& b,
            sangi::element_t<MA> omega = static_cast<sangi::element_t<MA>>(1.5),  // Relaxation factor
            sangi::element_t<MA> tolerance = std::numeric_limits<sangi::element_t<MA>>::epsilon() * 10,
            typename Matrix<sangi::element_t<MA>>::size_type max_iterations = 1000) {
            using T = sangi::element_t<MA>;
            SANGI_STATIC_ASSERT_SQUARE(MA);
            SANGI_STATIC_ASSERT_MATRIX_ROWS_EQ_VECTOR_SIZE(MA, VB);

            // Dimension check
            if (a.rows() != a.cols()) {
                assert(false && "DimensionError: sor_method: matrix must be square");
                throw DimensionError("sor_method: matrix must be square");
            }

            if (a.rows() != b.size()) {
                assert(false && "DimensionError: sor_method: matrix and vector dimensions mismatch");
                throw DimensionError("sor_method: matrix and vector dimensions mismatch");
            }

            // Range check on the relaxation factor
            if (omega <= 0 || omega >= 2) {
                throw std::invalid_argument("sor_method: relaxation parameter must be in (0,2)");
            }

            const auto n = a.rows();
            Vector<T> x(n, T(0));   // Set the initial solution to zero
            Vector<T> x_old(n);     // Previous solution

            // Zero check on the diagonal elements
            for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                if (std::abs(a(i, i)) < std::numeric_limits<T>::epsilon()) {
                    throw MathError("sor_method: zero diagonal element");
                }
            }

            // Iteration
            for (typename Matrix<T>::size_type iter = 0; iter < max_iterations; ++iter) {
                // Save the previous solution
                std::copy(x.begin(), x.end(), x_old.begin());

                // Compute the new solution
                for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                    T sum1 = 0;  // Contribution of already-updated elements
                    T sum2 = 0;  // Contribution of not-yet-updated elements

                    // Elements below i (already updated)
                    for (typename Matrix<T>::size_type j = 0; j < i; ++j) {
                        sum1 += a(i, j) * x[j];
                    }

                    // Elements above i (not yet updated)
                    for (typename Matrix<T>::size_type j = i + 1; j < n; ++j) {
                        sum2 += a(i, j) * x_old[j];
                    }

                    // SOR update formula
                    T x_gs = (b[i] - sum1 - sum2) / a(i, i);  // Gauss-Seidel solution
                    x[i] = x_old[i] + omega * (x_gs - x_old[i]);
                }

                // Convergence check
                T diff_norm = 0;
                for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                    diff_norm += (x[i] - x_old[i]) * (x[i] - x_old[i]);
                }
                diff_norm = std::sqrt(diff_norm);

                if (diff_norm < tolerance) {
                    break;
                }
            }

            return x;
        }

        //-----------------------------------------------------------------------------
        // Least-squares solvers
        //-----------------------------------------------------------------------------

        // Ordinary least-squares solver
        template<typename MA, typename VB>
            requires sangi::BaseMatrixLike<MA>
                  && sangi::BaseVectorLike<VB>
                  && std::same_as<sangi::element_t<MA>, sangi::element_t<VB>>
        Vector<sangi::element_t<MA>> least_squares(const MA& a, const VB& b,
            SolverType solver_type = SolverType::Auto) {
            using T = sangi::element_t<MA>;
            SANGI_STATIC_ASSERT_MATRIX_ROWS_EQ_VECTOR_SIZE(MA, VB);

            // Dimension check
            if (a.rows() != b.size()) {
                assert(false && "DimensionError: least_squares: matrix and vector dimensions mismatch");
                throw DimensionError("least_squares: matrix and vector dimensions mismatch");
            }

            const auto m = a.rows();
            const auto n = a.cols();

            // Solver selection logic
            if (solver_type == SolverType::Auto) {
                if (m >= n) {  // Overdetermined or exactly determined
                    // QR decomposition is numerically stable
                    return qr_solve(a, b);
                }
                else {  // Underdetermined system
                    // SVD decomposition provides the minimum-norm solution
                    return svd_solve(a, b);
                }
            }

            // Explicit solver selection
            switch (solver_type) {
            case SolverType::QR:
                return qr_solve(a, b);
            case SolverType::SVD:
                return svd_solve(a, b);
            default:
                // QR or SVD is typically used for least-squares problems
                return qr_solve(a, b);
            }
        }

        // String-based least-squares solver interface
        template<typename MA, typename VB>
            requires sangi::BaseMatrixLike<MA>
                  && sangi::BaseVectorLike<VB>
                  && std::same_as<sangi::element_t<MA>, sangi::element_t<VB>>
        Vector<sangi::element_t<MA>> least_squares(const MA& a, const VB& b, const std::string& solver_name) {
            return least_squares(a, b, parse_solver_type(solver_name));
        }

        // Weighted least-squares solver
        template<typename MA, typename VB, typename VW>
            requires sangi::BaseMatrixLike<MA>
                  && sangi::BaseVectorLike<VB>
                  && sangi::BaseVectorLike<VW>
                  && std::same_as<sangi::element_t<MA>, sangi::element_t<VB>>
                  && std::same_as<sangi::element_t<MA>, sangi::element_t<VW>>
        Vector<sangi::element_t<MA>> weighted_least_squares(const MA& a, const VB& b,
            const VW& weights) {
            using T = sangi::element_t<MA>;
            SANGI_STATIC_ASSERT_MATRIX_ROWS_EQ_VECTOR_SIZE(MA, VB);
            SANGI_STATIC_ASSERT_MATRIX_ROWS_EQ_VECTOR_SIZE(MA, VW);

            // Dimension check
            if (a.rows() != b.size() || a.rows() != weights.size()) {
                assert(false && "DimensionError: weighted_least_squares: dimensions mismatch");
                throw DimensionError("weighted_least_squares: dimensions mismatch");
            }

            const auto m = a.rows();
            const auto n = a.cols();

            // Transform the weight matrix (diagonal) and the data
            Matrix<T> a_weighted(m, n);
            Vector<T> b_weighted(m);

            for (typename Matrix<T>::size_type i = 0; i < m; ++i) {
                // Zero-weight check (for numerical stabilization)
                const T weight = std::max(weights[i], std::numeric_limits<T>::epsilon());
                const T sqrt_weight = std::sqrt(weight);

                b_weighted[i] = b[i] * sqrt_weight;

                for (typename Matrix<T>::size_type j = 0; j < n; ++j) {
                    a_weighted(i, j) = a(i, j) * sqrt_weight;
                }
            }

            // Solve as an ordinary least-squares problem
            return least_squares(a_weighted, b_weighted);
        }

        // Regularized least-squares solver (ridge regression)
        template<typename MA, typename VB>
            requires sangi::BaseMatrixLike<MA>
                  && sangi::BaseVectorLike<VB>
                  && std::same_as<sangi::element_t<MA>, sangi::element_t<VB>>
        Vector<sangi::element_t<MA>> regularized_least_squares(const MA& a, const VB& b,
            sangi::element_t<MA> lambda = static_cast<sangi::element_t<MA>>(0.1)) {
            using T = sangi::element_t<MA>;
            SANGI_STATIC_ASSERT_MATRIX_ROWS_EQ_VECTOR_SIZE(MA, VB);

            // Dimension check
            if (a.rows() != b.size()) {
                assert(false && "DimensionError: regularized_least_squares: dimensions mismatch");
                throw DimensionError("regularized_least_squares: dimensions mismatch");
            }

            const auto m = a.rows();
            const auto n = a.cols();

            // Normal equations: (A^T * A + lambda * I) * x = A^T * b

            // Compute A^T * A
            Matrix<T> ata(n, n);
            for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                for (typename Matrix<T>::size_type j = 0; j < n; ++j) {
                    ata(i, j) = 0;
                    for (typename Matrix<T>::size_type k = 0; k < m; ++k) {
                        ata(i, j) += a(k, i) * a(k, j);
                    }
                }
            }

            // Add the regularization term to the diagonal
            for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                ata(i, i) += lambda;
            }

            // Compute A^T * b
            Vector<T> atb(n);
            for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                atb[i] = 0;
                for (typename Matrix<T>::size_type j = 0; j < m; ++j) {
                    atb[i] += a(j, i) * b[j];
                }
            }

            // Solve the normal equations
            return cholesky_solve(ata, atb);
        }

        //-----------------------------------------------------------------------------
        // Solvers for special matrices
        //-----------------------------------------------------------------------------

        // Tridiagonal matrix solver
        template<typename VA, typename VB, typename VC, typename VD>
            requires sangi::BaseVectorLike<VA>
                  && sangi::BaseVectorLike<VB>
                  && sangi::BaseVectorLike<VC>
                  && sangi::BaseVectorLike<VD>
                  && std::same_as<sangi::element_t<VA>, sangi::element_t<VB>>
                  && std::same_as<sangi::element_t<VA>, sangi::element_t<VC>>
                  && std::same_as<sangi::element_t<VA>, sangi::element_t<VD>>
        Vector<sangi::element_t<VB>> tridiagonal_solve(const VA& a, const VB& b,
            const VC& c, const VD& d) {
            using T = sangi::element_t<VB>;
            // a: sub-diagonal (size = n-1)
            // b: main diagonal (size = n)
            // c: super-diagonal (size = n-1)
            // d: right-hand-side vector (size = n)
            SANGI_STATIC_ASSERT_VECTOR_SIZE_EQ_OTHER_MINUS_ONE(VA, VB);   // a.size == b.size - 1
            SANGI_STATIC_ASSERT_VECTOR_SIZE_EQ_OTHER_MINUS_ONE(VC, VB);   // c.size == b.size - 1
            SANGI_STATIC_ASSERT_VECTOR_SIZES_EQ(VB, VD);                  // b.size == d.size

            const auto n = b.size();

            // Dimension check
            if (a.size() != n - 1 || c.size() != n - 1 || d.size() != n) {
                assert(false && "DimensionError: tridiagonal_solve: dimensions mismatch");
                throw DimensionError("tridiagonal_solve: dimensions mismatch");
            }

            // Implementation of the Thomas algorithm
            Vector<T> c_prime(n - 1);
            Vector<T> d_prime(n);
            Vector<T> x(n);

            // Forward elimination
            c_prime[0] = c[0] / b[0];
            d_prime[0] = d[0] / b[0];

            for (typename Vector<T>::size_type i = 1; i < n - 1; ++i) {
                const T m = T(1) / (b[i] - a[i - 1] * c_prime[i - 1]);
                c_prime[i] = c[i] * m;
                d_prime[i] = (d[i] - a[i - 1] * d_prime[i - 1]) * m;
            }

            d_prime[n - 1] = (d[n - 1] - a[n - 2] * d_prime[n - 2]) /
                (b[n - 1] - a[n - 2] * c_prime[n - 2]);

            // Back substitution
            x[n - 1] = d_prime[n - 1];

            for (typename Vector<T>::size_type i = n - 1; i-- > 0; ) {
                x[i] = d_prime[i] - c_prime[i] * x[i + 1];
            }

            return x;
        }

        // Pentadiagonal matrix solver (arises in finite differences for 4th-order derivatives)
        //
        // | b0 c0 d0  0  0 ... | | x0 |   | f0 |
        // | a0 b1 c1 d1  0 ... | | x1 |   | f1 |
        // | e0 a1 b2 c2 d2 ... | | x2 | = | f2 |
        // |  0 e1 a2 b3 c3 ... | | x3 |   | f3 |
        // | ...                 | | .. |   | .. |
        //
        // e[i]: diagonal -2, a[i]: diagonal -1, b[i]: diagonal, c[i]: diagonal +1, d[i]: diagonal +2
        //
        // Algorithm: LU decomposition with bandwidth 2 (pentadiagonal-specialized Gauss elimination)
        // Complexity: O(n), memory: O(n)
        template<typename T>
        Vector<T> pentadiagonal_solve(
                const BaseVector<T>& e, const BaseVector<T>& a, const BaseVector<T>& b,
                const BaseVector<T>& c, const BaseVector<T>& d, const BaseVector<T>& f) {
            const auto n = b.size();
            if (n == 0) return {};

            if (n == 1) {
                Vector<T> x(1);
                x[0] = f[0] / b[0];
                return x;
            }

            if (n == 2) {
                // 2x2 system: e, d are empty
                Vector<T> x(2);
                T det = b[0] * b[1] - c[0] * a[0];
                x[0] = (b[1] * f[0] - c[0] * f[1]) / det;
                x[1] = (b[0] * f[1] - a[0] * f[0]) / det;
                return x;
            }

            // Dimension check
            if (e.size() != n - 2 || a.size() != n - 1 ||
                c.size() != n - 1 || d.size() != n - 2 || f.size() != n) {
                assert(false && "DimensionError: pentadiagonal_solve: dimensions mismatch ");
                throw DimensionError("pentadiagonal_solve: dimensions mismatch "
                    "(e:" + std::to_string(e.size()) + " a:" + std::to_string(a.size()) +
                    " b:" + std::to_string(b.size()) + " c:" + std::to_string(c.size()) +
                    " d:" + std::to_string(d.size()) + " f:" + std::to_string(f.size()) + ")");
            }

            // Work arrays (upper-triangular part of the LU decomposition + transformed RHS)
            std::vector<T> al(n - 1);  // sub-diagonal of L
            std::vector<T> be(n);      // diagonal of U
            std::vector<T> ga(n - 1);  // super-diagonal +1 of U
            std::vector<T> de(n - 2);  // super-diagonal +2 of U
            std::vector<T> rhs(n);     // transformed right-hand side

            // Initialization
            be[0] = b[0];
            ga[0] = c[0];
            de[0] = d[0];
            rhs[0] = f[0];

            // i = 1: e[0] is involved
            T mu = a[0] / be[0];
            be[1] = b[1] - mu * ga[0];
            ga[1] = c[1] - mu * de[0];
            if (n > 3) de[1] = d[1];  // de[1] stays equal to d[1] (mu * 0 = 0)
            rhs[1] = f[1] - mu * rhs[0];
            al[0] = mu;

            // i = 2 .. n-1: forward elimination (eliminating 2 rows)
            for (size_t i = 2; i < n; ++i) {
                // Eliminate diagonal -2 (e[i-2]) of row i
                T mu1 = e[i - 2] / be[i - 2];
                // Updated a'[i-1], b'[i], c'[i], d'[i], f'[i]
                T a_new = a[i - 1] - mu1 * ga[i - 2];
                T b_new = b[i] - mu1 * ((i - 2 < n - 2) ? de[i - 2] : T(0));
                T f_new = f[i] - mu1 * rhs[i - 2];

                // Eliminate diagonal -1 (a_new) of row i
                T mu2 = a_new / be[i - 1];
                be[i] = b_new - mu2 * ga[i - 1];
                if (i < n - 1) {
                    T c_i = (i < c.size()) ? c[i] : T(0);
                    ga[i] = c_i - mu2 * ((i - 1 < n - 2) ? de[i - 1] : T(0));
                }
                if (i < n - 2) {
                    de[i] = d[i];
                }
                rhs[i] = f_new - mu2 * rhs[i - 1];
                if (i < n - 1) al[i - 1] = mu2;
            }

            // Back substitution
            Vector<T> x(n);
            x[n - 1] = rhs[n - 1] / be[n - 1];

            if (n >= 2) {
                x[n - 2] = (rhs[n - 2] - ga[n - 2] * x[n - 1]) / be[n - 2];
            }

            if (n >= 3) {
                for (size_t i = n - 3; ; --i) {
                    x[i] = (rhs[i] - ga[i] * x[i + 1] - de[i] * x[i + 2]) / be[i];
                    if (i == 0) break;
                }
            }

            return x;
        }

        // Band matrix solver
/**
 * @brief Improved band matrix solver
 * @tparam T Element type
 * @param a Band matrix
 * @param b Right-hand-side vector
 * @param kl Lower bandwidth (number of non-zero elements below the diagonal)
 * @param ku Upper bandwidth (number of non-zero elements above the diagonal)
 * @return Solution vector
 */
        template<typename MA, typename VB>
            requires sangi::BaseMatrixLike<MA>
                  && sangi::BaseVectorLike<VB>
                  && std::same_as<sangi::element_t<MA>, sangi::element_t<VB>>
        Vector<sangi::element_t<MA>> band_solve(const MA& a, const VB& b,
            typename Matrix<sangi::element_t<MA>>::size_type kl,
            typename Matrix<sangi::element_t<MA>>::size_type ku) {
            using T = sangi::element_t<MA>;
            SANGI_STATIC_ASSERT_SQUARE(MA);
            SANGI_STATIC_ASSERT_MATRIX_ROWS_EQ_VECTOR_SIZE(MA, VB);

            const auto n = a.rows();

            // Dimension check
            if (n != a.cols() || n != b.size()) {
                assert(false && "DimensionError: band_solve: dimensions mismatch");
                throw DimensionError("band_solve: dimensions mismatch");
            }

            // Special handling for tridiagonal matrices (Thomas algorithm)
            if (kl == 1 && ku == 1) {
                Vector<T> diag_lower(n - 1);
                Vector<T> diag_main(n);
                Vector<T> diag_upper(n - 1);

                for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                    diag_main[i] = a(i, i);
                    if (i < n - 1) {
                        diag_upper[i] = a(i, i + 1);
                        diag_lower[i] = a(i + 1, i);
                    }
                }

                // Add scaling to improve stability
                T scale = T(0);
                for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                    scale = std::max(scale, std::abs(diag_main[i]));
                    if (i < n - 1) {
                        scale = std::max(scale, std::abs(diag_upper[i]));
                        scale = std::max(scale, std::abs(diag_lower[i]));
                    }
                }

                // Special handling when the scaling factor is close to zero
                if (scale < std::numeric_limits<T>::epsilon()) {
                    return Vector<T>(n, T(0));  // Return an all-zero solution
                }

                // Compute with the scaled diagonal elements
                Vector<T> scaled_lower(n - 1);
                Vector<T> scaled_main(n);
                Vector<T> scaled_upper(n - 1);
                Vector<T> scaled_b(n);

                for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                    scaled_main[i] = diag_main[i] / scale;
                    scaled_b[i] = b[i] / scale;
                    if (i < n - 1) {
                        scaled_upper[i] = diag_upper[i] / scale;
                        scaled_lower[i] = diag_lower[i] / scale;
                    }
                }

                // Improved version of the Thomas algorithm (better numerical stability)
                return tridiagonal_solve(scaled_lower, scaled_main, scaled_upper, scaled_b);
            }

            // LU decomposition specialized for band matrices (partial pivoting)
            // Partial pivoting expands the bandwidth of U up to kl+ku (fill-in)
            Matrix<T> lu(n, n, T(0));

            // Copy the original band matrix
            for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                const auto j_start = (i > kl) ? i - kl : typename Matrix<T>::size_type(0);
                const auto j_end = std::min(n, i + ku + 1);
                for (typename Matrix<T>::size_type j = j_start; j < j_end; ++j) {
                    lu(i, j) = a(i, j);
                }
            }

            // Pivot record: pivots[k] = the row swapped at step k
            std::vector<typename Matrix<T>::size_type> pivots(n);
            for (typename Matrix<T>::size_type i = 0; i < n; ++i)
                pivots[i] = i;

            // LU decomposition (partial pivoting)
            for (typename Matrix<T>::size_type k = 0; k < n - 1; ++k) {
                // Find the largest pivot in column k from row k onward (within bandwidth kl)
                typename Matrix<T>::size_type pivot_row = k;
                T max_val = std::abs(lu(k, k));
                const auto search_end = std::min(n, k + kl + 1);

                for (typename Matrix<T>::size_type i = k + 1; i < search_end; ++i) {
                    T abs_val = std::abs(lu(i, k));
                    if (abs_val > max_val) {
                        max_val = abs_val;
                        pivot_row = i;
                    }
                }

                pivots[k] = pivot_row;

                // Row swap (within the band is sufficient, but swap all columns to account for fill-in)
                if (pivot_row != k) {
                    const auto swap_start = (k > kl) ? k - kl : typename Matrix<T>::size_type(0);
                    const auto swap_end = std::min(n, k + kl + ku + 1);
                    for (typename Matrix<T>::size_type j = swap_start; j < swap_end; ++j)
                        std::swap(lu(k, j), lu(pivot_row, j));
                }

                // If the pivot is nearly zero, the matrix is singular
                if (std::abs(lu(k, k)) < std::numeric_limits<T>::epsilon() * T(100)) {
                    continue;
                }

                // Elimination (accounting for fill-in: bandwidth of U is at most kl+ku)
                for (typename Matrix<T>::size_type i = k + 1; i < search_end; ++i) {
                    lu(i, k) /= lu(k, k);
                    const auto j_end = std::min(n, k + kl + ku + 1);
                    for (typename Matrix<T>::size_type j = k + 1; j < j_end; ++j) {
                        lu(i, j) -= lu(i, k) * lu(k, j);
                    }
                }
            }

            // Compute the solution: Ly = Pb, Ux = y
            Vector<T> x = b;

            // Apply the pivot permutation to b in order
            for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                if (pivots[i] != i)
                    std::swap(x[i], x[pivots[i]]);
            }

            // Forward substitution (bandwidth of L is kl)
            for (typename Matrix<T>::size_type i = 1; i < n; ++i) {
                const auto j_start = (i > kl) ? i - kl : typename Matrix<T>::size_type(0);
                for (typename Matrix<T>::size_type j = j_start; j < i; ++j) {
                    x[i] -= lu(i, j) * x[j];
                }
            }

            // Back substitution (bandwidth of U is kl+ku: accounting for fill-in)
            for (typename Matrix<T>::size_type i = n; i-- > 0; ) {
                const auto j_end = std::min(n, i + kl + ku + 1);
                for (typename Matrix<T>::size_type j = i + 1; j < j_end; ++j) {
                    x[i] -= lu(i, j) * x[j];
                }
                x[i] /= lu(i, i);
            }

            return x;
        }

        //-----------------------------------------------------------------------------
        // Sparse matrix solvers
        //-----------------------------------------------------------------------------

        // Sparse matrix solver in CSR format
/**
 * @brief Improved sparse matrix solver in CSR format
 * @tparam T Element type
 * @param values Values of the non-zero elements in CSR format
 * @param row_ptr Row pointer array
 * @param col_idx Column index array
 * @param b Right-hand-side vector
 * @param tolerance Convergence threshold
 * @param max_iterations Maximum number of iterations
 * @return Solution vector
 */
        // Forward declaration: BiCGSTAB method (delegated from sparse_solve)
        template<typename T>
        Vector<T> sparse_bicgstab(const std::vector<T>& values,
            const std::vector<typename Matrix<T>::size_type>& row_ptr,
            const std::vector<typename Matrix<T>::size_type>& col_idx,
            const BaseVector<T>& b,
            T tolerance = std::numeric_limits<T>::epsilon() * 10,
            typename Matrix<T>::size_type max_iterations = 0);

        template<typename T>
        Vector<T> sparse_solve(const std::vector<T>& values,
            const std::vector<typename Matrix<T>::size_type>& row_ptr,
            const std::vector<typename Matrix<T>::size_type>& col_idx,
            const BaseVector<T>& b,
            T tolerance = std::numeric_limits<T>::epsilon() * T(100),
            typename Matrix<T>::size_type max_iterations = 0) {
            // General-purpose sparse matrix solver: delegates to the BiCGSTAB method (also handles non-symmetric matrices)
            return sparse_bicgstab(values, row_ptr, col_idx, b, tolerance, max_iterations);
        }

        // Conjugate gradient method for symmetric sparse matrices in CSR format
        template<typename T>
        Vector<T> sparse_cg(const std::vector<T>& values,
            const std::vector<typename Matrix<T>::size_type>& row_ptr,
            const std::vector<typename Matrix<T>::size_type>& col_idx,
            const BaseVector<T>& b,
            T tolerance = std::numeric_limits<T>::epsilon() * 10,
            typename Matrix<T>::size_type max_iterations = 0) {
            // Conjugate gradient method for sparse matrices
            const auto n = b.size();

            // Dimension check
            if (row_ptr.size() != n + 1) {
                assert(false && "DimensionError: sparse_cg: row_ptr size mismatch");
                throw DimensionError("sparse_cg: row_ptr size mismatch");
            }

            // Default maximum number of iterations
            if (max_iterations == 0) {
                max_iterations = 10 * n;
            }

            Vector<T> x(n, T(0));     // Set the initial solution to zero

            // Compute the norm of the right-hand-side vector
            const T b_norm = std::sqrt(dot(b, b));

            // Zero-vector check
            if (b_norm < std::numeric_limits<T>::epsilon()) {
                return x;  // If the right-hand side is zero, the zero vector is the solution
            }

            // Initial residual r = b - A*x = b (x is zero)
            Vector<T> r = b;
            Vector<T> p = r;          // Initial search direction
            Vector<T> ap(n);          // A*p buffer (allocated once outside the loop)

            T r_norm_squared = dot(r, r);   // Residual norm (squared)
            // Relative tolerance (squared): ‖r‖² ≤ (tolerance * ‖b‖)²
            const T rel_tol = tolerance * tolerance * b_norm * b_norm;

            // Iteration
            for (typename Matrix<T>::size_type iter = 0; iter < max_iterations; ++iter) {
                // Convergence check
                if (r_norm_squared <= rel_tol) {
                    break;
                }

                // Compute A*p (CSR format, ap already allocated outside the loop)
                ap.assign(n, T(0));
                for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                    for (typename Matrix<T>::size_type j = row_ptr[i]; j < row_ptr[i + 1]; ++j) {
                        ap[i] += values[j] * p[col_idx[j]];
                    }
                }

                // Step size α = r^T*r / (p^T*A*p)
                const T p_ap = dot(p, ap);

                // When p^T*A*p is close to zero (numerically unstable)
                if (std::abs(p_ap) < std::numeric_limits<T>::epsilon()) {
                    break;
                }

                const T alpha = r_norm_squared / p_ap;

                // Update the solution and residual (fused function)
                axpy(alpha, p, x);
                axpy(-alpha, ap, r);

                // New residual norm (squared)
                const T r_next_norm_squared = dot(r, r);

                // Update the direction vector: p = r + beta * p
                if (r_norm_squared < std::numeric_limits<T>::epsilon()) break;  // Avoid division by zero
                const T beta = r_next_norm_squared / r_norm_squared;
                axpby(T(1), r, beta, p);

                // Update the residual norm
                r_norm_squared = r_next_norm_squared;
            }

            return x;
        }

        // BiCGSTAB method for sparse matrices
        template<typename T>
        Vector<T> sparse_bicgstab(const std::vector<T>& values,
            const std::vector<typename Matrix<T>::size_type>& row_ptr,
            const std::vector<typename Matrix<T>::size_type>& col_idx,
            const BaseVector<T>& b,
            T tolerance,
            typename Matrix<T>::size_type max_iterations) {
            // BiCGSTAB method (stabilized biconjugate gradient method)
            const auto n = b.size();

            // Dimension check
            if (row_ptr.size() != n + 1) {
                assert(false && "DimensionError: sparse_bicgstab: row_ptr size mismatch");
                throw DimensionError("sparse_bicgstab: row_ptr size mismatch");
            }

            // Default maximum number of iterations
            if (max_iterations == 0) {
                max_iterations = 10 * n;
            }

            Vector<T> x(n, T(0));     // Set the initial solution to zero

            // Function for the CSR-format matrix-vector product
            auto csr_mv = [&](const Vector<T>& v, Vector<T>& result) {
                result.assign(n, T(0));
                for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                    for (typename Matrix<T>::size_type j = row_ptr[i]; j < row_ptr[i + 1]; ++j) {
                        result[i] += values[j] * v[col_idx[j]];
                    }
                }
                };

            // Initial residual r = b - A*x = b (x is zero)
            Vector<T> r = b;
            Vector<T> r_hat = r;  // Shadow factor

            // Compute the norm of the right-hand-side vector
            const T b_norm = std::sqrt(dot(b, b));

            // Zero-vector check
            if (b_norm < std::numeric_limits<T>::epsilon()) {
                return x;  // If the right-hand side is zero, the zero vector is the solution
            }

            const T rel_tol = tolerance * b_norm; // Relative tolerance

            // BiCGSTAB initialization
            T rho_prev = 1;
            T alpha = 1;
            T omega = 1;

            Vector<T> p(n, T(0));
            Vector<T> v(n, T(0));
            Vector<T> s(n);
            Vector<T> t(n);

            T rho = dot(r_hat, r);

            // Iteration
            for (typename Matrix<T>::size_type iter = 0; iter < max_iterations; ++iter) {
                // Convergence check
                const T r_norm = std::sqrt(dot(r, r));
                if (r_norm <= rel_tol) {
                    break;
                }

                // Compute ρ
                rho = dot(r_hat, r);

                // When ρ is close to zero (numerically unstable)
                if (std::abs(rho) < std::numeric_limits<T>::epsilon()) {
                    break;
                }

                const T beta = (rho / rho_prev) * (alpha / omega);

                // Update p: p = r + beta * (p - omega * v)
                {
                    auto* pp = p.data();
                    const auto* rp = r.data();
                    const auto* vp = v.data();
                    for (typename Matrix<T>::size_type i = 0; i < n; ++i)
                        pp[i] = rp[i] + beta * (pp[i] - omega * vp[i]);
                }

                // v = A*p
                csr_mv(p, v);

                // Compute α
                const T r_hat_v = dot(r_hat, v);

                // When r_hat_v is close to zero (numerically unstable)
                if (std::abs(r_hat_v) < std::numeric_limits<T>::epsilon()) {
                    break;
                }

                alpha = rho / r_hat_v;

                // s = r - α*v
                axpby(T(1), r, -alpha, v, s);

                // Convergence check
                const T s_norm = std::sqrt(dot(s, s));
                if (s_norm <= rel_tol) {
                    // x += α*p
                    axpy(alpha, p, x);
                    break;
                }

                // t = A*s
                csr_mv(s, t);

                // Compute ω
                const T t_t = dot(t, t);

                // When t_t is close to zero (numerically unstable)
                if (std::abs(t_t) < std::numeric_limits<T>::epsilon()) {
                    break;
                }

                omega = dot(t, s) / t_t;

                // Update x, r (fused loop)
                {
                    auto* xp = x.data();
                    const auto* pp = p.data();
                    const auto* sp = s.data();
                    auto* rp = r.data();
                    const auto* tp = t.data();
                    for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                        xp[i] += alpha * pp[i] + omega * sp[i];
                        rp[i] = sp[i] - omega * tp[i];
                    }
                }

                // When ω is close to zero (numerically unstable)
                if (std::abs(omega) < std::numeric_limits<T>::epsilon()) {
                    break;
                }

                // Update ρ_prev
                rho_prev = rho;
            }

            return x;
        }

        // ================================================================
        // Iterative Refinement
        // ================================================================

        /// Obtain a high-precision solution via LU decomposition + iterative refinement
        ///
        /// See: based on Okumura "Encyclopedia of Algorithms in C" invr.c (within regress.c),
        ///       templatized and integrated into the sangi API
        ///
        /// Principle:
        ///   1. Obtain an initial solution x₀ via LU decomposition
        ///   2. Compute the residual r = b - Ax with high precision
        ///   3. Use the LU decomposition to obtain d = A⁻¹r (no additional decomposition needed)
        ///   4. Update the solution by x ← x + d
        ///   5. Repeat until the residual norm is sufficiently small
        ///
        /// Effective for matrices with large condition numbers. Even when computing in double,
        /// it improves to near machine precision provided the residual can be computed accurately.
        ///
        /// @param A Coefficient matrix (n×n)
        /// @param b Right-hand-side vector
        /// @param maxIter Maximum number of iterations (typically 2-5 is enough)
        /// @param tol Relative tolerance for the residual norm
        /// @return Improved solution vector
        template<typename T>
        Vector<T> iterativeRefinement(
            const BaseMatrix<T>& A, const BaseVector<T>& b,
            size_t maxIter = 10,
            T tol = std::numeric_limits<T>::epsilon() * T(10))
        {
            if (A.rows() != A.cols()) {
                assert(false && "DimensionError: iterativeRefinement: matrix must be square");
                throw DimensionError("iterativeRefinement: matrix must be square");
            }
            if (A.rows() != b.size()) {
                assert(false && "DimensionError: iterativeRefinement: matrix and vector dimensions mismatch");
                throw DimensionError("iterativeRefinement: matrix and vector dimensions mismatch");
            }

            const auto n = A.rows();
            if (n == 0) return Vector<T>();

            // LU decomposition (performed once and reused across all iterations)
            auto [lu, pivots] = lu_decomposition(A);

            // Initial solution
            Vector<T> x = condest_detail::lu_solve_factored(lu, pivots, b);

            T bNorm = T(0);
            for (size_t i = 0; i < n; ++i) bNorm += std::abs(b[i]);
            if (bNorm == T(0)) return x;

            for (size_t iter = 0; iter < maxIter; ++iter) {
                // Compute the residual r = b - Ax
                // This is the key to precision: compute it as accurately as possible
                Vector<T> r(n);
                for (size_t i = 0; i < n; ++i) {
                    T s = b[i];
                    for (size_t j = 0; j < n; ++j) {
                        s -= A(i, j) * x[j];
                    }
                    r[i] = s;
                }

                // Check the residual norm
                T rNorm = T(0);
                for (size_t i = 0; i < n; ++i) rNorm += std::abs(r[i]);

                if (rNorm <= tol * bNorm) break;

                // Correction vector d = A⁻¹r (using the same LU decomposition)
                Vector<T> d = condest_detail::lu_solve_factored(lu, pivots, r);

                // Update the solution: x ← x + d
                for (size_t i = 0; i < n; ++i) {
                    x[i] += d[i];
                }
            }

            return x;
        }

        /// Result of iterative refinement (with detailed information)
        template<typename T>
        struct RefinementResult {
            Vector<T> solution;      // Improved solution
            T residualNorm;          // Final residual norm
            T conditionEstimate;     // Estimate of the condition number
            size_t iterations;       // Number of iterations performed
            bool converged;          // Whether it converged
        };

        /// Iterative refinement (with detailed information)
        template<typename T>
        RefinementResult<T> iterativeRefinementDetailed(
            const BaseMatrix<T>& A, const BaseVector<T>& b,
            size_t maxIter = 10,
            T tol = std::numeric_limits<T>::epsilon() * T(10))
        {
            if (A.rows() != A.cols()) {
                assert(false && "DimensionError: iterativeRefinementDetailed: matrix must be square");
                throw DimensionError("iterativeRefinementDetailed: matrix must be square");
            }
            if (A.rows() != b.size()) {
                assert(false && "DimensionError: iterativeRefinementDetailed: dimensions mismatch");
                throw DimensionError("iterativeRefinementDetailed: dimensions mismatch");
            }

            const auto n = A.rows();
            if (n == 0) return {Vector<T>(), T(0), T(0), 0, true};

            auto [lu, pivots] = lu_decomposition(A);
            Vector<T> x = condest_detail::lu_solve_factored(lu, pivots, b);

            T bNorm = T(0);
            for (size_t i = 0; i < n; ++i) bNorm += std::abs(b[i]);
            if (bNorm == T(0)) return {x, T(0), T(0), 0, true};

            T rNorm = T(0);
            size_t iter = 0;
            bool converged = false;

            for (iter = 0; iter < maxIter; ++iter) {
                Vector<T> r(n);
                for (size_t i = 0; i < n; ++i) {
                    T s = b[i];
                    for (size_t j = 0; j < n; ++j) s -= A(i, j) * x[j];
                    r[i] = s;
                }

                rNorm = T(0);
                for (size_t i = 0; i < n; ++i) rNorm += std::abs(r[i]);

                if (rNorm <= tol * bNorm) {
                    converged = true;
                    break;
                }

                Vector<T> d = condest_detail::lu_solve_factored(lu, pivots, r);
                for (size_t i = 0; i < n; ++i) x[i] += d[i];
            }

            // Condition number estimation
            T condEst = estimate_inv_norm1(lu, pivots);
            T norm1A = T(0);
            for (size_t j = 0; j < n; ++j) {
                T colSum = T(0);
                for (size_t i = 0; i < n; ++i)
                    colSum += std::abs(A(i, j));
                if (colSum > norm1A) norm1A = colSum;
            }
            condEst *= norm1A;

            return {std::move(x), rNorm / bNorm, condEst, iter, converged};
        }

    } // namespace algorithms
} // namespace sangi

#endif // SANGI_SOLVERS_HPP