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

// solvers.hpp
//
// 線形方程式系のソルバー
//
// このファイルでは、様々な線形方程式系のソルバーを実装します。
// 直接法および反復法のソルバーを含みます。
//
// 主な機能:
// - 直接法ソルバー（LU, Cholesky, QR, SVD）
// - 反復法ソルバー（共役勾配法, Jacobi法, Gauss-Seidel法）
// - 特殊行列向けソルバー（三重対角行列、帯行列など）

#ifndef CALX_SOLVERS_HPP
#define CALX_SOLVERS_HPP

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

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

namespace calx {
    namespace algorithms {

        // 前方宣言
        template<typename T>
        Vector<T> conjugate_gradient(const Matrix<T>& a, const Vector<T>& b,
            T tolerance,
            typename Matrix<T>::size_type max_iterations);

        //-----------------------------------------------------------------------------
        // ソルバー選択のためのユーティリティ
        //-----------------------------------------------------------------------------

        // ソルバータイプの列挙型
        enum class SolverType {
            Auto,     // 最適なソルバーを自動選択
            LU,       // LU分解ソルバー
            Cholesky, // Cholesky分解ソルバー
            QR,       // QR分解ソルバー
            SVD,      // 特異値分解ソルバー
            LDL,      // LDL分解ソルバー
            Iterative // 反復法ソルバー
        };

        // 文字列からソルバータイプへの変換
        inline SolverType parse_solver_type(const std::string& solver_name) {
            if (solver_name == "auto") return SolverType::Auto;
            if (solver_name == "lu") return SolverType::LU;
            if (solver_name == "cholesky") return SolverType::Cholesky;
            if (solver_name == "qr") return SolverType::QR;
            if (solver_name == "svd") return SolverType::SVD;
            if (solver_name == "ldl") return SolverType::LDL;
            if (solver_name == "iterative") return SolverType::Iterative;

            // デフォルトは自動選択
            return SolverType::Auto;
        }

        //-----------------------------------------------------------------------------
        // 直接法ソルバー
        //-----------------------------------------------------------------------------

        // 汎用ソルバーインターフェース (Ax = b)
        template<typename T>
        Vector<T> solve(const Matrix<T>& a, const Vector<T>& b,
            SolverType solver_type = SolverType::Auto) {
            // 次元チェック
            if (a.rows() != b.size()) {
                throw DimensionError("solve: matrix and vector dimensions mismatch");
            }

            // ソルバー選択ロジック
            if (solver_type == SolverType::Auto) {
                if (a.is_square()) {
                    // 対称性チェック
                    bool is_symmetric = true;
                    for (typename Matrix<T>::size_type i = 0; i < a.rows() && is_symmetric; ++i) {
                        for (typename Matrix<T>::size_type j = i + 1; j < a.cols() && is_symmetric; ++j) {
                            if (std::abs(a(i, j) - a(j, i)) > std::numeric_limits<T>::epsilon() * 10) {
                                is_symmetric = false;
                            }
                        }
                    }

                    if (is_symmetric) {
                        // 正定値チェック（簡易的な方法として対角要素が正かどうかをチェック）
                        bool is_positive_definite = true;
                        for (typename Matrix<T>::size_type i = 0; i < a.rows() && is_positive_definite; ++i) {
                            if (a(i, i) <= 0) {
                                is_positive_definite = false;
                            }
                        }

                        if (is_positive_definite) {
                            // 対称正定値行列にはCholesky分解が最適
                            return cholesky_solve(a, b);
                        }
                        else {
                            // 対称行列にはLDL分解が適している
                            return ldl_solve(a, b);
                        }
                    }
                    else {
                        // 非対称行列にはLU分解が一般的
                        return lu_solve(a, b);
                    }
                }
                else {
                    // 非正方行列には最小二乗解（QRまたはSVD）
                    return qr_solve(a, b);
                }
            }

            // 明示的なソルバー選択
            switch (solver_type) {
            case SolverType::LU:
                return lu_solve(a, b);
            case SolverType::Cholesky:
                return cholesky_solve(a, b);
            case SolverType::QR:
                return qr_solve(a, b);
            case SolverType::SVD:
                return svd_solve(a, b);
            case SolverType::LDL:
                return ldl_solve(a, b);
            case SolverType::Iterative:
                return conjugate_gradient(a, b);
            default:
                // 到達不能
                throw std::runtime_error("solve: unknown solver type");
            }
        }

        // 文字列ベースのソルバー選択インターフェース
        template<typename T>
        Vector<T> solve(const Matrix<T>& a, const Vector<T>& b, const std::string& solver_name) {
            return solve(a, b, parse_solver_type(solver_name));
        }

        // 複数右辺ベクトルのソルバー (AX = B)
        template<typename T>
        Matrix<T> solve_multi(const Matrix<T>& a, const Matrix<T>& b,
            SolverType solver_type = SolverType::Auto) {
            // 次元チェック
            if (a.rows() != b.rows()) {
                throw DimensionError("solve_multi: matrix dimensions mismatch");
            }

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

#if CALX_HAS_MKL
            if constexpr (std::is_same_v<T, float> || std::is_same_v<T, double>) {
                if (a.is_square() && (solver_type == SolverType::Auto || solver_type == SolverType::LU)) {
                    // LAPACKE_gesv: 複数右辺を一括で解く
                    Matrix<T> Acopy = a;
                    Matrix<T> x = b;
                    std::vector<lapack_int> ipiv(n);
                    lapack_int info;
                    if constexpr (std::is_same_v<T, float>)
                        info = LAPACKE_sgesv(LAPACK_ROW_MAJOR, (lapack_int)n, (lapack_int)nrhs,
                                             Acopy.data(), (lapack_int)n, ipiv.data(),
                                             x.data(), (lapack_int)nrhs);
                    else
                        info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, (lapack_int)n, (lapack_int)nrhs,
                                             Acopy.data(), (lapack_int)n, ipiv.data(),
                                             x.data(), (lapack_int)nrhs);
                    if (info > 0)
                        throw MathError("solve_multi: singular matrix detected");
                    return x;
                }
            }
#endif

            // 結果行列の初期化
            Matrix<T> x(n, nrhs);

            // 各右辺ベクトルに対してソルバーを適用
            Vector<T> bi(m);
            Vector<T> xi(n);

            for (typename Matrix<T>::size_type j = 0; j < nrhs; ++j) {
                // 右辺ベクトルの抽出
                for (typename Matrix<T>::size_type i = 0; i < m; ++i) {
                    bi[i] = b(i, j);
                }

                // 方程式を解く
                xi = solve(a, bi, solver_type);

                // 結果の格納
                for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                    x(i, j) = xi[i];
                }
            }

            return x;
        }

        // 文字列ベースの複数右辺ソルバーインターフェース
        template<typename T>
        Matrix<T> solve_multi(const Matrix<T>& a, const Matrix<T>& b, const std::string& solver_name) {
            return solve_multi(a, b, parse_solver_type(solver_name));
        }

        //-----------------------------------------------------------------------------
        // 反復法ソルバー
        //-----------------------------------------------------------------------------

        // 共役勾配法
        template<typename T>
        Vector<T> conjugate_gradient(const Matrix<T>& a, const Vector<T>& b,
            T tolerance = std::numeric_limits<T>::epsilon() * 10,
            typename Matrix<T>::size_type max_iterations = 1000) {
            // 次元チェック
            if (a.rows() != a.cols()) {
                throw DimensionError("conjugate_gradient: matrix must be square");
            }

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

            const auto n = a.rows();
            Vector<T> x(n, T(0));     // 初期解をゼロに設定
            Vector<T> r = b;          // 初期残差 r = b - A*x = b
            Vector<T> p = r;          // 初期探索方向
            Vector<T> ap(n);          // A*p バッファ (ループ外で1回確保)

            T r_norm = dot(r, r);     // 残差ノルム（二乗）
            const T b_norm = dot(b, b); // 右辺ベクトルノルム（二乗）

            // ゼロベクトルチェック
            if (b_norm < std::numeric_limits<T>::epsilon()) {
                return x;  // 右辺がゼロならゼロベクトルが解
            }

            // 相対許容誤差（二乗）: ‖r‖² ≤ (tolerance * ‖b‖)²
            const T rel_tol = tolerance * tolerance * b_norm;

            // 反復計算
            for (typename Matrix<T>::size_type iter = 0; iter < max_iterations; ++iter) {
                // 収束チェック
                if (r_norm <= rel_tol) {
                    break;
                }

                // A*p を計算 (ap はループ外で確保済み)
                for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                    T s = T(0);
                    for (typename Matrix<T>::size_type j = 0; j < n; ++j)
                        s += a(i, j) * p[j];
                    ap[i] = s;
                }

                // ステップサイズα = r^T*r / (p^T*A*p)
                const T pAp = dot(p, ap);
                if (std::abs(pAp) < std::numeric_limits<T>::epsilon()) break;
                const T alpha = r_norm / pAp;

                // 解と残差の更新 (融合ループ)
                axpy(alpha, p, x);
                axpy(-alpha, ap, r);

                // 新しい残差ノルム
                const T r_next_norm = dot(r, r);

                // 方向ベクトルの更新: p = r + beta * p
                if (r_norm < std::numeric_limits<T>::epsilon()) break;  // ゼロ除算回避
                const T beta = r_next_norm / r_norm;
                axpby(T(1), r, beta, p);

                // 残差ノルムの更新
                r_norm = r_next_norm;
            }

            return x;
        }

        // Jacobi法（ヤコビ法）
        template<typename T>
        Vector<T> jacobi_method(const Matrix<T>& a, const Vector<T>& b,
            T tolerance = std::numeric_limits<T>::epsilon() * 10,
            typename Matrix<T>::size_type max_iterations = 1000) {
            // 次元チェック
            if (a.rows() != a.cols()) {
                throw DimensionError("jacobi_method: matrix must be square");
            }

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

            const auto n = a.rows();
            Vector<T> x(n, T(0));    // 初期解をゼロに設定
            Vector<T> x_new(n);      // 新しい解

            // 対角成分のゼロチェック
            for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                if (std::abs(a(i, i)) < std::numeric_limits<T>::epsilon()) {
                    throw MathError("jacobi_method: zero diagonal element");
                }
            }

            // 反復計算
            for (typename Matrix<T>::size_type iter = 0; iter < max_iterations; ++iter) {
                // 新しい解の計算
                for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                    T sum = b[i];
                    for (typename Matrix<T>::size_type j = 0; j < n; ++j) {
                        if (i != j) {
                            sum -= a(i, j) * x[j];
                        }
                    }
                    x_new[i] = sum / a(i, i);
                }

                // 収束チェック
                T diff_norm = 0;
                for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                    diff_norm += (x_new[i] - x[i]) * (x_new[i] - x[i]);
                }
                diff_norm = std::sqrt(diff_norm);

                if (diff_norm < tolerance) {
                    break;
                }

                // 解の更新
                std::swap(x, x_new);
            }

            return x;
        }

        // Gauss-Seidel法
        template<typename T>
        Vector<T> gauss_seidel_method(const Matrix<T>& a, const Vector<T>& b,
            T tolerance = std::numeric_limits<T>::epsilon() * 10,
            typename Matrix<T>::size_type max_iterations = 1000) {
            // 次元チェック
            if (a.rows() != a.cols()) {
                throw DimensionError("gauss_seidel_method: matrix must be square");
            }

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

            const auto n = a.rows();
            Vector<T> x(n, T(0));   // 初期解をゼロに設定
            Vector<T> x_old(n);     // 前回の解

            // 対角成分のゼロチェック
            for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                if (std::abs(a(i, i)) < std::numeric_limits<T>::epsilon()) {
                    throw MathError("gauss_seidel_method: zero diagonal element");
                }
            }

            // 反復計算
            for (typename Matrix<T>::size_type iter = 0; iter < max_iterations; ++iter) {
                // 前回の解を保存
                std::copy(x.begin(), x.end(), x_old.begin());

                // 新しい解の計算
                for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                    T sum = b[i];

                    // i未満の要素（すでに更新済み）
                    for (typename Matrix<T>::size_type j = 0; j < i; ++j) {
                        sum -= a(i, j) * x[j];
                    }

                    // i超の要素（まだ更新されていない）
                    for (typename Matrix<T>::size_type j = i + 1; j < n; ++j) {
                        sum -= a(i, j) * x_old[j];
                    }

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

                // 収束チェック
                T diff_norm = 0;
                for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                    diff_norm += (x[i] - x_old[i]) * (x[i] - x_old[i]);
                }
                diff_norm = std::sqrt(diff_norm);

                if (diff_norm < tolerance) {
                    break;
                }
            }

            return x;
        }

        // SOR法（逐次過緩和法）
        template<typename T>
        Vector<T> sor_method(const Matrix<T>& a, const Vector<T>& b,
            T omega = static_cast<T>(1.5),  // 緩和係数
            T tolerance = std::numeric_limits<T>::epsilon() * 10,
            typename Matrix<T>::size_type max_iterations = 1000) {
            // 次元チェック
            if (a.rows() != a.cols()) {
                throw DimensionError("sor_method: matrix must be square");
            }

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

            // 緩和係数の範囲チェック
            if (omega <= 0 || omega >= 2) {
                throw std::invalid_argument("sor_method: relaxation parameter must be in (0,2)");
            }

            const auto n = a.rows();
            Vector<T> x(n, T(0));   // 初期解をゼロに設定
            Vector<T> x_old(n);     // 前回の解

            // 対角成分のゼロチェック
            for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                if (std::abs(a(i, i)) < std::numeric_limits<T>::epsilon()) {
                    throw MathError("sor_method: zero diagonal element");
                }
            }

            // 反復計算
            for (typename Matrix<T>::size_type iter = 0; iter < max_iterations; ++iter) {
                // 前回の解を保存
                std::copy(x.begin(), x.end(), x_old.begin());

                // 新しい解の計算
                for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                    T sum1 = 0;  // 更新済みの要素の寄与
                    T sum2 = 0;  // 未更新の要素の寄与

                    // i未満の要素（すでに更新済み）
                    for (typename Matrix<T>::size_type j = 0; j < i; ++j) {
                        sum1 += a(i, j) * x[j];
                    }

                    // i超の要素（まだ更新されていない）
                    for (typename Matrix<T>::size_type j = i + 1; j < n; ++j) {
                        sum2 += a(i, j) * x_old[j];
                    }

                    // SOR更新式
                    T x_gs = (b[i] - sum1 - sum2) / a(i, i);  // Gauss-Seidel解
                    x[i] = x_old[i] + omega * (x_gs - x_old[i]);
                }

                // 収束チェック
                T diff_norm = 0;
                for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                    diff_norm += (x[i] - x_old[i]) * (x[i] - x_old[i]);
                }
                diff_norm = std::sqrt(diff_norm);

                if (diff_norm < tolerance) {
                    break;
                }
            }

            return x;
        }

        //-----------------------------------------------------------------------------
        // 最小二乗法ソルバー
        //-----------------------------------------------------------------------------

        // 通常の最小二乗法ソルバー
        template<typename T>
        Vector<T> least_squares(const Matrix<T>& a, const Vector<T>& b,
            SolverType solver_type = SolverType::Auto) {
            // 次元チェック
            if (a.rows() != b.size()) {
                throw DimensionError("least_squares: matrix and vector dimensions mismatch");
            }

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

            // ソルバー選択ロジック
            if (solver_type == SolverType::Auto) {
                if (m >= n) {  // 過剰決定または正確な決定
                    // QR分解が数値的に安定
                    return qr_solve(a, b);
                }
                else {  // 劣決定系
                    // SVD分解が最小ノルム解を提供
                    return svd_solve(a, b);
                }
            }

            // 明示的なソルバー選択
            switch (solver_type) {
            case SolverType::QR:
                return qr_solve(a, b);
            case SolverType::SVD:
                return svd_solve(a, b);
            default:
                // 最小二乗問題には通常QRまたはSVDが使用される
                return qr_solve(a, b);
            }
        }

        // 文字列ベースの最小二乗法ソルバーインターフェース
        template<typename T>
        Vector<T> least_squares(const Matrix<T>& a, const Vector<T>& b, const std::string& solver_name) {
            return least_squares(a, b, parse_solver_type(solver_name));
        }

        // 重み付き最小二乗法ソルバー
        template<typename T>
        Vector<T> weighted_least_squares(const Matrix<T>& a, const Vector<T>& b,
            const Vector<T>& weights) {
            // 次元チェック
            if (a.rows() != b.size() || a.rows() != weights.size()) {
                throw DimensionError("weighted_least_squares: dimensions mismatch");
            }

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

            // 重み行列（対角行列）とデータの変換
            Matrix<T> a_weighted(m, n);
            Vector<T> b_weighted(m);

            for (typename Matrix<T>::size_type i = 0; i < m; ++i) {
                // ゼロ重みのチェック（数値的に安定化のため）
                const T weight = std::max(weights[i], std::numeric_limits<T>::epsilon());
                const T sqrt_weight = std::sqrt(weight);

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

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

            // 通常の最小二乗問題として解く
            return least_squares(a_weighted, b_weighted);
        }

        // 正則化最小二乗法ソルバー（リッジ回帰）
        template<typename T>
        Vector<T> regularized_least_squares(const Matrix<T>& a, const Vector<T>& b,
            T lambda = static_cast<T>(0.1)) {
            // 次元チェック
            if (a.rows() != b.size()) {
                throw DimensionError("regularized_least_squares: dimensions mismatch");
            }

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

            // 正規方程式: (A^T * A + lambda * I) * x = A^T * b

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

            // 対角項に正則化項を追加
            for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                ata(i, i) += lambda;
            }

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

            // 正規方程式を解く
            return cholesky_solve(ata, atb);
        }

        //-----------------------------------------------------------------------------
        // 特殊な行列のためのソルバー
        //-----------------------------------------------------------------------------

        // 三重対角行列ソルバー
        template<typename T>
        Vector<T> tridiagonal_solve(const Vector<T>& a, const Vector<T>& b,
            const Vector<T>& c, const Vector<T>& d) {
            // a: 対角成分の下
            // b: 対角成分
            // c: 対角成分の上
            // d: 右辺ベクトル

            const auto n = b.size();

            // 次元チェック
            if (a.size() != n - 1 || c.size() != n - 1 || d.size() != n) {
                throw DimensionError("tridiagonal_solve: dimensions mismatch");
            }

            // トーマスアルゴリズムの実装
            Vector<T> c_prime(n - 1);
            Vector<T> d_prime(n);
            Vector<T> x(n);

            // 前進消去
            c_prime[0] = c[0] / b[0];
            d_prime[0] = d[0] / b[0];

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

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

            // 後退代入
            x[n - 1] = d_prime[n - 1];

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

            return x;
        }

        // 五重対角行列ソルバー (4階微分の有限差分で出現)
        //
        // | b0 c0 d0  0  0 ... | | x0 |   | f0 |
        // | a0 b1 c1 d1  0 ... | | x1 |   | f1 |
        // | e0 a1 b2 c2 d2 ... | | x2 | = | f2 |
        // |  0 e1 a2 b3 c3 ... | | x3 |   | f3 |
        // | ...                 | | .. |   | .. |
        //
        // e[i]: 対角 -2, a[i]: 対角 -1, b[i]: 対角, c[i]: 対角 +1, d[i]: 対角 +2
        //
        // アルゴリズム: 帯幅 2 の LU 分解 (Gauss 消去の五重対角特化版)
        // 計算量: O(n), メモリ: O(n)
        template<typename T>
        Vector<T> pentadiagonal_solve(
                const Vector<T>& e, const Vector<T>& a, const Vector<T>& b,
                const Vector<T>& c, const Vector<T>& d, const Vector<T>& f) {
            const auto n = b.size();
            if (n == 0) return {};

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

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

            // 次元チェック
            if (e.size() != n - 2 || a.size() != n - 1 ||
                c.size() != n - 1 || d.size() != n - 2 || f.size() != n) {
                throw DimensionError("pentadiagonal_solve: dimensions mismatch "
                    "(e:" + std::to_string(e.size()) + " a:" + std::to_string(a.size()) +
                    " b:" + std::to_string(b.size()) + " c:" + std::to_string(c.size()) +
                    " d:" + std::to_string(d.size()) + " f:" + std::to_string(f.size()) + ")");
            }

            // 作業配列 (LU 分解の上三角部分 + 変換された RHS)
            std::vector<T> al(n - 1);  // L の副対角
            std::vector<T> be(n);      // U の対角
            std::vector<T> ga(n - 1);  // U の超対角 +1
            std::vector<T> de(n - 2);  // U の超対角 +2
            std::vector<T> rhs(n);     // 変換された右辺

            // 初期化
            be[0] = b[0];
            ga[0] = c[0];
            de[0] = d[0];
            rhs[0] = f[0];

            // i = 1: e[0] が関与
            T mu = a[0] / be[0];
            be[1] = b[1] - mu * ga[0];
            ga[1] = c[1] - mu * de[0];
            if (n > 3) de[1] = d[1];  // de[1] は d[1] のまま (mu * 0 = 0)
            rhs[1] = f[1] - mu * rhs[0];
            al[0] = mu;

            // i = 2 .. n-1: 前進消去 (2行分の消去)
            for (size_t i = 2; i < n; ++i) {
                // 行 i の対角 -2 (e[i-2]) を消去
                T mu1 = e[i - 2] / be[i - 2];
                // 更新後の a'[i-1], b'[i], c'[i], d'[i], f'[i]
                T a_new = a[i - 1] - mu1 * ga[i - 2];
                T b_new = b[i] - mu1 * ((i - 2 < n - 2) ? de[i - 2] : T(0));
                T f_new = f[i] - mu1 * rhs[i - 2];

                // 行 i の対角 -1 (a_new) を消去
                T mu2 = a_new / be[i - 1];
                be[i] = b_new - mu2 * ga[i - 1];
                if (i < n - 1) {
                    T c_i = (i < c.size()) ? c[i] : T(0);
                    ga[i] = c_i - mu2 * ((i - 1 < n - 2) ? de[i - 1] : T(0));
                }
                if (i < n - 2) {
                    de[i] = d[i];
                }
                rhs[i] = f_new - mu2 * rhs[i - 1];
                if (i < n - 1) al[i - 1] = mu2;
            }

            // 後退代入
            Vector<T> x(n);
            x[n - 1] = rhs[n - 1] / be[n - 1];

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

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

            return x;
        }

        // 帯行列ソルバー
/**
 * @brief 改良版帯行列ソルバー
 * @tparam T 要素の型
 * @param a 帯行列
 * @param b 右辺ベクトル
 * @param kl 下帯域幅（対角成分より下の非ゼロ要素の数）
 * @param ku 上帯域幅（対角成分より上の非ゼロ要素の数）
 * @return 解ベクトル
 */
        template<typename T>
        Vector<T> band_solve(const Matrix<T>& a, const Vector<T>& b,
            typename Matrix<T>::size_type kl,
            typename Matrix<T>::size_type ku) {

            const auto n = a.rows();

            // 次元チェック
            if (n != a.cols() || n != b.size()) {
                throw DimensionError("band_solve: dimensions mismatch");
            }

            // 三重対角行列の特別処理（トーマスアルゴリズム）
            if (kl == 1 && ku == 1) {
                Vector<T> diag_lower(n - 1);
                Vector<T> diag_main(n);
                Vector<T> diag_upper(n - 1);

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

                // 安定性向上のためにスケーリング処理を追加
                T scale = T(0);
                for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                    scale = std::max(scale, std::abs(diag_main[i]));
                    if (i < n - 1) {
                        scale = std::max(scale, std::abs(diag_upper[i]));
                        scale = std::max(scale, std::abs(diag_lower[i]));
                    }
                }

                // スケーリング係数がゼロに近い場合の特別処理
                if (scale < std::numeric_limits<T>::epsilon()) {
                    return Vector<T>(n, T(0));  // すべてゼロの解を返す
                }

                // スケーリングされた対角要素で計算
                Vector<T> scaled_lower(n - 1);
                Vector<T> scaled_main(n);
                Vector<T> scaled_upper(n - 1);
                Vector<T> scaled_b(n);

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

                // トーマスアルゴリズムの改良版（数値的安定性向上）
                return tridiagonal_solve(scaled_lower, scaled_main, scaled_upper, scaled_b);
            }

            // 帯行列専用の LU 分解（部分ピボット選択）
            // 部分ピボットにより U の帯幅は最大 kl+ku に拡大する（fill-in）
            Matrix<T> lu(n, n, T(0));

            // 元の帯行列をコピー
            for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                const auto j_start = (i > kl) ? i - kl : typename Matrix<T>::size_type(0);
                const auto j_end = std::min(n, i + ku + 1);
                for (typename Matrix<T>::size_type j = j_start; j < j_end; ++j) {
                    lu(i, j) = a(i, j);
                }
            }

            // ピボット記録: pivots[k] = ステップ k で交換した行
            std::vector<typename Matrix<T>::size_type> pivots(n);
            for (typename Matrix<T>::size_type i = 0; i < n; ++i)
                pivots[i] = i;

            // LU 分解（部分ピボット）
            for (typename Matrix<T>::size_type k = 0; k < n - 1; ++k) {
                // 列 k の k 行目以降（帯域 kl 以内）で最大ピボットを探す
                typename Matrix<T>::size_type pivot_row = k;
                T max_val = std::abs(lu(k, k));
                const auto search_end = std::min(n, k + kl + 1);

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

                pivots[k] = pivot_row;

                // 行交換 (帯域内のみで十分だが、fill-in を考慮して全列交換)
                if (pivot_row != k) {
                    const auto swap_start = (k > kl) ? k - kl : typename Matrix<T>::size_type(0);
                    const auto swap_end = std::min(n, k + kl + ku + 1);
                    for (typename Matrix<T>::size_type j = swap_start; j < swap_end; ++j)
                        std::swap(lu(k, j), lu(pivot_row, j));
                }

                // ピボットがほぼゼロなら特異行列
                if (std::abs(lu(k, k)) < std::numeric_limits<T>::epsilon() * T(100)) {
                    continue;
                }

                // 消去（fill-in を考慮: U の帯幅は最大 kl+ku）
                for (typename Matrix<T>::size_type i = k + 1; i < search_end; ++i) {
                    lu(i, k) /= lu(k, k);
                    const auto j_end = std::min(n, k + kl + ku + 1);
                    for (typename Matrix<T>::size_type j = k + 1; j < j_end; ++j) {
                        lu(i, j) -= lu(i, k) * lu(k, j);
                    }
                }
            }

            // 解の計算: Ly = Pb, Ux = y
            Vector<T> x = b;

            // ピボット置換を b に順次適用
            for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                if (pivots[i] != i)
                    std::swap(x[i], x[pivots[i]]);
            }

            // 前進代入 (L の帯幅は kl)
            for (typename Matrix<T>::size_type i = 1; i < n; ++i) {
                const auto j_start = (i > kl) ? i - kl : typename Matrix<T>::size_type(0);
                for (typename Matrix<T>::size_type j = j_start; j < i; ++j) {
                    x[i] -= lu(i, j) * x[j];
                }
            }

            // 後退代入 (U の帯幅は kl+ku: fill-in を考慮)
            for (typename Matrix<T>::size_type i = n; i-- > 0; ) {
                const auto j_end = std::min(n, i + kl + ku + 1);
                for (typename Matrix<T>::size_type j = i + 1; j < j_end; ++j) {
                    x[i] -= lu(i, j) * x[j];
                }
                x[i] /= lu(i, i);
            }

            return x;
        }

        //-----------------------------------------------------------------------------
        // 疎行列のソルバー
        //-----------------------------------------------------------------------------

        // CSRフォーマットの疎行列ソルバー
/**
 * @brief 改良版CSR形式の疎行列ソルバー
 * @tparam T 要素の型
 * @param values CSR形式の非ゼロ要素の値
 * @param row_ptr 行ポインタ配列
 * @param col_idx 列インデックス配列
 * @param b 右辺ベクトル
 * @param tolerance 収束判定閾値
 * @param max_iterations 最大反復回数
 * @return 解ベクトル
 */
        // 前方宣言: BiCGSTAB法（sparse_solveから委譲）
        template<typename T>
        Vector<T> sparse_bicgstab(const std::vector<T>& values,
            const std::vector<typename Matrix<T>::size_type>& row_ptr,
            const std::vector<typename Matrix<T>::size_type>& col_idx,
            const Vector<T>& b,
            T tolerance,
            typename Matrix<T>::size_type max_iterations);

        template<typename T>
        Vector<T> sparse_solve(const std::vector<T>& values,
            const std::vector<typename Matrix<T>::size_type>& row_ptr,
            const std::vector<typename Matrix<T>::size_type>& col_idx,
            const Vector<T>& b,
            T tolerance = std::numeric_limits<T>::epsilon() * T(100),
            typename Matrix<T>::size_type max_iterations = 0) {

            // 汎用疎行列ソルバー: BiCGSTAB法に委譲（非対称行列にも対応）
            return sparse_bicgstab(values, row_ptr, col_idx, b, tolerance, max_iterations);
        }

        // CSR形式の対称疎行列に対する共役勾配法
        template<typename T>
        Vector<T> sparse_cg(const std::vector<T>& values,
            const std::vector<typename Matrix<T>::size_type>& row_ptr,
            const std::vector<typename Matrix<T>::size_type>& col_idx,
            const Vector<T>& b,
            T tolerance = std::numeric_limits<T>::epsilon() * 10,
            typename Matrix<T>::size_type max_iterations = 0) {
            // 疎行列に対する共役勾配法
            const auto n = b.size();

            // 次元チェック
            if (row_ptr.size() != n + 1) {
                throw DimensionError("sparse_cg: row_ptr size mismatch");
            }

            // デフォルトの最大反復回数
            if (max_iterations == 0) {
                max_iterations = 10 * n;
            }

            Vector<T> x(n, T(0));     // 初期解をゼロに設定

            // 右辺ベクトルのノルム計算
            const T b_norm = std::sqrt(dot(b, b));

            // ゼロベクトルチェック
            if (b_norm < std::numeric_limits<T>::epsilon()) {
                return x;  // 右辺がゼロならゼロベクトルが解
            }

            // 初期残差 r = b - A*x = b（xはゼロ）
            Vector<T> r = b;
            Vector<T> p = r;          // 初期探索方向
            Vector<T> ap(n);          // A*p バッファ (ループ外で1回確保)

            T r_norm_squared = dot(r, r);   // 残差ノルム（二乗）
            // 相対許容誤差（二乗）: ‖r‖² ≤ (tolerance * ‖b‖)²
            const T rel_tol = tolerance * tolerance * b_norm * b_norm;

            // 反復計算
            for (typename Matrix<T>::size_type iter = 0; iter < max_iterations; ++iter) {
                // 収束チェック
                if (r_norm_squared <= rel_tol) {
                    break;
                }

                // A*p を計算（CSR形式, ap はループ外で確保済み）
                ap.assign(n, T(0));
                for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                    for (typename Matrix<T>::size_type j = row_ptr[i]; j < row_ptr[i + 1]; ++j) {
                        ap[i] += values[j] * p[col_idx[j]];
                    }
                }

                // ステップサイズα = r^T*r / (p^T*A*p)
                const T p_ap = dot(p, ap);

                // p^T*A*p がゼロに近い場合（数値的不安定）
                if (std::abs(p_ap) < std::numeric_limits<T>::epsilon()) {
                    break;
                }

                const T alpha = r_norm_squared / p_ap;

                // 解と残差の更新 (融合関数)
                axpy(alpha, p, x);
                axpy(-alpha, ap, r);

                // 新しい残差ノルム（二乗）
                const T r_next_norm_squared = dot(r, r);

                // 方向ベクトルの更新: p = r + beta * p
                if (r_norm_squared < std::numeric_limits<T>::epsilon()) break;  // ゼロ除算回避
                const T beta = r_next_norm_squared / r_norm_squared;
                axpby(T(1), r, beta, p);

                // 残差ノルムの更新
                r_norm_squared = r_next_norm_squared;
            }

            return x;
        }

        // 疎行列に対するBiCGSTAB法
        template<typename T>
        Vector<T> sparse_bicgstab(const std::vector<T>& values,
            const std::vector<typename Matrix<T>::size_type>& row_ptr,
            const std::vector<typename Matrix<T>::size_type>& col_idx,
            const Vector<T>& b,
            T tolerance = std::numeric_limits<T>::epsilon() * 10,
            typename Matrix<T>::size_type max_iterations = 0) {
            // BiCGSTAB法（安定化双共役勾配法）
            const auto n = b.size();

            // 次元チェック
            if (row_ptr.size() != n + 1) {
                throw DimensionError("sparse_bicgstab: row_ptr size mismatch");
            }

            // デフォルトの最大反復回数
            if (max_iterations == 0) {
                max_iterations = 10 * n;
            }

            Vector<T> x(n, T(0));     // 初期解をゼロに設定

            // CSR形式の行列-ベクトル積の関数
            auto csr_mv = [&](const Vector<T>& v, Vector<T>& result) {
                result.assign(n, T(0));
                for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                    for (typename Matrix<T>::size_type j = row_ptr[i]; j < row_ptr[i + 1]; ++j) {
                        result[i] += values[j] * v[col_idx[j]];
                    }
                }
                };

            // 初期残差 r = b - A*x = b（xはゼロ）
            Vector<T> r = b;
            Vector<T> r_hat = r;  // 影響因子

            // 右辺ベクトルのノルム計算
            const T b_norm = std::sqrt(dot(b, b));

            // ゼロベクトルチェック
            if (b_norm < std::numeric_limits<T>::epsilon()) {
                return x;  // 右辺がゼロならゼロベクトルが解
            }

            const T rel_tol = tolerance * b_norm; // 相対許容誤差

            // BiCGSTAB初期化
            T rho_prev = 1;
            T alpha = 1;
            T omega = 1;

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

            T rho = dot(r_hat, r);

            // 反復計算
            for (typename Matrix<T>::size_type iter = 0; iter < max_iterations; ++iter) {
                // 収束チェック
                const T r_norm = std::sqrt(dot(r, r));
                if (r_norm <= rel_tol) {
                    break;
                }

                // ρの計算
                rho = dot(r_hat, r);

                // ρがゼロに近い場合（数値的不安定）
                if (std::abs(rho) < std::numeric_limits<T>::epsilon()) {
                    break;
                }

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

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

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

                // αの計算
                const T r_hat_v = dot(r_hat, v);

                // r_hat_vがゼロに近い場合（数値的不安定）
                if (std::abs(r_hat_v) < std::numeric_limits<T>::epsilon()) {
                    break;
                }

                alpha = rho / r_hat_v;

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

                // 収束チェック
                const T s_norm = std::sqrt(dot(s, s));
                if (s_norm <= rel_tol) {
                    // x += α*p
                    axpy(alpha, p, x);
                    break;
                }

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

                // ωの計算
                const T t_t = dot(t, t);

                // t_tがゼロに近い場合（数値的不安定）
                if (std::abs(t_t) < std::numeric_limits<T>::epsilon()) {
                    break;
                }

                omega = dot(t, s) / t_t;

                // x, rの更新 (融合ループ)
                {
                    auto* xp = x.data();
                    const auto* pp = p.data();
                    const auto* sp = s.data();
                    auto* rp = r.data();
                    const auto* tp = t.data();
                    for (typename Matrix<T>::size_type i = 0; i < n; ++i) {
                        xp[i] += alpha * pp[i] + omega * sp[i];
                        rp[i] = sp[i] - omega * tp[i];
                    }
                }

                // ωがゼロに近い場合（数値的不安定）
                if (std::abs(omega) < std::numeric_limits<T>::epsilon()) {
                    break;
                }

                // ρ_prev の更新
                rho_prev = rho;
            }

            return x;
        }

        // ================================================================
        // 反復改善法 (Iterative Refinement)
        // ================================================================

        /// LU 分解 + 反復改善法で高精度な解を求める
        ///
        /// 参考: 奥村「C言語によるアルゴリズム事典」invr.c (regress.c 内) を元に
        ///       テンプレート化・calx API 統合
        ///
        /// 原理:
        ///   1. LU 分解で初期解 x₀ を求める
        ///   2. 残差 r = b - Ax を高精度に計算
        ///   3. LU 分解を使って d = A⁻¹r を求める (追加の分解不要)
        ///   4. x ← x + d で解を更新
        ///   5. 残差ノルムが十分小さくなるまで繰り返す
        ///
        /// 条件数が大きい行列で有効。double で計算しても
        /// 残差を正確に計算できれば機械精度近くまで改善される。
        ///
        /// @param A 係数行列 (n×n)
        /// @param b 右辺ベクトル
        /// @param maxIter 最大反復回数 (通常 2-5 回で十分)
        /// @param tol 残差ノルムの相対許容誤差
        /// @return 改善された解ベクトル
        template<typename T>
        Vector<T> iterativeRefinement(
            const Matrix<T>& A, const Vector<T>& b,
            size_t maxIter = 10,
            T tol = std::numeric_limits<T>::epsilon() * T(10))
        {
            if (A.rows() != A.cols()) {
                throw DimensionError("iterativeRefinement: matrix must be square");
            }
            if (A.rows() != b.size()) {
                throw DimensionError("iterativeRefinement: matrix and vector dimensions mismatch");
            }

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

            // LU 分解 (1回だけ行い、全反復で再利用)
            auto [lu, pivots] = lu_decomposition(A);

            // 初期解
            Vector<T> x = condest_detail::lu_solve_factored(lu, pivots, b);

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

            for (size_t iter = 0; iter < maxIter; ++iter) {
                // 残差 r = b - Ax を計算
                // ここが精度の鍵: 可能な限り正確に計算する
                Vector<T> r(n);
                for (size_t i = 0; i < n; ++i) {
                    T s = b[i];
                    for (size_t j = 0; j < n; ++j) {
                        s -= A(i, j) * x[j];
                    }
                    r[i] = s;
                }

                // 残差のノルムをチェック
                T rNorm = T(0);
                for (size_t i = 0; i < n; ++i) rNorm += std::abs(r[i]);

                if (rNorm <= tol * bNorm) break;

                // 補正ベクトル d = A⁻¹r (同じ LU 分解を使用)
                Vector<T> d = condest_detail::lu_solve_factored(lu, pivots, r);

                // 解を更新: x ← x + d
                for (size_t i = 0; i < n; ++i) {
                    x[i] += d[i];
                }
            }

            return x;
        }

        /// 反復改善法の結果 (詳細情報付き)
        template<typename T>
        struct RefinementResult {
            Vector<T> solution;      // 改善された解
            T residualNorm;          // 最終残差ノルム
            T conditionEstimate;     // 条件数の推定値
            size_t iterations;       // 実行した反復回数
            bool converged;          // 収束したか
        };

        /// 反復改善法 (詳細情報付き)
        template<typename T>
        RefinementResult<T> iterativeRefinementDetailed(
            const Matrix<T>& A, const Vector<T>& b,
            size_t maxIter = 10,
            T tol = std::numeric_limits<T>::epsilon() * T(10))
        {
            if (A.rows() != A.cols()) {
                throw DimensionError("iterativeRefinementDetailed: matrix must be square");
            }
            if (A.rows() != b.size()) {
                throw DimensionError("iterativeRefinementDetailed: dimensions mismatch");
            }

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

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

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

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

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

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

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

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

            // 条件数推定
            T condEst = estimate_inv_norm1(lu, pivots);
            T norm1A = T(0);
            for (size_t j = 0; j < n; ++j) {
                T colSum = T(0);
                for (size_t i = 0; i < n; ++i)
                    colSum += std::abs(A(i, j));
                if (colSum > norm1A) norm1A = colSum;
            }
            condEst *= norm1A;

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

    } // namespace algorithms
} // namespace calx

#endif // CALX_SOLVERS_HPP