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

/**
 * @file iterative_solvers.hpp
 * @brief 反復法ソルバー API 拡充
 *
 * 追加ソルバー:
 *   - MINRES (対称不定行列用)
 *   - FGMRES (可変前処理対応の Flexible GMRES)
 *
 * 統一 API:
 *   - IterativeSolver<T> : ビルダーパターンによるソルバー構成
 *   - SolverLog<T>       : 収束履歴の記録
 *   - solve()            : 行列性質に基づく自動ソルバー選択
 */

#ifndef CALX_ITERATIVE_SOLVERS_HPP
#define CALX_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>

namespace calx {

namespace sparse_algorithms {

// ============================================================================
// SolverLog — 収束履歴の記録
// ============================================================================

template<typename T>
struct SolverLog {
    std::vector<T> residual_history;  ///< 各反復の残差ノルム
    std::vector<T> time_history;      ///< (将来用) 各反復の累積時間

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

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

    /// 収束率 (最後の数反復の幾何平均)
    [[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 — 対称不定行列用
// ============================================================================

/**
 * @brief MINRES (Minimum Residual Method)
 *
 * 対称 (不定可) 疎行列 A に対して Ax = b を解く。
 * CG が正定値を要求するのに対し、MINRES は不定行列にも使える。
 *
 * 実装: Lanczos 三重対角化 + GMRES 式 Givens QR
 * Lanczos で対称行列を三重対角化し、GMRES と同じ Givens 回転手順で
 * Hessenberg 上三角分解を行う。対称性から Hessenberg = 三重対角となり
 * 各列の回転が 1 回で済む (GMRES の j 回に対して定数回)。
 */
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 方式を三重対角に特化)
    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 基底
    std::vector<Vector<T>> Q(max_iter + 1);
    Q[0] = r / r_norm;

    // 三重対角行列を Hessenberg 形式で保持 (列ごとに 2 要素: H[j][j], H[j+1][j])
    // ただし Givens 回転適用前の値を保持するため、前の列も必要
    // → GMRES と同じ構造を使う (H は (m+1)×m だが三重対角なので帯幅 2)
    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 ステップ (= 対称行列の Arnoldi ステップ)
        Vector<T> w = A * Q[j];

        // Modified Gram-Schmidt (対称性から j-1 と j のみ非ゼロ)
        if (j > 0) {
            H[j - 1][j] = dot(Q[j - 1], w);  // = beta_j (三重対角の上側)
            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}) {
            // Givens 回転を適用してから終了
            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];

        // 以前の Givens 回転を H 列に適用 (三重対角 → G_{j-2} が fill-in を生むため最大 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 番目の Givens 回転
        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;
        }
    }

    // 後退代入: 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];
        }
    }

    // 解の更新: 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 (可変前処理対応)
// ============================================================================

/**
 * @brief Flexible GMRES (可変前処理対応)
 *
 * 通常の GMRES では各反復で同一の前処理が必要だが、
 * FGMRES では反復ごとに異なる前処理を適用できる。
 * 前処理付き基底ベクトルを別途保持するため、メモリは 2 倍。
 *
 * 参考: 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 基底 V と前処理済み基底 Z
        std::vector<Vector<T>> V(m + 1);
        std::vector<Vector<T>> Z(m);
        V[0] = r / r_norm;

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

        // Givens 回転パラメータ
        std::vector<T> cs(m, T{0}), sn(m, T{0});

        // 変換済み右辺
        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 (前処理済み基底を保持)
            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}) {
                // 以前の Givens 回転を適用
                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];

            // 以前の Givens 回転を H 列に適用
            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 番目の Givens 回転
            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;
            }
        }

        // 後退代入
        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 (前処理済み基底を使用)
        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};
}

// ============================================================================
// ログ付き CG / BiCGSTAB / GMRES — 既存ソルバーのラッパー
// ============================================================================

/**
 * @brief ログ付き 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 ログ付き 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)};
}

// ============================================================================
// SolverType 列挙
// ============================================================================

enum class SolverType {
    CG,         ///< 共役勾配法 (SPD)
    BiCG,       ///< BiCG (非対称)
    BiCGSTAB,   ///< BiCGSTAB (非対称, 安定)
    GMRES,      ///< GMRES(m) (非対称)
    FGMRES,     ///< Flexible GMRES (可変前処理)
    MINRES,     ///< MINRES (対称不定)
    Auto        ///< 自動選択
};

enum class PreconditionerType {
    None,       ///< 前処理なし
    Jacobi,     ///< Jacobi (対角)
    SSOR,       ///< SSOR
    ILU0,       ///< ILU(0)
    IC0,        ///< IC(0) (SPD 用)
    Custom      ///< カスタム
};

// ============================================================================
// IterativeSolver — ビルダーパターンによる統一 API
// ============================================================================

/**
 * @brief 反復法ソルバーの統一インタフェース
 *
 * 使用例:
 * @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) {}

    /// ソルバーの種類を設定
    IterativeSolver& solver(SolverType type) {
        solver_type_ = type;
        return *this;
    }

    /// 組み込み前処理を設定
    IterativeSolver& preconditioner(PreconditionerType type) {
        precond_type_ = type;
        return *this;
    }

    /// カスタム前処理を設定
    IterativeSolver& preconditioner(Preconditioner<T> p) {
        precond_type_ = PreconditionerType::Custom;
        custom_precond_ = std::move(p);
        return *this;
    }

    /// SSOR の緩和係数を設定
    IterativeSolver& ssorOmega(T omega) {
        ssor_omega_ = omega;
        return *this;
    }

    /// 許容誤差を設定
    IterativeSolver& tolerance(T abs_tol, T rel_tol) {
        criteria_.set_absolute_tolerance(abs_tol);
        criteria_.set_relative_tolerance(rel_tol);
        return *this;
    }

    /// 最大反復回数を設定
    IterativeSolver& maxIterations(std::size_t n) {
        criteria_.set_max_iterations(n);
        return *this;
    }

    /// GMRES リスタートパラメータを設定
    IterativeSolver& restart(std::size_t m) {
        restart_ = m;
        return *this;
    }

    /// 初期推定値を設定
    IterativeSolver& initialGuess(Vector<T> x0) {
        x0_ = std::move(x0);
        return *this;
    }

    /// 収束履歴の記録を有効化
    IterativeSolver& enableLog() {
        log_enabled_ = true;
        return *this;
    }

    /// 収束履歴を取得 (solve() 後に呼ぶ)
    [[nodiscard]] const SolverLog<T>& log() const { return log_; }

    /// ソルバーを実行
    [[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);

        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::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::Custom:
            return custom_precond_;
        case PreconditionerType::None:
        default:
            return {};
        }
    }

    /// 行列の性質から自動的にソルバーを選択
    SolverType autoSelect() const {
        // 対称判定 (サンプリング)
        bool symmetric = isApproximatelySymmetric();
        if (symmetric) {
            // SPD の可能性 → CG を試行 (不定なら MINRES にフォールバック)
            return SolverType::CG;
        }
        // 非対称 → BiCGSTAB (一般的に安定)
        return SolverType::BiCGSTAB;
    }

    /// 対称性の簡易判定 (対角近傍のサンプリング)
    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() — 自動選択の便利関数
// ============================================================================

/**
 * @brief 疎行列の連立一次方程式を自動的に適切なソルバーで解く
 *
 * 行列の性質を判定し、対称なら CG、非対称なら BiCGSTAB を選択。
 * ILU(0) 前処理を自動適用。
 */
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 calx

#endif // CALX_ITERATIVE_SOLVERS_HPP