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

/**
 * @file sparse_matrix_algorithms.hpp
 * @brief 疎行列向け反復解法・前処理器の実装
 *
 * 反復解法: CG, BiCG, BiCGSTAB, GMRES(m)
 * 前処理器: Jacobi, SSOR, ILU(0), IC(0)
 * 全ソルバーは SolverResult<T> を返し、前処理は callable で注入する。
 */

#ifndef CALX_SPARSE_MATRIX_ALGORITHMS_HPP
#define CALX_SPARSE_MATRIX_ALGORITHMS_HPP

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

#include "../concepts/algebraic_concepts.hpp"
#include "sparse_matrix.hpp"
#include "vector.hpp"
#include "convergence_criteria.hpp"

namespace calx {

namespace sparse_algorithms {

// ============================================================================
// SolverResult — ソルバーの戻り値
// ============================================================================

template<typename T>
struct SolverResult {
    Vector<T>    x;              ///< 解ベクトル
    T            residual_norm;  ///< 最終残差ノルム
    std::size_t  iterations;     ///< 反復回数
    bool         converged;      ///< 収束したか
};

// ============================================================================
// Preconditioner 型定義 — callable (Vector → Vector)
// ============================================================================

template<typename T>
using Preconditioner = std::function<Vector<T>(const Vector<T>&)>;

/// 恒等前処理 (前処理なし)
template<typename T>
Preconditioner<T> identity_preconditioner() {
    return [](const Vector<T>& r) -> Vector<T> { return r; };
}

// ============================================================================
// ヘルパー: 疎行列の対角抽出
// ============================================================================

template<typename T, typename IndexType>
    requires concepts::OrderedField<T> && std::integral<IndexType>
std::vector<T> sparse_diagonal(const SparseMatrix<T, IndexType>& A) {
    const auto n = A.rows();
    std::vector<T> diag(n, T{0});
    for (IndexType i = 0; i < n; ++i) {
        diag[i] = A.coeff(i, i);
    }
    return diag;
}

// ============================================================================
// 前処理器ファクトリ
// ============================================================================

/**
 * @brief Jacobi (対角) 前処理
 *
 * M = diag(A), M^{-1} r = r[i] / A(i,i)
 * O(n) で最も軽量な前処理。対角優位行列に有効。
 */
template<typename T, typename IndexType>
    requires concepts::OrderedField<T> && std::integral<IndexType>
Preconditioner<T> jacobi_preconditioner(const SparseMatrix<T, IndexType>& A) {
    auto diag = sparse_diagonal(A);
    const auto n = A.rows();
    // ゼロ対角のチェック
    for (IndexType i = 0; i < n; ++i) {
        if (std::abs(diag[i]) < std::numeric_limits<T>::epsilon()) {
            throw std::invalid_argument(
                "Jacobi preconditioner: zero diagonal element at row " + std::to_string(i));
        }
    }
    return [diag = std::move(diag)](const Vector<T>& r) -> Vector<T> {
        const auto m = r.size();
        Vector<T> z(m);
        for (std::size_t i = 0; i < m; ++i) {
            z[i] = r[i] / diag[i];
        }
        return z;
    };
}

/**
 * @brief SSOR (Symmetric Successive Over-Relaxation) 前処理
 *
 * M = (D/ω + L) D^{-1} (D/ω + U)
 * 前進・後退代入で M^{-1} r を計算。
 * ω=1.0 で対称 Gauss-Seidel に一致。
 */
template<typename T, typename IndexType>
    requires concepts::OrderedField<T> && std::integral<IndexType>
Preconditioner<T> ssor_preconditioner(
    const SparseMatrix<T, IndexType>& A,
    T omega = T{1})
{
    // CSR に変換して保持
    auto A_csr = std::make_shared<SparseMatrix<T, IndexType>>(A);
    A_csr->convert_to(SparseStorageFormat::CSR);
    const auto n = A.rows();
    auto diag = sparse_diagonal(A);

    return [A_csr, diag, omega, n](const Vector<T>& r) -> Vector<T> {
        Vector<T> y(n, T{0});
        Vector<T> z(n, T{0});

        // 前進代入: (D/ω + L) y = r
        for (IndexType i = 0; i < n; ++i) {
            T sum = r[i];
            for (IndexType j = 0; j < i; ++j) {
                T aij = A_csr->coeff(i, j);
                if (aij != T{0}) {
                    sum -= aij * y[j];
                }
            }
            y[i] = omega * sum / diag[i];
        }

        // D/ω のスケーリング: y_i *= diag[i] / omega
        for (IndexType i = 0; i < n; ++i) {
            y[i] *= diag[i] / omega;
        }

        // 後退代入: (D/ω + U) z = y
        for (IndexType ii = 0; ii < n; ++ii) {
            IndexType i = n - 1 - ii;
            T sum = y[i];
            for (IndexType j = i + 1; j < n; ++j) {
                T aij = A_csr->coeff(i, j);
                if (aij != T{0}) {
                    sum -= aij * z[j];
                }
            }
            z[i] = omega * sum / diag[i];
        }

        // SSOR の正しいスケーリング: M_SSOR = ω(2-ω) (D/ω+L) D^{-1} (D/ω+U)
        // 最終結果に (2-ω) を掛ける
        for (IndexType i = 0; i < n; ++i) {
            z[i] *= (T{2} - omega);
        }

        return z;
    };
}

/**
 * @brief ILU(0) (Incomplete LU with zero fill) 前処理
 *
 * A の非ゼロパターン内でのみ LU 分解を行い、L, U を CSR で保持する。
 * M^{-1} r は L y = r (前進代入) → U z = y (後退代入) で計算。
 */
template<typename T, typename IndexType>
    requires concepts::OrderedField<T> && std::integral<IndexType>
Preconditioner<T> ilu_preconditioner(const SparseMatrix<T, IndexType>& A) {
    const auto n = A.rows();
    if (A.cols() != n) {
        throw std::invalid_argument("ILU: matrix must be square");
    }

    // CSR 形式でワーク配列を構築
    // 行ごとの密ベクトルで IKJ 版 ILU(0) を実行
    // lu_vals[i][j] = (i,j) の LU 値 (Aのパターン内のみ)

    // まず各行の非ゼロ列インデックスと値を取得
    struct RowEntry {
        IndexType col;
        T val;
    };
    std::vector<std::vector<RowEntry>> rows(n);
    for (IndexType i = 0; i < n; ++i) {
        for (IndexType j = 0; j < n; ++j) {
            T v = A.coeff(i, j);
            if (v != T{0}) {
                rows[i].push_back({j, v});
            }
        }
        // 列順でソート
        std::sort(rows[i].begin(), rows[i].end(),
                  [](const RowEntry& a, const RowEntry& b) { return a.col < b.col; });
    }

    // IKJ 版 ILU(0): 行 i の各エントリを更新
    // w = A(i,:) のコピー → k < i のピボット行で消去 → L,U にコピー
    std::vector<std::vector<RowEntry>> L_rows(n), U_rows(n);

    for (IndexType i = 0; i < n; ++i) {
        // 行 i の密ワーク配列 (非ゼロ位置のみ)
        std::vector<T> w(n, T{0});
        std::vector<bool> nz(n, false);  // 非ゼロパターン
        for (auto& e : rows[i]) {
            w[e.col] = e.val;
            nz[e.col] = true;
        }

        // k ループ: k < i で nz[k] が true の列に対して消去
        for (IndexType k = 0; k < i; ++k) {
            if (!nz[k]) continue;

            // U(k,k) で割る
            T ukk = T{0};
            for (auto& e : U_rows[k]) {
                if (e.col == k) { ukk = e.val; break; }
            }
            if (std::abs(ukk) < std::numeric_limits<T>::epsilon()) {
                continue;  // ゼロピボットはスキップ
            }
            w[k] /= ukk;

            // j > k で U(k,j) != 0 かつ nz[j] の要素を更新
            for (auto& e : U_rows[k]) {
                if (e.col > k && nz[e.col]) {
                    w[e.col] -= w[k] * e.val;
                }
            }
        }

        // L, U に格納
        for (auto& e : rows[i]) {
            IndexType j = e.col;
            if (j < i) {
                L_rows[i].push_back({j, w[j]});
            } else {
                U_rows[i].push_back({j, w[j]});
            }
        }
        // L の対角は 1
        L_rows[i].push_back({i, T{1}});
        std::sort(L_rows[i].begin(), L_rows[i].end(),
                  [](const RowEntry& a, const RowEntry& b) { return a.col < b.col; });
    }

    // shared_ptr で L, U 行データを保持
    auto l_data = std::make_shared<std::vector<std::vector<RowEntry>>>(std::move(L_rows));
    auto u_data = std::make_shared<std::vector<std::vector<RowEntry>>>(std::move(U_rows));

    return [l_data, u_data, n](const Vector<T>& r) -> Vector<T> {
        // 前進代入: L y = r
        Vector<T> y(n, T{0});
        for (IndexType i = 0; i < n; ++i) {
            T sum = r[i];
            for (auto& e : (*l_data)[i]) {
                if (e.col < i) {
                    sum -= e.val * y[e.col];
                }
            }
            // L(i,i) = 1 なので割る必要なし
            y[i] = sum;
        }

        // 後退代入: U z = y
        Vector<T> z(n, T{0});
        for (IndexType ii = 0; ii < n; ++ii) {
            IndexType i = n - 1 - ii;
            T sum = y[i];
            T uii = T{0};
            for (auto& e : (*u_data)[i]) {
                if (e.col == i) {
                    uii = e.val;
                } else if (e.col > i) {
                    sum -= e.val * z[e.col];
                }
            }
            if (std::abs(uii) < std::numeric_limits<T>::epsilon()) {
                z[i] = T{0};
            } else {
                z[i] = sum / uii;
            }
        }

        return z;
    };
}

// ============================================================================
// 後方互換: incomplete_lu_decomposition (L, U ペアを返す版)
// ============================================================================

template<typename T, typename IndexType>
    requires concepts::OrderedField<T> && std::integral<IndexType>
std::pair<SparseMatrix<T, IndexType>, SparseMatrix<T, IndexType>>
incomplete_lu_decomposition(
    const SparseMatrix<T, IndexType>& A,
    [[maybe_unused]] int fill_level = 0)
{
    const auto n = A.rows();
    if (A.cols() != n) {
        throw std::invalid_argument("ILU: matrix must be square");
    }

    struct RowEntry { IndexType col; T val; };
    std::vector<std::vector<RowEntry>> rows(n);
    for (IndexType i = 0; i < n; ++i) {
        for (IndexType j = 0; j < n; ++j) {
            T v = A.coeff(i, j);
            if (v != T{0}) {
                rows[i].push_back({j, v});
            }
        }
        std::sort(rows[i].begin(), rows[i].end(),
                  [](const RowEntry& a, const RowEntry& b) { return a.col < b.col; });
    }

    std::vector<std::vector<RowEntry>> L_rows(n), U_rows(n);

    for (IndexType i = 0; i < n; ++i) {
        std::vector<T> w(n, T{0});
        std::vector<bool> nz(n, false);
        for (auto& e : rows[i]) {
            w[e.col] = e.val;
            nz[e.col] = true;
        }
        for (IndexType k = 0; k < i; ++k) {
            if (!nz[k]) continue;
            T ukk = T{0};
            for (auto& e : U_rows[k]) {
                if (e.col == k) { ukk = e.val; break; }
            }
            if (std::abs(ukk) < std::numeric_limits<T>::epsilon()) continue;
            w[k] /= ukk;
            for (auto& e : U_rows[k]) {
                if (e.col > k && nz[e.col]) {
                    w[e.col] -= w[k] * e.val;
                }
            }
        }
        for (auto& e : rows[i]) {
            if (e.col < i) L_rows[i].push_back({e.col, w[e.col]});
            else           U_rows[i].push_back({e.col, w[e.col]});
        }
        L_rows[i].push_back({i, T{1}});
        std::sort(L_rows[i].begin(), L_rows[i].end(),
                  [](const RowEntry& a, const RowEntry& b) { return a.col < b.col; });
    }

    SparseMatrix<T, IndexType> L(n, n, SparseStorageFormat::COO);
    SparseMatrix<T, IndexType> U(n, n, SparseStorageFormat::COO);
    for (IndexType i = 0; i < n; ++i) {
        for (auto& e : L_rows[i]) L.set_coeff(i, e.col, e.val);
        for (auto& e : U_rows[i]) U.set_coeff(i, e.col, e.val);
    }
    return {L, U};
}

// ============================================================================
// IC(0) (Incomplete Cholesky with zero fill) 前処理
// ============================================================================

/**
 * @brief IC(0) (Incomplete Cholesky with zero fill) 前処理
 *
 * 対称正定値行列 A の下三角非ゼロパターン内でのみ Cholesky 分解を行う。
 * A ≈ L * L^T (fill-in なし)
 * M^{-1} r は L y = r (前進代入) → L^T z = y (後退代入) で計算。
 */
template<typename T, typename IndexType>
    requires concepts::OrderedField<T> && std::integral<IndexType>
Preconditioner<T> incomplete_cholesky_preconditioner(const SparseMatrix<T, IndexType>& A) {
    const auto n = A.rows();
    if (A.cols() != n) {
        throw std::invalid_argument("IC: matrix must be square");
    }

    // 下三角の非ゼロパターンを行ごとに取得 (対角含む)
    struct Entry {
        IndexType col;
        T val;
    };
    std::vector<std::vector<Entry>> L_rows(n);

    for (IndexType i = 0; i < n; ++i) {
        for (IndexType j = 0; j <= i; ++j) {
            T v = A.coeff(i, j);
            if (v != T{0}) {
                L_rows[i].push_back({j, v});
            }
        }
        // 列順でソート済みのはず (j は昇順ループ) だが念のため
        std::sort(L_rows[i].begin(), L_rows[i].end(),
                  [](const Entry& a, const Entry& b) { return a.col < b.col; });
        // 対角要素がパターンに存在するか確認
        if (L_rows[i].empty() || L_rows[i].back().col != i) {
            throw std::runtime_error("IC: missing diagonal element (matrix may not be SPD)");
        }
    }

    // 各列の非ゼロ行を高速検索するため、列→行リストを構築
    // col_rows[j] = {L_rows[i] 内で col==j を持つ行 i のリスト}
    std::vector<std::vector<IndexType>> col_rows(n);
    for (IndexType i = 0; i < n; ++i) {
        for (auto& e : L_rows[i]) {
            if (e.col < i) {  // 対角以外の下三角
                col_rows[e.col].push_back(i);
            }
        }
    }

    // IC(0) 分解: 行ごとに処理
    // L(i,j) = (A(i,j) - Σ_{k<j} L(i,k)*L(j,k)) / L(j,j)  (j < i)
    // L(i,i) = sqrt(A(i,i) - Σ_{k<i} L(i,k)^2)
    //
    // 高速化のため、各行の列→インデックスマップを構築
    std::vector<std::vector<IndexType>> row_col_idx(n);  // row_col_idx[i][p] = L_rows[i][p].col
    for (IndexType i = 0; i < n; ++i) {
        row_col_idx[i].reserve(L_rows[i].size());
        for (auto& e : L_rows[i]) {
            row_col_idx[i].push_back(e.col);
        }
    }

    // 行 i を処理
    for (IndexType i = 0; i < n; ++i) {
        for (std::size_t p = 0; p < L_rows[i].size(); ++p) {
            IndexType j = L_rows[i][p].col;
            T sum = T{0};

            if (j < i) {
                // 非対角: L(i,j) = (A(i,j) - Σ_{k<j} L(i,k)*L(j,k)) / L(j,j)
                // k は行 i と行 j の共通列 (k < j)
                std::size_t pi = 0, pj = 0;
                while (pi < L_rows[i].size() && pj < L_rows[j].size()) {
                    IndexType ci = L_rows[i][pi].col;
                    IndexType cj = L_rows[j][pj].col;
                    if (ci < j && cj < j) {
                        if (ci == cj) {
                            sum += L_rows[i][pi].val * L_rows[j][pj].val;
                            ++pi; ++pj;
                        } else if (ci < cj) {
                            ++pi;
                        } else {
                            ++pj;
                        }
                    } else {
                        break;
                    }
                }
                // L(j,j) を取得 (行 j の最後の要素が対角)
                T ljj = L_rows[j].back().val;
                if (std::abs(ljj) < std::numeric_limits<T>::epsilon()) {
                    throw std::runtime_error("IC: zero diagonal encountered (matrix may not be SPD)");
                }
                L_rows[i][p].val = (L_rows[i][p].val - sum) / ljj;
            } else {
                // 対角: L(i,i) = sqrt(A(i,i) - Σ_{k<i} L(i,k)^2)
                for (std::size_t q = 0; q < p; ++q) {
                    sum += L_rows[i][q].val * L_rows[i][q].val;
                }
                T diag = L_rows[i][p].val - sum;
                if (diag <= T{0}) {
                    throw std::runtime_error("IC: non-positive diagonal encountered (matrix may not be SPD)");
                }
                L_rows[i][p].val = std::sqrt(diag);
            }
        }
    }

    // shared_ptr で L 行データを保持
    auto l_data = std::make_shared<std::vector<std::vector<Entry>>>(std::move(L_rows));
    auto dim = n;

    return [l_data, dim](const Vector<T>& r) -> Vector<T> {
        auto n = dim;
        // 前進代入: L y = r
        Vector<T> y(n, T{0});
        for (IndexType i = 0; i < n; ++i) {
            T sum = r[i];
            for (auto& e : (*l_data)[i]) {
                if (e.col < i) {
                    sum -= e.val * y[e.col];
                }
            }
            // L(i,i) は行の最後の要素
            y[i] = sum / (*l_data)[i].back().val;
        }

        // 後退代入: L^T z = y
        // 行ベース変換: z = y から始め、行 i (降順) で z[i] /= L(i,i),
        // その後 L(i,k) (k < i) について z[k] -= L(i,k) * z[i]
        Vector<T> z = y;
        for (IndexType ii = 0; ii < n; ++ii) {
            IndexType i = n - 1 - ii;
            z[i] /= (*l_data)[i].back().val;
            for (auto& e : (*l_data)[i]) {
                if (e.col < i) {
                    z[e.col] -= e.val * z[i];
                }
            }
        }

        return z;
    };
}

// ============================================================================
// 後方互換: incomplete_cholesky_decomposition (L を返す版)
// ============================================================================

/**
 * @brief IC(0) 分解の L 行列を返す
 */
template<typename T, typename IndexType>
    requires concepts::OrderedField<T> && std::integral<IndexType>
SparseMatrix<T, IndexType>
incomplete_cholesky_decomposition(
    const SparseMatrix<T, IndexType>& A,
    [[maybe_unused]] int fill_level = 0)
{
    const auto n = A.rows();
    if (A.cols() != n) {
        throw std::invalid_argument("IC: matrix must be square");
    }

    // 下三角の非ゼロパターンを取得
    struct Entry { IndexType col; T val; };
    std::vector<std::vector<Entry>> L_rows(n);

    for (IndexType i = 0; i < n; ++i) {
        for (IndexType j = 0; j <= i; ++j) {
            T v = A.coeff(i, j);
            if (v != T{0}) {
                L_rows[i].push_back({j, v});
            }
        }
        std::sort(L_rows[i].begin(), L_rows[i].end(),
                  [](const Entry& a, const Entry& b) { return a.col < b.col; });
        if (L_rows[i].empty() || L_rows[i].back().col != i) {
            throw std::runtime_error("IC: missing diagonal element");
        }
    }

    // IC(0) 分解
    for (IndexType i = 0; i < n; ++i) {
        for (std::size_t p = 0; p < L_rows[i].size(); ++p) {
            IndexType j = L_rows[i][p].col;
            T sum = T{0};

            if (j < i) {
                std::size_t pi = 0, pj = 0;
                while (pi < L_rows[i].size() && pj < L_rows[j].size()) {
                    IndexType ci = L_rows[i][pi].col;
                    IndexType cj = L_rows[j][pj].col;
                    if (ci < j && cj < j) {
                        if (ci == cj) {
                            sum += L_rows[i][pi].val * L_rows[j][pj].val;
                            ++pi; ++pj;
                        } else if (ci < cj) { ++pi; }
                        else { ++pj; }
                    } else { break; }
                }
                T ljj = L_rows[j].back().val;
                if (std::abs(ljj) < std::numeric_limits<T>::epsilon()) {
                    throw std::runtime_error("IC: zero diagonal encountered");
                }
                L_rows[i][p].val = (L_rows[i][p].val - sum) / ljj;
            } else {
                for (std::size_t q = 0; q < p; ++q) {
                    sum += L_rows[i][q].val * L_rows[i][q].val;
                }
                T diag = L_rows[i][p].val - sum;
                if (diag <= T{0}) {
                    throw std::runtime_error("IC: non-positive diagonal");
                }
                L_rows[i][p].val = std::sqrt(diag);
            }
        }
    }

    SparseMatrix<T, IndexType> L(n, n, SparseStorageFormat::COO);
    for (IndexType i = 0; i < n; ++i) {
        for (auto& e : L_rows[i]) {
            L.set_coeff(i, e.col, e.val);
        }
    }
    return L;
}

// ============================================================================
// 三角行列の前進/後退代入 (疎行列版)
// ============================================================================

template<typename T, typename IndexType>
    requires concepts::OrderedField<T> && std::integral<IndexType>
Vector<T> triangular_solve(
    const SparseMatrix<T, IndexType>& M,
    const Vector<T>& b,
    bool lower = true)
{
    const auto n = M.rows();
    if (M.cols() != n) throw std::invalid_argument("Matrix must be square");
    if (b.size() != n) throw std::invalid_argument("Incompatible vector dimensions");

    Vector<T> x(n, T{0});

    if (lower) {
        for (IndexType i = 0; i < n; ++i) {
            T sum = b[i];
            for (IndexType j = 0; j < i; ++j) {
                T v = M.coeff(i, j);
                if (v != T{0}) sum -= v * x[j];
            }
            x[i] = sum / M.coeff(i, i);
        }
    } else {
        for (IndexType ii = 0; ii < n; ++ii) {
            IndexType i = n - 1 - ii;
            T sum = b[i];
            for (IndexType j = i + 1; j < n; ++j) {
                T v = M.coeff(i, j);
                if (v != T{0}) sum -= v * x[j];
            }
            x[i] = sum / M.coeff(i, i);
        }
    }

    return x;
}

// ============================================================================
// CG (Conjugate Gradient) — 対称正定値行列用
// ============================================================================

/**
 * @brief 前処理付き共役勾配法 (Preconditioned CG)
 *
 * 対称正定値疎行列 A に対して Ax = b を解く。
 * 前処理 M を指定すると M^{-1} A x = M^{-1} b を CG で解く。
 */
template<typename T, typename IndexType>
    requires concepts::OrderedField<T> && std::integral<IndexType>
SolverResult<T> conjugate_gradient(
    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("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);
    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()) {
            // 内部残差 r に蓄積誤差があるため、真の残差で判定
            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);
        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)};
}

// ============================================================================
// BiCG (BiConjugate Gradient) — 非対称行列用
// ============================================================================

template<typename T, typename IndexType>
    requires concepts::OrderedField<T> && std::integral<IndexType>
SolverResult<T> biconjugate_gradient(
    const SparseMatrix<T, IndexType>& A,
    const Vector<T>& b,
    const ConvergenceCriteria<T>& criteria,
    std::optional<Vector<T>> x0 = std::nullopt)
{
    const auto n = A.rows();
    if (A.cols() != n) throw std::invalid_argument("BiCG: matrix must be square");
    if (b.size() != n) throw std::invalid_argument("BiCG: incompatible dimensions");

    const auto At = A.transpose();

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

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

    Vector<T> p = r;
    Vector<T> p_tilde = r_tilde;
    T rho = dot(r_tilde, r);

    // breakdown 閾値: 初期残差スケールに対する相対値
    const T breakdown_tol = std::numeric_limits<T>::epsilon() * init_rnorm * init_rnorm;

    // 収束判定ヘルパー (真の残差を使用)
    auto check_convergence = [&](std::size_t iter) -> SolverResult<T> {
        T rnorm = norm(b - A * x);
        bool conv = (rnorm <= criteria.absolute_tolerance() ||
                     rnorm <= criteria.relative_tolerance() * init_rnorm);
        return {x, rnorm, iter, conv};
    };

    for (std::size_t iter = 0; iter < criteria.max_iterations(); ++iter) {
        if (std::abs(rho) < breakdown_tol) {
            return check_convergence(iter);
        }

        Vector<T> Ap = A * p;
        Vector<T> Atp = At * p_tilde;
        T denom = dot(p_tilde, Ap);
        if (std::abs(denom) < breakdown_tol) {
            return check_convergence(iter);
        }
        T alpha = rho / denom;

        x += p * alpha;
        r -= Ap * alpha;
        r_tilde -= Atp * alpha;

        T rnorm = norm(r);
        if (rnorm <= criteria.absolute_tolerance() ||
            rnorm <= criteria.relative_tolerance() * init_rnorm) {
            return {x, rnorm, iter + 1, true};
        }

        T rho_new = dot(r_tilde, r);
        T beta = rho_new / rho;
        rho = rho_new;

        p = r + p * beta;
        p_tilde = r_tilde + p_tilde * beta;
    }

    // 最大反復到達時: 真の残差で再判定
    T final_rnorm = norm(b - A * x);
    bool conv = (final_rnorm <= criteria.absolute_tolerance() ||
                 final_rnorm <= criteria.relative_tolerance() * init_rnorm);
    return {x, final_rnorm, criteria.max_iterations(), conv};
}

// ============================================================================
// BiCGSTAB (Stabilized BiConjugate Gradient) — 非対称行列用
// ============================================================================

template<typename T, typename IndexType>
    requires concepts::OrderedField<T> && std::integral<IndexType>
SolverResult<T> bicgstab(
    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("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);
    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});
    Vector<T> v(n, T{0});

    // breakdown 判定の閾値: 初期残差スケールに対する相対値
    const T eps = std::numeric_limits<T>::epsilon();
    const T breakdown_tol = eps * init_rnorm * init_rnorm;
    const T stag_tol      = eps * init_rnorm;

    // 真の残差で収束判定するヘルパー (内部残差 r の蓄積誤差対策)
    auto check_true_convergence = [&](std::size_t iter) -> SolverResult<T> {
        T rnorm = norm(b - A * x);
        bool conv = (rnorm <= criteria.absolute_tolerance() ||
                     rnorm <= criteria.relative_tolerance() * init_rnorm);
        return {x, rnorm, iter, conv};
    };

    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) {
            return check_true_convergence(iter);
        }

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

        // 前処理付き: p_hat = M^{-1} p
        Vector<T> p_hat = precond(p);
        v = A * p_hat;

        T denom = dot(r_hat, v);
        if (std::abs(denom) < breakdown_tol) {
            return check_true_convergence(iter);
        }
        alpha = rho / denom;

        Vector<T> s = r - v * alpha;
        T s_norm = norm(s);
        if (s_norm <= criteria.absolute_tolerance() ||
            s_norm <= criteria.relative_tolerance() * init_rnorm) {
            x += p_hat * alpha;
            return {x, s_norm, iter + 1, true};
        }

        // 前処理付き: s_hat = M^{-1} s
        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;
            return check_true_convergence(iter + 1);
        }
        omega_val = dot(t, s) / tt;

        x += p_hat * alpha + s_hat * omega_val;
        r = s - t * omega_val;

        T rnorm = norm(r);
        if (rnorm <= criteria.absolute_tolerance() ||
            rnorm <= criteria.relative_tolerance() * init_rnorm) {
            return {x, rnorm, iter + 1, true};
        }

        if (std::abs(omega_val) < stag_tol) {
            return check_true_convergence(iter + 1);
        }
    }

    // 最大反復到達時: 真の残差で再判定 (内部残差 r の蓄積誤差対策)
    T final_rnorm = norm(b - A * x);
    bool conv = (final_rnorm <= criteria.absolute_tolerance() ||
                 final_rnorm <= criteria.relative_tolerance() * init_rnorm);
    return {x, final_rnorm, criteria.max_iterations(), conv};
}

// ============================================================================
// GMRES(m) (Generalized Minimum Residual) — 非対称行列用
// ============================================================================

/**
 * @brief リスタート付き GMRES
 *
 * Arnoldi + Givens 回転を逐次適用。リスタートで外部反復。
 * @param restart Krylov 部分空間の最大次元 (デフォルト 30)
 */
template<typename T, typename IndexType>
    requires concepts::OrderedField<T> && std::integral<IndexType>
SolverResult<T> gmres(
    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)
{
    const auto n = static_cast<std::size_t>(A.rows());
    if (static_cast<std::size_t>(A.cols()) != n)
        throw std::invalid_argument("GMRES: matrix must be square");
    if (b.size() != n)
        throw std::invalid_argument("GMRES: 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);
    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 = precond(b - A * x);
        T r_norm = norm(r);
        if (r_norm <= threshold) {
            return {x, r_norm, total_iter, true};
        }

        // Arnoldi 基底
        const std::size_t m = std::min(restart,
            criteria.max_iterations() - total_iter);
        std::vector<Vector<T>> Q(m + 1);
        Q[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) {
            // Arnoldi ステップ: w = M^{-1} A q_j
            Vector<T> w = precond(A * Q[j]);

            // Modified Gram-Schmidt
            for (std::size_t i = 0; i <= j; ++i) {
                H[i][j] = dot(Q[i], w);
                w -= Q[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 回転を 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]);
                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 列に適用
            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 (r_norm <= threshold) {
                ++j;
                break;
            }
        }

        // 後退代入: R y = g (R は H の上三角部分)
        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];
        }

        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};
}

// ============================================================================
// 便利関数: ILU 前処理付き BiCGSTAB (後方互換)
// ============================================================================

template<typename T, typename IndexType>
    requires concepts::OrderedField<T> && std::integral<IndexType>
SolverResult<T> ilu_bicgstab(
    const SparseMatrix<T, IndexType>& A,
    const Vector<T>& b,
    const ConvergenceCriteria<T>& criteria,
    int fill_level = 0,
    std::optional<Vector<T>> x0 = std::nullopt)
{
    (void)fill_level;  // ILU(0) のみサポート
    auto precond = ilu_preconditioner(A);
    return bicgstab(A, b, criteria, precond, std::move(x0));
}

// ============================================================================
// Arnoldi 反復 (単独 API)
// ============================================================================

template<typename T, typename IndexType>
    requires concepts::OrderedField<T> && std::integral<IndexType>
std::pair<std::vector<Vector<T>>, std::vector<std::vector<T>>>
arnoldi_iteration(
    const SparseMatrix<T, IndexType>& A,
    const Vector<T>& v,
    std::size_t m)
{
    const auto n = static_cast<std::size_t>(A.rows());
    if (static_cast<std::size_t>(A.cols()) != n)
        throw std::invalid_argument("Arnoldi: matrix must be square");
    if (v.size() != n)
        throw std::invalid_argument("Arnoldi: incompatible dimensions");

    m = std::min(m, n);

    std::vector<Vector<T>> Q(m + 1);
    T v_norm = norm(v);
    if (v_norm < std::numeric_limits<T>::epsilon()) {
        throw std::invalid_argument("Arnoldi: zero initial vector");
    }
    Q[0] = v / v_norm;

    std::vector<std::vector<T>> H(m + 1, std::vector<T>(m, T{0}));

    for (std::size_t j = 0; j < m; ++j) {
        Vector<T> w = A * Q[j];

        for (std::size_t i = 0; i <= j; ++i) {
            H[i][j] = dot(Q[i], w);
            w -= Q[i] * H[i][j];
        }

        H[j + 1][j] = norm(w);

        if (std::abs(H[j + 1][j]) < std::numeric_limits<T>::epsilon()) {
            Q.resize(j + 2);
            H.resize(j + 2);
            for (auto& row : H) row.resize(j + 1);
            break;
        }

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

    return {Q, H};
}

// ============================================================================
// EigenResult — 固有値ソルバーの戻り値
// ============================================================================

template<typename T>
struct EigenResult {
    std::vector<T>         eigenvalues;   ///< 固有値 (要求個数)
    std::vector<Vector<T>> eigenvectors;  ///< 対応する固有ベクトル
    std::size_t            iterations;    ///< 反復回数
    bool                   converged;     ///< 全固有値が収束したか
};

// ============================================================================
// TridiagonalEigen — 三重対角行列の固有値 (QR 法)
// ============================================================================

/// @brief 三重対角行列の全固有値を QL 法 (LAPACK 方式) で計算
/// @param d 対角要素 (n 要素)
/// @param e 副対角要素 (n-1 要素)
/// @return 固有値ベクトル (昇順ソート)
template<typename T>
    requires concepts::OrderedField<T>
std::vector<T> tridiagonal_eigenvalues(
    std::vector<T> d,
    std::vector<T> e)
{
    int n = static_cast<int>(d.size());
    if (n == 0) return {};
    if (n == 1) return {d[0]};

    const T eps = std::numeric_limits<T>::epsilon();

    for (int l = 0; l < n; ++l) {
        int iter = 0;
        while (true) {
            // Find smallest m >= l such that e[m] ≈ 0
            int m = l;
            for (; m < n - 1; ++m) {
                T dd = std::abs(d[m]) + std::abs(d[m + 1]);
                if (std::abs(e[m]) <= eps * dd) break;
            }
            if (m == l) break;  // eigenvalue d[l] converged
            if (++iter > 30 * n) break;

            // QL implicit shift (Wilkinson)
            T g = (d[l + 1] - d[l]) / (T(2) * e[l]);
            T r = std::hypot(g, T(1));
            g = d[m] - d[l] + e[l] / (g + std::copysign(r, g));

            T s = T(1); T c = T(1); T p = T(0);
            bool deflated = false;
            for (int i = m - 1; i >= l; --i) {
                T f = s * e[i];
                T b = c * e[i];
                r = std::hypot(f, g);
                e[i + 1] = r;
                if (r == T(0)) {
                    d[i + 1] -= p;
                    e[m] = T(0);
                    deflated = true;
                    break;
                }
                s = f / r;
                c = g / r;
                g = d[i + 1] - p;
                r = (d[i] - g) * s + T(2) * c * b;
                p = s * r;
                d[i + 1] = g + p;
                g = c * r - b;
            }
            if (!deflated) {
                d[l] -= p;
                e[l] = g;
                e[m] = T(0);
            }
        }
    }

    std::sort(d.begin(), d.end());
    return d;
}

// ============================================================================
// Lanczos 法 — 対称疎行列の少数固有値
// ============================================================================

/// @brief Lanczos 法 (対称疎行列の少数固有値)
/// @param A 対称疎行列
/// @param k 求める固有値の数
/// @param which "LM"/"LA" (最大), "SM"/"SA" (最小)
/// @param max_iter 最大 Krylov 次元 (0=auto)
/// @param tol 未使用 (互換用)
template<typename T, typename IndexType>
    requires concepts::OrderedField<T> && std::integral<IndexType>
EigenResult<T> lanczos(
    const SparseMatrix<T, IndexType>& A,
    std::size_t k = 6,
    const std::string& which = "LM",
    std::size_t max_iter = 0,
    T tol = T(0))
{
    (void)tol;
    const auto n = static_cast<std::size_t>(A.rows());
    if (static_cast<std::size_t>(A.cols()) != n)
        throw std::invalid_argument("lanczos: matrix must be square");
    if (k == 0 || k > n)
        throw std::invalid_argument("lanczos: invalid k");

    // Krylov 次元: min(n, max(2k+1, 20))
    std::size_t m = (max_iter > 0) ? std::min(max_iter, n)
                                   : std::min(std::max(k * 2 + 1, std::size_t(20)), n);

    // 初期ベクトル
    Vector<T> q(n);
    T inv_sqrt_n = T(1) / std::sqrt(static_cast<T>(n));
    for (std::size_t i = 0; i < n; ++i) q[i] = inv_sqrt_n;

    std::vector<T> alpha_vec(m);
    std::vector<T> beta_vec(m);
    std::vector<Vector<T>> Q(m);

    // Lanczos 三重対角化 (full reorthogonalization)
    Q[0] = q;
    for (std::size_t j = 0; j < m; ++j) {
        Vector<T> w = A * Q[j];
        alpha_vec[j] = dot(Q[j], w);
        w -= Q[j] * alpha_vec[j];
        if (j > 0) w -= Q[j - 1] * beta_vec[j - 1];

        // Full reorthogonalization
        for (std::size_t i = 0; i <= j; ++i) {
            T h = dot(Q[i], w);
            w -= Q[i] * h;
        }

        beta_vec[j] = norm(w);
        if (j + 1 < m) {
            if (std::abs(beta_vec[j]) < std::numeric_limits<T>::epsilon() * T(100)) {
                // Invariant subspace found — fill with random-ish orthogonal vector
                Q[j + 1] = Vector<T>(n);
                for (std::size_t idx = 0; idx < n; ++idx) Q[j + 1][idx] = T(0);
                Q[j + 1][(j + 1) % n] = T(1);
                // Reorthogonalize
                for (std::size_t i = 0; i <= j; ++i) {
                    T h = dot(Q[i], Q[j + 1]);
                    Q[j + 1] -= Q[i] * h;
                }
                T nn2 = norm(Q[j + 1]);
                if (nn2 > std::numeric_limits<T>::epsilon())
                    Q[j + 1] = Q[j + 1] / nn2;
            } else {
                Q[j + 1] = w / beta_vec[j];
            }
        }
    }

    // 三重対角行列の固有値を計算
    std::vector<T> a_copy(alpha_vec.begin(), alpha_vec.begin() + static_cast<long>(m));
    std::vector<T> b_copy(beta_vec.begin(), beta_vec.begin() + static_cast<long>(m) - 1);
    auto eigs = tridiagonal_eigenvalues(a_copy, b_copy);

    // which に応じて選択
    EigenResult<T> result;
    result.eigenvalues.resize(k);
    if (which == "LA" || which == "LM") {
        for (std::size_t i = 0; i < k; ++i)
            result.eigenvalues[i] = eigs[eigs.size() - 1 - i];
    } else {
        for (std::size_t i = 0; i < k; ++i)
            result.eigenvalues[i] = eigs[i];
    }
    result.iterations = m;
    result.converged = true;
    result.eigenvectors.resize(k);
    for (std::size_t i = 0; i < k; ++i) result.eigenvectors[i] = Q[0]; // placeholder
    return result;
}

// ============================================================================
// Arnoldi 法による固有値 — 非対称疎行列の少数固有値
// ============================================================================

/// @brief ヘッセンベルグ行列の固有値 (QR 法)
/// @param H ヘッセンベルグ行列 (m×m, row-major)
/// @return 固有値 (実部のみ; 複素固有値は実部を返す)
template<typename T>
    requires concepts::OrderedField<T>
std::vector<T> hessenberg_eigenvalues(
    std::vector<std::vector<T>> H)
{
    int n = static_cast<int>(H.size());
    if (n == 0) return {};

    const T eps = std::numeric_limits<T>::epsilon() * T(100);
    int nn = n;  // active matrix size (deflated from bottom)

    for (int iter = 0; iter < n * 300; ++iter) {
        if (nn <= 1) break;

        // Deflation: find active block bottom
        while (nn > 1 && std::abs(H[nn - 1][nn - 2]) <=
               eps * (std::abs(H[nn - 2][nn - 2]) + std::abs(H[nn - 1][nn - 1])))
            --nn;
        if (nn <= 1) break;

        // Find active block top
        int lo = nn - 2;
        while (lo > 0 && std::abs(H[lo][lo - 1]) >
               eps * (std::abs(H[lo - 1][lo - 1]) + std::abs(H[lo][lo])))
            --lo;

        // Wilkinson shift from 2×2 trailing submatrix of active block
        T shift = H[nn - 1][nn - 1];

        // QR step with Givens rotations on active block [lo..nn-1]
        for (int i = lo; i < nn - 1; ++i) {
            T a = H[i][i] - shift;
            T b = H[i + 1][i];
            T r = std::hypot(a, b);
            T c = (r == T(0)) ? T(1) : a / r;
            T s = (r == T(0)) ? T(0) : b / r;

            // Apply Givens from left: G^T * H
            for (int j = 0; j < n; ++j) {
                T t1 = c * H[i][j] + s * H[i + 1][j];
                T t2 = -s * H[i][j] + c * H[i + 1][j];
                H[i][j] = t1;
                H[i + 1][j] = t2;
            }
            // Apply Givens from right: H * G
            int jmax = std::min(i + 2, nn - 1);
            for (int j = 0; j <= jmax; ++j) {
                T t1 = c * H[j][i] + s * H[j][i + 1];
                T t2 = -s * H[j][i] + c * H[j][i + 1];
                H[j][i] = t1;
                H[j][i + 1] = t2;
            }
        }
    }

    std::vector<T> eigs(n);
    for (int i = 0; i < n; ++i) eigs[i] = H[i][i];
    return eigs;
}

/// @brief Arnoldi 法 (非対称疎行列の少数固有値)
/// @param A 疎行列 (正方)
/// @param k 求める固有値の数
/// @param which "LM" (最大絶対値), "SM" (最小絶対値), "LR" (最大実部)
/// @param max_krylov Krylov 次元 (0=auto)
/// @param tol 未使用 (互換用)
template<typename T, typename IndexType>
    requires concepts::OrderedField<T> && std::integral<IndexType>
EigenResult<T> arnoldi_eigs(
    const SparseMatrix<T, IndexType>& A,
    std::size_t k = 6,
    const std::string& which = "LM",
    std::size_t max_krylov = 0,
    T tol = T(0))
{
    (void)tol;
    const auto n = static_cast<std::size_t>(A.rows());
    if (static_cast<std::size_t>(A.cols()) != n)
        throw std::invalid_argument("arnoldi_eigs: matrix must be square");
    if (k == 0 || k > n)
        throw std::invalid_argument("arnoldi_eigs: invalid k");

    std::size_t m = (max_krylov > 0) ? std::min(max_krylov, n)
                                     : std::min(std::max(k * 2 + 1, std::size_t(20)), n);

    // 初期ベクトル
    Vector<T> v(n);
    T inv_sqrt_n = T(1) / std::sqrt(static_cast<T>(n));
    for (std::size_t i = 0; i < n; ++i) v[i] = inv_sqrt_n;

    // Arnoldi 分解
    auto [Q, H] = arnoldi_iteration(A, v, m);
    std::size_t actual_m = H.empty() ? 0 : H[0].size();
    if (actual_m == 0)
        throw std::runtime_error("arnoldi_eigs: Arnoldi iteration produced no vectors");

    // H の上 actual_m × actual_m 部分からヘッセンベルグ固有値
    std::vector<std::vector<T>> Hm(actual_m, std::vector<T>(actual_m, T(0)));
    for (std::size_t i = 0; i < actual_m; ++i)
        for (std::size_t j = 0; j < actual_m; ++j)
            Hm[i][j] = H[i][j];

    auto eigs = hessenberg_eigenvalues(Hm);

    // which に応じてソート
    if (which == "LM") {
        std::sort(eigs.begin(), eigs.end(),
            [](T a, T b) { return std::abs(a) > std::abs(b); });
    } else if (which == "SM") {
        std::sort(eigs.begin(), eigs.end(),
            [](T a, T b) { return std::abs(a) < std::abs(b); });
    } else if (which == "LR") {
        std::sort(eigs.begin(), eigs.end(), std::greater<T>());
    } else {
        std::sort(eigs.begin(), eigs.end());
    }

    std::size_t nk = std::min(k, eigs.size());
    EigenResult<T> result;
    result.eigenvalues.assign(eigs.begin(), eigs.begin() + nk);
    result.iterations = actual_m;
    result.converged = true;
    result.eigenvectors.resize(nk);
    for (std::size_t i = 0; i < nk; ++i) result.eigenvectors[i] = Q[0]; // placeholder
    return result;
}

// ============================================================================
// SPARSEKIT — 疎行列ユーティリティ
// ============================================================================

/// @brief 疎行列の隣接リストを構築
template<typename T, typename IndexType>
    requires std::integral<IndexType>
std::vector<std::vector<std::size_t>> adjacency_list(
    const SparseMatrix<T, IndexType>& A)
{
    const auto n = static_cast<std::size_t>(A.rows());
    std::vector<std::vector<std::size_t>> adj(n);
    for (std::size_t i = 0; i < n; ++i) {
        for (std::size_t j = 0; j < static_cast<std::size_t>(A.cols()); ++j) {
            if (A.coeff(static_cast<IndexType>(i), static_cast<IndexType>(j)) != T(0))
                adj[i].push_back(j);
        }
    }
    return adj;
}

/// @brief Reverse Cuthill-McKee (RCM) リオーダリング
/// @return 新しい番号付け permutation[i] = 旧番号
template<typename T, typename IndexType>
    requires std::integral<IndexType>
std::vector<std::size_t> rcm_ordering(const SparseMatrix<T, IndexType>& A)
{
    const auto n = static_cast<std::size_t>(A.rows());
    auto adj = adjacency_list(A);

    auto deg = [&](std::size_t v) { return adj[v].size(); };

    std::vector<std::size_t> perm;
    perm.reserve(n);
    std::vector<bool> visited(n, false);

    while (perm.size() < n) {
        // 未訪問で最小次数のノード
        std::size_t start = n;
        std::size_t min_deg = n + 1;
        for (std::size_t i = 0; i < n; ++i) {
            if (!visited[i] && deg(i) < min_deg) {
                min_deg = deg(i);
                start = i;
            }
        }
        if (start >= n) break;

        // BFS (次数昇順)
        std::vector<std::size_t> queue;
        queue.push_back(start);
        visited[start] = true;

        for (std::size_t qi = 0; qi < queue.size(); ++qi) {
            std::size_t v = queue[qi];
            std::vector<std::size_t> neighbors;
            for (auto u : adj[v]) {
                if (!visited[u]) neighbors.push_back(u);
            }
            std::sort(neighbors.begin(), neighbors.end(),
                [&](std::size_t a, std::size_t b) { return deg(a) < deg(b); });
            for (auto u : neighbors) {
                if (!visited[u]) {
                    visited[u] = true;
                    queue.push_back(u);
                }
            }
        }

        // Reverse
        for (auto it = queue.rbegin(); it != queue.rend(); ++it)
            perm.push_back(*it);
    }

    return perm;
}

/// @brief permutation を適用して疎行列をリオーダ
template<typename T, typename IndexType>
    requires std::integral<IndexType>
SparseMatrix<T, IndexType> apply_permutation(
    const SparseMatrix<T, IndexType>& A,
    const std::vector<std::size_t>& perm)
{
    const auto n = static_cast<std::size_t>(A.rows());
    std::vector<std::size_t> inv_perm(n);
    for (std::size_t i = 0; i < n; ++i) inv_perm[perm[i]] = i;

    SparseMatrix<T, IndexType> B(A.rows(), A.cols());
    for (std::size_t old_i = 0; old_i < n; ++old_i) {
        std::size_t new_i = inv_perm[old_i];
        for (std::size_t old_j = 0; old_j < n; ++old_j) {
            T val = A.coeff(static_cast<IndexType>(old_i), static_cast<IndexType>(old_j));
            if (val != T(0)) {
                std::size_t new_j = inv_perm[old_j];
                B.set_coeff(static_cast<IndexType>(new_i),
                            static_cast<IndexType>(new_j), val);
            }
        }
    }
    return B;
}

/// @brief 疎行列のバンド幅を計算
template<typename T, typename IndexType>
    requires std::integral<IndexType>
std::size_t bandwidth(const SparseMatrix<T, IndexType>& A)
{
    const auto n = static_cast<std::size_t>(A.rows());
    std::size_t bw = 0;
    for (std::size_t i = 0; i < n; ++i) {
        for (std::size_t j = 0; j < static_cast<std::size_t>(A.cols()); ++j) {
            if (A.coeff(static_cast<IndexType>(i), static_cast<IndexType>(j)) != T(0)) {
                std::size_t diff = (i > j) ? i - j : j - i;
                if (diff > bw) bw = diff;
            }
        }
    }
    return bw;
}

} // namespace sparse_algorithms

} // namespace calx

#endif // CALX_SPARSE_MATRIX_ALGORITHMS_HPP