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

/**
 * @file iterative_solvers.hpp
 * @brief Iterative solver API extension
 *
 * Additional solvers:
 *   - MINRES (for symmetric indefinite matrices)
 *   - FGMRES (Flexible GMRES with variable preconditioning)
 *
 * Unified API:
 *   - IterativeSolver<T> : builder-pattern solver configuration
 *   - SolverLog<T>       : convergence history recording
 *   - solve()            : automatic solver selection based on matrix properties
 */

#ifndef SANGI_ITERATIVE_SOLVERS_HPP
#define SANGI_ITERATIVE_SOLVERS_HPP

#include <math/core/sparse_matrix_algorithms.hpp>
#include <math/core/sparse_matrix.hpp>
#include <math/core/vector.hpp>
#include <math/core/convergence_criteria.hpp>

#include <optional>
#include <functional>
#include <vector>
#include <cmath>
#include <limits>
#include <stdexcept>
#include <string>
#include <utility>
#include <random>

namespace sangi {

namespace sparse_algorithms {

// ============================================================================
// SolverLog — convergence history recording
// ============================================================================

template<typename T>
struct SolverLog {
    std::vector<T> residual_history;  ///< residual norm at each iteration
    std::vector<T> time_history;      ///< (future use) cumulative time at each iteration

    void clear() {
        residual_history.clear();
        time_history.clear();
    }

    void record(T rnorm) {
        residual_history.push_back(rnorm);
    }

    /// Convergence rate (geometric mean over the last few iterations)
    [[nodiscard]] T convergence_rate(std::size_t window = 5) const {
        auto n = residual_history.size();
        if (n < 2) return T(1);
        auto w = std::min(window, n - 1);
        T ratio = residual_history[n - 1] / residual_history[n - 1 - w];
        return std::pow(std::abs(ratio), T(1) / static_cast<T>(w));
    }
};

// ============================================================================
// MINRES — for symmetric indefinite matrices
// ============================================================================

/**
 * @brief MINRES (Minimum Residual Method)
 *
 * Solves Ax = b for a symmetric (possibly indefinite) sparse matrix A.
 * Whereas CG requires positive definiteness, MINRES also works for indefinite matrices.
 *
 * Implementation: Lanczos tridiagonalization + GMRES-style Givens QR
 * Lanczos tridiagonalizes the symmetric matrix, and the same Givens rotation procedure
 * as GMRES performs the Hessenberg upper-triangular factorization. By symmetry,
 * Hessenberg = tridiagonal, so each column requires only a single rotation
 * (a constant number, versus j times for GMRES).
 */
template<typename T, typename IndexType>
    requires concepts::OrderedField<T> && std::integral<IndexType>
SolverResult<T> minres(
    const SparseMatrix<T, IndexType>& A,
    const Vector<T>& b,
    const ConvergenceCriteria<T>& criteria,
    Preconditioner<T> precond = {},
    std::optional<Vector<T>> x0 = std::nullopt,
    SolverLog<T>* log = nullptr)
{
    const auto n = static_cast<std::size_t>(A.rows());
    if (static_cast<std::size_t>(A.cols()) != n)
        throw std::invalid_argument("MINRES: matrix must be square");
    if (b.size() != n)
        throw std::invalid_argument("MINRES: incompatible dimensions");

    (void)precond;

    Vector<T> x = x0.value_or(Vector<T>(n, T{0}));

    const T b_norm = norm(b);
    const T threshold = criteria.absolute_tolerance() +
                        criteria.relative_tolerance() * b_norm;

    // Lanczos + Givens QR (GMRES approach specialized to tridiagonal)
    const std::size_t max_iter = std::min(criteria.max_iterations(), n);

    Vector<T> r = b - A * x;
    T r_norm = norm(r);
    if (log) log->record(r_norm);
    if (r_norm <= threshold) {
        return {x, r_norm, 0, true};
    }

    // Lanczos basis
    std::vector<Vector<T>> Q(max_iter + 1);
    Q[0] = r / r_norm;

    // Store the tridiagonal matrix in Hessenberg form (2 elements per column: H[j][j], H[j+1][j])
    // However, the previous column is also needed to retain the values before Givens rotation
    // → use the same structure as GMRES (H is (m+1)×m, but band width 2 since tridiagonal)
    std::vector<std::vector<T>> H(max_iter + 1, std::vector<T>(max_iter, T{0}));

    std::vector<T> cs(max_iter, T{0}), sn(max_iter, T{0});
    std::vector<T> g(max_iter + 1, T{0});
    g[0] = r_norm;

    Vector<T> v_prev(n, T{0});
    std::size_t j = 0;
    for (; j < max_iter; ++j) {
        // Lanczos step (= Arnoldi step for a symmetric matrix)
        Vector<T> w = A * Q[j];

        // Modified Gram-Schmidt (by symmetry, only j-1 and j are nonzero)
        if (j > 0) {
            H[j - 1][j] = dot(Q[j - 1], w);  // = beta_j (upper part of the tridiagonal)
            w -= Q[j - 1] * H[j - 1][j];
        }
        H[j][j] = dot(Q[j], w);  // = alpha_j
        w -= Q[j] * H[j][j];

        H[j + 1][j] = norm(w);  // = beta_{j+1}

        // Happy breakdown
        if (std::abs(H[j + 1][j]) < std::numeric_limits<T>::epsilon() * T{100}) {
            // Apply the Givens rotations, then finish
            for (std::size_t i = 0; i < j; ++i) {
                T tmp = cs[i] * H[i][j] + sn[i] * H[i + 1][j];
                H[i + 1][j] = -sn[i] * H[i][j] + cs[i] * H[i + 1][j];
                H[i][j] = tmp;
            }
            T rr = std::sqrt(H[j][j] * H[j][j] + H[j + 1][j] * H[j + 1][j]);
            if (rr > T{0}) {
                cs[j] = H[j][j] / rr;
                sn[j] = H[j + 1][j] / rr;
                H[j][j] = rr;
                H[j + 1][j] = T{0};
                T g_new = -sn[j] * g[j];
                g[j] = cs[j] * g[j];
                g[j + 1] = g_new;
            }
            ++j;
            break;
        }

        Q[j + 1] = w / H[j + 1][j];

        // Apply previous Givens rotations to the H column (tridiagonal → G_{j-2} generates fill-in, so at most 2)
        for (std::size_t i = (j >= 2 ? j - 2 : 0); i < j; ++i) {
            T tmp = cs[i] * H[i][j] + sn[i] * H[i + 1][j];
            H[i + 1][j] = -sn[i] * H[i][j] + cs[i] * H[i + 1][j];
            H[i][j] = tmp;
        }

        // j-th Givens rotation
        T rr = std::sqrt(H[j][j] * H[j][j] + H[j + 1][j] * H[j + 1][j]);
        cs[j] = H[j][j] / rr;
        sn[j] = H[j + 1][j] / rr;
        H[j][j] = rr;
        H[j + 1][j] = T{0};

        T g_new = -sn[j] * g[j];
        g[j] = cs[j] * g[j];
        g[j + 1] = g_new;

        r_norm = std::abs(g[j + 1]);
        if (log) log->record(r_norm);

        if (r_norm <= threshold) {
            ++j;
            break;
        }
    }

    // Back substitution: R y = g
    const std::size_t dim = j;
    std::vector<T> y(dim, T{0});
    for (std::size_t ii = 0; ii < dim; ++ii) {
        std::size_t idx = dim - 1 - ii;
        y[idx] = g[idx];
        for (std::size_t k = idx + 1; k < dim; ++k) {
            y[idx] -= H[idx][k] * y[k];
        }
        if (std::abs(H[idx][idx]) > std::numeric_limits<T>::epsilon()) {
            y[idx] /= H[idx][idx];
        }
    }

    // Solution update: x += Q * y
    for (std::size_t k = 0; k < dim; ++k) {
        x += Q[k] * y[k];
    }

    T final_rnorm = norm(b - A * x);
    return {x, final_rnorm, dim, final_rnorm <= threshold};
}

// ============================================================================
// FGMRES — Flexible GMRES (with variable preconditioning)
// ============================================================================

/**
 * @brief Flexible GMRES (with variable preconditioning)
 *
 * Whereas ordinary GMRES requires the same preconditioner at every iteration,
 * FGMRES can apply a different preconditioner at each iteration.
 * Because the preconditioned basis vectors are stored separately, memory usage is doubled.
 *
 * Reference: Saad (1993), "A Flexible Inner-Outer Preconditioned GMRES Algorithm"
 */
template<typename T, typename IndexType>
    requires concepts::OrderedField<T> && std::integral<IndexType>
SolverResult<T> fgmres(
    const SparseMatrix<T, IndexType>& A,
    const Vector<T>& b,
    const ConvergenceCriteria<T>& criteria,
    std::size_t restart = 30,
    Preconditioner<T> precond = {},
    std::optional<Vector<T>> x0 = std::nullopt,
    SolverLog<T>* log = nullptr)
{
    const auto n = static_cast<std::size_t>(A.rows());
    if (static_cast<std::size_t>(A.cols()) != n)
        throw std::invalid_argument("FGMRES: matrix must be square");
    if (b.size() != n)
        throw std::invalid_argument("FGMRES: incompatible dimensions");

    if (!precond) precond = identity_preconditioner<T>();
    if (restart == 0 || restart > n) restart = n;

    Vector<T> x = x0.value_or(Vector<T>(n, T{0}));

    const T b_norm = norm(b);
    if (b_norm <= criteria.absolute_tolerance()) {
        return {x, b_norm, 0, true};
    }
    const T threshold = criteria.absolute_tolerance() +
                        criteria.relative_tolerance() * b_norm;

    std::size_t total_iter = 0;

    while (total_iter < criteria.max_iterations()) {
        Vector<T> r = b - A * x;
        T r_norm = norm(r);
        if (log) log->record(r_norm);
        if (r_norm <= threshold) {
            return {x, r_norm, total_iter, true};
        }

        const std::size_t m = std::min(restart,
            criteria.max_iterations() - total_iter);

        // Arnoldi basis V and preconditioned basis Z
        std::vector<Vector<T>> V(m + 1);
        std::vector<Vector<T>> Z(m);
        V[0] = r / r_norm;

        // Upper Hessenberg matrix (m+1 x m)
        std::vector<std::vector<T>> H(m + 1, std::vector<T>(m, T{0}));

        // Givens rotation parameters
        std::vector<T> cs(m, T{0}), sn(m, T{0});

        // Transformed right-hand side
        std::vector<T> g(m + 1, T{0});
        g[0] = r_norm;

        std::size_t j = 0;
        for (; j < m; ++j) {
            // FGMRES: z_j = M^{-1} v_j (retain the preconditioned basis)
            Z[j] = precond(V[j]);
            Vector<T> w = A * Z[j];

            // Modified Gram-Schmidt
            for (std::size_t i = 0; i <= j; ++i) {
                H[i][j] = dot(V[i], w);
                w -= V[i] * H[i][j];
            }
            H[j + 1][j] = norm(w);

            // Happy breakdown
            if (std::abs(H[j + 1][j]) < std::numeric_limits<T>::epsilon() * T{100}) {
                // Apply previous Givens rotations
                for (std::size_t i = 0; i < j; ++i) {
                    T tmp = cs[i] * H[i][j] + sn[i] * H[i + 1][j];
                    H[i + 1][j] = -sn[i] * H[i][j] + cs[i] * H[i + 1][j];
                    H[i][j] = tmp;
                }
                T rr = std::sqrt(H[j][j] * H[j][j] + H[j + 1][j] * H[j + 1][j]);
                if (rr > T{0}) {
                    cs[j] = H[j][j] / rr;
                    sn[j] = H[j + 1][j] / rr;
                    H[j][j] = rr;
                    H[j + 1][j] = T{0};
                    T g_new = -sn[j] * g[j];
                    g[j] = cs[j] * g[j];
                    g[j + 1] = g_new;
                }
                ++j;
                break;
            }

            V[j + 1] = w / H[j + 1][j];

            // Apply previous Givens rotations to the H column
            for (std::size_t i = 0; i < j; ++i) {
                T tmp = cs[i] * H[i][j] + sn[i] * H[i + 1][j];
                H[i + 1][j] = -sn[i] * H[i][j] + cs[i] * H[i + 1][j];
                H[i][j] = tmp;
            }

            // j-th Givens rotation
            T rr = std::sqrt(H[j][j] * H[j][j] + H[j + 1][j] * H[j + 1][j]);
            cs[j] = H[j][j] / rr;
            sn[j] = H[j + 1][j] / rr;
            H[j][j] = rr;
            H[j + 1][j] = T{0};

            T g_new = -sn[j] * g[j];
            g[j] = cs[j] * g[j];
            g[j + 1] = g_new;

            r_norm = std::abs(g[j + 1]);
            ++total_iter;
            if (log) log->record(r_norm);

            if (r_norm <= threshold) {
                ++j;
                break;
            }
        }

        // Back substitution
        const std::size_t dim = j;
        std::vector<T> y(dim, T{0});
        for (std::size_t ii = 0; ii < dim; ++ii) {
            std::size_t idx = dim - 1 - ii;
            y[idx] = g[idx];
            for (std::size_t k = idx + 1; k < dim; ++k) {
                y[idx] -= H[idx][k] * y[k];
            }
            if (std::abs(H[idx][idx]) > std::numeric_limits<T>::epsilon()) {
                y[idx] /= H[idx][idx];
            }
        }

        // FGMRES: x += Z * y (using the preconditioned basis)
        for (std::size_t k = 0; k < dim; ++k) {
            x += Z[k] * y[k];
        }

        if (r_norm <= threshold) {
            return {x, r_norm, total_iter, true};
        }
    }

    T final_rnorm = norm(b - A * x);
    return {x, final_rnorm, total_iter, false};
}

// ============================================================================
// Logged CG / BiCGSTAB / GMRES — wrappers around existing solvers
// ============================================================================

/**
 * @brief Logged CG
 */
template<typename T, typename IndexType>
    requires concepts::OrderedField<T> && std::integral<IndexType>
SolverResult<T> conjugate_gradient_logged(
    const SparseMatrix<T, IndexType>& A,
    const Vector<T>& b,
    const ConvergenceCriteria<T>& criteria,
    SolverLog<T>& log,
    Preconditioner<T> precond = {},
    std::optional<Vector<T>> x0 = std::nullopt)
{
    const auto n = A.rows();
    if (A.cols() != n) throw std::invalid_argument("CG: matrix must be square");
    if (b.size() != n) throw std::invalid_argument("CG: incompatible dimensions");

    if (!precond) precond = identity_preconditioner<T>();

    Vector<T> x = x0.value_or(Vector<T>(n, T{0}));
    Vector<T> r = b - A * x;

    const T init_rnorm = norm(r);
    log.record(init_rnorm);
    if (init_rnorm <= criteria.absolute_tolerance()) {
        return {x, init_rnorm, 0, true};
    }

    Vector<T> z = precond(r);
    Vector<T> p = z;
    T rz = dot(r, z);

    auto is_converged = [&](T rnorm) {
        return rnorm <= criteria.absolute_tolerance() ||
               rnorm <= criteria.relative_tolerance() * init_rnorm;
    };

    for (std::size_t iter = 0; iter < criteria.max_iterations(); ++iter) {
        Vector<T> Ap = A * p;
        T pAp = dot(p, Ap);
        if (std::abs(pAp) < std::numeric_limits<T>::epsilon()) {
            T rnorm = norm(b - A * x);
            return {x, rnorm, iter, is_converged(rnorm)};
        }
        T alpha = rz / pAp;
        x += p * alpha;
        r -= Ap * alpha;

        T rnorm = norm(r);
        log.record(rnorm);
        if (is_converged(rnorm)) {
            return {x, rnorm, iter + 1, true};
        }

        z = precond(r);
        T rz_new = dot(r, z);
        T beta = rz_new / rz;
        rz = rz_new;
        p = z + p * beta;
    }

    T final_rnorm = norm(b - A * x);
    return {x, final_rnorm, criteria.max_iterations(), is_converged(final_rnorm)};
}

/**
 * @brief Logged BiCGSTAB
 */
template<typename T, typename IndexType>
    requires concepts::OrderedField<T> && std::integral<IndexType>
SolverResult<T> bicgstab_logged(
    const SparseMatrix<T, IndexType>& A,
    const Vector<T>& b,
    const ConvergenceCriteria<T>& criteria,
    SolverLog<T>& log,
    Preconditioner<T> precond = {},
    std::optional<Vector<T>> x0 = std::nullopt)
{
    const auto n = A.rows();
    if (A.cols() != n) throw std::invalid_argument("BiCGSTAB: matrix must be square");
    if (b.size() != n) throw std::invalid_argument("BiCGSTAB: incompatible dimensions");

    if (!precond) precond = identity_preconditioner<T>();

    Vector<T> x = x0.value_or(Vector<T>(n, T{0}));
    Vector<T> r = b - A * x;
    Vector<T> r_hat = r;

    const T init_rnorm = norm(r);
    log.record(init_rnorm);
    if (init_rnorm <= criteria.absolute_tolerance()) {
        return {x, init_rnorm, 0, true};
    }

    T rho = T{1}, alpha = T{1}, omega_val = T{1};
    Vector<T> p(n, T{0}), v(n, T{0});
    const T eps = std::numeric_limits<T>::epsilon();
    const T breakdown_tol = eps * init_rnorm * init_rnorm;

    auto is_converged = [&](T rnorm) {
        return rnorm <= criteria.absolute_tolerance() ||
               rnorm <= criteria.relative_tolerance() * init_rnorm;
    };

    for (std::size_t iter = 0; iter < criteria.max_iterations(); ++iter) {
        T rho_new = dot(r_hat, r);
        if (std::abs(rho_new) < breakdown_tol) {
            T rnorm = norm(b - A * x);
            return {x, rnorm, iter, is_converged(rnorm)};
        }

        if (iter == 0) {
            p = r;
        } else {
            T beta = (rho_new / rho) * (alpha / omega_val);
            p = r + (p - v * omega_val) * beta;
        }
        rho = rho_new;

        Vector<T> p_hat = precond(p);
        v = A * p_hat;
        T denom = dot(r_hat, v);
        if (std::abs(denom) < breakdown_tol) {
            T rnorm = norm(b - A * x);
            return {x, rnorm, iter, is_converged(rnorm)};
        }
        alpha = rho / denom;

        Vector<T> s = r - v * alpha;
        T s_norm = norm(s);
        if (is_converged(s_norm)) {
            x += p_hat * alpha;
            log.record(s_norm);
            return {x, s_norm, iter + 1, true};
        }

        Vector<T> s_hat = precond(s);
        Vector<T> t = A * s_hat;
        T tt = dot(t, t);
        if (std::abs(tt) < breakdown_tol) {
            x += p_hat * alpha;
            T rnorm = norm(b - A * x);
            return {x, rnorm, iter + 1, is_converged(rnorm)};
        }
        omega_val = dot(t, s) / tt;
        x += p_hat * alpha + s_hat * omega_val;
        r = s - t * omega_val;

        T rnorm = norm(r);
        log.record(rnorm);
        if (is_converged(rnorm)) {
            return {x, rnorm, iter + 1, true};
        }

        if (std::abs(omega_val) < eps * init_rnorm) {
            T true_rnorm = norm(b - A * x);
            return {x, true_rnorm, iter + 1, is_converged(true_rnorm)};
        }
    }

    T final_rnorm = norm(b - A * x);
    return {x, final_rnorm, criteria.max_iterations(), is_converged(final_rnorm)};
}

// ============================================================================
// IDR(s) — Induced Dimension Reduction (nonsymmetric, transpose-free)
// ============================================================================

/**
 * @brief IDR(s) (Sonneveld & van Gijzen 2008)
 *
 * A transpose-free Krylov method that solves Ax=b for a nonsymmetric sparse matrix A.
 * Using an s-dimensional shadow space P, it successively pushes the residual into the
 * nested subspaces G_j (= (I−ω_jA)(G_{j−1}∩P^⊥)).
 * It generalizes BiCGSTAB; raising s yields convergence closer to GMRES, while
 * matrix-vector products and memory grow only in proportion to s. Supports right preconditioning.
 *
 * @param s dimension of the shadow space (default 4). Larger gives faster convergence but more memory.
 */
template<typename T, typename IndexType>
    requires concepts::OrderedField<T> && std::integral<IndexType>
SolverResult<T> idrs(
    const SparseMatrix<T, IndexType>& A,
    const Vector<T>& b,
    const ConvergenceCriteria<T>& criteria,
    std::size_t s = 4,
    Preconditioner<T> precond = {},
    std::optional<Vector<T>> x0 = std::nullopt)
{
    const auto n = A.rows();
    if (A.cols() != n) throw std::invalid_argument("IDR(s): matrix must be square");
    if (b.size() != n) throw std::invalid_argument("IDR(s): incompatible dimensions");
    if (s < 1) s = 1;
    if (static_cast<std::size_t>(n) < s) s = static_cast<std::size_t>(n);
    if (!precond) precond = identity_preconditioner<T>();

    Vector<T> x = x0.value_or(Vector<T>(n, T{0}));
    Vector<T> r = b - A * x;
    const T init_rnorm = norm(r);
    if (init_rnorm <= criteria.absolute_tolerance())
        return {x, init_rnorm, 0, true};

    auto is_converged = [&](T rn) {
        return rn <= criteria.absolute_tolerance() ||
               rn <= criteria.relative_tolerance() * init_rnorm;
    };
    const T eps = std::numeric_limits<T>::epsilon();

    // shadow space P (s vectors): fixed-seed random numbers + modified Gram-Schmidt orthonormalization
    std::vector<Vector<T>> P(s, Vector<T>(n, T{0}));
    std::mt19937_64 rng(0x1D125ULL);
    std::uniform_real_distribution<double> dist(-1.0, 1.0);
    for (std::size_t k = 0; k < s; ++k) {
        for (std::size_t i = 0; i < static_cast<std::size_t>(n); ++i)
            P[k][i] = static_cast<T>(dist(rng));
        for (std::size_t j = 0; j < k; ++j)
            P[k] -= P[j] * dot(P[j], P[k]);
        T nrm = norm(P[k]);
        if (nrm > eps) P[k] *= (T(1) / nrm);
        else P[k][k % static_cast<std::size_t>(n)] = T(1);
    }

    std::vector<Vector<T>> G(s, Vector<T>(n, T{0}));
    std::vector<Vector<T>> U(s, Vector<T>(n, T{0}));
    std::vector<T> M(s * s, T{0});
    for (std::size_t i = 0; i < s; ++i) M[i * s + i] = T(1);
    auto Mref = [&](std::size_t i, std::size_t j) -> T& { return M[i * s + j]; };
    T om = T(1);

    std::size_t iter = 0;
    while (iter < criteria.max_iterations()) {
        std::vector<T> f(s);
        for (std::size_t i = 0; i < s; ++i) f[i] = dot(P[i], r);

        for (std::size_t k = 0; k < s; ++k) {
            // Solve the small system M(k:s,k:s) c = f(k:s) by Gaussian elimination with partial pivoting
            const std::size_t m = s - k;
            std::vector<T> As(m * m), c(m);
            for (std::size_t a = 0; a < m; ++a) {
                c[a] = f[k + a];
                for (std::size_t bb = 0; bb < m; ++bb) As[a * m + bb] = Mref(k + a, k + bb);
            }
            for (std::size_t col = 0; col < m; ++col) {
                std::size_t pv = col;
                for (std::size_t rr = col + 1; rr < m; ++rr)
                    if (std::abs(As[rr * m + col]) > std::abs(As[pv * m + col])) pv = rr;
                if (pv != col) {
                    for (std::size_t cc = 0; cc < m; ++cc) std::swap(As[col * m + cc], As[pv * m + cc]);
                    std::swap(c[col], c[pv]);
                }
                T dia = As[col * m + col];
                if (std::abs(dia) < eps) dia = (dia >= T(0) ? T(1) : T(-1)) * eps;
                for (std::size_t rr = col + 1; rr < m; ++rr) {
                    T fac = As[rr * m + col] / dia;
                    for (std::size_t cc = col; cc < m; ++cc) As[rr * m + cc] -= fac * As[col * m + cc];
                    c[rr] -= fac * c[col];
                }
            }
            for (std::size_t ii = m; ii-- > 0;) {
                T sum = c[ii];
                for (std::size_t cc = ii + 1; cc < m; ++cc) sum -= As[ii * m + cc] * c[cc];
                c[ii] = sum / As[ii * m + ii];
            }

            // v = r − Σ_{j≥k} c_j G_j ;  precondition and form U_k, G_k
            Vector<T> v = r;
            for (std::size_t j = k; j < s; ++j) v -= G[j] * c[j - k];
            Vector<T> vp = precond(v);
            Vector<T> uk = vp * om;
            for (std::size_t j = k; j < s; ++j) uk += U[j] * c[j - k];
            U[k] = uk;
            G[k] = A * U[k];

            // bi-orthogonalize against P[0..k-1]
            for (std::size_t i = 0; i < k; ++i) {
                T mii = Mref(i, i);
                if (std::abs(mii) < eps) continue;
                T alpha = dot(P[i], G[k]) / mii;
                G[k] -= G[i] * alpha;
                U[k] -= U[i] * alpha;
            }
            for (std::size_t i = k; i < s; ++i) Mref(i, k) = dot(P[i], G[k]);

            T mkk = Mref(k, k);
            if (std::abs(mkk) < eps) mkk = (mkk >= T(0) ? T(1) : T(-1)) * eps;
            T beta = f[k] / mkk;
            r -= G[k] * beta;
            x += U[k] * beta;

            T rnorm = norm(r);
            if (is_converged(rnorm)) return {x, rnorm, iter + 1, true};

            for (std::size_t i = k + 1; i < s; ++i) f[i] -= beta * Mref(i, k);
        }

        // Move to the next G space (omega step)
        Vector<T> vp = precond(r);
        Vector<T> t = A * vp;
        T tt = dot(t, t);
        if (std::abs(tt) < eps) {
            T rnorm = norm(b - A * x);
            return {x, rnorm, iter + 1, is_converged(rnorm)};
        }
        om = dot(t, r) / tt;
        x += vp * om;
        r -= t * om;
        ++iter;

        T rnorm = norm(r);
        if (is_converged(rnorm)) return {x, rnorm, iter, true};
    }

    T final_rnorm = norm(b - A * x);
    return {x, final_rnorm, criteria.max_iterations(), is_converged(final_rnorm)};
}

// ============================================================================
// QMR — Quasi-Minimal Residual (nonsymmetric, three-term recurrence + quasi-minimization)
// ============================================================================

/**
 * @brief QMR (Freund & Nachtigal 1991, without look-ahead)
 *
 * Solves Ax=b for a nonsymmetric sparse matrix A. On the basis constructed by
 * nonsymmetric Lanczos biorthogonalization, it minimizes the pseudo-norm of the
 * Hessenberg system (quasi-minimal), so the residual oscillates less and converges
 * more smoothly than BiCG. Since it uses A^T, transpose() is formed once.
 *
 * Preconditioning is applied as left preconditioning assuming a symmetric preconditioner
 * (Jacobi/SSOR/IC etc.). If you want to use a nonsymmetric ILU, FGMRES / BiCGSTAB is recommended.
 */
template<typename T, typename IndexType>
    requires concepts::OrderedField<T> && std::integral<IndexType>
SolverResult<T> qmr(
    const SparseMatrix<T, IndexType>& A,
    const Vector<T>& b,
    const ConvergenceCriteria<T>& criteria,
    Preconditioner<T> precond = {},
    std::optional<Vector<T>> x0 = std::nullopt)
{
    const auto n = A.rows();
    if (A.cols() != n) throw std::invalid_argument("QMR: matrix must be square");
    if (b.size() != n) throw std::invalid_argument("QMR: incompatible dimensions");
    if (!precond) precond = identity_preconditioner<T>();

    const SparseMatrix<T, IndexType> At = A.transpose();

    Vector<T> x = x0.value_or(Vector<T>(n, T{0}));
    Vector<T> r = b - A * x;
    const T init_rnorm = norm(r);
    if (init_rnorm <= criteria.absolute_tolerance())
        return {x, init_rnorm, 0, true};

    auto is_converged = [&](T rn) {
        return rn <= criteria.absolute_tolerance() ||
               rn <= criteria.relative_tolerance() * init_rnorm;
    };
    const T eps = std::numeric_limits<T>::epsilon();
    const T tiny = eps * init_rnorm;

    Vector<T> v_tilde = r;
    Vector<T> y = precond(v_tilde);      // M1^{-1} v~
    T rho = norm(y);
    Vector<T> w_tilde = r;
    Vector<T> z = w_tilde;               // M2^{-T} = I
    T xi = norm(z);

    T gamma = T(1), eta = T(-1), theta = T(0), epsilon_v = T(1);
    Vector<T> p(n, T{0}), q(n, T{0}), d(n, T{0}), svec(n, T{0}), v(n, T{0}), w(n, T{0});

    for (std::size_t iter = 1; iter <= criteria.max_iterations(); ++iter) {
        if (std::abs(rho) < tiny || std::abs(xi) < tiny) {
            T rn = norm(b - A * x);
            return {x, rn, iter - 1, is_converged(rn)};
        }
        v = v_tilde * (T(1) / rho);
        y *= (T(1) / rho);
        w = w_tilde * (T(1) / xi);
        z *= (T(1) / xi);

        T delta = dot(z, y);
        if (std::abs(delta) < tiny) {
            T rn = norm(b - A * x);
            return {x, rn, iter - 1, is_converged(rn)};
        }

        Vector<T> y_tilde = y;            // M2^{-1} = I
        Vector<T> z_tilde = precond(z);   // M1^{-T} = M1^{-1} (symmetric preconditioner)

        if (iter == 1) {
            p = y_tilde;
            q = z_tilde;
        } else {
            p = y_tilde - p * (xi * delta / epsilon_v);
            q = z_tilde - q * (rho * delta / epsilon_v);
        }
        Vector<T> p_tilde = A * p;
        epsilon_v = dot(q, p_tilde);
        if (std::abs(epsilon_v) < tiny) {
            T rn = norm(b - A * x);
            return {x, rn, iter - 1, is_converged(rn)};
        }
        T beta = epsilon_v / delta;
        if (std::abs(beta) < tiny) {
            T rn = norm(b - A * x);
            return {x, rn, iter - 1, is_converged(rn)};
        }

        v_tilde = p_tilde - v * beta;
        y = precond(v_tilde);
        T rho_prev = rho;
        rho = norm(y);

        w_tilde = (At * q) - w * beta;
        z = w_tilde;                      // M2^{-T} = I
        xi = norm(z);

        T gamma_prev = gamma;
        T theta_prev = theta;
        theta = rho / (gamma_prev * std::abs(beta));
        gamma = T(1) / std::sqrt(T(1) + theta * theta);
        if (std::abs(gamma) < eps) {
            T rn = norm(b - A * x);
            return {x, rn, iter, is_converged(rn)};
        }
        eta = -eta * rho_prev * (gamma * gamma) / (beta * gamma_prev * gamma_prev);

        if (iter == 1) {
            d = p * eta;
            svec = p_tilde * eta;
        } else {
            T sc = (theta_prev * gamma) * (theta_prev * gamma);
            d = p * eta + d * sc;
            svec = p_tilde * eta + svec * sc;
        }
        x += d;
        r -= svec;

        T rnorm = norm(r);
        if (is_converged(rnorm)) return {x, rnorm, iter, true};
    }

    T final_rnorm = norm(b - A * x);
    return {x, final_rnorm, criteria.max_iterations(), is_converged(final_rnorm)};
}

// ============================================================================
// SolverType enumeration
// ============================================================================

enum class SolverType {
    CG,         ///< Conjugate gradient method (SPD)
    BiCG,       ///< BiCG (nonsymmetric)
    BiCGSTAB,   ///< BiCGSTAB (nonsymmetric, stable)
    GMRES,      ///< GMRES(m) (nonsymmetric)
    FGMRES,     ///< Flexible GMRES (variable preconditioning)
    MINRES,     ///< MINRES (symmetric indefinite)
    IDR,        ///< IDR(s) (nonsymmetric, transpose-free)
    QMR,        ///< QMR (nonsymmetric, quasi-minimal residual)
    Auto        ///< automatic selection
};

enum class PreconditionerType {
    None,       ///< no preconditioning
    Jacobi,     ///< Jacobi (diagonal)
    SSOR,       ///< SSOR
    ILU0,       ///< ILU(0)
    IC0,        ///< IC(0) (for SPD)
    AMG,        ///< algebraic multigrid (aggregation type, for SPD sparse matrices)
    Custom      ///< custom
};

// ============================================================================
// IterativeSolver — unified API via the builder pattern
// ============================================================================

/**
 * @brief Unified interface for iterative solvers
 *
 * Usage example:
 * @code
 * auto result = IterativeSolver<double>(A, b)
 *     .solver(SolverType::BiCGSTAB)
 *     .preconditioner(PreconditionerType::ILU0)
 *     .tolerance(1e-10, 1e-8)
 *     .maxIterations(500)
 *     .enableLog()
 *     .solve();
 * @endcode
 */
template<typename T, typename IndexType = std::size_t>
    requires concepts::OrderedField<T> && std::integral<IndexType>
class IterativeSolver {
public:
    IterativeSolver(const SparseMatrix<T, IndexType>& A, const Vector<T>& b)
        : A_(A), b_(b) {}

    /// Set the solver type
    IterativeSolver& solver(SolverType type) {
        solver_type_ = type;
        return *this;
    }

    /// Set a built-in preconditioner
    IterativeSolver& preconditioner(PreconditionerType type) {
        precond_type_ = type;
        return *this;
    }

    /// Set a custom preconditioner
    IterativeSolver& preconditioner(Preconditioner<T> p) {
        precond_type_ = PreconditionerType::Custom;
        custom_precond_ = std::move(p);
        return *this;
    }

    /// Set the SSOR relaxation factor
    IterativeSolver& ssorOmega(T omega) {
        ssor_omega_ = omega;
        return *this;
    }

    /// Set the tolerances
    IterativeSolver& tolerance(T abs_tol, T rel_tol) {
        criteria_.set_absolute_tolerance(abs_tol);
        criteria_.set_relative_tolerance(rel_tol);
        return *this;
    }

    /// Set the maximum number of iterations
    IterativeSolver& maxIterations(std::size_t n) {
        criteria_.set_max_iterations(n);
        return *this;
    }

    /// Set the GMRES restart parameter
    IterativeSolver& restart(std::size_t m) {
        restart_ = m;
        return *this;
    }

    /// Set the IDR(s) shadow space dimension s
    IterativeSolver& idrShadowSize(std::size_t s) {
        idr_s_ = s;
        return *this;
    }

    /// Set the initial guess
    IterativeSolver& initialGuess(Vector<T> x0) {
        x0_ = std::move(x0);
        return *this;
    }

    /// Enable convergence history recording
    IterativeSolver& enableLog() {
        log_enabled_ = true;
        return *this;
    }

    /// Get the convergence history (call after solve())
    [[nodiscard]] const SolverLog<T>& log() const { return log_; }

    /// Run the solver
    [[nodiscard]] SolverResult<T> solve() {
        auto precond = buildPreconditioner();
        auto type = (solver_type_ == SolverType::Auto) ? autoSelect() : solver_type_;

        if (log_enabled_) log_.clear();
        SolverLog<T>* log_ptr = log_enabled_ ? &log_ : nullptr;

        switch (type) {
        case SolverType::CG:
            if (log_ptr) {
                return conjugate_gradient_logged(A_, b_, criteria_, *log_ptr, precond, x0_);
            }
            return conjugate_gradient(A_, b_, criteria_, precond, x0_);

        case SolverType::BiCG:
            return biconjugate_gradient(A_, b_, criteria_, x0_);

        case SolverType::BiCGSTAB:
            if (log_ptr) {
                return bicgstab_logged(A_, b_, criteria_, *log_ptr, precond, x0_);
            }
            return bicgstab(A_, b_, criteria_, precond, x0_);

        case SolverType::GMRES:
            return gmres(A_, b_, criteria_, restart_, precond, x0_);

        case SolverType::FGMRES:
            return fgmres(A_, b_, criteria_, restart_, precond, x0_, log_ptr);

        case SolverType::MINRES:
            return minres(A_, b_, criteria_, precond, x0_, log_ptr);

        case SolverType::IDR:
            return idrs(A_, b_, criteria_, idr_s_, precond, x0_);

        case SolverType::QMR:
            return qmr(A_, b_, criteria_, precond, x0_);

        default:
            return conjugate_gradient(A_, b_, criteria_, precond, x0_);
        }
    }

private:
    const SparseMatrix<T, IndexType>& A_;
    const Vector<T>& b_;
    SolverType solver_type_ = SolverType::Auto;
    PreconditionerType precond_type_ = PreconditionerType::None;
    Preconditioner<T> custom_precond_;
    ConvergenceCriteria<T> criteria_;
    std::size_t restart_ = 30;
    std::size_t idr_s_ = 4;
    std::optional<Vector<T>> x0_;
    T ssor_omega_ = T{1};
    bool log_enabled_ = false;
    SolverLog<T> log_;

    Preconditioner<T> buildPreconditioner() {
        switch (precond_type_) {
        case PreconditionerType::Jacobi:
            return jacobi_preconditioner(A_);
        case PreconditionerType::SSOR:
            return ssor_preconditioner(A_, ssor_omega_);
        case PreconditionerType::ILU0:
            return ilu_preconditioner(A_);
        case PreconditionerType::IC0:
            return incomplete_cholesky_preconditioner(A_);
        case PreconditionerType::AMG:
            return amg_preconditioner(A_);
        case PreconditionerType::Custom:
            return custom_precond_;
        case PreconditionerType::None:
        default:
            return {};
        }
    }

    /// Automatically select a solver based on matrix properties
    SolverType autoSelect() const {
        // Symmetry test (sampling)
        bool symmetric = isApproximatelySymmetric();
        if (symmetric) {
            // Possibly SPD → try CG (fall back to MINRES if indefinite)
            return SolverType::CG;
        }
        // Nonsymmetric → BiCGSTAB (generally stable)
        return SolverType::BiCGSTAB;
    }

    /// Quick symmetry test (sampling near the diagonal)
    bool isApproximatelySymmetric() const {
        const auto n = A_.rows();
        const std::size_t samples = std::min(static_cast<std::size_t>(n),
                                              std::size_t(50));
        const T eps = std::numeric_limits<T>::epsilon() * T(1000);
        for (std::size_t k = 0; k < samples; ++k) {
            std::size_t i = k * static_cast<std::size_t>(n) / samples;
            for (std::size_t j = i + 1; j < std::min(i + 5, static_cast<std::size_t>(n)); ++j) {
                T aij = A_.coeff(static_cast<IndexType>(i), static_cast<IndexType>(j));
                T aji = A_.coeff(static_cast<IndexType>(j), static_cast<IndexType>(i));
                T scale = std::max(std::abs(aij), std::abs(aji));
                if (scale > T(0) && std::abs(aij - aji) > eps * scale) {
                    return false;
                }
            }
        }
        return true;
    }
};

// ============================================================================
// solve() — convenience function for automatic selection
// ============================================================================

/**
 * @brief Automatically solve a sparse linear system with an appropriate solver
 *
 * Determines the matrix properties and selects CG if symmetric, BiCGSTAB if nonsymmetric.
 * Automatically applies ILU(0) preconditioning.
 */
template<typename T, typename IndexType>
    requires concepts::OrderedField<T> && std::integral<IndexType>
SolverResult<T> solve(
    const SparseMatrix<T, IndexType>& A,
    const Vector<T>& b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>())
{
    return IterativeSolver<T, IndexType>(A, b)
        .preconditioner(PreconditionerType::ILU0)
        .solver(SolverType::Auto)
        .tolerance(criteria.absolute_tolerance(), criteria.relative_tolerance())
        .maxIterations(criteria.max_iterations())
        .solve();
}

} // namespace sparse_algorithms

} // namespace sangi

#endif // SANGI_ITERATIVE_SOLVERS_HPP