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

/**
 * @file sparse_matrix.hpp
 * @brief Sparse matrix class implementation
 * @author MKL Algebra Library
 */

#ifndef SANGI_SPARSE_MATRIX_HPP
#define SANGI_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 sangi {

/**
 * @brief Sparse matrix storage formats
 */
enum class SparseStorageFormat {
    COO,  ///< Coordinate format
    CSR,  ///< Compressed Sparse Row
    CSC   ///< Compressed Sparse Column
};

template<typename T, typename IndexType>
    requires concepts::AdditiveAbelianGroup<T> && std::integral<IndexType>
class SparseMatrix;

// Trait used to exclude SparseMatrix from the scalar-multiply overloads. Without
// it, evaluating Field<SparseMatrix> (which probes `a*b`) re-enters overload
// resolution for SparseMatrix*SparseMatrix, whose scalar-operator candidates ask
// Scalar<SparseMatrix> = Field<SparseMatrix> || ... — a circular constraint that
// MSVC reports as C7608. The trait is a plain (non-recursive) type check, so when
// placed first in a `requires` conjunction it short-circuits before Scalar<> runs.
template<typename U> struct is_sparse_matrix : std::false_type {};
template<typename T, typename IndexType>
struct is_sparse_matrix<SparseMatrix<T, IndexType>> : std::true_type {};
template<typename U>
inline constexpr bool is_sparse_matrix_v =
    is_sparse_matrix<std::remove_cv_t<std::remove_reference_t<U>>>::value;

/**
 * @brief Sparse matrix class
 *
 * Provides an efficient representation and operations for sparse matrices
 * (matrices with many zero elements).
 *
 * @tparam T Element type (must satisfy the additive abelian group requirements)
 * @tparam IndexType Index type (default is 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:
    // Internal representation: implementation varies by storage format
    SparseStorageFormat storage_format_;
    IndexType rows_;
    IndexType cols_;

    // COO format data
    std::vector<triplet_type> coo_data_;

    // CSR format data
    std::vector<T> csr_values_;
    std::vector<IndexType> csr_col_indices_;
    std::vector<IndexType> csr_row_ptr_;

    // CSC format data
    std::vector<T> csc_values_;
    std::vector<IndexType> csc_row_indices_;
    std::vector<IndexType> csc_col_ptr_;

    // Whether the current storage format is valid
    bool coo_valid_ = false;
    bool csr_valid_ = false;
    bool csc_valid_ = false;

    // Threshold for numerical types (used for zero detection).
    // Changed to non-static member initialization so that it works even with
    // non-constexpr types (such as sangi::Float multi-precision types).
    // Built-in types are folded as constants, while Float performs a single
    // per-instance Float(0) initialization.
    T epsilon_ = std::is_floating_point_v<T> ? static_cast<T>(1e-10) : T{};

public:
    /**
     * @brief Default constructor
     */
    SparseMatrix()
        : storage_format_(SparseStorageFormat::COO), rows_(0), cols_(0), coo_valid_(true) {}

    /**
     * @brief Size-specifying constructor
     *
     * @param rows Number of rows
     * @param cols Number of columns
     * @param format Storage format (default is 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 Initializer list constructor
     *
     * @param list List of triplets (row, col, value)
     * @param rows Number of rows
     * @param cols Number of columns
     * @param format Storage format (default is 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 Copy constructor
     */
    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 Move constructor
     */
    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 Copy assignment operator
     */
    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 Move assignment operator
     */
    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 Set the epsilon value used for zero detection
     *
     * @param epsilon New epsilon value
     */
    void set_epsilon(T epsilon) {
        epsilon_ = epsilon;
    }

    /**
     * @brief Get the current epsilon value
     *
     * @return The current epsilon value
     */
    T get_epsilon() const {
        return epsilon_;
    }

    /**
     * @brief Get the number of rows
     *
     * @return Number of rows
     */
    IndexType rows() const {
        return rows_;
    }

    /**
     * @brief Get the number of columns
     *
     * @return Number of columns
     */
    IndexType cols() const {
        return cols_;
    }

    /**
     * @brief Get the number of non-zero elements
     *
     * @return Number of non-zero elements
     */
    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 Get the current storage format
     *
     * @return Storage format
     */
    SparseStorageFormat storage_format() const {
        return storage_format_;
    }

    /// Read-only access to the internal CSC arrays (for direct solvers such as SparseLU).
    /// Calls ensure_csc_valid() internally, so the initial conversion runs even from 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 Change the storage format
     *
     * @param format New storage format
     */
    void convert_to(SparseStorageFormat format) {
        if (format == storage_format_) {
            return;  // Already in the same format
        }

        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 Get the element at the specified position
     *
     * @param i Row index
     * @param j Column index
     * @return Element value
     */
    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 Set a value (overwrites existing value or inserts a new one)
     *
     * @param i Row index
     * @param j Column index
     * @param value Value to set
     */
    void set_coeff(IndexType i, IndexType j, const T& value) {
        validate_indices(i, j);

        // Do not store values close to zero
        if (std::abs(value) <= epsilon_) {
            // Remove the existing value
            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 Remove an element
     *
     * @param i Row index
     * @param j Column index
     */
    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;
        }

        // Invalidate the other formats as well
        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 Get a row vector
     *
     * @param i Row index
     * @return Row vector
     */
    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 Get a column vector
     *
     * @param j Column index
     * @return Column vector
     */
    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 Matrix transpose
     *
     * @return Transposed matrix
     */
    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: {
                // (M is CSR.) Mᵀ in CSC is just M's CSR arrays reinterpreted:
                //   values  : same (row-major) order
                //   row idx : M's column indices  (= Mᵀ's row indices)
                //   col_ptr : M's row_ptr          (= Mᵀ's column boundaries)
                // col_ptr therefore has size rows_+1 = result.cols_+1 (correct even
                // for rectangular matrices; the old count-based version mis-sized it).
                result.csc_values_ = csr_values_;
                result.csc_row_indices_ = csr_col_indices_;
                result.csc_col_ptr_ = csr_row_ptr_;
                result.csc_valid_ = true;
                // The constructor pre-filled an (empty) CSR for `result`; discard it
                // and rebuild CSR from the freshly populated CSC.
                result.csr_valid_ = false;
                result.coo_valid_ = false;
                result.storage_format_ = SparseStorageFormat::CSC;
                result.convert_to(SparseStorageFormat::CSR);
                break;
            }
            case SparseStorageFormat::CSC: {
                // (M is CSC.) Mᵀ in CSR is M's CSC arrays reinterpreted symmetrically.
                result.csr_values_ = csc_values_;
                result.csr_col_indices_ = csc_row_indices_;
                result.csr_row_ptr_ = csc_col_ptr_;
                result.csr_valid_ = true;
                result.csc_valid_ = false;
                result.coo_valid_ = false;
                result.storage_format_ = SparseStorageFormat::CSR;
                result.convert_to(SparseStorageFormat::CSC);
                break;
            }
        }

        return result;
    }

    /**
     * @brief Element-wise scalar multiplication
     *
     * @param scalar Scalar value
     * @return Matrix scaled by the scalar
     */
    template<typename ScalarType>
        requires (!is_sparse_matrix_v<ScalarType>) && 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 Scalar division of the matrix
     *
     * @param scalar Scalar value
     * @return Matrix divided by the scalar
     */
    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 Matrix addition
     *
     * @param other Matrix to add
     * @return Result of the addition
     */
    SparseMatrix operator+(const SparseMatrix& other) const {
        if (rows_ != other.rows_ || cols_ != other.cols_) {
            throw std::invalid_argument("Matrix dimensions must match");
        }

        // Addition is simplest in COO format
        ensure_coo_valid();
        other.ensure_coo_valid();

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

        // Copy the existing data
        result.coo_data_ = coo_data_;

        // Add or accumulate elements from the other matrix
        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);
            }
        }

        // Remove elements that became zero
        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;

        // Convert to another format if necessary
        if (storage_format_ != SparseStorageFormat::COO) {
            result.convert_to(storage_format_);
        }

        return result;
    }

    /**
     * @brief Matrix subtraction
     *
     * @param other Matrix to subtract
     * @return Result of the subtraction
     */
    SparseMatrix operator-(const SparseMatrix& other) const {
        if (rows_ != other.rows_ || cols_ != other.cols_) {
            throw std::invalid_argument("Matrix dimensions must match");
        }

        // Subtraction is simplest in COO format
        ensure_coo_valid();
        other.ensure_coo_valid();

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

        // Copy the existing data
        result.coo_data_ = coo_data_;

        // Subtract elements of the other matrix
        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);
            }
        }

        // Remove elements that became zero
        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;

        // Convert to another format if necessary
        if (storage_format_ != SparseStorageFormat::COO) {
            result.convert_to(storage_format_);
        }

        return result;
    }

    /**
     * @brief Matrix-vector product
     *
     * @param vec Right-hand side vector
     * @return Resulting vector
     */
    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: sparse matrix x dense matrix -> dense matrix (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 Matrix-matrix product
     *
     * @param other Right-hand side matrix
     * @return Resulting matrix
     */
    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);

        // Convert to the optimal formats
        ensure_csr_valid();  // Left operand uses CSR
        other.ensure_csc_valid();  // Right operand uses 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;

                // Efficient sparse-sparse multiplication
                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;

        // Convert to another format if necessary
        if (storage_format_ != SparseStorageFormat::COO) {
            result.convert_to(storage_format_);
        }

        return result;
    }

    /**
     * @brief Element access operator for the matrix
     *
     * @param i Row index
     * @param j Column index
     * @return Reference (proxy) to the element
     */
    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) {}

        // Implicit conversion operator
        operator T() const {
            return matrix_.coeff(i_, j_);
        }

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

        // Compound assignment operators
        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 Element access operator (non-const)
     *
     * @param i Row index
     * @param j Column index
     * @return Element proxy
     */
    ElementProxy operator()(IndexType i, IndexType j) {
        validate_indices(i, j);
        return ElementProxy(*this, i, j);
    }

    /**
     * @brief Element access operator (const)
     *
     * @param i Row index
     * @param j Column index
     * @return Element value
     */
    T operator()(IndexType i, IndexType j) const {
        return coeff(i, j);
    }

    /**
     * @brief Get a block (submatrix) of the matrix
     *
     * @param start_row Starting row index
     * @param start_col Starting column index
     * @param block_rows Number of rows in the block
     * @param block_cols Number of columns in the block
     * @return Submatrix
     */
    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;

        // Convert to another format if necessary
        if (storage_format_ != SparseStorageFormat::COO) {
            result.convert_to(storage_format_);
        }

        return result;
    }

    /**
     * @brief Create a diagonal matrix
     *
     * @param diag Vector of diagonal entries
     * @return Diagonal matrix
     */
    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 Create an identity matrix
     *
     * @param n Matrix size
     * @return Identity matrix
     */
    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 Create a tridiagonal matrix
     *
     * @param diag Vector of diagonal entries
     * @param lower Vector of subdiagonal entries
     * @param upper Vector of superdiagonal entries
     * @return Tridiagonal matrix
     */
    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);

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

        // Subdiagonal entries
        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]);
            }
        }

        // Superdiagonal entries
        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 Create a zero matrix
     *
     * @param rows Number of rows
     * @param cols Number of columns
     * @return Zero matrix
     */
    static SparseMatrix zero(IndexType rows, IndexType cols) {
        return SparseMatrix(rows, cols);
    }

private:
    /**
     * @brief Validate the indices
     *
     * @param i Row index
     * @param j Column index
     */
    void validate_indices(IndexType i, IndexType j) const {
        if (i >= rows_ || j >= cols_) {
            throw std::out_of_range("Matrix indices out of bounds");
        }
    }

    /**
     * @brief Validate the row index
     *
     * @param i Row index
     */
    void validate_row_index(IndexType i) const {
        if (i >= rows_) {
            throw std::out_of_range("Row index out of bounds");
        }
    }

    /**
     * @brief Validate the column index
     *
     * @param j Column index
     */
    void validate_col_index(IndexType j) const {
        if (j >= cols_) {
            throw std::out_of_range("Column index out of bounds");
        }
    }

    /**
     * @brief Ensure that COO format data is valid
     */
    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 Ensure that CSR format data is valid
     */
    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 Ensure that CSC format data is valid
     */
    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 conversion
     */
    void convert_coo_to_csr() {
        if (!coo_valid_ || coo_data_.empty()) {
            // Empty matrix case
            csr_values_.clear();
            csr_col_indices_.clear();
            csr_row_ptr_.resize(rows_ + 1, 0);
            csr_valid_ = true;
            return;
        }

        // Sort COO data in row-major order
        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));
            });

        // Initialize the CSR data structures
        const size_type nnz = coo_data_.size();
        csr_values_.resize(nnz);
        csr_col_indices_.resize(nnz);
        csr_row_ptr_.resize(rows_ + 1, 0);

        // Count non-zero elements in each row
        for (const auto& [i, j, val] : coo_data_) {
            csr_row_ptr_[i + 1]++;
        }

        // Compute the prefix sum to obtain row pointers
        for (IndexType i = 0; i < rows_; ++i) {
            csr_row_ptr_[i + 1] += csr_row_ptr_[i];
        }

        // Fill in the values and column indices
        for (const auto& [i, j, val] : coo_data_) {
            const IndexType pos = csr_row_ptr_[i]++;
            csr_values_[pos] = val;
            csr_col_indices_[pos] = j;
        }

        // Fix up the row pointers (they were mutated during prefix sum)
        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 conversion
     */
    void convert_coo_to_csc() {
        if (!coo_valid_ || coo_data_.empty()) {
            // Empty matrix case
            csc_values_.clear();
            csc_row_indices_.clear();
            csc_col_ptr_.resize(cols_ + 1, 0);
            csc_valid_ = true;
            return;
        }

        // Sort COO data in column-major order
        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));
            });

        // Initialize the CSC data structures
        const size_type nnz = coo_data_.size();
        csc_values_.resize(nnz);
        csc_row_indices_.resize(nnz);
        csc_col_ptr_.resize(cols_ + 1, 0);

        // Count non-zero elements in each column
        for (const auto& [i, j, val] : coo_data_) {
            csc_col_ptr_[j + 1]++;
        }

        // Compute the prefix sum to obtain column pointers
        for (IndexType j = 0; j < cols_; ++j) {
            csc_col_ptr_[j + 1] += csc_col_ptr_[j];
        }

        // Fill in the values and row indices
        for (const auto& [i, j, val] : coo_data_) {
            const IndexType pos = csc_col_ptr_[j]++;
            csc_values_[pos] = val;
            csc_row_indices_[pos] = i;
        }

        // Fix up the column pointers (they were mutated during prefix sum)
        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 conversion
     */
    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 conversion
     */
    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 conversion
     */
    void convert_csr_to_csc() {
        if (!csr_valid_) {
            throw std::runtime_error("CSR data is not valid");
        }

        // First convert to COO, then convert COO to CSC
        convert_csr_to_coo();
        convert_coo_to_csc();
    }

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

        // First convert to COO, then convert COO to CSR
        convert_csc_to_coo();
        convert_coo_to_csr();
    }

    /**
     * @brief Convert an initializer list to CSR format
     *
     * @param list List of triplets
     */
    void convert_to_csr(const std::initializer_list<triplet_type>& list) {
        std::vector<triplet_type> sorted_list(list);

        // Sort in row-major order
        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));
            });

        // Initialize the CSR data structures
        const size_type nnz = sorted_list.size();
        csr_values_.resize(nnz);
        csr_col_indices_.resize(nnz);
        csr_row_ptr_.resize(rows_ + 1, 0);

        // Count non-zero elements in each row
        for (const auto& [i, j, val] : sorted_list) {
            csr_row_ptr_[i + 1]++;
        }

        // Compute the prefix sum to obtain row pointers
        for (IndexType i = 0; i < rows_; ++i) {
            csr_row_ptr_[i + 1] += csr_row_ptr_[i];
        }

        // Fill in the values and column indices
        for (const auto& [i, j, val] : sorted_list) {
            const IndexType pos = csr_row_ptr_[i]++;
            csr_values_[pos] = val;
            csr_col_indices_[pos] = j;
        }

        // Fix up the row pointers
        for (IndexType i = rows_; i > 0; --i) {
            csr_row_ptr_[i] = csr_row_ptr_[i - 1];
        }
        csr_row_ptr_[0] = 0;
    }

    /**
     * @brief Convert an initializer list to CSC format
     *
     * @param list List of triplets
     */
    void convert_to_csc(const std::initializer_list<triplet_type>& list) {
        std::vector<triplet_type> sorted_list(list);

        // Sort in column-major order
        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));
            });

        // Initialize the CSC data structures
        const size_type nnz = sorted_list.size();
        csc_values_.resize(nnz);
        csc_row_indices_.resize(nnz);
        csc_col_ptr_.resize(cols_ + 1, 0);

        // Count non-zero elements in each column
        for (const auto& [i, j, val] : sorted_list) {
            csc_col_ptr_[j + 1]++;
        }

        // Compute the prefix sum to obtain column pointers
        for (IndexType j = 0; j < cols_; ++j) {
            csc_col_ptr_[j + 1] += csc_col_ptr_[j];
        }

        // Fill in the values and row indices
        for (const auto& [i, j, val] : sorted_list) {
            const IndexType pos = csc_col_ptr_[j]++;
            csc_values_[pos] = val;
            csc_row_indices_[pos] = i;
        }

        // Fix up the column pointers
        for (IndexType j = cols_; j > 0; --j) {
            csc_col_ptr_[j] = csc_col_ptr_[j - 1];
        }
        csc_col_ptr_[0] = 0;
    }

    /**
     * @brief Element lookup in COO format
     */
    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 Element lookup in CSR format
     */
    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];
            }
            // Assumes column indices are sorted
            if (csr_col_indices_[k] > j) {
                break;
            }
        }
        return T{0};
    }

    /**
     * @brief Element lookup in CSC format
     */
    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];
            }
            // Assumes row indices are sorted
            if (csc_row_indices_[k] > i) {
                break;
            }
        }
        return T{0};
    }

    /**
     * @brief Element assignment in COO format
     */
    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;
            }
        }
        // If no existing element is found, append a new one
        coo_data_.emplace_back(i, j, value);
    }

    /**
     * @brief Element assignment in CSR format
     */
    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) {
                // Update the existing element
                csr_values_[k] = value;
                return;
            }
            if (csr_col_indices_[k] > j) {
                // Insert a new element
                csr_values_.insert(csr_values_.begin() + k, value);
                csr_col_indices_.insert(csr_col_indices_.begin() + k, j);
                // Update row pointers
                for (IndexType m = i + 1; m <= rows_; ++m) {
                    ++csr_row_ptr_[m];
                }
                return;
            }
        }

        // Append at the end of the row
        csr_values_.insert(csr_values_.begin() + end, value);
        csr_col_indices_.insert(csr_col_indices_.begin() + end, j);
        // Update row pointers
        for (IndexType m = i + 1; m <= rows_; ++m) {
            ++csr_row_ptr_[m];
        }
    }

    /**
     * @brief Element assignment in CSC format
     */
    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) {
                // Update the existing element
                csc_values_[k] = value;
                return;
            }
            if (csc_row_indices_[k] > i) {
                // Insert a new element
                csc_values_.insert(csc_values_.begin() + k, value);
                csc_row_indices_.insert(csc_row_indices_.begin() + k, i);
                // Update column pointers
                for (IndexType m = j + 1; m <= cols_; ++m) {
                    ++csc_col_ptr_[m];
                }
                return;
            }
        }

        // Append at the end of the column
        csc_values_.insert(csc_values_.begin() + end, value);
        csc_row_indices_.insert(csc_row_indices_.begin() + end, i);
        // Update column pointers
        for (IndexType m = j + 1; m <= cols_; ++m) {
            ++csc_col_ptr_[m];
        }
    }

    /**
     * @brief Element removal in COO format
     */
    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 Element removal in CSR format
     */
    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) {
                // Remove the element
                csr_values_.erase(csr_values_.begin() + k);
                csr_col_indices_.erase(csr_col_indices_.begin() + k);
                // Update row pointers
                for (IndexType m = i + 1; m <= rows_; ++m) {
                    --csr_row_ptr_[m];
                }
                return;
            }
        }
    }

    /**
     * @brief Element removal in CSC format
     */
    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) {
                // Remove the element
                csc_values_.erase(csc_values_.begin() + k);
                csc_row_indices_.erase(csc_row_indices_.begin() + k);
                // Update column pointers
                for (IndexType m = j + 1; m <= cols_; ++m) {
                    --csc_col_ptr_[m];
                }
                return;
            }
        }
    }
};

/**
 * @brief Free-standing transpose function
 *
 * @tparam T Element type
 * @tparam IndexType Index type
 * @param matrix Matrix to transpose
 * @return Transposed matrix
 */
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 Free-standing scalar multiplication
 *
 * @tparam T Element type
 * @tparam ScalarType Scalar type
 * @tparam IndexType Index type
 * @param scalar Scalar value
 * @param matrix Matrix
 * @return Matrix scaled by the scalar
 */
template<typename T, typename ScalarType, typename IndexType>
    requires (!is_sparse_matrix_v<ScalarType>) && concepts::Scalar<ScalarType>
          && concepts::AdditiveAbelianGroup<T> && std::integral<IndexType>
SparseMatrix<T, IndexType> operator*(const ScalarType& scalar, const SparseMatrix<T, IndexType>& matrix) {
    return matrix * scalar;
}

/**
 * @brief Output stream operator overload
 *
 * @tparam T Element type
 * @tparam IndexType Index type
 * @param os Output stream
 * @param matrix Matrix to output
 * @return The updated output stream
 */
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";

    // Output only the non-zero elements
    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-aware specializations (implement as needed)
#ifdef SANGI_USE_MKL
// MKL-specific specializations and extensions go here
#endif

} // namespace sangi

#endif // SANGI_SPARSE_MATRIX_HPP