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

/**
 * @file sparse_matrix.hpp
 * @brief 疎行列クラスの実装
 * @author MKL Algebra Library
 */

#ifndef CALX_SPARSE_MATRIX_HPP
#define CALX_SPARSE_MATRIX_HPP

#include <vector>
#include <utility>
#include <algorithm>
#include <stdexcept>
#include <type_traits>
#include <initializer_list>
#include <iosfwd>
#include <iomanip>

#include "../concepts/algebraic_concepts.hpp"
#include "common.hpp"
#include "traits.hpp"
#include "vector.hpp"

namespace calx {

/**
 * @brief 疎行列の格納形式
 */
enum class SparseStorageFormat {
    COO,  ///< 座標形式 (Coordinate format)
    CSR,  ///< 圧縮行格納形式 (Compressed Sparse Row)
    CSC   ///< 圧縮列格納形式 (Compressed Sparse Column)
};

/**
 * @brief 疎行列クラス
 * 
 * 疎行列（ゼロ要素が多い行列）の効率的な表現と操作を提供します。
 * 
 * @tparam T 要素の型（加法アーベル群の要件を満たすこと）
 * @tparam IndexType インデックスの型（デフォルトはsize_t）
 */
template<typename T, typename IndexType = std::size_t>
    requires concepts::AdditiveAbelianGroup<T> && std::integral<IndexType>
class SparseMatrix {
public:
    using value_type = T;
    using index_type = IndexType;
    using size_type = std::size_t;
    using triplet_type = std::tuple<IndexType, IndexType, T>; // (row, col, value)

private:
    // 内部表現: 格納形式によって異なる実装
    SparseStorageFormat storage_format_;
    IndexType rows_;
    IndexType cols_;

    // COO形式データ
    std::vector<triplet_type> coo_data_;

    // CSR形式データ
    std::vector<T> csr_values_;
    std::vector<IndexType> csr_col_indices_;
    std::vector<IndexType> csr_row_ptr_;

    // CSC形式データ
    std::vector<T> csc_values_;
    std::vector<IndexType> csc_row_indices_;
    std::vector<IndexType> csc_col_ptr_;

    // 現在の格納形式が有効かどうか
    bool coo_valid_ = false;
    bool csr_valid_ = false;
    bool csc_valid_ = false;

    // 数値型に対する閾値 (ゼロ判定用)
    static constexpr auto default_epsilon = 
        std::is_floating_point_v<T> ? static_cast<T>(1e-10) : static_cast<T>(0);
    T epsilon_ = default_epsilon;

public:
    /**
     * @brief デフォルトコンストラクタ
     */
    SparseMatrix() 
        : storage_format_(SparseStorageFormat::COO), rows_(0), cols_(0), coo_valid_(true) {}

    /**
     * @brief サイズ指定コンストラクタ
     * 
     * @param rows 行数
     * @param cols 列数
     * @param format 格納形式（デフォルトはCOO）
     */
    SparseMatrix(IndexType rows, IndexType cols, 
                 SparseStorageFormat format = SparseStorageFormat::COO)
        : storage_format_(format), rows_(rows), cols_(cols) 
    {
        switch (storage_format_) {
            case SparseStorageFormat::COO:
                coo_valid_ = true;
                break;
            case SparseStorageFormat::CSR:
                csr_row_ptr_.resize(rows + 1, 0);
                csr_valid_ = true;
                break;
            case SparseStorageFormat::CSC:
                csc_col_ptr_.resize(cols + 1, 0);
                csc_valid_ = true;
                break;
        }
    }

    /**
     * @brief 初期化リストコンストラクタ
     * 
     * @param list 三つ組み (row, col, value) のリスト
     * @param rows 行数
     * @param cols 列数
     * @param format 格納形式（デフォルトはCOO）
     */
    SparseMatrix(std::initializer_list<triplet_type> list, 
                 IndexType rows, IndexType cols,
                 SparseStorageFormat format = SparseStorageFormat::COO)
        : storage_format_(format), rows_(rows), cols_(cols)
    {
        for (const auto& [i, j, val] : list) {
            if (i >= rows || j >= cols) {
                throw std::out_of_range("Index out of bounds");
            }
        }

        switch (storage_format_) {
            case SparseStorageFormat::COO:
                coo_data_ = list;
                coo_valid_ = true;
                break;
            case SparseStorageFormat::CSR:
                convert_to_csr(list);
                csr_valid_ = true;
                break;
            case SparseStorageFormat::CSC:
                convert_to_csc(list);
                csc_valid_ = true;
                break;
        }
    }

    /**
     * @brief コピーコンストラクタ
     */
    SparseMatrix(const SparseMatrix& other)
        : storage_format_(other.storage_format_),
          rows_(other.rows_),
          cols_(other.cols_),
          coo_data_(other.coo_data_),
          csr_values_(other.csr_values_),
          csr_col_indices_(other.csr_col_indices_),
          csr_row_ptr_(other.csr_row_ptr_),
          csc_values_(other.csc_values_),
          csc_row_indices_(other.csc_row_indices_),
          csc_col_ptr_(other.csc_col_ptr_),
          coo_valid_(other.coo_valid_),
          csr_valid_(other.csr_valid_),
          csc_valid_(other.csc_valid_),
          epsilon_(other.epsilon_)
    {}

    /**
     * @brief ムーブコンストラクタ
     */
    SparseMatrix(SparseMatrix&& other) noexcept
        : storage_format_(other.storage_format_),
          rows_(other.rows_),
          cols_(other.cols_),
          coo_data_(std::move(other.coo_data_)),
          csr_values_(std::move(other.csr_values_)),
          csr_col_indices_(std::move(other.csr_col_indices_)),
          csr_row_ptr_(std::move(other.csr_row_ptr_)),
          csc_values_(std::move(other.csc_values_)),
          csc_row_indices_(std::move(other.csc_row_indices_)),
          csc_col_ptr_(std::move(other.csc_col_ptr_)),
          coo_valid_(other.coo_valid_),
          csr_valid_(other.csr_valid_),
          csc_valid_(other.csc_valid_),
          epsilon_(other.epsilon_)
    {
        other.rows_ = 0;
        other.cols_ = 0;
        other.coo_valid_ = false;
        other.csr_valid_ = false;
        other.csc_valid_ = false;
    }

    /**
     * @brief コピー代入演算子
     */
    SparseMatrix& operator=(const SparseMatrix& other) {
        if (this != &other) {
            storage_format_ = other.storage_format_;
            rows_ = other.rows_;
            cols_ = other.cols_;
            coo_data_ = other.coo_data_;
            csr_values_ = other.csr_values_;
            csr_col_indices_ = other.csr_col_indices_;
            csr_row_ptr_ = other.csr_row_ptr_;
            csc_values_ = other.csc_values_;
            csc_row_indices_ = other.csc_row_indices_;
            csc_col_ptr_ = other.csc_col_ptr_;
            coo_valid_ = other.coo_valid_;
            csr_valid_ = other.csr_valid_;
            csc_valid_ = other.csc_valid_;
            epsilon_ = other.epsilon_;
        }
        return *this;
    }

    /**
     * @brief ムーブ代入演算子
     */
    SparseMatrix& operator=(SparseMatrix&& other) noexcept {
        if (this != &other) {
            storage_format_ = other.storage_format_;
            rows_ = other.rows_;
            cols_ = other.cols_;
            coo_data_ = std::move(other.coo_data_);
            csr_values_ = std::move(other.csr_values_);
            csr_col_indices_ = std::move(other.csr_col_indices_);
            csr_row_ptr_ = std::move(other.csr_row_ptr_);
            csc_values_ = std::move(other.csc_values_);
            csc_row_indices_ = std::move(other.csc_row_indices_);
            csc_col_ptr_ = std::move(other.csc_col_ptr_);
            coo_valid_ = other.coo_valid_;
            csr_valid_ = other.csr_valid_;
            csc_valid_ = other.csc_valid_;
            epsilon_ = other.epsilon_;

            other.rows_ = 0;
            other.cols_ = 0;
            other.coo_valid_ = false;
            other.csr_valid_ = false;
            other.csc_valid_ = false;
        }
        return *this;
    }

    /**
     * @brief ゼロ判定用イプシロン値の設定
     * 
     * @param epsilon 新しいイプシロン値
     */
    void set_epsilon(T epsilon) {
        epsilon_ = epsilon;
    }

    /**
     * @brief 現在のイプシロン値を取得
     * 
     * @return 現在のイプシロン値
     */
    T get_epsilon() const {
        return epsilon_;
    }

    /**
     * @brief 行数を取得
     * 
     * @return 行数
     */
    IndexType rows() const {
        return rows_;
    }

    /**
     * @brief 列数を取得
     * 
     * @return 列数
     */
    IndexType cols() const {
        return cols_;
    }

    /**
     * @brief 非ゼロ要素数を取得
     * 
     * @return 非ゼロ要素数
     */
    size_type nnz() const {
        if (coo_valid_) {
            return coo_data_.size();
        } else if (csr_valid_) {
            return csr_values_.size();
        } else if (csc_valid_) {
            return csc_values_.size();
        }
        return 0;
    }

    /**
     * @brief 現在の格納形式を取得
     * 
     * @return 格納形式
     */
    SparseStorageFormat storage_format() const {
        return storage_format_;
    }

    /// CSC 内部配列への読み取りアクセス (SparseLU 等の直接ソルバー用)
    /// ensure_csc_valid() を内部で呼ぶので const でも初回変換が走る
    const std::vector<T>& csc_values() const { ensure_csc_valid(); return csc_values_; }
    const std::vector<IndexType>& csc_row_indices() const { ensure_csc_valid(); return csc_row_indices_; }
    const std::vector<IndexType>& csc_col_ptr() const { ensure_csc_valid(); return csc_col_ptr_; }

    /**
     * @brief 格納形式を変更
     *
     * @param format 新しい格納形式
     */
    void convert_to(SparseStorageFormat format) {
        if (format == storage_format_) {
            return;  // 既に同じ形式
        }

        switch (format) {
            case SparseStorageFormat::COO:
                ensure_coo_valid();
                break;
            case SparseStorageFormat::CSR:
                ensure_csr_valid();
                break;
            case SparseStorageFormat::CSC:
                ensure_csc_valid();
                break;
        }

        storage_format_ = format;
    }

    /**
     * @brief 指定した位置の要素を取得
     * 
     * @param i 行インデックス
     * @param j 列インデックス
     * @return 要素の値
     */
    T coeff(IndexType i, IndexType j) const {
        validate_indices(i, j);

        switch (storage_format_) {
            case SparseStorageFormat::COO:
                return coeff_coo(i, j);
            case SparseStorageFormat::CSR:
                return coeff_csr(i, j);
            case SparseStorageFormat::CSC:
                return coeff_csc(i, j);
            default:
                throw std::runtime_error("Invalid storage format");
        }
    }

    /**
     * @brief 値の設定（既存の値を上書き、または新規追加）
     * 
     * @param i 行インデックス
     * @param j 列インデックス
     * @param value 設定する値
     */
    void set_coeff(IndexType i, IndexType j, const T& value) {
        validate_indices(i, j);

        // ゼロに近い値は格納しない
        if (std::abs(value) <= epsilon_) {
            // 既存の値を削除
            remove_coeff(i, j);
            return;
        }

        switch (storage_format_) {
            case SparseStorageFormat::COO:
                set_coeff_coo(i, j, value);
                break;
            case SparseStorageFormat::CSR:
                set_coeff_csr(i, j, value);
                break;
            case SparseStorageFormat::CSC:
                set_coeff_csc(i, j, value);
                break;
        }
    }

    /**
     * @brief 要素の削除
     * 
     * @param i 行インデックス
     * @param j 列インデックス
     */
    void remove_coeff(IndexType i, IndexType j) {
        validate_indices(i, j);

        switch (storage_format_) {
            case SparseStorageFormat::COO:
                remove_coeff_coo(i, j);
                break;
            case SparseStorageFormat::CSR:
                remove_coeff_csr(i, j);
                break;
            case SparseStorageFormat::CSC:
                remove_coeff_csc(i, j);
                break;
        }

        // 他の形式も無効化
        if (storage_format_ != SparseStorageFormat::COO) {
            coo_valid_ = false;
        }
        if (storage_format_ != SparseStorageFormat::CSR) {
            csr_valid_ = false;
        }
        if (storage_format_ != SparseStorageFormat::CSC) {
            csc_valid_ = false;
        }
    }

    /**
     * @brief 行ベクトルの取得
     * 
     * @param i 行インデックス
     * @return 行ベクトル
     */
    Vector<T> row(IndexType i) const {
        validate_row_index(i);

        Vector<T> result(cols_);

        switch (storage_format_) {
            case SparseStorageFormat::COO: {
                for (const auto& [row, col, val] : coo_data_) {
                    if (row == i) {
                        result[col] = val;
                    }
                }
                break;
            }
            case SparseStorageFormat::CSR: {
                const IndexType start = csr_row_ptr_[i];
                const IndexType end = csr_row_ptr_[i + 1];
                for (IndexType k = start; k < end; ++k) {
                    result[csr_col_indices_[k]] = csr_values_[k];
                }
                break;
            }
            case SparseStorageFormat::CSC: {
                for (IndexType j = 0; j < cols_; ++j) {
                    const IndexType start = csc_col_ptr_[j];
                    const IndexType end = csc_col_ptr_[j + 1];
                    for (IndexType k = start; k < end; ++k) {
                        if (csc_row_indices_[k] == i) {
                            result[j] = csc_values_[k];
                            break;
                        }
                    }
                }
                break;
            }
        }

        return result;
    }

    /**
     * @brief 列ベクトルの取得
     * 
     * @param j 列インデックス
     * @return 列ベクトル
     */
    Vector<T> col(IndexType j) const {
        validate_col_index(j);

        Vector<T> result(rows_);

        switch (storage_format_) {
            case SparseStorageFormat::COO: {
                for (const auto& [row, col, val] : coo_data_) {
                    if (col == j) {
                        result[row] = val;
                    }
                }
                break;
            }
            case SparseStorageFormat::CSR: {
                for (IndexType i = 0; i < rows_; ++i) {
                    const IndexType start = csr_row_ptr_[i];
                    const IndexType end = csr_row_ptr_[i + 1];
                    for (IndexType k = start; k < end; ++k) {
                        if (csr_col_indices_[k] == j) {
                            result[i] = csr_values_[k];
                            break;
                        }
                    }
                }
                break;
            }
            case SparseStorageFormat::CSC: {
                const IndexType start = csc_col_ptr_[j];
                const IndexType end = csc_col_ptr_[j + 1];
                for (IndexType k = start; k < end; ++k) {
                    result[csc_row_indices_[k]] = csc_values_[k];
                }
                break;
            }
        }

        return result;
    }

    /**
     * @brief 行列の転置
     * 
     * @return 転置行列
     */
    SparseMatrix transpose() const {
        SparseMatrix result(cols_, rows_, storage_format_);

        switch (storage_format_) {
            case SparseStorageFormat::COO: {
                result.coo_data_.reserve(coo_data_.size());
                for (const auto& [i, j, val] : coo_data_) {
                    result.coo_data_.emplace_back(j, i, val);
                }
                result.coo_valid_ = true;
                break;
            }
            case SparseStorageFormat::CSR: {
                // CSRの転置はCSCとして生成し、CSRに変換
                result.csc_values_ = csr_values_;
                result.csc_row_indices_ = csr_col_indices_;
                result.csc_col_ptr_ = std::vector<IndexType>(cols_ + 1, 0);

                // CSCのcol_ptrを計算（転置後の行ポインタに相当）
                for (IndexType i = 0; i < csr_values_.size(); ++i) {
                    result.csc_col_ptr_[csr_col_indices_[i] + 1]++;
                }
                
                for (IndexType i = 0; i < cols_; ++i) {
                    result.csc_col_ptr_[i + 1] += result.csc_col_ptr_[i];
                }

                result.csc_valid_ = true;
                result.convert_to(SparseStorageFormat::CSR);
                break;
            }
            case SparseStorageFormat::CSC: {
                // CSCの転置はCSRとして生成し、CSCに変換
                result.csr_values_ = csc_values_;
                result.csr_col_indices_ = csc_row_indices_;
                result.csr_row_ptr_ = std::vector<IndexType>(rows_ + 1, 0);

                // CSRのrow_ptrを計算（転置後の列ポインタに相当）
                for (IndexType i = 0; i < csc_values_.size(); ++i) {
                    result.csr_row_ptr_[csc_row_indices_[i] + 1]++;
                }
                
                for (IndexType i = 0; i < rows_; ++i) {
                    result.csr_row_ptr_[i + 1] += result.csr_row_ptr_[i];
                }

                result.csr_valid_ = true;
                result.convert_to(SparseStorageFormat::CSC);
                break;
            }
        }

        return result;
    }

    /**
     * @brief 行列の要素ごとのスカラー倍
     * 
     * @param scalar スカラー値
     * @return スカラー倍された行列
     */
    template<typename ScalarType>
        requires concepts::Scalar<ScalarType>
    SparseMatrix operator*(const ScalarType& scalar) const {
        SparseMatrix result(*this);

        switch (storage_format_) {
            case SparseStorageFormat::COO: {
                for (auto& [i, j, val] : result.coo_data_) {
                    val *= scalar;
                }
                break;
            }
            case SparseStorageFormat::CSR: {
                for (auto& val : result.csr_values_) {
                    val *= scalar;
                }
                break;
            }
            case SparseStorageFormat::CSC: {
                for (auto& val : result.csc_values_) {
                    val *= scalar;
                }
                break;
            }
        }

        return result;
    }

    /**
     * @brief 行列のスカラー除算
     * 
     * @param scalar スカラー値
     * @return スカラー除算された行列
     */
    template<typename ScalarType>
        requires concepts::Scalar<ScalarType>
    SparseMatrix operator/(const ScalarType& scalar) const {
        if (scalar == ScalarType{0}) {
            throw std::invalid_argument("Division by zero");
        }

        SparseMatrix result(*this);

        switch (storage_format_) {
            case SparseStorageFormat::COO: {
                for (auto& [i, j, val] : result.coo_data_) {
                    val /= scalar;
                }
                break;
            }
            case SparseStorageFormat::CSR: {
                for (auto& val : result.csr_values_) {
                    val /= scalar;
                }
                break;
            }
            case SparseStorageFormat::CSC: {
                for (auto& val : result.csc_values_) {
                    val /= scalar;
                }
                break;
            }
        }

        return result;
    }

    /**
     * @brief 行列の加算
     * 
     * @param other 加算する行列
     * @return 加算結果
     */
    SparseMatrix operator+(const SparseMatrix& other) const {
        if (rows_ != other.rows_ || cols_ != other.cols_) {
            throw std::invalid_argument("Matrix dimensions must match");
        }

        // 加算はCOO形式で行うのが最も単純
        ensure_coo_valid();
        other.ensure_coo_valid();

        SparseMatrix result(rows_, cols_, SparseStorageFormat::COO);

        // 既存のデータをコピー
        result.coo_data_ = coo_data_;

        // 他方の行列の要素を追加または加算
        for (const auto& [i, j, val] : other.coo_data_) {
            bool found = false;
            for (auto& [ri, rj, rval] : result.coo_data_) {
                if (ri == i && rj == j) {
                    rval += val;
                    found = true;
                    break;
                }
            }
            if (!found) {
                result.coo_data_.emplace_back(i, j, val);
            }
        }

        // ゼロになった要素を除去
        result.coo_data_.erase(
            std::remove_if(result.coo_data_.begin(), result.coo_data_.end(),
                [this](const auto& t) { 
                    return std::abs(std::get<2>(t)) <= epsilon_; 
                }),
            result.coo_data_.end()
        );

        result.coo_valid_ = true;
        
        // 必要に応じて他の形式に変換
        if (storage_format_ != SparseStorageFormat::COO) {
            result.convert_to(storage_format_);
        }

        return result;
    }

    /**
     * @brief 行列の減算
     * 
     * @param other 減算する行列
     * @return 減算結果
     */
    SparseMatrix operator-(const SparseMatrix& other) const {
        if (rows_ != other.rows_ || cols_ != other.cols_) {
            throw std::invalid_argument("Matrix dimensions must match");
        }

        // 減算はCOO形式で行うのが最も単純
        ensure_coo_valid();
        other.ensure_coo_valid();

        SparseMatrix result(rows_, cols_, SparseStorageFormat::COO);

        // 既存のデータをコピー
        result.coo_data_ = coo_data_;

        // 他方の行列の要素を減算
        for (const auto& [i, j, val] : other.coo_data_) {
            bool found = false;
            for (auto& [ri, rj, rval] : result.coo_data_) {
                if (ri == i && rj == j) {
                    rval -= val;
                    found = true;
                    break;
                }
            }
            if (!found) {
                result.coo_data_.emplace_back(i, j, -val);
            }
        }

        // ゼロになった要素を除去
        result.coo_data_.erase(
            std::remove_if(result.coo_data_.begin(), result.coo_data_.end(),
                [this](const auto& t) { 
                    return std::abs(std::get<2>(t)) <= epsilon_; 
                }),
            result.coo_data_.end()
        );

        result.coo_valid_ = true;
        
        // 必要に応じて他の形式に変換
        if (storage_format_ != SparseStorageFormat::COO) {
            result.convert_to(storage_format_);
        }

        return result;
    }

    /**
     * @brief 行列とベクトルの積
     * 
     * @param vec 右側のベクトル
     * @return 結果ベクトル
     */
    Vector<T> operator*(const Vector<T>& vec) const {
        if (cols_ != vec.size()) {
            throw std::invalid_argument("Matrix and vector dimensions must match");
        }

        Vector<T> result(rows_, T{0});

        switch (storage_format_) {
            case SparseStorageFormat::COO: {
                for (const auto& [i, j, val] : coo_data_) {
                    result[i] += val * vec[j];
                }
                break;
            }
            case SparseStorageFormat::CSR: {
                for (IndexType i = 0; i < rows_; ++i) {
                    const IndexType start = csr_row_ptr_[i];
                    const IndexType end = csr_row_ptr_[i + 1];
                    for (IndexType k = start; k < end; ++k) {
                        result[i] += csr_values_[k] * vec[csr_col_indices_[k]];
                    }
                }
                break;
            }
            case SparseStorageFormat::CSC: {
                for (IndexType j = 0; j < cols_; ++j) {
                    const T vec_j = vec[j];
                    const IndexType start = csc_col_ptr_[j];
                    const IndexType end = csc_col_ptr_[j + 1];
                    for (IndexType k = start; k < end; ++k) {
                        result[csc_row_indices_[k]] += csc_values_[k] * vec_j;
                    }
                }
                break;
            }
        }

        return result;
    }

    /// SpMM: 疎行列 × 密行列 → 密行列 (GSL-7)
    Matrix<T> multiply(const Matrix<T>& dense) const {
        if (cols_ != dense.rows()) {
            throw std::invalid_argument("SpMM: dimension mismatch");
        }
        Matrix<T> result(rows_, dense.cols(), T{0});
        switch (storage_format_) {
            case SparseStorageFormat::COO:
                for (const auto& [i, j, val] : coo_data_)
                    for (IndexType c = 0; c < dense.cols(); ++c)
                        result(i, c) += val * dense(j, c);
                break;
            case SparseStorageFormat::CSR:
                for (IndexType i = 0; i < rows_; ++i) {
                    const IndexType start = csr_row_ptr_[i];
                    const IndexType end = csr_row_ptr_[i + 1];
                    for (IndexType k = start; k < end; ++k)
                        for (IndexType c = 0; c < dense.cols(); ++c)
                            result(i, c) += csr_values_[k] * dense(csr_col_indices_[k], c);
                }
                break;
            case SparseStorageFormat::CSC:
                for (IndexType j = 0; j < cols_; ++j) {
                    const IndexType start = csc_col_ptr_[j];
                    const IndexType end = csc_col_ptr_[j + 1];
                    for (IndexType k = start; k < end; ++k)
                        for (IndexType c = 0; c < dense.cols(); ++c)
                            result(csc_row_indices_[k], c) += csc_values_[k] * dense(j, c);
                }
                break;
        }
        return result;
    }

    /**
     * @brief 行列同士の積
     * 
     * @param other 右側の行列
     * @return 結果行列
     */
    SparseMatrix operator*(const SparseMatrix& other) const {
        if (cols_ != other.rows_) {
            throw std::invalid_argument("Matrix dimensions must match for multiplication");
        }

        SparseMatrix result(rows_, other.cols_, SparseStorageFormat::COO);

        // 最適な形式に変換
        ensure_csr_valid();  // 左側はCSR
        other.ensure_csc_valid();  // 右側はCSC

        for (IndexType i = 0; i < rows_; ++i) {
            const IndexType row_start = csr_row_ptr_[i];
            const IndexType row_end = csr_row_ptr_[i + 1];

            for (IndexType j = 0; j < other.cols_; ++j) {
                const IndexType col_start = other.csc_col_ptr_[j];
                const IndexType col_end = other.csc_col_ptr_[j + 1];

                T sum{0};
                IndexType k1 = row_start;
                IndexType k2 = col_start;

                // 疎行列同士の効率的な乗算
                while (k1 < row_end && k2 < col_end) {
                    IndexType col1 = csr_col_indices_[k1];
                    IndexType row2 = other.csc_row_indices_[k2];

                    if (col1 < row2) {
                        ++k1;
                    } else if (col1 > row2) {
                        ++k2;
                    } else {
                        sum += csr_values_[k1] * other.csc_values_[k2];
                        ++k1;
                        ++k2;
                    }
                }

                if (std::abs(sum) > epsilon_) {
                    result.coo_data_.emplace_back(i, j, sum);
                }
            }
        }

        result.coo_valid_ = true;

        // 必要に応じて他の形式に変換
        if (storage_format_ != SparseStorageFormat::COO) {
            result.convert_to(storage_format_);
        }

        return result;
    }

    /**
     * @brief 行列の要素アクセス演算子
     * 
     * @param i 行インデックス
     * @param j 列インデックス
     * @return 要素への参照（プロキシ）
     */
    class ElementProxy {
    private:
        SparseMatrix<T, IndexType>& matrix_;
        IndexType i_;
        IndexType j_;

    public:
        ElementProxy(SparseMatrix<T, IndexType>& matrix, IndexType i, IndexType j)
            : matrix_(matrix), i_(i), j_(j) {}

        // 暗黙の変換演算子
        operator T() const {
            return matrix_.coeff(i_, j_);
        }

        // 代入演算子
        ElementProxy& operator=(const T& value) {
            matrix_.set_coeff(i_, j_, value);
            return *this;
        }

        // 複合代入演算子
        ElementProxy& operator+=(const T& value) {
            matrix_.set_coeff(i_, j_, matrix_.coeff(i_, j_) + value);
            return *this;
        }

        ElementProxy& operator-=(const T& value) {
            matrix_.set_coeff(i_, j_, matrix_.coeff(i_, j_) - value);
            return *this;
        }

        ElementProxy& operator*=(const T& value) {
            matrix_.set_coeff(i_, j_, matrix_.coeff(i_, j_) * value);
            return *this;
        }

        ElementProxy& operator/=(const T& value) {
            matrix_.set_coeff(i_, j_, matrix_.coeff(i_, j_) / value);
            return *this;
        }
    };

    /**
     * @brief 行列の要素アクセス演算子（非const版）
     * 
     * @param i 行インデックス
     * @param j 列インデックス
     * @return 要素プロキシ
     */
    ElementProxy operator()(IndexType i, IndexType j) {
        validate_indices(i, j);
        return ElementProxy(*this, i, j);
    }

    /**
     * @brief 行列の要素アクセス演算子（const版）
     * 
     * @param i 行インデックス
     * @param j 列インデックス
     * @return 要素の値
     */
    T operator()(IndexType i, IndexType j) const {
        return coeff(i, j);
    }

    /**
     * @brief 行列のブロック（部分行列）を取得
     * 
     * @param start_row 開始行インデックス
     * @param start_col 開始列インデックス
     * @param block_rows ブロックの行数
     * @param block_cols ブロックの列数
     * @return 部分行列
     */
    SparseMatrix block(IndexType start_row, IndexType start_col, 
                       IndexType block_rows, IndexType block_cols) const {
        if (start_row + block_rows > rows_ || start_col + block_cols > cols_) {
            throw std::out_of_range("Block indices out of bounds");
        }

        ensure_coo_valid();

        SparseMatrix result(block_rows, block_cols, SparseStorageFormat::COO);

        for (const auto& [i, j, val] : coo_data_) {
            if (i >= start_row && i < start_row + block_rows &&
                j >= start_col && j < start_col + block_cols) {
                result.coo_data_.emplace_back(i - start_row, j - start_col, val);
            }
        }

        result.coo_valid_ = true;

        // 必要に応じて他の形式に変換
        if (storage_format_ != SparseStorageFormat::COO) {
            result.convert_to(storage_format_);
        }

        return result;
    }

    /**
     * @brief 対角行列の生成
     * 
     * @param diag 対角成分のベクトル
     * @return 対角行列
     */
    static SparseMatrix diagonal(const Vector<T>& diag) {
        const IndexType n = diag.size();
        SparseMatrix result(n, n, SparseStorageFormat::COO);

        for (IndexType i = 0; i < n; ++i) {
            if (std::abs(diag[i]) > result.epsilon_) {
                result.coo_data_.emplace_back(i, i, diag[i]);
            }
        }

        result.coo_valid_ = true;
        return result;
    }

    /**
     * @brief 単位行列の生成
     * 
     * @param n 行列サイズ
     * @return 単位行列
     */
    static SparseMatrix identity(IndexType n) {
        SparseMatrix result(n, n, SparseStorageFormat::COO);
        result.coo_data_.reserve(n);

        for (IndexType i = 0; i < n; ++i) {
            result.coo_data_.emplace_back(i, i, T{1});
        }

        result.coo_valid_ = true;
        return result;
    }

    /**
     * @brief 三重対角行列の生成
     * 
     * @param diag 対角成分のベクトル
     * @param lower 下対角成分のベクトル
     * @param upper 上対角成分のベクトル
     * @return 三重対角行列
     */
    static SparseMatrix tridiagonal(const Vector<T>& diag, 
                                 const Vector<T>& lower,
                                 const Vector<T>& upper) {
        const IndexType n = diag.size();
        
        if (lower.size() != n - 1 || upper.size() != n - 1) {
            throw std::invalid_argument("Diagonal and off-diagonal sizes must match");
        }
        
        SparseMatrix result(n, n, SparseStorageFormat::COO);
        
        // 対角成分
        for (IndexType i = 0; i < n; ++i) {
            if (std::abs(diag[i]) > result.epsilon_) {
                result.coo_data_.emplace_back(i, i, diag[i]);
            }
        }
        
        // 下対角成分
        for (IndexType i = 0; i < n - 1; ++i) {
            if (std::abs(lower[i]) > result.epsilon_) {
                result.coo_data_.emplace_back(i + 1, i, lower[i]);
            }
        }
        
        // 上対角成分
        for (IndexType i = 0; i < n - 1; ++i) {
            if (std::abs(upper[i]) > result.epsilon_) {
                result.coo_data_.emplace_back(i, i + 1, upper[i]);
            }
        }
        
        result.coo_valid_ = true;
        return result;
    }

    /**
     * @brief ゼロ行列の生成
     * 
     * @param rows 行数
     * @param cols 列数
     * @return ゼロ行列
     */
    static SparseMatrix zero(IndexType rows, IndexType cols) {
        return SparseMatrix(rows, cols);
    }

private:
    /**
     * @brief インデックスの妥当性を検証
     * 
     * @param i 行インデックス
     * @param j 列インデックス
     */
    void validate_indices(IndexType i, IndexType j) const {
        if (i >= rows_ || j >= cols_) {
            throw std::out_of_range("Matrix indices out of bounds");
        }
    }

    /**
     * @brief 行インデックスの妥当性を検証
     * 
     * @param i 行インデックス
     */
    void validate_row_index(IndexType i) const {
        if (i >= rows_) {
            throw std::out_of_range("Row index out of bounds");
        }
    }

    /**
     * @brief 列インデックスの妥当性を検証
     * 
     * @param j 列インデックス
     */
    void validate_col_index(IndexType j) const {
        if (j >= cols_) {
            throw std::out_of_range("Column index out of bounds");
        }
    }

    /**
     * @brief COO形式のデータが有効であることを確認
     */
    void ensure_coo_valid() const {
        if (coo_valid_) {
            return;
        }

        auto& self = const_cast<SparseMatrix&>(*this);

        if (csr_valid_) {
            self.convert_csr_to_coo();
        } else if (csc_valid_) {
            self.convert_csc_to_coo();
        } else {
            throw std::runtime_error("No valid storage format");
        }
    }

    /**
     * @brief CSR形式のデータが有効であることを確認
     */
    void ensure_csr_valid() const {
        if (csr_valid_) {
            return;
        }

        auto& self = const_cast<SparseMatrix&>(*this);

        if (coo_valid_) {
            self.convert_coo_to_csr();
        } else if (csc_valid_) {
            self.convert_csc_to_csr();
        } else {
            throw std::runtime_error("No valid storage format");
        }
    }

    /**
     * @brief CSC形式のデータが有効であることを確認
     */
    void ensure_csc_valid() const {
        if (csc_valid_) {
            return;
        }

        auto& self = const_cast<SparseMatrix&>(*this);

        if (coo_valid_) {
            self.convert_coo_to_csc();
        } else if (csr_valid_) {
            self.convert_csr_to_csc();
        } else {
            throw std::runtime_error("No valid storage format");
        }
    }

    /**
     * @brief COO→CSR変換
     */
    void convert_coo_to_csr() {
        if (!coo_valid_ || coo_data_.empty()) {
            // 空の行列の場合
            csr_values_.clear();
            csr_col_indices_.clear();
            csr_row_ptr_.resize(rows_ + 1, 0);
            csr_valid_ = true;
            return;
        }

        // COOデータをソート（行優先）
        std::sort(coo_data_.begin(), coo_data_.end(),
            [](const triplet_type& a, const triplet_type& b) {
                return std::tie(std::get<0>(a), std::get<1>(a)) < 
                       std::tie(std::get<0>(b), std::get<1>(b));
            });

        // CSRデータ構造を初期化
        const size_type nnz = coo_data_.size();
        csr_values_.resize(nnz);
        csr_col_indices_.resize(nnz);
        csr_row_ptr_.resize(rows_ + 1, 0);

        // 各行の非ゼロ要素数をカウント
        for (const auto& [i, j, val] : coo_data_) {
            csr_row_ptr_[i + 1]++;
        }

        // 累積和を計算して行ポインタを得る
        for (IndexType i = 0; i < rows_; ++i) {
            csr_row_ptr_[i + 1] += csr_row_ptr_[i];
        }

        // 値と列インデックスを設定
        for (const auto& [i, j, val] : coo_data_) {
            const IndexType pos = csr_row_ptr_[i]++;
            csr_values_[pos] = val;
            csr_col_indices_[pos] = j;
        }

        // 行ポインタを修正（累積和の計算中に変更されたため）
        for (IndexType i = rows_; i > 0; --i) {
            csr_row_ptr_[i] = csr_row_ptr_[i - 1];
        }
        csr_row_ptr_[0] = 0;

        csr_valid_ = true;
    }

    /**
     * @brief COO→CSC変換
     */
    void convert_coo_to_csc() {
        if (!coo_valid_ || coo_data_.empty()) {
            // 空の行列の場合
            csc_values_.clear();
            csc_row_indices_.clear();
            csc_col_ptr_.resize(cols_ + 1, 0);
            csc_valid_ = true;
            return;
        }

        // COOデータをソート（列優先）
        std::sort(coo_data_.begin(), coo_data_.end(),
            [](const triplet_type& a, const triplet_type& b) {
                return std::tie(std::get<1>(a), std::get<0>(a)) < 
                       std::tie(std::get<1>(b), std::get<0>(b));
            });

        // CSCデータ構造を初期化
        const size_type nnz = coo_data_.size();
        csc_values_.resize(nnz);
        csc_row_indices_.resize(nnz);
        csc_col_ptr_.resize(cols_ + 1, 0);

        // 各列の非ゼロ要素数をカウント
        for (const auto& [i, j, val] : coo_data_) {
            csc_col_ptr_[j + 1]++;
        }

        // 累積和を計算して列ポインタを得る
        for (IndexType j = 0; j < cols_; ++j) {
            csc_col_ptr_[j + 1] += csc_col_ptr_[j];
        }

        // 値と行インデックスを設定
        for (const auto& [i, j, val] : coo_data_) {
            const IndexType pos = csc_col_ptr_[j]++;
            csc_values_[pos] = val;
            csc_row_indices_[pos] = i;
        }

        // 列ポインタを修正（累積和の計算中に変更されたため）
        for (IndexType j = cols_; j > 0; --j) {
            csc_col_ptr_[j] = csc_col_ptr_[j - 1];
        }
        csc_col_ptr_[0] = 0;

        csc_valid_ = true;
    }

    /**
     * @brief CSR→COO変換
     */
    void convert_csr_to_coo() {
        if (!csr_valid_) {
            throw std::runtime_error("CSR data is not valid");
        }

        const size_type nnz = csr_values_.size();
        coo_data_.resize(nnz);

        size_type idx = 0;
        for (IndexType i = 0; i < rows_; ++i) {
            for (IndexType j = csr_row_ptr_[i]; j < csr_row_ptr_[i + 1]; ++j) {
                coo_data_[idx++] = std::make_tuple(i, csr_col_indices_[j], csr_values_[j]);
            }
        }

        coo_valid_ = true;
    }

    /**
     * @brief CSC→COO変換
     */
    void convert_csc_to_coo() {
        if (!csc_valid_) {
            throw std::runtime_error("CSC data is not valid");
        }

        const size_type nnz = csc_values_.size();
        coo_data_.resize(nnz);

        size_type idx = 0;
        for (IndexType j = 0; j < cols_; ++j) {
            for (IndexType i = csc_col_ptr_[j]; i < csc_col_ptr_[j + 1]; ++i) {
                coo_data_[idx++] = std::make_tuple(csc_row_indices_[i], j, csc_values_[i]);
            }
        }

        coo_valid_ = true;
    }

    /**
     * @brief CSR→CSC変換
     */
    void convert_csr_to_csc() {
        if (!csr_valid_) {
            throw std::runtime_error("CSR data is not valid");
        }

        // 一旦COOに変換して、そこからCSCに変換する
        convert_csr_to_coo();
        convert_coo_to_csc();
    }

    /**
     * @brief CSC→CSR変換
     */
    void convert_csc_to_csr() {
        if (!csc_valid_) {
            throw std::runtime_error("CSC data is not valid");
        }

        // 一旦COOに変換して、そこからCSRに変換する
        convert_csc_to_coo();
        convert_coo_to_csr();
    }

    /**
     * @brief 初期化リストからCSR形式に変換
     * 
     * @param list 三つ組みのリスト
     */
    void convert_to_csr(const std::initializer_list<triplet_type>& list) {
        std::vector<triplet_type> sorted_list(list);
        
        // 行優先でソート
        std::sort(sorted_list.begin(), sorted_list.end(),
            [](const triplet_type& a, const triplet_type& b) {
                return std::tie(std::get<0>(a), std::get<1>(a)) < 
                       std::tie(std::get<0>(b), std::get<1>(b));
            });

        // CSRデータ構造を初期化
        const size_type nnz = sorted_list.size();
        csr_values_.resize(nnz);
        csr_col_indices_.resize(nnz);
        csr_row_ptr_.resize(rows_ + 1, 0);

        // 各行の非ゼロ要素数をカウント
        for (const auto& [i, j, val] : sorted_list) {
            csr_row_ptr_[i + 1]++;
        }

        // 累積和を計算して行ポインタを得る
        for (IndexType i = 0; i < rows_; ++i) {
            csr_row_ptr_[i + 1] += csr_row_ptr_[i];
        }

        // 値と列インデックスを設定
        for (const auto& [i, j, val] : sorted_list) {
            const IndexType pos = csr_row_ptr_[i]++;
            csr_values_[pos] = val;
            csr_col_indices_[pos] = j;
        }

        // 行ポインタを修正
        for (IndexType i = rows_; i > 0; --i) {
            csr_row_ptr_[i] = csr_row_ptr_[i - 1];
        }
        csr_row_ptr_[0] = 0;
    }

    /**
     * @brief 初期化リストからCSC形式に変換
     * 
     * @param list 三つ組みのリスト
     */
    void convert_to_csc(const std::initializer_list<triplet_type>& list) {
        std::vector<triplet_type> sorted_list(list);
        
        // 列優先でソート
        std::sort(sorted_list.begin(), sorted_list.end(),
            [](const triplet_type& a, const triplet_type& b) {
                return std::tie(std::get<1>(a), std::get<0>(a)) < 
                       std::tie(std::get<1>(b), std::get<0>(b));
            });

        // CSCデータ構造を初期化
        const size_type nnz = sorted_list.size();
        csc_values_.resize(nnz);
        csc_row_indices_.resize(nnz);
        csc_col_ptr_.resize(cols_ + 1, 0);

        // 各列の非ゼロ要素数をカウント
        for (const auto& [i, j, val] : sorted_list) {
            csc_col_ptr_[j + 1]++;
        }

        // 累積和を計算して列ポインタを得る
        for (IndexType j = 0; j < cols_; ++j) {
            csc_col_ptr_[j + 1] += csc_col_ptr_[j];
        }

        // 値と行インデックスを設定
        for (const auto& [i, j, val] : sorted_list) {
            const IndexType pos = csc_col_ptr_[j]++;
            csc_values_[pos] = val;
            csc_row_indices_[pos] = i;
        }

        // 列ポインタを修正
        for (IndexType j = cols_; j > 0; --j) {
            csc_col_ptr_[j] = csc_col_ptr_[j - 1];
        }
        csc_col_ptr_[0] = 0;
    }

    /**
     * @brief COO形式での要素取得
     */
    T coeff_coo(IndexType i, IndexType j) const {
        for (const auto& [row, col, val] : coo_data_) {
            if (row == i && col == j) {
                return val;
            }
        }
        return T{0};
    }

    /**
     * @brief CSR形式での要素取得
     */
    T coeff_csr(IndexType i, IndexType j) const {
        const IndexType start = csr_row_ptr_[i];
        const IndexType end = csr_row_ptr_[i + 1];
        
        for (IndexType k = start; k < end; ++k) {
            if (csr_col_indices_[k] == j) {
                return csr_values_[k];
            }
            // 列インデックスはソートされているという前提
            if (csr_col_indices_[k] > j) {
                break;
            }
        }
        return T{0};
    }

    /**
     * @brief CSC形式での要素取得
     */
    T coeff_csc(IndexType i, IndexType j) const {
        const IndexType start = csc_col_ptr_[j];
        const IndexType end = csc_col_ptr_[j + 1];
        
        for (IndexType k = start; k < end; ++k) {
            if (csc_row_indices_[k] == i) {
                return csc_values_[k];
            }
            // 行インデックスはソートされているという前提
            if (csc_row_indices_[k] > i) {
                break;
            }
        }
        return T{0};
    }

    /**
     * @brief COO形式での値設定
     */
    void set_coeff_coo(IndexType i, IndexType j, const T& value) {
        for (auto& [row, col, val] : coo_data_) {
            if (row == i && col == j) {
                val = value;
                return;
            }
        }
        // 既存の要素が見つからない場合は追加
        coo_data_.emplace_back(i, j, value);
    }

    /**
     * @brief CSR形式での値設定
     */
    void set_coeff_csr(IndexType i, IndexType j, const T& value) {
        const IndexType start = csr_row_ptr_[i];
        const IndexType end = csr_row_ptr_[i + 1];
        
        for (IndexType k = start; k < end; ++k) {
            if (csr_col_indices_[k] == j) {
                // 既存の要素を更新
                csr_values_[k] = value;
                return;
            }
            if (csr_col_indices_[k] > j) {
                // 新しい要素を挿入
                csr_values_.insert(csr_values_.begin() + k, value);
                csr_col_indices_.insert(csr_col_indices_.begin() + k, j);
                // 行ポインタを更新
                for (IndexType m = i + 1; m <= rows_; ++m) {
                    ++csr_row_ptr_[m];
                }
                return;
            }
        }
        
        // 行の最後に追加
        csr_values_.insert(csr_values_.begin() + end, value);
        csr_col_indices_.insert(csr_col_indices_.begin() + end, j);
        // 行ポインタを更新
        for (IndexType m = i + 1; m <= rows_; ++m) {
            ++csr_row_ptr_[m];
        }
    }

    /**
     * @brief CSC形式での値設定
     */
    void set_coeff_csc(IndexType i, IndexType j, const T& value) {
        const IndexType start = csc_col_ptr_[j];
        const IndexType end = csc_col_ptr_[j + 1];
        
        for (IndexType k = start; k < end; ++k) {
            if (csc_row_indices_[k] == i) {
                // 既存の要素を更新
                csc_values_[k] = value;
                return;
            }
            if (csc_row_indices_[k] > i) {
                // 新しい要素を挿入
                csc_values_.insert(csc_values_.begin() + k, value);
                csc_row_indices_.insert(csc_row_indices_.begin() + k, i);
                // 列ポインタを更新
                for (IndexType m = j + 1; m <= cols_; ++m) {
                    ++csc_col_ptr_[m];
                }
                return;
            }
        }
        
        // 列の最後に追加
        csc_values_.insert(csc_values_.begin() + end, value);
        csc_row_indices_.insert(csc_row_indices_.begin() + end, i);
        // 列ポインタを更新
        for (IndexType m = j + 1; m <= cols_; ++m) {
            ++csc_col_ptr_[m];
        }
    }

    /**
     * @brief COO形式での要素削除
     */
    void remove_coeff_coo(IndexType i, IndexType j) {
        auto it = std::find_if(coo_data_.begin(), coo_data_.end(),
            [i, j](const triplet_type& t) {
                return std::get<0>(t) == i && std::get<1>(t) == j;
            });
        
        if (it != coo_data_.end()) {
            coo_data_.erase(it);
        }
    }

    /**
     * @brief CSR形式での要素削除
     */
    void remove_coeff_csr(IndexType i, IndexType j) {
        const IndexType start = csr_row_ptr_[i];
        const IndexType end = csr_row_ptr_[i + 1];
        
        for (IndexType k = start; k < end; ++k) {
            if (csr_col_indices_[k] == j) {
                // 要素を削除
                csr_values_.erase(csr_values_.begin() + k);
                csr_col_indices_.erase(csr_col_indices_.begin() + k);
                // 行ポインタを更新
                for (IndexType m = i + 1; m <= rows_; ++m) {
                    --csr_row_ptr_[m];
                }
                return;
            }
        }
    }

    /**
     * @brief CSC形式での要素削除
     */
    void remove_coeff_csc(IndexType i, IndexType j) {
        const IndexType start = csc_col_ptr_[j];
        const IndexType end = csc_col_ptr_[j + 1];
        
        for (IndexType k = start; k < end; ++k) {
            if (csc_row_indices_[k] == i) {
                // 要素を削除
                csc_values_.erase(csc_values_.begin() + k);
                csc_row_indices_.erase(csc_row_indices_.begin() + k);
                // 列ポインタを更新
                for (IndexType m = j + 1; m <= cols_; ++m) {
                    --csc_col_ptr_[m];
                }
                return;
            }
        }
    }
};

/**
 * @brief フリースタンディングの転置関数
 * 
 * @tparam T 要素の型
 * @tparam IndexType インデックスの型
 * @param matrix 転置する行列
 * @return 転置行列
 */
template<typename T, typename IndexType>
    requires concepts::AdditiveAbelianGroup<T> && std::integral<IndexType>
SparseMatrix<T, IndexType> transpose(const SparseMatrix<T, IndexType>& matrix) {
    return matrix.transpose();
}

/**
 * @brief フリースタンディングのスカラー乗算
 * 
 * @tparam T 要素の型
 * @tparam ScalarType スカラーの型
 * @tparam IndexType インデックスの型
 * @param scalar スカラー値
 * @param matrix 行列
 * @return スカラー倍された行列
 */
template<typename T, typename ScalarType, typename IndexType>
    requires concepts::AdditiveAbelianGroup<T> && concepts::Scalar<ScalarType> && std::integral<IndexType>
SparseMatrix<T, IndexType> operator*(const ScalarType& scalar, const SparseMatrix<T, IndexType>& matrix) {
    return matrix * scalar;
}

/**
 * @brief 出力ストリーム演算子のオーバーロード
 * 
 * @tparam T 要素の型
 * @tparam IndexType インデックスの型
 * @param os 出力ストリーム
 * @param matrix 出力する行列
 * @return 更新された出力ストリーム
 */
template<typename T, typename IndexType>
    requires concepts::AdditiveAbelianGroup<T> && std::integral<IndexType>
std::ostream& operator<<(std::ostream& os, const SparseMatrix<T, IndexType>& matrix) {
    os << "SparseMatrix " << matrix.rows() << "x" << matrix.cols() 
       << ", nnz=" << matrix.nnz() << "\n";

    // 非ゼロ要素だけを出力
    for (IndexType i = 0; i < matrix.rows(); ++i) {
        for (IndexType j = 0; j < matrix.cols(); ++j) {
            T val = matrix(i, j);
            if (val != T{0}) {
                os << "(" << i << "," << j << "): " << val << "\n";
            }
        }
    }
    return os;
}

// MKL対応の特殊化（必要に応じて実装）
#ifdef CALX_USE_MKL
// MKL用の特殊化や拡張をここに実装
#endif

} // namespace calx

#endif // CALX_SPARSE_MATRIX_HPP