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

// matrix.hpp
// Matrix class definitions
//
// This file defines matrix classes for the MKL algebra library.
// It contains implementations of dynamic-size and fixed-size matrices.

#ifndef SANGI_MATRIX_HPP
#define SANGI_MATRIX_HPP

#include <cassert>
#include <math/core/common.hpp>
#include <math/core/traits.hpp>
#include <math/core/generic_math.hpp>
#include <math/core/computation_policy.hpp>
#include <math/core/vector.hpp>
#include <math/core/BaseMatrix.hpp>

#include <algorithm>
#include <array>
#include <bit>
#include <initializer_list>
#include <memory>
#include <optional>
#include <numeric>
#include <vector>
#include <iomanip>
#include <iostream>
#include <sstream>

// GCC -O2 emits an array-bounds false positive when expanding
// PacketTraits::load/store, so suppress it for the whole file.
#ifdef __GNUC__
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Warray-bounds"
#endif

namespace sangi {

    // Forward declarations (referenced from operator^)
    template<typename T> class Matrix;
    template<typename T>
    [[nodiscard]] Matrix<T> power(const Matrix<T>& base, unsigned int exponent);

    //-----------------------------------------------------------------------------
    // Matrix storage definitions (implemented inside the detail namespace)
    //-----------------------------------------------------------------------------
    namespace detail {

        // Dynamic-size matrix storage (shared_ptr<T[]> based, compatible with the sangi Vector<T>)
        //
        // The internal storage is held as std::shared_ptr<T[]>, with size_/capacity_
        // tracking used and allocated element counts separately (reproducing std::vector semantics).
        // Using shared_ptr allows future integration with BaseMatrix<T>'s shared_storage_.
        template<typename T>
        class DynamicMatrixStorage {
        public:
            using value_type = T;
            using size_type = std::size_t;
            using reference = T&;
            using const_reference = const T&;
            using pointer = T*;
            using const_pointer = const T*;
            using iterator = T*;
            using const_iterator = const T*;

            // Default constructor
            DynamicMatrixStorage() noexcept
                : rows_(0), cols_(0), data_(), size_(0), capacity_(0) {}

            // Size-specifying constructor (initialized with default value = T())
            DynamicMatrixStorage(size_type rows, size_type cols)
                : rows_(rows), cols_(cols)
                , data_(rows * cols > 0 ? std::shared_ptr<T[]>(new T[rows * cols]()) : std::shared_ptr<T[]>())
                , size_(rows * cols), capacity_(rows * cols) {}

            // Constructor specifying size and initial value
            DynamicMatrixStorage(size_type rows, size_type cols, const T& value)
                : rows_(rows), cols_(cols)
                , data_(rows * cols > 0 ? std::shared_ptr<T[]>(new T[rows * cols]) : std::shared_ptr<T[]>())
                , size_(rows * cols), capacity_(rows * cols) {
                for (size_type i = 0; i < size_; ++i) data_[i] = value;
            }

            // Copy constructor (deep copy)
            DynamicMatrixStorage(const DynamicMatrixStorage& other)
                : rows_(other.rows_), cols_(other.cols_)
                , data_(other.size_ > 0 ? std::shared_ptr<T[]>(new T[other.size_]) : std::shared_ptr<T[]>())
                , size_(other.size_), capacity_(other.size_) {
                for (size_type i = 0; i < size_; ++i) data_[i] = other.data_[i];
            }

            // Move constructor
            DynamicMatrixStorage(DynamicMatrixStorage&& other) noexcept
                : rows_(other.rows_), cols_(other.cols_)
                , data_(std::move(other.data_))
                , size_(other.size_), capacity_(other.capacity_) {
                other.rows_ = 0;
                other.cols_ = 0;
                other.size_ = 0;
                other.capacity_ = 0;
            }

            // Copy assignment (deep copy)
            DynamicMatrixStorage& operator=(const DynamicMatrixStorage& other) {
                if (this == &other) return *this;
                if (other.size_ > 0) {
                    std::shared_ptr<T[]> new_data(new T[other.size_]);
                    for (size_type i = 0; i < other.size_; ++i) new_data[i] = other.data_[i];
                    data_ = std::move(new_data);
                } else {
                    data_.reset();
                }
                rows_ = other.rows_;
                cols_ = other.cols_;
                size_ = other.size_;
                capacity_ = other.size_;
                return *this;
            }

            // Move assignment
            DynamicMatrixStorage& operator=(DynamicMatrixStorage&& other) noexcept {
                if (this == &other) return *this;
                rows_ = other.rows_;
                cols_ = other.cols_;
                data_ = std::move(other.data_);
                size_ = other.size_;
                capacity_ = other.capacity_;
                other.rows_ = 0;
                other.cols_ = 0;
                other.size_ = 0;
                other.capacity_ = 0;
                return *this;
            }

            // Destructor
            ~DynamicMatrixStorage() = default;

            // Element access in (row, col) form — bounds checked via assert in Debug builds
            reference operator()(size_type row, size_type col) {
                assert(row < rows_ && "Matrix::operator(): row out of range");
                assert(col < cols_ && "Matrix::operator(): col out of range");
                return data_[row * cols_ + col];
            }

            const_reference operator()(size_type row, size_type col) const {
                assert(row < rows_ && "Matrix::operator(): row out of range");
                assert(col < cols_ && "Matrix::operator(): col out of range");
                return data_[row * cols_ + col];
            }

            // Range-checked element access
            reference at(size_type row, size_type col) {
                check_index(row, col);
                return data_[row * cols_ + col];
            }

            const_reference at(size_type row, size_type col) const {
                check_index(row, col);
                return data_[row * cols_ + col];
            }

            // Data pointer access
            pointer data() noexcept {
                return data_.get();
            }

            const_pointer data() const noexcept {
                return data_.get();
            }

            // Iterators (raw T* — satisfies the random_access_iterator concept)
            iterator begin() noexcept { return data_.get(); }
            const_iterator begin() const noexcept { return data_.get(); }
            const_iterator cbegin() const noexcept { return data_.get(); }
            iterator end() noexcept { return data_.get() + size_; }
            const_iterator end() const noexcept { return data_.get() + size_; }
            const_iterator cend() const noexcept { return data_.get() + size_; }

            // Capacity
            bool empty() const noexcept { return size_ == 0; }

            size_type size() const noexcept { return size_; }

            size_type rows() const noexcept { return rows_; }

            size_type cols() const noexcept { return cols_; }

            size_type capacity() const noexcept { return capacity_; }

            // Memory allocation (capacity expansion only, size unchanged; std::vector::reserve compatible)
            void reserve(size_type new_cap) {
                if (new_cap <= capacity_) return;
                std::shared_ptr<T[]> new_data(new T[new_cap]());
                for (size_type i = 0; i < size_; ++i) new_data[i] = data_[i];
                data_ = std::move(new_data);
                capacity_ = new_cap;
            }

            // Shrink (adjust capacity down to size; std::vector::shrink_to_fit compatible)
            void shrink_to_fit() {
                if (capacity_ == size_) return;
                if (size_ == 0) {
                    data_.reset();
                    capacity_ = 0;
                    return;
                }
                std::shared_ptr<T[]> new_data(new T[size_]);
                for (size_type i = 0; i < size_; ++i) new_data[i] = data_[i];
                data_ = std::move(new_data);
                capacity_ = size_;
            }

            // Resize (new elements initialized with T())
            void resize(size_type rows, size_type cols) {
                const size_type new_size = rows * cols;
                if (new_size > capacity_) {
                    std::shared_ptr<T[]> new_data(new T[new_size]());
                    for (size_type i = 0; i < size_; ++i) new_data[i] = data_[i];
                    data_ = std::move(new_data);
                    capacity_ = new_size;
                } else if (new_size > size_) {
                    for (size_type i = size_; i < new_size; ++i) data_[i] = T();
                }
                rows_ = rows;
                cols_ = cols;
                size_ = new_size;
            }

            void resize(size_type rows, size_type cols, const value_type& value) {
                const size_type new_size = rows * cols;
                if (new_size > capacity_) {
                    std::shared_ptr<T[]> new_data(new T[new_size]);
                    for (size_type i = 0; i < size_; ++i) new_data[i] = data_[i];
                    for (size_type i = size_; i < new_size; ++i) new_data[i] = value;
                    data_ = std::move(new_data);
                    capacity_ = new_size;
                } else if (new_size > size_) {
                    for (size_type i = size_; i < new_size; ++i) data_[i] = value;
                }
                rows_ = rows;
                cols_ = cols;
                size_ = new_size;
            }

            // Clear (capacity preserved; std::vector::clear compatible)
            void clear() noexcept {
                rows_ = 0;
                cols_ = 0;
                size_ = 0;
            }

            // Change row/column counts
            void reshape(size_type rows, size_type cols) {
                if (rows * cols != rows_ * cols_) {
                    throw DimensionError("Matrix reshape: total element count must be preserved");
                }
                rows_ = rows;
                cols_ = cols;
            }

            // Whether the matrix is square
            bool is_square() const noexcept {
                return rows_ == cols_;
            }

            // Obtain the shared_ptr (for BaseMatrix integration and for sharing lifetime with views)
            const std::shared_ptr<T[]>& shared_data() const noexcept {
                return data_;
            }

        private:
            // Index check
            void check_index(size_type row, size_type col) const {
                if (row >= rows_) {
                    throw IndexError("Matrix row index out of range: " + std::to_string(row) +
                        " >= " + std::to_string(rows_));
                }
                if (col >= cols_) {
                    throw IndexError("Matrix column index out of range: " + std::to_string(col) +
                        " >= " + std::to_string(cols_));
                }
            }

            size_type rows_;
            size_type cols_;
            std::shared_ptr<T[]> data_;
            size_type size_;
            size_type capacity_;
        };

        // Fixed-size matrix storage (std::array based)
        template<typename T, std::size_t Rows, std::size_t Cols>
        class StaticMatrixStorage {
        public:
            using value_type = T;
            using size_type = std::size_t;
            using reference = T&;
            using const_reference = const T&;
            using pointer = T*;
            using const_pointer = const T*;
            using iterator = pointer;
            using const_iterator = const_pointer;

            // Default constructor (zero-initialized)
            StaticMatrixStorage() : data_{} {}

            // Constructor initializing with a value
            explicit StaticMatrixStorage(const T& value) {
                std::fill(begin(), end(), value);
            }

            // Copy/move constructors and assignments
            StaticMatrixStorage(const StaticMatrixStorage&) = default;
            StaticMatrixStorage(StaticMatrixStorage&&) noexcept = default;
            StaticMatrixStorage& operator=(const StaticMatrixStorage&) = default;
            StaticMatrixStorage& operator=(StaticMatrixStorage&&) noexcept = default;

            // Destructor
            ~StaticMatrixStorage() = default;

            // Element access in (row, col) form — bounds checked via assert in Debug builds
            reference operator()(size_type row, size_type col) {
                assert(row < Rows && "StaticMatrix::operator(): row out of range");
                assert(col < Cols && "StaticMatrix::operator(): col out of range");
                return data_[row * Cols + col];
            }

            const_reference operator()(size_type row, size_type col) const {
                assert(row < Rows && "StaticMatrix::operator(): row out of range");
                assert(col < Cols && "StaticMatrix::operator(): col out of range");
                return data_[row * Cols + col];
            }

            // Range-checked element access
            reference at(size_type row, size_type col) {
                check_index(row, col);
                return data_[row * Cols + col];
            }

            const_reference at(size_type row, size_type col) const {
                check_index(row, col);
                return data_[row * Cols + col];
            }

            // Data pointer access
            pointer data() noexcept {
                return data_.data();
            }

            const_pointer data() const noexcept {
                return data_.data();
            }

            // Iterators
            iterator begin() noexcept {
                return data_.data();
            }

            const_iterator begin() const noexcept {
                return data_.data();
            }

            const_iterator cbegin() const noexcept {
                return data_.data();
            }

            iterator end() noexcept {
                return data_.data() + Rows * Cols;
            }

            const_iterator end() const noexcept {
                return data_.data() + Rows * Cols;
            }

            const_iterator cend() const noexcept {
                return data_.data() + Rows * Cols;
            }

            // Capacity
            constexpr bool empty() const noexcept {
                return size() == 0;
            }

            constexpr size_type size() const noexcept {
                return Rows * Cols;
            }

            constexpr size_type rows() const noexcept {
                return Rows;
            }

            constexpr size_type cols() const noexcept {
                return Cols;
            }

            // Whether the matrix is square
            constexpr bool is_square() const noexcept {
                return Rows == Cols;
            }

        private:
            // Index check
            void check_index(size_type row, size_type col) const {
                if (row >= Rows) {
                    throw IndexError("StaticMatrix row index out of range: " + std::to_string(row) +
                        " >= " + std::to_string(Rows));
                }
                if (col >= Cols) {
                    throw IndexError("StaticMatrix column index out of range: " + std::to_string(col) +
                        " >= " + std::to_string(Cols));
                }
            }

            std::array<T, Rows* Cols> data_;
        };

        // Aligned matrix storage (MKL-optimized)
        template<typename T>
        class AlignedMatrixStorage {
        public:
            using value_type = T;
            using size_type = std::size_t;
            using reference = T&;
            using const_reference = const T&;
            using pointer = T*;
            using const_pointer = const T*;
            using iterator = pointer;
            using const_iterator = const_pointer;

            // Alignment value
            static constexpr size_type alignment = memory::simd_alignment();

            // Default constructor
            AlignedMatrixStorage() : data_(nullptr), rows_(0), cols_(0), capacity_(0) {}

            // Size-specifying constructor
            AlignedMatrixStorage(size_type rows, size_type cols)
                : data_(nullptr), rows_(rows), cols_(cols), capacity_(rows* cols) {
                if (capacity_ > 0) {
                    data_ = memory::aligned_alloc<T>(capacity_, alignment);
                    std::uninitialized_default_construct_n(data_, capacity_);
                }
            }

            // Constructor specifying size and initial value
            AlignedMatrixStorage(size_type rows, size_type cols, const T& value)
                : data_(nullptr), rows_(rows), cols_(cols), capacity_(rows* cols) {
                if (capacity_ > 0) {
                    data_ = memory::aligned_alloc<T>(capacity_, alignment);
                    std::uninitialized_fill_n(data_, capacity_, value);
                }
            }

            // Copy constructor
            AlignedMatrixStorage(const AlignedMatrixStorage& other)
                : data_(nullptr), rows_(other.rows_), cols_(other.cols_), capacity_(other.capacity_) {
                if (capacity_ > 0) {
                    data_ = memory::aligned_alloc<T>(capacity_, alignment);
                    std::uninitialized_copy_n(other.data_, capacity_, data_);
                }
            }

            // Move constructor
            AlignedMatrixStorage(AlignedMatrixStorage&& other) noexcept
                : data_(other.data_), rows_(other.rows_), cols_(other.cols_), capacity_(other.capacity_) {
                other.data_ = nullptr;
                other.rows_ = 0;
                other.cols_ = 0;
                other.capacity_ = 0;
            }

            // Copy assignment operator
            AlignedMatrixStorage& operator=(const AlignedMatrixStorage& other) {
                if (this == &other) return *this;

                // Discard the old data
                clear();
                if (data_) {
                    memory::aligned_free(data_);
                }

                // Allocate the new data
                rows_ = other.rows_;
                cols_ = other.cols_;
                capacity_ = other.capacity_;

                if (capacity_ > 0) {
                    data_ = memory::aligned_alloc<T>(capacity_, alignment);
                    std::uninitialized_copy_n(other.data_, capacity_, data_);
                }
                else {
                    data_ = nullptr;
                }

                return *this;
            }

            // Move assignment operator
            AlignedMatrixStorage& operator=(AlignedMatrixStorage&& other) noexcept {
                if (this == &other) return *this;

                // Discard the old data
                clear();
                if (data_) {
                    memory::aligned_free(data_);
                }

                // Move the other side's data
                data_ = other.data_;
                rows_ = other.rows_;
                cols_ = other.cols_;
                capacity_ = other.capacity_;

                other.data_ = nullptr;
                other.rows_ = 0;
                other.cols_ = 0;
                other.capacity_ = 0;

                return *this;
            }

            // Destructor
            ~AlignedMatrixStorage() {
                clear();
                if (data_) {
                    memory::aligned_free(data_);
                }
            }

            // Element access in (row, col) form
            reference operator()(size_type row, size_type col) {
                return data_[row * cols_ + col];
            }

            const_reference operator()(size_type row, size_type col) const {
                return data_[row * cols_ + col];
            }

            // Range-checked element access
            reference at(size_type row, size_type col) {
                check_index(row, col);
                return data_[row * cols_ + col];
            }

            const_reference at(size_type row, size_type col) const {
                check_index(row, col);
                return data_[row * cols_ + col];
            }

            // Data pointer access
            pointer data() noexcept {
                return data_;
            }

            const_pointer data() const noexcept {
                return data_;
            }

            // Iterators
            iterator begin() noexcept {
                return data_;
            }

            const_iterator begin() const noexcept {
                return data_;
            }

            const_iterator cbegin() const noexcept {
                return data_;
            }

            iterator end() noexcept {
                return data_ + rows_ * cols_;
            }

            const_iterator end() const noexcept {
                return data_ + rows_ * cols_;
            }

            const_iterator cend() const noexcept {
                return data_ + rows_ * cols_;
            }

            // Capacity
            bool empty() const noexcept {
                return rows_ == 0 || cols_ == 0;
            }

            size_type size() const noexcept {
                return rows_ * cols_;
            }

            size_type rows() const noexcept {
                return rows_;
            }

            size_type cols() const noexcept {
                return cols_;
            }

            size_type capacity() const noexcept {
                return capacity_;
            }

            // Memory allocation
            void reserve(size_type new_cap) {
                if (new_cap <= capacity_) return;

                // Allocate new memory
                pointer new_data = memory::aligned_alloc<T>(new_cap, alignment);

                // Move the old data
                if (size() > 0) {
                    std::uninitialized_move_n(data_, size(), new_data);
                    // Destroy the old objects
                    std::destroy_n(data_, size());
                }

                // Free the old memory
                if (data_) {
                    memory::aligned_free(data_);
                }

                // Install the new memory
                data_ = new_data;
                capacity_ = new_cap;
            }

            // Clear
            void clear() noexcept {
                if (data_) {
                    std::destroy_n(data_, size());
                }
                rows_ = 0;
                cols_ = 0;
            }

            // Resize
            void resize(size_type rows, size_type cols) {
                size_type new_size = rows * cols;

                if (new_size > capacity_) {
                    // Grow
                    reserve(new_size);
                }

                if (new_size > size()) {
                    // Append elements
                    std::uninitialized_default_construct_n(data_ + size(), new_size - size());
                }
                else if (new_size < size()) {
                    // Remove elements
                    std::destroy_n(data_ + new_size, size() - new_size);
                }

                rows_ = rows;
                cols_ = cols;
            }

            void resize(size_type rows, size_type cols, const value_type& value) {
                size_type new_size = rows * cols;

                if (new_size > capacity_) {
                    // Grow
                    reserve(new_size);
                }

                if (new_size > size()) {
                    // Append elements
                    std::uninitialized_fill_n(data_ + size(), new_size - size(), value);
                }
                else if (new_size < size()) {
                    // Remove elements
                    std::destroy_n(data_ + new_size, size() - new_size);
                }

                rows_ = rows;
                cols_ = cols;
            }

            // Change row/column counts
            void reshape(size_type rows, size_type cols) {
                if (rows * cols != rows_ * cols_) {
                    throw DimensionError("Matrix reshape: total element count must be preserved");
                }
                rows_ = rows;
                cols_ = cols;
            }

            // Whether the matrix is square
            bool is_square() const noexcept {
                return rows_ == cols_;
            }

        private:
            // Index check
            void check_index(size_type row, size_type col) const {
                if (row >= rows_) {
                    throw IndexError("AlignedMatrix row index out of range: " + std::to_string(row) +
                        " >= " + std::to_string(rows_));
                }
                if (col >= cols_) {
                    throw IndexError("AlignedMatrix column index out of range: " + std::to_string(col) +
                        " >= " + std::to_string(cols_));
                }
            }

            pointer data_;
            size_type rows_;
            size_type cols_;
            size_type capacity_;
        };

        // Matrix view implementation
        template<typename MatrixType>
        class MatrixView {
        public:
            // Obtain type definitions from the host type
            using value_type = typename std::remove_reference_t<MatrixType>::value_type;
            using size_type = std::size_t;
            using reference = value_type&;
            using const_reference = const value_type&;

            // Reference type to the underlying matrix
            using matrix_reference = MatrixType&;

            // Create a view of the entire matrix
            explicit MatrixView(matrix_reference matrix)
                : matrix_(matrix), row_offset_(0), col_offset_(0),
                num_rows_(matrix.rows()), num_cols_(matrix.cols()) {
            }

            // Create a partial view of the matrix
            MatrixView(matrix_reference matrix,
                size_type row_offset, size_type col_offset,
                size_type num_rows, size_type num_cols)
                : matrix_(matrix), row_offset_(row_offset), col_offset_(col_offset),
                num_rows_(num_rows), num_cols_(num_cols) {

                // Range check
                if (row_offset + num_rows > matrix.rows() || col_offset + num_cols > matrix.cols()) {
                    throw IndexError("MatrixView: view range exceeds matrix bounds");
                }
            }

            // Copy constructor
            MatrixView(const MatrixView&) = default;

            // Copy assignment is forbidden (copying the view itself is not allowed)
            MatrixView& operator=(const MatrixView&) = delete;

            // Assignment that copies matrix contents
            template<typename OtherMatrixType>
            MatrixView& operator=(const OtherMatrixType& other) {
                assign(other);
                return *this;
            }

            // Element access
            reference operator()(size_type row, size_type col) {
                return matrix_(row_offset_ + row, col_offset_ + col);
            }

            const_reference operator()(size_type row, size_type col) const {
                return matrix_(row_offset_ + row, col_offset_ + col);
            }

            // Range-checked element access
            reference at(size_type row, size_type col) {
                if (row >= num_rows_ || col >= num_cols_) {
                    throw IndexError("MatrixView::at: index out of range");
                }
                return matrix_(row_offset_ + row, col_offset_ + col);
            }

            const_reference at(size_type row, size_type col) const {
                if (row >= num_rows_ || col >= num_cols_) {
                    throw IndexError("MatrixView::at: index out of range");
                }
                return matrix_(row_offset_ + row, col_offset_ + col);
            }

            // Capacity related
            bool empty() const noexcept {
                return num_rows_ == 0 || num_cols_ == 0;
            }

            size_type rows() const noexcept {
                return num_rows_;
            }

            size_type cols() const noexcept {
                return num_cols_;
            }

            size_type size() const noexcept {
                return num_rows_ * num_cols_;
            }

            // Offset accessors
            size_type row_offset() const noexcept {
                return row_offset_;
            }

            size_type col_offset() const noexcept {
                return col_offset_;
            }

            // Obtain the underlying matrix
            matrix_reference parent() noexcept {
                return matrix_;
            }

            const matrix_reference parent() const noexcept {
                return matrix_;
            }

            // Create a subview
            MatrixView subview(size_type row_offset, size_type col_offset,
                size_type num_rows, size_type num_cols) const {
                if (row_offset + num_rows > num_rows_ || col_offset + num_cols > num_cols_) {
                    throw IndexError("MatrixView::subview: subview range exceeds view bounds");
                }

                return MatrixView(matrix_,
                    row_offset_ + row_offset,
                    col_offset_ + col_offset,
                    num_rows, num_cols);
            }

            // Retrieve a row
            template<typename VectorType>
            void get_row(size_type row_idx, VectorType& dest) const {
                if (row_idx >= num_rows_) {
                    throw IndexError("MatrixView::get_row: row index out of range");
                }

                if (dest.size() != num_cols_) {
                    dest.resize(num_cols_);
                }

                for (size_type j = 0; j < num_cols_; ++j) {
                    dest[j] = matrix_(row_offset_ + row_idx, col_offset_ + j);
                }
            }

            // Retrieve a column
            template<typename VectorType>
            void get_column(size_type col_idx, VectorType& dest) const {
                if (col_idx >= num_cols_) {
                    throw IndexError("MatrixView::get_column: column index out of range");
                }

                if (dest.size() != num_rows_) {
                    dest.resize(num_rows_);
                }

                for (size_type i = 0; i < num_rows_; ++i) {
                    dest[i] = matrix_(row_offset_ + i, col_offset_ + col_idx);
                }
            }

            // Set a row
            template<typename VectorType>
            void set_row(size_type row_idx, const VectorType& src) {
                if (row_idx >= num_rows_) {
                    throw IndexError("MatrixView::set_row: row index out of range");
                }

                if (src.size() != num_cols_) {
                    throw DimensionError("MatrixView::set_row: vector size mismatch");
                }

                for (size_type j = 0; j < num_cols_; ++j) {
                    matrix_(row_offset_ + row_idx, col_offset_ + j) = src[j];
                }
            }

            // Set a column
            template<typename VectorType>
            void set_column(size_type col_idx, const VectorType& src) {
                if (col_idx >= num_cols_) {
                    throw IndexError("MatrixView::set_column: column index out of range");
                }

                if (src.size() != num_rows_) {
                    throw DimensionError("MatrixView::set_column: vector size mismatch");
                }

                for (size_type i = 0; i < num_rows_; ++i) {
                    matrix_(row_offset_ + i, col_offset_ + col_idx) = src[i];
                }
            }

            // Copy a value
            void fill(const value_type& value) {
                for (size_type i = 0; i < num_rows_; ++i) {
                    for (size_type j = 0; j < num_cols_; ++j) {
                        matrix_(row_offset_ + i, col_offset_ + j) = value;
                    }
                }
            }

            // Zero-clear
            void zero() {
                fill(numeric_traits<value_type>::zero());
            }

            // Assignment from another matrix
            template<typename OtherMatrixType>
            void assign(const OtherMatrixType& src) {
                if (src.rows() != num_rows_ || src.cols() != num_cols_) {
                    throw DimensionError("MatrixView::assign: matrix size mismatch");
                }

                for (size_type i = 0; i < num_rows_; ++i) {
                    for (size_type j = 0; j < num_cols_; ++j) {
                        matrix_(row_offset_ + i, col_offset_ + j) = src(i, j);
                    }
                }
            }

            // Whether the matrix is square
            bool is_square() const noexcept {
                return num_rows_ == num_cols_;
            }

        private:
            matrix_reference matrix_;  // Reference to the underlying matrix
            size_type row_offset_;     // Row offset
            size_type col_offset_;     // Column offset
            size_type num_rows_;       // Row count of the view
            size_type num_cols_;       // Column count of the view
        };

        // Read-only matrix view class
        template<typename MatrixType>
        class ConstMatrixView {
        public:
            // Obtain type definitions from the host type
            using value_type = typename std::remove_reference_t<MatrixType>::value_type;
            using size_type = std::size_t;
            using reference = const value_type&;  // const reference
            using const_reference = const value_type&;

            // Reference type to the underlying matrix (always const)
            using matrix_reference = const MatrixType&;

            // Create a view of the entire matrix
            explicit ConstMatrixView(matrix_reference matrix)
                : matrix_(matrix), row_offset_(0), col_offset_(0),
                num_rows_(matrix.rows()), num_cols_(matrix.cols()) {
            }

            // Create a partial view of the matrix
            ConstMatrixView(matrix_reference matrix,
                size_type row_offset, size_type col_offset,
                size_type num_rows, size_type num_cols)
                : matrix_(matrix), row_offset_(row_offset), col_offset_(col_offset),
                num_rows_(num_rows), num_cols_(num_cols) {

                // Range check
                if (row_offset + num_rows > matrix.rows() || col_offset + num_cols > matrix.cols()) {
                    throw IndexError("ConstMatrixView: view range exceeds matrix bounds");
                }
            }

            // Converting constructor from MatrixView
            explicit ConstMatrixView(const MatrixView<MatrixType>& view)
                : matrix_(view.parent()),
                row_offset_(view.row_offset()),
                col_offset_(view.col_offset()),
                num_rows_(view.rows()),
                num_cols_(view.cols()) {
            }

            // Copy constructor
            ConstMatrixView(const ConstMatrixView&) = default;

            // Copy assignment is forbidden (copying the view itself is not allowed)
            ConstMatrixView& operator=(const ConstMatrixView&) = delete;

            // Element access (always const)
            const_reference operator()(size_type row, size_type col) const {
                return matrix_(row_offset_ + row, col_offset_ + col);
            }

            // Range-checked element access (always const)
            const_reference at(size_type row, size_type col) const {
                if (row >= num_rows_ || col >= num_cols_) {
                    throw IndexError("ConstMatrixView::at: index out of range");
                }
                return matrix_(row_offset_ + row, col_offset_ + col);
            }

            // Capacity related
            bool empty() const noexcept {
                return num_rows_ == 0 || num_cols_ == 0;
            }

            size_type rows() const noexcept {
                return num_rows_;
            }

            size_type cols() const noexcept {
                return num_cols_;
            }

            size_type size() const noexcept {
                return num_rows_ * num_cols_;
            }

            // Offset accessors
            size_type row_offset() const noexcept {
                return row_offset_;
            }

            size_type col_offset() const noexcept {
                return col_offset_;
            }

            // Obtain the underlying matrix (always const)
            matrix_reference parent() const noexcept {
                return matrix_;
            }

            // Create a subview
            ConstMatrixView subview(size_type row_offset, size_type col_offset,
                size_type num_rows, size_type num_cols) const {
                if (row_offset + num_rows > num_rows_ || col_offset + num_cols > num_cols_) {
                    throw IndexError("ConstMatrixView::subview: subview range exceeds view bounds");
                }

                return ConstMatrixView(matrix_,
                    row_offset_ + row_offset,
                    col_offset_ + col_offset,
                    num_rows, num_cols);
            }

            // Retrieve a row
            template<typename VectorType>
            void get_row(size_type row_idx, VectorType& dest) const {
                if (row_idx >= num_rows_) {
                    throw IndexError("ConstMatrixView::get_row: row index out of range");
                }

                if (dest.size() != num_cols_) {
                    dest.resize(num_cols_);
                }

                for (size_type j = 0; j < num_cols_; ++j) {
                    dest[j] = matrix_(row_offset_ + row_idx, col_offset_ + j);
                }
            }

            // Retrieve a column
            template<typename VectorType>
            void get_column(size_type col_idx, VectorType& dest) const {
                if (col_idx >= num_cols_) {
                    throw IndexError("ConstMatrixView::get_column: column index out of range");
                }

                if (dest.size() != num_rows_) {
                    dest.resize(num_rows_);
                }

                for (size_type i = 0; i < num_rows_; ++i) {
                    dest[i] = matrix_(row_offset_ + i, col_offset_ + col_idx);
                }
            }

            // Whether the matrix is square
            bool is_square() const noexcept {
                return num_rows_ == num_cols_;
            }

        private:
            matrix_reference matrix_;  // Reference to the underlying matrix (const)
            size_type row_offset_;     // Row offset
            size_type col_offset_;     // Column offset
            size_type num_rows_;       // Row count of the view
            size_type num_cols_;       // Column count of the view
        };

        // Helper functions for creating views
        template<typename MatrixType>
        MatrixView<MatrixType> make_matrix_view(MatrixType& matrix) {
            return MatrixView<MatrixType>(matrix);
        }

        template<typename MatrixType>
        MatrixView<MatrixType> make_matrix_view(MatrixType& matrix,
            std::size_t row_offset, std::size_t col_offset,
            std::size_t num_rows, std::size_t num_cols) {
            return MatrixView<MatrixType>(matrix, row_offset, col_offset, num_rows, num_cols);
        }

        template<typename MatrixType>
        ConstMatrixView<MatrixType> make_const_matrix_view(const MatrixType& matrix) {
            return ConstMatrixView<MatrixType>(matrix);
        }

        template<typename MatrixType>
        ConstMatrixView<MatrixType> make_const_matrix_view(const MatrixType& matrix,
            std::size_t row_offset, std::size_t col_offset,
            std::size_t num_rows, std::size_t num_cols) {
            return ConstMatrixView<MatrixType>(matrix, row_offset, col_offset, num_rows, num_cols);
        }

        // Matrix base class
        template<typename T, typename StoragePolicy>
        class MatrixBase {
        public:
            using value_type = T;
            using size_type = std::size_t;
            using reference = T&;
            using const_reference = const T&;
            using iterator = typename StoragePolicy::iterator;
            using const_iterator = typename StoragePolicy::const_iterator;

            // Default constructor
            MatrixBase() = default;

            // Size-specifying constructor
            explicit MatrixBase(size_type rows, size_type cols) : storage_(rows, cols) {}

            // Constructor specifying size and initial value
            MatrixBase(size_type rows, size_type cols, const T& value) : storage_(rows, cols, value) {}

            // Copy/move constructors and assignments
            MatrixBase(const MatrixBase&) = default;
            MatrixBase(MatrixBase&&) noexcept = default;
            MatrixBase& operator=(const MatrixBase&) = default;
            MatrixBase& operator=(MatrixBase&&) noexcept = default;

            // Destructor
            ~MatrixBase() = default;

            // Element access
            reference operator()(size_type row, size_type col) {
                return storage_(row, col);
            }

            const_reference operator()(size_type row, size_type col) const {
                return storage_(row, col);
            }

            // Range-checked element access
            reference at(size_type row, size_type col) {
                return storage_.at(row, col);
            }

            const_reference at(size_type row, size_type col) const {
                return storage_.at(row, col);
            }

            // Iterators
            iterator begin() noexcept {
                return storage_.begin();
            }

            const_iterator begin() const noexcept {
                return storage_.begin();
            }

            const_iterator cbegin() const noexcept {
                return storage_.cbegin();
            }

            iterator end() noexcept {
                return storage_.end();
            }

            const_iterator end() const noexcept {
                return storage_.end();
            }

            const_iterator cend() const noexcept {
                return storage_.cend();
            }

            // Capacity
            bool empty() const noexcept {
                return storage_.empty();
            }

            size_type rows() const noexcept {
                return storage_.rows();
            }

            size_type cols() const noexcept {
                return storage_.cols();
            }

            size_type size() const noexcept {
                return storage_.size();
            }

            // Resize
            void resize(size_type rows, size_type cols) {
                storage_.resize(rows, cols);
            }

            void resize(size_type rows, size_type cols, const value_type& value) {
                storage_.resize(rows, cols, value);
            }

            // Clear
            void clear() noexcept {
                storage_.clear();
            }

            // Change row/column counts
            void reshape(size_type rows, size_type cols) {
                storage_.reshape(rows, cols);
            }

            // Set value for all elements
            void fill(const T& value) {
                for (size_type i = 0; i < rows(); ++i) {
                    for (size_type j = 0; j < cols(); ++j) {
                        storage_(i, j) = value;
                    }
                }
            }

            // Zero-clear
            void zero() {
                fill(numeric_traits<T>::zero());
            }

            // Set to identity matrix
            void identity() {
                if (!is_square()) {
                    throw DimensionError("identity(): matrix must be square");
                }

                zero();

                const size_type n = rows();
                for (size_type i = 0; i < n; ++i) {
                    storage_(i, i) = numeric_traits<T>::one();
                }
            }

            // Retrieve a row
            template<typename VectorType>
            void get_row(size_type row_idx, VectorType& dest) const {
                if (row_idx >= rows()) {
                    throw IndexError("get_row: row index out of range");
                }

                if (dest.size() != cols()) {
                    dest.resize(cols());
                }

                for (size_type j = 0; j < cols(); ++j) {
                    dest[j] = storage_(row_idx, j);
                }
            }

            // Retrieve a column
            template<typename VectorType>
            void get_column(size_type col_idx, VectorType& dest) const {
                if (col_idx >= cols()) {
                    throw IndexError("get_column: column index out of range");
                }

                if (dest.size() != rows()) {
                    dest.resize(rows());
                }

                for (size_type i = 0; i < rows(); ++i) {
                    dest[i] = storage_(i, col_idx);
                }
            }

            // Set a row
            template<typename VectorType>
            void set_row(size_type row_idx, const VectorType& src) {
                if (row_idx >= rows()) {
                    throw IndexError("set_row: row index out of range");
                }

                if (src.size() != cols()) {
                    throw DimensionError("set_row: vector size mismatch");
                }

                for (size_type j = 0; j < cols(); ++j) {
                    storage_(row_idx, j) = src[j];
                }
            }

            // Set a column
            template<typename VectorType>
            void set_column(size_type col_idx, const VectorType& src) {
                if (col_idx >= cols()) {
                    throw IndexError("set_column: column index out of range");
                }

                if (src.size() != rows()) {
                    throw DimensionError("set_column: vector size mismatch");
                }

                for (size_type i = 0; i < rows(); ++i) {
                    storage_(i, col_idx) = src[i];
                }
            }

            // Obtain the storage (for internal implementation use)
            StoragePolicy& storage() noexcept {
                return storage_;
            }

            const StoragePolicy& storage() const noexcept {
                return storage_;
            }

            // Raw pointer into the internal array (row-major contiguous memory)
            T* data() noexcept {
                return storage_.data();
            }

            const T* data() const noexcept {
                return storage_.data();
            }

            // Whether the matrix is square
            bool is_square() const noexcept {
                return storage_.is_square();
            }

            // Obtain a submatrix view
            MatrixView<MatrixBase> submatrix(size_type row_offset, size_type col_offset,
                size_type num_rows, size_type num_cols) {
                return MatrixView<MatrixBase>(*this, row_offset, col_offset, num_rows, num_cols);
            }

            ConstMatrixView<MatrixBase> submatrix(size_type row_offset, size_type col_offset,
                size_type num_rows, size_type num_cols) const {
                return ConstMatrixView<MatrixBase>(*this, row_offset, col_offset, num_rows, num_cols);
            }

            // Obtain a row view
            MatrixView<MatrixBase> row(size_type row_idx) {
                if (row_idx >= rows()) {
                    throw IndexError("row: row index out of range");
                }
                return MatrixView<MatrixBase>(*this, row_idx, 0, 1, cols());
            }

            ConstMatrixView<MatrixBase> row(size_type row_idx) const {
                if (row_idx >= rows()) {
                    throw IndexError("row: row index out of range");
                }
                return ConstMatrixView<MatrixBase>(*this, row_idx, 0, 1, cols());
            }

            // Obtain a column view
            MatrixView<MatrixBase> column(size_type col_idx) {
                if (col_idx >= cols()) {
                    throw IndexError("column: column index out of range");
                }
                return MatrixView<MatrixBase>(*this, 0, col_idx, rows(), 1);
            }

            ConstMatrixView<MatrixBase> column(size_type col_idx) const {
                if (col_idx >= cols()) {
                    throw IndexError("column: column index out of range");
                }
                return ConstMatrixView<MatrixBase>(*this, 0, col_idx, rows(), 1);
            }

            // Block operations API (Eigen-style)

            /// block(i, j, p, q): alias for submatrix
            MatrixView<MatrixBase> block(size_type i, size_type j, size_type p, size_type q) {
                return submatrix(i, j, p, q);
            }
            ConstMatrixView<MatrixBase> block(size_type i, size_type j, size_type p, size_type q) const {
                return submatrix(i, j, p, q);
            }

            /// Top-left corner: p rows and q columns starting at (0,0)
            MatrixView<MatrixBase> topLeftCorner(size_type p, size_type q) {
                return submatrix(0, 0, p, q);
            }
            ConstMatrixView<MatrixBase> topLeftCorner(size_type p, size_type q) const {
                return submatrix(0, 0, p, q);
            }

            /// Top-right corner
            MatrixView<MatrixBase> topRightCorner(size_type p, size_type q) {
                return submatrix(0, cols() - q, p, q);
            }
            ConstMatrixView<MatrixBase> topRightCorner(size_type p, size_type q) const {
                return submatrix(0, cols() - q, p, q);
            }

            /// Bottom-left corner
            MatrixView<MatrixBase> bottomLeftCorner(size_type p, size_type q) {
                return submatrix(rows() - p, 0, p, q);
            }
            ConstMatrixView<MatrixBase> bottomLeftCorner(size_type p, size_type q) const {
                return submatrix(rows() - p, 0, p, q);
            }

            /// Bottom-right corner
            MatrixView<MatrixBase> bottomRightCorner(size_type p, size_type q) {
                return submatrix(rows() - p, cols() - q, p, q);
            }
            ConstMatrixView<MatrixBase> bottomRightCorner(size_type p, size_type q) const {
                return submatrix(rows() - p, cols() - q, p, q);
            }

            /// Top n rows
            MatrixView<MatrixBase> topRows(size_type n) {
                return submatrix(0, 0, n, cols());
            }
            ConstMatrixView<MatrixBase> topRows(size_type n) const {
                return submatrix(0, 0, n, cols());
            }

            /// Bottom n rows
            MatrixView<MatrixBase> bottomRows(size_type n) {
                return submatrix(rows() - n, 0, n, cols());
            }
            ConstMatrixView<MatrixBase> bottomRows(size_type n) const {
                return submatrix(rows() - n, 0, n, cols());
            }

            /// Left n columns
            MatrixView<MatrixBase> leftCols(size_type n) {
                return submatrix(0, 0, rows(), n);
            }
            ConstMatrixView<MatrixBase> leftCols(size_type n) const {
                return submatrix(0, 0, rows(), n);
            }

            /// Right n columns
            MatrixView<MatrixBase> rightCols(size_type n) {
                return submatrix(0, cols() - n, rows(), n);
            }
            ConstMatrixView<MatrixBase> rightCols(size_type n) const {
                return submatrix(0, cols() - n, rows(), n);
            }

        protected:
            StoragePolicy storage_;
        };

    } // namespace detail

    //-----------------------------------------------------------------------------
    // Triangular / symmetric matrix views
    //-----------------------------------------------------------------------------

    /// Kind of triangular matrix
    enum class TriangularMode { Upper, Lower, UnitUpper, UnitLower, StrictlyUpper, StrictlyLower };

    /**
     * @brief TriangularView — triangular matrix view (Eigen-compatible)
     *
     * References the upper/lower triangle of a matrix without copying and performs triangular operations.
     * - solve(b): solves the triangular linear system Ax=b by forward/back substitution
     * - operator*(v): triangular matrix-vector product (skips zero entries)
     * - toDense(): produces a dense matrix that copies only the triangular part
     */
    template<typename MatrixType, TriangularMode Mode>
    class TriangularView {
    public:
        using T = typename std::remove_cvref_t<MatrixType>::value_type;
        using size_type = std::size_t;

        static constexpr bool is_upper = (Mode == TriangularMode::Upper ||
                                          Mode == TriangularMode::UnitUpper ||
                                          Mode == TriangularMode::StrictlyUpper);
        static constexpr bool is_unit  = (Mode == TriangularMode::UnitUpper ||
                                          Mode == TriangularMode::UnitLower);
        static constexpr bool is_strict = (Mode == TriangularMode::StrictlyUpper ||
                                           Mode == TriangularMode::StrictlyLower);

        explicit TriangularView(MatrixType& mat) : mat_(mat) {
            if (mat_.rows() != mat_.cols())
                throw DimensionError("TriangularView: matrix must be square");
        }

        size_type rows() const { return mat_.rows(); }
        size_type cols() const { return mat_.cols(); }

        /// Element access (zero outside the triangular region; unit diagonal returns 1)
        T operator()(size_type i, size_type j) const {
            if constexpr (is_upper) {
                if (i > j) return T(0);
                if (is_strict && i == j) return T(0);
                if (is_unit && i == j) return T(1);
            } else {
                if (i < j) return T(0);
                if (is_strict && i == j) return T(0);
                if (is_unit && i == j) return T(1);
            }
            return mat_(i, j);
        }

        /// Solve the triangular linear system: this * x = b -> x
        /// AVX2 FMA dot product + row-wise back/forward substitution
        Vector<T> solve(const Vector<T>& b) const {
            const auto n = rows();
            if (b.size() != n)
                throw DimensionError("TriangularView::solve: dimension mismatch");
            Vector<T> x = b;
            T* __restrict xp = x.data();
            const T* __restrict A = mat_.data();
            using PT = PacketTraits<T>;
            constexpr size_type W = PT::size;  // AVX2 double: 4

            if constexpr (is_upper) {
                for (size_type ii = 0; ii < n; ++ii) {
                    size_type i = n - 1 - ii;
                    const T* __restrict row = A + i * n;
                    size_type j = i + 1;
                    size_type len = n - j;

                    if constexpr (W > 1) {
                        auto acc0 = PT::set1(T(0));
                        auto acc1 = PT::set1(T(0));
                        size_type end2w = j + (len & ~(2 * W - 1));
                        for (; j < end2w; j += 2 * W) {
                            acc0 = PT::fmadd(PT::load(row + j),     PT::load(xp + j),     acc0);
                            acc1 = PT::fmadd(PT::load(row + j + W), PT::load(xp + j + W), acc1);
                        }
                        size_type endw = i + 1 + (len & ~(W - 1));
                        for (; j < endw; j += W)
                            acc0 = PT::fmadd(PT::load(row + j), PT::load(xp + j), acc0);
                        T sum = PT::reduce_add(PT::add(acc0, acc1));
                        for (; j < n; ++j) sum += row[j] * xp[j];
                        xp[i] -= sum;
                    } else {
                        T sum = T(0);
                        for (; j < n; ++j) sum += row[j] * xp[j];
                        xp[i] -= sum;
                    }
                    if constexpr (!is_unit && !is_strict)
                        xp[i] /= row[i];
                }
            } else {
                for (size_type i = 0; i < n; ++i) {
                    const T* __restrict row = A + i * n;
                    size_type j = 0;

                    if constexpr (W > 1) {
                        auto acc0 = PT::set1(T(0));
                        auto acc1 = PT::set1(T(0));
                        size_type end2w = (i / (2 * W)) * (2 * W);
                        for (; j < end2w; j += 2 * W) {
                            acc0 = PT::fmadd(PT::load(row + j),     PT::load(xp + j),     acc0);
                            acc1 = PT::fmadd(PT::load(row + j + W), PT::load(xp + j + W), acc1);
                        }
                        size_type endw = (i / W) * W;
                        for (; j < endw; j += W)
                            acc0 = PT::fmadd(PT::load(row + j), PT::load(xp + j), acc0);
                        T sum = PT::reduce_add(PT::add(acc0, acc1));
                        for (; j < i; ++j) sum += row[j] * xp[j];
                        xp[i] -= sum;
                    } else {
                        T sum = T(0);
                        for (; j < i; ++j) sum += row[j] * xp[j];
                        xp[i] -= sum;
                    }
                    if constexpr (!is_unit && !is_strict)
                        xp[i] /= row[i];
                }
            }
            return x;
        }

        /// Matrix-vector product: y = this * x (triangular part only)
        Vector<T> operator*(const Vector<T>& x) const {
            const auto n = rows();
            if (x.size() != n)
                throw DimensionError("TriangularView::operator*: dimension mismatch");
            Vector<T> y(n, T(0));
            for (size_type i = 0; i < n; ++i) {
                size_type jstart = is_upper ? i : 0;
                size_type jend   = is_upper ? n : i + 1;
                if constexpr (is_strict) {
                    if (is_upper) ++jstart; else --jend;
                }
                for (size_type j = jstart; j < jend; ++j) {
                    if constexpr (is_unit) {
                        y[i] += (i == j ? T(1) : mat_(i, j)) * x[j];
                    } else {
                        y[i] += mat_(i, j) * x[j];
                    }
                }
            }
            return y;
        }

        /// Produce a dense matrix containing only the triangular part
        Matrix<T> toDense() const {
            const auto n = rows();
            Matrix<T> result(n, n, T(0));
            for (size_type i = 0; i < n; ++i)
                for (size_type j = 0; j < n; ++j)
                    result(i, j) = (*this)(i, j);
            return result;
        }

        /// Write the triangular part back into the underlying matrix (e.g., for constructing symmetric matrices)
        template<typename SrcMatrix>
        void setFrom(const SrcMatrix& src) requires (!std::is_const_v<MatrixType>) {
            const auto n = rows();
            for (size_type i = 0; i < n; ++i) {
                size_type jstart = is_upper ? i : 0;
                size_type jend   = is_upper ? n : i + 1;
                for (size_type j = jstart; j < jend; ++j)
                    mat_(i, j) = src(i, j);
            }
        }

    private:
        MatrixType& mat_;
    };

    /**
     * @brief SelfAdjointView — symmetric matrix view (Eigen-compatible)
     *
     * References the full symmetric matrix using only the upper or lower triangle.
     * Elements on the unstored side are returned via transposition.
     */
    template<typename MatrixType, TriangularMode Mode>
    class SelfAdjointView {
        static_assert(Mode == TriangularMode::Upper || Mode == TriangularMode::Lower,
                      "SelfAdjointView requires Upper or Lower mode");
    public:
        using T = typename std::remove_cvref_t<MatrixType>::value_type;
        using size_type = std::size_t;

        static constexpr bool from_upper = (Mode == TriangularMode::Upper);

        explicit SelfAdjointView(MatrixType& mat) : mat_(mat) {
            if (mat_.rows() != mat_.cols())
                throw DimensionError("SelfAdjointView: matrix must be square");
        }

        size_type rows() const { return mat_.rows(); }
        size_type cols() const { return mat_.cols(); }

        /// Element access (the unstored side returns the transposed entry)
        T operator()(size_type i, size_type j) const {
            if constexpr (from_upper) {
                return (i <= j) ? mat_(i, j) : mat_(j, i);
            } else {
                return (i >= j) ? mat_(i, j) : mat_(j, i);
            }
        }

        /// Symmetric matrix-vector product: y = this * x
        /// Simultaneous AVX2 FMA dot product + axpy
        Vector<T> operator*(const Vector<T>& x) const {
            const auto n = rows();
            if (x.size() != n)
                throw DimensionError("SelfAdjointView::operator*: dimension mismatch");
            Vector<T> y(n, T(0));
            const T* __restrict A = mat_.data();
            const T* __restrict xp = x.data();
            T* __restrict yp = y.data();
            using PT = PacketTraits<T>;
            constexpr size_type W = PT::size;

            if constexpr (from_upper) {
                for (size_type i = 0; i < n; ++i) {
                    const T* __restrict row = A + i * n;
                    T xi = xp[i];
                    size_type j = i + 1;
                    size_type len = n - j;

                    if constexpr (W > 1) {
                        auto vsum = PT::set1(T(0));
                        auto vxi = PT::set1(xi);
                        size_type endw = j + (len & ~(W - 1));
                        for (; j < endw; j += W) {
                            auto va = PT::load(row + j);
                            vsum = PT::fmadd(va, PT::load(xp + j), vsum);
                            PT::store(yp + j, PT::fmadd(va, vxi, PT::load(yp + j)));
                        }
                        T sum = row[i] * xi + PT::reduce_add(vsum);
                        for (; j < n; ++j) {
                            T aij = row[j];
                            sum += aij * xp[j];
                            yp[j] += aij * xi;
                        }
                        yp[i] += sum;
                    } else {
                        T sum = row[i] * xi;
                        for (; j < n; ++j) {
                            T aij = row[j];
                            sum += aij * xp[j];
                            yp[j] += aij * xi;
                        }
                        yp[i] += sum;
                    }
                }
            } else {
                for (size_type i = 0; i < n; ++i) {
                    const T* __restrict row = A + i * n;
                    T xi = xp[i];
                    size_type j = 0;

                    if constexpr (W > 1) {
                        auto vsum = PT::set1(T(0));
                        auto vxi = PT::set1(xi);
                        size_type endw = (i / W) * W;
                        for (; j < endw; j += W) {
                            auto va = PT::load(row + j);
                            vsum = PT::fmadd(va, PT::load(xp + j), vsum);
                            PT::store(yp + j, PT::fmadd(va, vxi, PT::load(yp + j)));
                        }
                        T sum = row[i] * xi + PT::reduce_add(vsum);
                        for (; j < i; ++j) {
                            T aij = row[j];
                            sum += aij * xp[j];
                            yp[j] += aij * xi;
                        }
                        yp[i] += sum;
                    } else {
                        T sum = row[i] * xi;
                        for (; j < i; ++j) {
                            T aij = row[j];
                            sum += aij * xp[j];
                            yp[j] += aij * xi;
                        }
                        yp[i] += sum;
                    }
                }
            }
            return y;
        }

        /// Produce a dense copy of the symmetric matrix
        Matrix<T> toDense() const {
            const auto n = rows();
            Matrix<T> result(n, n);
            for (size_type i = 0; i < n; ++i)
                for (size_type j = 0; j < n; ++j)
                    result(i, j) = (*this)(i, j);
            return result;
        }

    private:
        MatrixType& mat_;
    };

    //-----------------------------------------------------------------------------
    // Matrix operation functions
    //-----------------------------------------------------------------------------

    // Matrix addition: C = A + B
    template<typename MatrixA, typename MatrixB, typename MatrixC,
        typename ComputePolicy = DefaultComputePolicy<typename MatrixA::value_type>>
        void add(const MatrixA& a, const MatrixB& b, MatrixC& c) {
        // Dimension check
        if (a.rows() != b.rows() || a.cols() != b.cols() ||
            a.rows() != c.rows() || a.cols() != c.cols()) {
            throw DimensionError("Matrix addition: dimension mismatch");
        }

        // Delegate the work to the compute policy
        ComputePolicy::add(a, b, c);
    }

    // Matrix subtraction: C = A - B
    template<typename MatrixA, typename MatrixB, typename MatrixC,
        typename ComputePolicy = DefaultComputePolicy<typename MatrixA::value_type>>
        void subtract(const MatrixA& a, const MatrixB& b, MatrixC& c) {
        using size_type = typename MatrixA::size_type;

        // Dimension check
        if (a.rows() != b.rows() || a.cols() != b.cols() ||
            a.rows() != c.rows() || a.cols() != c.cols()) {
            throw DimensionError("Matrix subtraction: dimension mismatch");
        }

        const size_type rows = a.rows();
        const size_type cols = a.cols();

        for (size_type i = 0; i < rows; ++i) {
            for (size_type j = 0; j < cols; ++j) {
                c(i, j) = a(i, j) - b(i, j);
            }
        }
    }

    // Matrix-scalar multiplication: B = alpha * A
    template<typename MatrixA, typename Scalar, typename MatrixB,
        typename ComputePolicy = DefaultComputePolicy<typename MatrixA::value_type>>
        void scale(const MatrixA& a, const Scalar& alpha, MatrixB& b) {
        using size_type = typename MatrixA::size_type;

        // Dimension check
        if (a.rows() != b.rows() || a.cols() != b.cols()) {
            throw DimensionError("Matrix scaling: dimension mismatch");
        }

        const size_type rows = a.rows();
        const size_type cols = a.cols();

        for (size_type i = 0; i < rows; ++i) {
            for (size_type j = 0; j < cols; ++j) {
                b(i, j) = a(i, j) * alpha;
            }
        }
    }

    // Matrix-matrix multiplication: C = A * B
    template<typename MatrixA, typename MatrixB, typename MatrixC,
        typename ComputePolicy = DefaultComputePolicy<typename MatrixA::value_type>>
        void multiply(const MatrixA& a, const MatrixB& b, MatrixC& c) {
        // Dimension check
        if (a.cols() != b.rows() ||
            c.rows() != a.rows() || c.cols() != b.cols()) {
            throw DimensionError("Matrix multiplication: dimension mismatch");
        }

        // Delegate the work to the compute policy
        ComputePolicy::multiply(a, b, c);
    }

    // fused GEMM: C = alpha * A * B + beta * C
    template<typename MatrixA, typename MatrixB, typename MatrixC,
        typename ComputePolicy = DefaultComputePolicy<typename MatrixA::value_type>>
        void gemm(const MatrixA& a, const MatrixB& b, MatrixC& c,
                  typename MatrixA::value_type alpha = typename MatrixA::value_type{1},
                  typename MatrixA::value_type beta = typename MatrixA::value_type{0}) {
        if (a.cols() != b.rows() ||
            c.rows() != a.rows() || c.cols() != b.cols()) {
            throw DimensionError("GEMM: dimension mismatch");
        }

        ComputePolicy::gemm(a, b, c, alpha, beta);
    }

    // B-transposed multiply: C = A * B^T  (B is not physically transposed)
    // A: M x K, B: N x K (pre-transpose), C: M x N
    // Uses the MKL TransB=true flag, or falls back to a naive implementation
    template<typename MatrixA, typename MatrixB, typename MatrixC,
        typename ComputePolicy = DefaultComputePolicy<typename MatrixA::value_type>>
    void multiply_trans_b(const MatrixA& a, const MatrixB& b, MatrixC& c) {
        if (a.cols() != b.cols() || c.rows() != a.rows() || c.cols() != b.rows()) {
            throw DimensionError("multiply_trans_b: dimension mismatch");
        }
        ComputePolicy::multiply_trans_b(a, b, c);
    }

    // A-transposed multiply: C = A^T * B  (A is not physically transposed)
    // A: K x M, B: K x N, C: M x N
    template<typename MatrixA, typename MatrixB, typename MatrixC,
        typename ComputePolicy = DefaultComputePolicy<typename MatrixA::value_type>>
    void multiply_trans_a(const MatrixA& a, const MatrixB& b, MatrixC& c) {
        if (a.rows() != b.rows() || c.rows() != a.cols() || c.cols() != b.cols()) {
            throw DimensionError("multiply_trans_a: dimension mismatch");
        }
        ComputePolicy::multiply_trans_a(a, b, c);
    }

    // B-transposed multiply + accumulate: C = alpha * A * B^T + beta * C
    // Intended to fuse Linear's bias broadcast + GEMM into a single cblas call
    template<typename MatrixA, typename MatrixB, typename MatrixC,
        typename ComputePolicy = DefaultComputePolicy<typename MatrixA::value_type>>
    void multiply_trans_b_accumulate(const MatrixA& a, const MatrixB& b, MatrixC& c,
                                      typename MatrixA::value_type alpha,
                                      typename MatrixA::value_type beta) {
        if (a.cols() != b.cols() || c.rows() != a.rows() || c.cols() != b.rows()) {
            throw DimensionError("multiply_trans_b_accumulate: dimension mismatch");
        }
        ComputePolicy::multiply_trans_b_accumulate(a, b, c, alpha, beta);
    }

    // Matrix-vector multiplication: y = A * x
    // Small matrices (m<=256): process two rows at once to reuse the x vector in cache
    // Large matrices: 4-way unrolled FMA
    template<typename MatrixA, typename VectorX, typename VectorY,
        typename ComputePolicy = DefaultComputePolicy<typename MatrixA::value_type>>
        void multiply_vector(const MatrixA& a, const VectorX& x, VectorY& y) {
        using T = typename MatrixA::value_type;
        using size_type = typename MatrixA::size_type;

        if (a.cols() != x.size() || y.size() != a.rows()) {
            throw DimensionError("Matrix-vector multiplication: dimension mismatch");
        }

        const size_type m = a.rows();
        const size_type n = a.cols();
        const T* a_data = a.data();
        const T* x_data = x.data();
        const size_type lda = a.cols();

        if constexpr (has_simd_packet_v<T>) {
            using PT = PacketTraits<T>;
            constexpr size_type PS = PT::size;
            constexpr size_type STRIDE = PS * 4;
            const size_type unroll_end = n - (n % STRIDE);
            const size_type vec_end = n - (n % PS);

            // Small-matrix path: process two rows at once (reuse the x vector in L1 cache)
            size_type i = 0;
            if (m <= 256) {
                for (; i + 1 < m; i += 2) {
                    const T* row0 = &a_data[i * lda];
                    const T* row1 = &a_data[(i + 1) * lda];
                    auto a0 = PT::set1(T{0}), a1 = a0, b0 = a0, b1 = a0;
                    auto c0 = a0, c1 = a0, d0 = a0, d1 = a0;
                    size_type j = 0;
                    for (; j < unroll_end; j += STRIDE) {
                        auto x0 = PT::load(&x_data[j]);
                        auto x1 = PT::load(&x_data[j + PS]);
                        auto x2 = PT::load(&x_data[j + PS * 2]);
                        auto x3 = PT::load(&x_data[j + PS * 3]);
                        a0 = PT::fmadd(PT::load(&row0[j]),          x0, a0);
                        a1 = PT::fmadd(PT::load(&row0[j + PS]),     x1, a1);
                        b0 = PT::fmadd(PT::load(&row0[j + PS * 2]), x2, b0);
                        b1 = PT::fmadd(PT::load(&row0[j + PS * 3]), x3, b1);
                        c0 = PT::fmadd(PT::load(&row1[j]),          x0, c0);
                        c1 = PT::fmadd(PT::load(&row1[j + PS]),     x1, c1);
                        d0 = PT::fmadd(PT::load(&row1[j + PS * 2]), x2, d0);
                        d1 = PT::fmadd(PT::load(&row1[j + PS * 3]), x3, d1);
                    }
                    auto sum0_p = PT::add(PT::add(a0, a1), PT::add(b0, b1));
                    auto sum1_p = PT::add(PT::add(c0, c1), PT::add(d0, d1));
                    // Add the tail into the packet accumulators (reduce_add is done only once at the end)
                    for (; j < vec_end; j += PS) {
                        auto xp = PT::load(&x_data[j]);
                        sum0_p = PT::fmadd(PT::load(&row0[j]), xp, sum0_p);
                        sum1_p = PT::fmadd(PT::load(&row1[j]), xp, sum1_p);
                    }
                    T s0 = PT::reduce_add(sum0_p);
                    T s1 = PT::reduce_add(sum1_p);
                    for (; j < n; ++j) {
                        s0 += row0[j] * x_data[j];
                        s1 += row1[j] * x_data[j];
                    }
                    y[i] = s0;
                    y[i + 1] = s1;
                }
            }

            // Remaining (odd row count or large matrix)
            for (; i < m; ++i) {
                const T* a_row = &a_data[i * lda];
                auto acc0 = PT::set1(T{0}), acc1 = acc0, acc2 = acc0, acc3 = acc0;
                size_type j = 0;
                for (; j < unroll_end; j += STRIDE) {
                    acc0 = PT::fmadd(PT::load(&a_row[j]),          PT::load(&x_data[j]),          acc0);
                    acc1 = PT::fmadd(PT::load(&a_row[j + PS]),     PT::load(&x_data[j + PS]),     acc1);
                    acc2 = PT::fmadd(PT::load(&a_row[j + PS * 2]), PT::load(&x_data[j + PS * 2]), acc2);
                    acc3 = PT::fmadd(PT::load(&a_row[j + PS * 3]), PT::load(&x_data[j + PS * 3]), acc3);
                }
                auto acc = PT::add(PT::add(acc0, acc1), PT::add(acc2, acc3));
                // Also add the tail into the packet accumulators
                for (; j < vec_end; j += PS)
                    acc = PT::fmadd(PT::load(&a_row[j]), PT::load(&x_data[j]), acc);
                T sum = PT::reduce_add(acc);
                for (; j < n; ++j)
                    sum += a_row[j] * x_data[j];
                y[i] = sum;
            }
        } else {
            for (size_type i = 0; i < m; ++i) {
                T sum = T{0};
                for (size_type j = 0; j < n; ++j)
                    sum += a(i, j) * x[j];
                y[i] = sum;
            }
        }
    }

    // Matrix transpose: B = A^T (delegates to computation_policy's optimized transpose)
    template<typename MatrixA, typename MatrixB,
        typename ComputePolicy = DefaultComputePolicy<typename MatrixA::value_type>>
        void transpose(const MatrixA& a, MatrixB& b) {
        ComputePolicy::transpose(a, b);
    }

    // Extract the diagonal of a matrix: diag = diag(A)
    template<typename MatrixA, typename VectorDiag>
    void extract_diagonal(const MatrixA& a, VectorDiag& diag) {
        using size_type = typename MatrixA::size_type;

        const size_type min_dim = std::min(a.rows(), a.cols());

        // Adjust size
        if (diag.size() != min_dim) {
            diag.resize(min_dim);
        }

        // Extract the diagonal entries
        for (size_type i = 0; i < min_dim; ++i) {
            diag[i] = a(i, i);
        }
    }

    // Construct a diagonal matrix: A = diag(v)
    template<typename VectorDiag, typename MatrixA>
    void diagonal_matrix(const VectorDiag& diag, MatrixA& a) {
        using size_type = typename MatrixA::size_type;

        const size_type n = diag.size();

        // Dimension check
        if (a.rows() != n || a.cols() != n) {
            throw DimensionError("diagonal_matrix: dimension mismatch");
        }

        // Zero-initialize
        a.zero();

        // Set the diagonal entries
        for (size_type i = 0; i < n; ++i) {
            a(i, i) = diag[i];
        }
    }

    // Trace computation: tr(A)
    template<typename MatrixA>
    auto trace(const MatrixA& a) -> typename MatrixA::value_type {
        using value_type = typename MatrixA::value_type;
        using size_type = typename MatrixA::size_type;

        if (!a.is_square()) {
            throw DimensionError("trace: matrix must be square");
        }

        const size_type n = a.rows();
        value_type result = numeric_traits<value_type>::zero();

        for (size_type i = 0; i < n; ++i) {
            result += a(i, i);
        }

        return result;
    }

    // Frobenius norm: ||A||_F
    template<typename MatrixA>
    auto frobenius_norm(const MatrixA& a) -> typename MatrixA::value_type {
        using value_type = typename MatrixA::value_type;
        using size_type = typename MatrixA::size_type;

        const size_type rows = a.rows();
        const size_type cols = a.cols();

        value_type sum_squares = numeric_traits<value_type>::zero();

        for (size_type i = 0; i < rows; ++i) {
            for (size_type j = 0; j < cols; ++j) {
                sum_squares += a(i, j) * a(i, j);
            }
        }

        if constexpr (std::is_floating_point_v<value_type>) {
            return static_cast<value_type>(std::sqrt(static_cast<double>(sum_squares)));
        } else {
            // Multi-precision Float etc.: sqrt directly on value_type without going through double (precision preserved)
            return detail::generic_sqrt(sum_squares);
        }
    }

    // Coefficient-wise reductions (Eigen parity). Constrained to matrix-like
    // types so vector overloads of the same names are never shadowed.
    template<typename MatrixA>
        requires requires(const MatrixA& m) { m.rows(); m.cols(); m(0, 0); }
    auto sum(const MatrixA& a) -> typename MatrixA::value_type {
        using value_type = typename MatrixA::value_type;
        using size_type = typename MatrixA::size_type;
        value_type result = numeric_traits<value_type>::zero();
        for (size_type i = 0; i < a.rows(); ++i)
            for (size_type j = 0; j < a.cols(); ++j)
                result += a(i, j);
        return result;
    }

    template<typename MatrixA>
        requires requires(const MatrixA& m) { m.rows(); m.cols(); m(0, 0); }
    auto mean(const MatrixA& a) -> typename MatrixA::value_type {
        using value_type = typename MatrixA::value_type;
        const auto count = a.rows() * a.cols();
        if (count == 0) throw DimensionError("mean: empty matrix");
        return sum(a) / static_cast<value_type>(count);
    }

    // minCoeff/maxCoeff with optional (row, col) of the extremum (real types only).
    template<typename MatrixA>
        requires requires(const MatrixA& m) { m.rows(); m.cols(); m(0, 0); }
    auto minCoeff(const MatrixA& a, typename MatrixA::size_type* row = nullptr,
                  typename MatrixA::size_type* col = nullptr)
        -> typename MatrixA::value_type {
        using value_type = typename MatrixA::value_type;
        using size_type = typename MatrixA::size_type;
        static_assert(!numeric_traits<value_type>::is_complex,
            "minCoeff requires an ordered (non-complex) scalar type");
        if (a.rows() == 0 || a.cols() == 0)
            throw DimensionError("minCoeff: empty matrix");
        value_type best = a(0, 0);
        size_type br = 0, bc = 0;
        for (size_type i = 0; i < a.rows(); ++i)
            for (size_type j = 0; j < a.cols(); ++j)
                if (a(i, j) < best) { best = a(i, j); br = i; bc = j; }
        if (row) *row = br;
        if (col) *col = bc;
        return best;
    }

    template<typename MatrixA>
        requires requires(const MatrixA& m) { m.rows(); m.cols(); m(0, 0); }
    auto maxCoeff(const MatrixA& a, typename MatrixA::size_type* row = nullptr,
                  typename MatrixA::size_type* col = nullptr)
        -> typename MatrixA::value_type {
        using value_type = typename MatrixA::value_type;
        using size_type = typename MatrixA::size_type;
        static_assert(!numeric_traits<value_type>::is_complex,
            "maxCoeff requires an ordered (non-complex) scalar type");
        if (a.rows() == 0 || a.cols() == 0)
            throw DimensionError("maxCoeff: empty matrix");
        value_type best = a(0, 0);
        size_type br = 0, bc = 0;
        for (size_type i = 0; i < a.rows(); ++i)
            for (size_type j = 0; j < a.cols(); ++j)
                if (a(i, j) > best) { best = a(i, j); br = i; bc = j; }
        if (row) *row = br;
        if (col) *col = bc;
        return best;
    }

    // L1 matrix norm (maximum column sum): max_j Sum_i |a(i,j)|
    template<typename MatrixA>
    auto matrix_norm_l1(const MatrixA& a) -> typename MatrixA::value_type {
        using value_type = typename MatrixA::value_type;
        using size_type = typename MatrixA::size_type;

        const size_type rows = a.rows();
        const size_type cols = a.cols();

        value_type max_col_sum = numeric_traits<value_type>::zero();

        for (size_type j = 0; j < cols; ++j) {
            value_type col_sum = numeric_traits<value_type>::zero();
            for (size_type i = 0; i < rows; ++i) {
                col_sum += detail::generic_abs(a(i, j));
            }
            if (col_sum > max_col_sum) max_col_sum = col_sum;
        }

        return max_col_sum;
    }

    // L-infinity matrix norm (maximum row sum): max_i Sum_j |a(i,j)|
    template<typename MatrixA>
    auto matrix_norm_linf(const MatrixA& a) -> typename MatrixA::value_type {
        using value_type = typename MatrixA::value_type;
        using size_type = typename MatrixA::size_type;

        const size_type rows = a.rows();
        const size_type cols = a.cols();

        value_type max_row_sum = numeric_traits<value_type>::zero();

        for (size_type i = 0; i < rows; ++i) {
            value_type row_sum = numeric_traits<value_type>::zero();
            for (size_type j = 0; j < cols; ++j) {
                row_sum += detail::generic_abs(a(i, j));
            }
            if (row_sum > max_row_sum) max_row_sum = row_sum;
        }

        return max_row_sum;
    }

    // Determinant computation (basic implementation, intended for small matrices)
    template<typename MatrixA>
    auto determinant(const MatrixA& a) -> typename MatrixA::value_type {
        using value_type = typename MatrixA::value_type;
        using size_type = typename MatrixA::size_type;

        if (!a.is_square()) {
            throw DimensionError("determinant: matrix must be square");
        }

        const size_type n = a.rows();

        // 1x1 matrix
        if (n == 1) {
            return a(0, 0);
        }

        // 2x2 matrix
        if (n == 2) {
            return a(0, 0) * a(1, 1) - a(0, 1) * a(1, 0);
        }

        // 3x3 matrix
        if (n == 3) {
            return a(0, 0) * (a(1, 1) * a(2, 2) - a(1, 2) * a(2, 1)) -
                a(0, 1) * (a(1, 0) * a(2, 2) - a(1, 2) * a(2, 0)) +
                a(0, 2) * (a(1, 0) * a(2, 1) - a(1, 1) * a(2, 0));
        }

        // 4x4 and larger: select algorithm based on type
        if constexpr (numeric_traits<value_type>::is_integer) {
            // Integer type: Bareiss algorithm (computes the determinant without division)
            // M[i][j] = (M[i][j]*M[k][k] - M[i][k]*M[k][j]) / prev_pivot
            // This division is guaranteed to be exact (by Sylvester's identity).
            Matrix<value_type> M(n, n);
            for (size_type i = 0; i < n; ++i)
                for (size_type j = 0; j < n; ++j)
                    M(i, j) = a(i, j);
            int sign = 1;
            value_type prev_pivot = numeric_traits<value_type>::one();

            for (size_type k = 0; k < n; ++k) {
                // Partial pivoting
                size_type pivot_row = k;
                for (size_type i = k + 1; i < n; ++i) {
                    if (M(i, k) != numeric_traits<value_type>::zero() &&
                        (M(pivot_row, k) == numeric_traits<value_type>::zero() ||
                         numeric_traits<value_type>::pivotBetter(M(i, k), M(pivot_row, k))))
                        pivot_row = i;
                }

                if (pivot_row != k) {
                    for (size_type j = 0; j < n; ++j)
                        std::swap(M(k, j), M(pivot_row, j));
                    sign = -sign;
                }

                if (M(k, k) == numeric_traits<value_type>::zero())
                    return numeric_traits<value_type>::zero();

                for (size_type i = k + 1; i < n; ++i) {
                    for (size_type j = k + 1; j < n; ++j) {
                        // The Bareiss division divides exactly by the Sylvester identity. We verify
                        // this by multiplying back (same convention as PolyMatrix). Non-divisibility / overflow corruption is caught in debug
                        // (because intermediate-product overflow of fixed-width integers silently produces wrong values). When an exact
                        // fail-closed determinant is required, also use numeric_bridge::certifyDeterminant(LU).
                        const value_type num = M(i, j) * M(k, k) - M(i, k) * M(k, j);
                        const value_type quot = num / prev_pivot;
                        assert(quot * prev_pivot == num &&
                               "matrix determinant (Bareiss): non-exact division "
                               "(integer overflow or non-exact ring)");
                        M(i, j) = quot;
                    }
                }
                // Zero out column k for rows i > k (not required for subsequent steps, but explicit)
                for (size_type i = k + 1; i < n; ++i)
                    M(i, k) = numeric_traits<value_type>::zero();

                prev_pivot = M(k, k);
            }
            return M(n - 1, n - 1) * value_type(sign);
        } else {
            // Floating-point / rational: compute determinant in O(n^3) via Gaussian elimination
            Matrix<value_type> lu(n, n);
            for (size_type i = 0; i < n; ++i)
                for (size_type j = 0; j < n; ++j)
                    lu(i, j) = a(i, j);
            value_type det = numeric_traits<value_type>::one();
            int sign = 1;

            for (size_type k = 0; k < n; ++k) {
                // Partial pivoting (uses pivotBetter to choose the best criterion for the type)
                size_type pivot_row = k;
                for (size_type i = k + 1; i < n; ++i) {
                    if (lu(i, k) != numeric_traits<value_type>::zero() &&
                        (lu(pivot_row, k) == numeric_traits<value_type>::zero() ||
                         numeric_traits<value_type>::pivotBetter(lu(i, k), lu(pivot_row, k))))
                        pivot_row = i;
                }

                if (pivot_row != k) {
                    for (size_type j = 0; j < n; ++j)
                        std::swap(lu(k, j), lu(pivot_row, j));
                    sign = -sign;
                }

                if (lu(k, k) == numeric_traits<value_type>::zero())
                    return numeric_traits<value_type>::zero(); // singular

                for (size_type i = k + 1; i < n; ++i) {
                    value_type factor = lu(i, k) / lu(k, k);
                    for (size_type j = k + 1; j < n; ++j)
                        lu(i, j) -= factor * lu(k, j);
                }
            }

            for (size_type i = 0; i < n; ++i)
                det *= lu(i, i);
            if (sign < 0) det = -det;
            return det;
        }
    }

    // Matrix inverse (Gauss-Jordan elimination)
    // Requires division, so it can only be used with element types that form a field
    template<typename T>
    Matrix<T> inverse(const Matrix<T>& a) {
        using size_type = typename Matrix<T>::size_type;

        if (!a.is_square()) {
            throw DimensionError("inverse: matrix must be square");
        }

        const size_type n = a.rows();

        if (n == 0) {
            throw DimensionError("inverse: matrix is empty");
        }

        // 1x1 matrix
        if (n == 1) {
            T val = a(0, 0);
            if (val == numeric_traits<T>::zero()) {
                throw std::runtime_error("inverse: singular matrix");
            }
            Matrix<T> result(1, 1);
            result(0, 0) = numeric_traits<T>::one() / val;
            return result;
        }

        // 2x2 matrix (direct formula)
        if (n == 2) {
            T det = a(0, 0) * a(1, 1) - a(0, 1) * a(1, 0);
            if (det == numeric_traits<T>::zero()) {
                throw std::runtime_error("inverse: singular matrix");
            }
            T inv_det = numeric_traits<T>::one() / det;
            Matrix<T> result(2, 2);
            result(0, 0) = a(1, 1) * inv_det;
            result(0, 1) = -a(0, 1) * inv_det;
            result(1, 0) = -a(1, 0) * inv_det;
            result(1, 1) = a(0, 0) * inv_det;
            return result;
        }

        // General: Gauss-Jordan elimination
        // [A | I] -> [I | A^{-1}]
        Matrix<T> aug(n, 2 * n);
        for (size_type i = 0; i < n; ++i) {
            for (size_type j = 0; j < n; ++j) {
                aug(i, j) = a(i, j);
            }
            aug(i, n + i) = numeric_traits<T>::one();
        }

        for (size_type col = 0; col < n; ++col) {
            // Partial pivoting (uses pivotBetter to choose the best criterion for the type)
            size_type pivot_row = col;
            for (size_type i = col + 1; i < n; ++i) {
                if (!(aug(i, col) == numeric_traits<T>::zero()) &&
                    (aug(pivot_row, col) == numeric_traits<T>::zero() ||
                     numeric_traits<T>::pivotBetter(aug(i, col), aug(pivot_row, col))))
                    pivot_row = i;
            }

            // Swap the pivot rows
            if (pivot_row != col) {
                for (size_type j = 0; j < 2 * n; ++j) {
                    T tmp = aug(col, j);
                    aug(col, j) = aug(pivot_row, j);
                    aug(pivot_row, j) = tmp;
                }
            }

            T pivot = aug(col, col);
            if (pivot == numeric_traits<T>::zero()) {
                throw std::runtime_error("inverse: singular matrix");
            }

            // Scale the pivot row
            T inv_pivot = numeric_traits<T>::one() / pivot;
            for (size_type j = 0; j < 2 * n; ++j) {
                aug(col, j) *= inv_pivot;
            }

            // Eliminate the other rows
            for (size_type i = 0; i < n; ++i) {
                if (i != col) {
                    T factor = aug(i, col);
                    for (size_type j = 0; j < 2 * n; ++j) {
                        aug(i, j) -= factor * aug(col, j);
                    }
                }
            }
        }

        // Extract the result
        Matrix<T> result(n, n);
        for (size_type i = 0; i < n; ++i) {
            for (size_type j = 0; j < n; ++j) {
                result(i, j) = aug(i, n + j);
            }
        }
        return result;
    }

    // Hadamard (element-wise) product of two matrices of identical shape.
    // The Tensor<T> overload lives in Tensor.hpp; this is the dense-matrix counterpart.
    template<typename T>
    [[nodiscard]] Matrix<T> hadamard(const Matrix<T>& a, const Matrix<T>& b) {
        if (a.rows() != b.rows() || a.cols() != b.cols()) {
            throw DimensionError("hadamard: shape mismatch");
        }
        using size_type = typename Matrix<T>::size_type;
        Matrix<T> result(a.rows(), a.cols());
        for (size_type i = 0; i < a.rows(); ++i)
            for (size_type j = 0; j < a.cols(); ++j)
                result(i, j) = a(i, j) * b(i, j);
        return result;
    }

    // log|det(A)| accumulated in the log domain to avoid overflow/underflow for
    // matrices whose determinant is astronomically large or small (e.g. the
    // covariance matrices behind Gaussian log-likelihoods). Uses LU with partial
    // pivoting and sums log|U_ii|. Throws if A is singular.
    template<typename MatrixA>
    [[nodiscard]] auto logdet(const MatrixA& a) -> typename MatrixA::value_type {
        using value_type = typename MatrixA::value_type;
        using size_type = typename MatrixA::size_type;
        static_assert(!numeric_traits<value_type>::is_integer,
                      "logdet requires a floating-point matrix");

        if (!a.is_square()) {
            throw DimensionError("logdet: matrix must be square");
        }

        const size_type n = a.rows();
        if (n == 0) return value_type(0); // det of the empty matrix is 1, log 1 = 0

        Matrix<value_type> lu(n, n);
        for (size_type i = 0; i < n; ++i)
            for (size_type j = 0; j < n; ++j)
                lu(i, j) = a(i, j);

        value_type acc = numeric_traits<value_type>::zero();
        for (size_type k = 0; k < n; ++k) {
            size_type pivot_row = k;
            for (size_type i = k + 1; i < n; ++i) {
                if (lu(i, k) != numeric_traits<value_type>::zero() &&
                    (lu(pivot_row, k) == numeric_traits<value_type>::zero() ||
                     numeric_traits<value_type>::pivotBetter(lu(i, k), lu(pivot_row, k))))
                    pivot_row = i;
            }
            if (pivot_row != k) {
                for (size_type j = 0; j < n; ++j)
                    std::swap(lu(k, j), lu(pivot_row, j));
            }
            if (lu(k, k) == numeric_traits<value_type>::zero()) {
                throw MathError("logdet: singular matrix (determinant is zero)");
            }
            if constexpr (std::is_floating_point_v<value_type>) {
                acc += static_cast<value_type>(
                    std::log(std::abs(static_cast<double>(lu(k, k)))));
            } else {
                // Multi-precision Float etc.: log|diag| kept as value_type (precision preserved)
                acc += detail::generic_log(detail::generic_abs(lu(k, k)));
            }
            for (size_type i = k + 1; i < n; ++i) {
                value_type factor = lu(i, k) / lu(k, k);
                for (size_type j = k + 1; j < n; ++j)
                    lu(i, j) -= factor * lu(k, j);
            }
        }
        return acc;
    }

    // Matrix inverse via textbook Gauss-Jordan elimination on [A | I] with partial
    // pivoting, applied uniformly for every size (no 1x1/2x2 shortcuts). Provided
    // as an explicitly named, didactic counterpart to inverse().
    template<typename T>
    [[nodiscard]] Matrix<T> inverse_gauss_jordan(const Matrix<T>& a) {
        using size_type = typename Matrix<T>::size_type;

        if (!a.is_square()) {
            throw DimensionError("inverse_gauss_jordan: matrix must be square");
        }

        const size_type n = a.rows();
        if (n == 0) {
            throw DimensionError("inverse_gauss_jordan: matrix is empty");
        }

        // [A | I]
        Matrix<T> aug(n, 2 * n);
        for (size_type i = 0; i < n; ++i) {
            for (size_type j = 0; j < n; ++j)
                aug(i, j) = a(i, j);
            aug(i, n + i) = numeric_traits<T>::one();
        }

        for (size_type col = 0; col < n; ++col) {
            // Partial pivoting
            size_type pivot_row = col;
            for (size_type i = col + 1; i < n; ++i) {
                if (!(aug(i, col) == numeric_traits<T>::zero()) &&
                    (aug(pivot_row, col) == numeric_traits<T>::zero() ||
                     numeric_traits<T>::pivotBetter(aug(i, col), aug(pivot_row, col))))
                    pivot_row = i;
            }
            if (pivot_row != col) {
                for (size_type j = 0; j < 2 * n; ++j)
                    std::swap(aug(col, j), aug(pivot_row, j));
            }

            T pivot = aug(col, col);
            if (pivot == numeric_traits<T>::zero()) {
                throw std::runtime_error("inverse_gauss_jordan: singular matrix");
            }

            // Normalize the pivot row, then eliminate the pivot column everywhere else.
            T inv_pivot = numeric_traits<T>::one() / pivot;
            for (size_type j = 0; j < 2 * n; ++j)
                aug(col, j) *= inv_pivot;

            for (size_type i = 0; i < n; ++i) {
                if (i == col) continue;
                T factor = aug(i, col);
                if (factor == numeric_traits<T>::zero()) continue;
                for (size_type j = 0; j < 2 * n; ++j)
                    aug(i, j) -= factor * aug(col, j);
            }
        }

        Matrix<T> result(n, n);
        for (size_type i = 0; i < n; ++i)
            for (size_type j = 0; j < n; ++j)
                result(i, j) = aug(i, n + j);
        return result;
    }

    //-----------------------------------------------------------------------------
    // NoAliasProxy: alias safety hint
    //-----------------------------------------------------------------------------
    //
    // Today operator* returns Matrix by value, so A = A*B is safe via a temporary object.
    // This API exists so that alias checks can be skipped if MatMulExpr is later evaluated lazily.
    // Usage: C.noalias() = A * B;  // declare that C does not overlap with A or B
    //
    template<typename MatrixType>
    struct NoAliasProxy {
        MatrixType& ref;

        // Assignment from a Matrix (copy / move)
        template<typename U>
        MatrixType& operator=(U&& rhs) {
            ref = std::forward<U>(rhs);
            return ref;
        }

        // Add-assign
        template<typename U>
        MatrixType& operator+=(U&& rhs) {
            ref += std::forward<U>(rhs);
            return ref;
        }

        // Subtract-assign
        template<typename U>
        MatrixType& operator-=(U&& rhs) {
            ref -= std::forward<U>(rhs);
            return ref;
        }
    };

    //-----------------------------------------------------------------------------
    // Main matrix classes
    //-----------------------------------------------------------------------------

    // Dynamic-size matrix (shares sangi::BaseMatrix<T> as a common parent with StaticMatrix and can be treated polymorphically)
    template<typename T>
    class Matrix : public BaseMatrix<T>,
                   public detail::MatrixBase<T, detail::DynamicMatrixStorage<T>>,
                   public MatExpr<Matrix<T>> {
    public:
        using Base = detail::MatrixBase<T, detail::DynamicMatrixStorage<T>>;
        using BMBase = BaseMatrix<T>;
        using value_type = typename Base::value_type;
        using size_type = typename Base::size_type;
        using reference = typename Base::reference;
        using const_reference = typename Base::const_reference;
        using iterator = typename Base::iterator;
        using const_iterator = typename Base::const_iterator;

        // Resolve name collisions between BaseMatrix and MatrixBase: prefer the existing MatrixBase API (backward compatibility)
        using Base::operator();
        using Base::at;
        using Base::data;
        using Base::begin;
        using Base::end;
        using Base::cbegin;
        using Base::cend;
        using Base::rows;
        using Base::cols;
        using Base::size;
        using Base::empty;
        using Base::block;
        using Base::is_square;

        // Constructor
        Matrix() : BMBase(), Base() {}

        // Constructor specifying row and column counts
        Matrix(size_type rows, size_type cols) : BMBase(), Base(rows, cols) {
            sync_basematrix();
        }

        // Constructor specifying row count, column count, and initial value
        Matrix(size_type rows, size_type cols, const T& value) : BMBase(), Base(rows, cols, value) {
            sync_basematrix();
        }

        // Constructor using an initializer list (2D)
        Matrix(std::initializer_list<std::initializer_list<T>> init) : BMBase(), Base() {
            size_type rows = init.size();
            size_type cols = rows > 0 ? init.begin()->size() : 0;

            this->Base::resize(rows, cols);

            size_type i = 0;
            for (const auto& row : init) {
                size_type j = 0;
                for (const auto& val : row) {
                    if (j < cols) {
                        (*this)(i, j) = val;
                        ++j;
                    }
                }
                ++i;
            }
            sync_basematrix();
        }

        // Copy/move constructors and assignments
        Matrix(const Matrix& other) : BMBase(), Base(static_cast<const Base&>(other)),
                                       MatExpr<Matrix<T>>() {
            sync_basematrix();
        }
        Matrix(Matrix&& other) noexcept : BMBase(), Base(static_cast<Base&&>(other)),
                                          MatExpr<Matrix<T>>() {
            sync_basematrix();
        }
        Matrix& operator=(const Matrix& other) {
            if (this != &other) {
                Base::operator=(static_cast<const Base&>(other));
                sync_basematrix();
            }
            return *this;
        }
        Matrix& operator=(Matrix&& other) noexcept {
            if (this != &other) {
                Base::operator=(static_cast<Base&&>(other));
                sync_basematrix();
            }
            return *this;
        }

        // SIMD packet load (linear index, row-major contiguous memory)
        auto packet(std::size_t idx) const {
            return PacketTraits<T>::load(this->data() + idx);
        }

        // Construction from an expression template
        template<typename E>
        Matrix(const MatExpr<E>& expr) : BMBase(), Base(expr.derived().rows(), expr.derived().cols()) {
            assignFromMatExpr(expr.derived());
            sync_basematrix();
        }

        // resize/clear/reshape also update the BaseMatrix view
        void resize(size_type r, size_type c) {
            Base::resize(r, c);
            sync_basematrix();
        }
        void resize(size_type r, size_type c, const value_type& v) {
            Base::resize(r, c, v);
            sync_basematrix();
        }
        void clear() noexcept {
            Base::clear();
            sync_basematrix();
        }
        void reshape(size_type r, size_type c) {
            Base::reshape(r, c);
            sync_basematrix();
        }

    public:
        // Assignment from an expression template
        template<typename E>
        Matrix& operator=(const MatExpr<E>& expr) {
            const E& e = expr.derived();
            this->resize(e.rows(), e.cols());
            assignFromMatExpr(e);
            return *this;
        }

    private:
        // Sync the BaseMatrix view from MatrixBase's storage_
        void sync_basematrix() noexcept {
            BMBase::set_view(this->storage_.shared_data(),
                             this->storage_.data(),
                             this->storage_.rows(),
                             this->storage_.cols(),
                             this->storage_.cols(),
                             1);
        }
        template<typename E>
        void assignFromMatExpr(const E& e) {
            const size_type total = this->rows() * this->cols();
            if constexpr (has_simd_packet_v<T>) {
                using PT = PacketTraits<T>;
                constexpr size_type PS = PT::size;
                constexpr size_type PS4 = PS * 4;
                T* __restrict dst = this->data();
                const size_type unroll_end = total - (total % PS4);
                const size_type vec_end = total - (total % PS);
                size_type i = 0;
                for (; i < unroll_end; i += PS4) {
                    PT::store(dst + i,          e.packet(i));
                    PT::store(dst + i + PS,     e.packet(i + PS));
                    PT::store(dst + i + PS * 2, e.packet(i + PS * 2));
                    PT::store(dst + i + PS * 3, e.packet(i + PS * 3));
                }
                for (; i < vec_end; i += PS)
                    PT::store(dst + i, e.packet(i));
                // Tail goes scalar (linear -> (row, col) conversion)
                for (; i < total; ++i)
                    dst[i] = static_cast<T>(e(i / this->cols(), i % this->cols()));
            } else {
                for (size_type i = 0; i < this->rows(); ++i)
                    for (size_type j = 0; j < this->cols(); ++j)
                        (*this)(i, j) = static_cast<T>(e(i, j));
            }
        }

    public:

        // Copy construction from another Matrix
        template<typename OtherT, typename OtherStorage>
        Matrix(const detail::MatrixBase<OtherT, OtherStorage>& other)
            : BMBase(), Base(other.rows(), other.cols()) {

            for (size_type i = 0; i < this->rows(); ++i) {
                for (size_type j = 0; j < this->cols(); ++j) {
                    (*this)(i, j) = static_cast<T>(other(i, j));
                }
            }
            sync_basematrix();
        }

        // Deep-copy ctor from a BaseMatrix<T> view.
        // Needed to make the `Matrix<T> R = a;` pattern (where a is a BaseMatrix<T> argument)
        // work in the body of LinAlg functions. Whether the caller passes Matrix<T> or StaticMatrix,
        // exactly one deep copy runs here (= only the unavoidable copy).
        // Note: Matrix<T> itself derives from BaseMatrix<T>, but the normal Matrix -> Matrix copy ctor
        // (Matrix(const Matrix&)) is selected by overload resolution (exact match). This base copy ctor
        // is used for other BaseMatrix-derived classes such as StaticMatrix or for BaseMatrix view arguments.
        Matrix(const BaseMatrix<T>& other)
            : BMBase(), Base(other.rows(), other.cols()) {
            for (size_type i = 0; i < this->rows(); ++i) {
                for (size_type j = 0; j < this->cols(); ++j) {
                    (*this)(i, j) = other(i, j);
                }
            }
            sync_basematrix();
        }

        // Copy construction from a MatrixView
        template<typename MatrixType>
        Matrix(const detail::MatrixView<MatrixType>& view)
            : BMBase(), Base(view.rows(), view.cols()) {

            for (size_type i = 0; i < this->rows(); ++i) {
                for (size_type j = 0; j < this->cols(); ++j) {
                    (*this)(i, j) = static_cast<T>(view(i, j));
                }
            }
            sync_basematrix();
        }

        // Copy construction from a ConstMatrixView
        template<typename MatrixType>
        Matrix(const detail::ConstMatrixView<MatrixType>& view)
            : BMBase(), Base(view.rows(), view.cols()) {

            for (size_type i = 0; i < this->rows(); ++i) {
                for (size_type j = 0; j < this->cols(); ++j) {
                    (*this)(i, j) = static_cast<T>(view(i, j));
                }
            }
            sync_basematrix();
        }

        // Destructor
        ~Matrix() = default;

        /// Triangular matrix view
        template<TriangularMode Mode>
        TriangularView<Matrix, Mode> triangularView() {
            return TriangularView<Matrix, Mode>(*this);
        }
        template<TriangularMode Mode>
        TriangularView<const Matrix, Mode> triangularView() const {
            return TriangularView<const Matrix, Mode>(const_cast<const Matrix&>(*this));
        }

        /// Symmetric (self-adjoint) matrix view
        template<TriangularMode Mode>
        SelfAdjointView<Matrix, Mode> selfAdjointView() {
            return SelfAdjointView<Matrix, Mode>(*this);
        }
        template<TriangularMode Mode>
        SelfAdjointView<const Matrix, Mode> selfAdjointView() const {
            return SelfAdjointView<const Matrix, Mode>(const_cast<const Matrix&>(*this));
        }

        // Add operator
        Matrix& operator+=(const Matrix& rhs) {
            if (this->rows() != rhs.rows() || this->cols() != rhs.cols()) {
                throw DimensionError("Matrix addition: dimension mismatch");
            }

            for (size_type i = 0; i < this->rows(); ++i) {
                for (size_type j = 0; j < this->cols(); ++j) {
                    (*this)(i, j) += rhs(i, j);
                }
            }

            return *this;
        }

        // Subtract operator
        Matrix& operator-=(const Matrix& rhs) {
            if (this->rows() != rhs.rows() || this->cols() != rhs.cols()) {
                throw DimensionError("Matrix subtraction: dimension mismatch");
            }

            for (size_type i = 0; i < this->rows(); ++i) {
                for (size_type j = 0; j < this->cols(); ++j) {
                    (*this)(i, j) -= rhs(i, j);
                }
            }

            return *this;
        }

        // Scalar multiplication operator
        Matrix& operator*=(const T& scalar) {
            for (size_type i = 0; i < this->rows(); ++i) {
                for (size_type j = 0; j < this->cols(); ++j) {
                    (*this)(i, j) *= scalar;
                }
            }

            return *this;
        }

        // Scalar division operator
        Matrix& operator/=(const T& scalar) {
            // Division-by-zero check
            if (scalar == T{ 0 }) {
                throw std::invalid_argument("Division by zero");
            }

            for (size_type i = 0; i < this->rows(); ++i) {
                for (size_type j = 0; j < this->cols(); ++j) {
                    (*this)(i, j) /= scalar;
                }
            }

            return *this;
        }

        // Compound assignment from an expression template
        template<typename E>
        Matrix& operator+=(const MatExpr<E>& rhs) {
            const E& r = rhs.derived();
            for (size_type i = 0; i < this->rows(); ++i)
                for (size_type j = 0; j < this->cols(); ++j)
                    (*this)(i, j) += static_cast<T>(r(i, j));
            return *this;
        }

        template<typename E>
        Matrix& operator-=(const MatExpr<E>& rhs) {
            const E& r = rhs.derived();
            for (size_type i = 0; i < this->rows(); ++i)
                for (size_type j = 0; j < this->cols(); ++j)
                    (*this)(i, j) -= static_cast<T>(r(i, j));
            return *this;
        }

        // Obtain the transposed matrix (allocates a new Matrix and returns it)
        // For large matrices this has the overhead of a heap allocation + zero-initialization;
        // use transpose_into() to write into a pre-allocated buffer when high performance is needed.
        Matrix transpose() const {
            Matrix result(this->cols(), this->rows());
            DefaultComputePolicy<T>::transpose(*this, result);
            return result;
        }

        // Transpose into an existing buffer (allocation-free)
        // dst must already be sized as (cols x rows)
        void transpose_into(Matrix& dst) const {
            DefaultComputePolicy<T>::transpose(*this, dst);
        }

        // Create the identity matrix (static method)
        static Matrix identity(size_type size) {
            Matrix result(size, size);

            // Initialize with zero
            result.zero();

            // Set the diagonal entries to 1
            for (size_type i = 0; i < size; ++i) {
                result(i, i) = numeric_traits<T>::one();
            }

            return result;
        }

        // Determinant computation
        T determinant() const {
            return sangi::determinant(*this);
        }

        // Frobenius norm computation
        T norm() const {
            return sangi::frobenius_norm(*this);
        }

        // L1 matrix norm (maximum column sum)
        T norm_l1() const {
            return sangi::matrix_norm_l1(*this);
        }

        // L-infinity matrix norm (maximum row sum)
        T norm_linf() const {
            return sangi::matrix_norm_linf(*this);
        }

        // Trace computation
        T trace() const {
            return sangi::trace(*this);
        }

        // Inverse computation
        [[nodiscard]] Matrix inverse() const {
            return sangi::inverse(*this);
        }

        /// operator^: A^'T' = transpose, A^'H' = Hermitian transpose, A^(-1) = inverse, A^n = power
        ///
        /// Examples:
        ///   auto At = A ^ 'T';    // transpose
        ///   auto Ah = A ^ 'H';    // Hermitian transpose (conjugate transpose)
        ///   auto Ainv = A ^ (-1); // inverse
        ///   auto A3 = A ^ 3;      // A*A*A
        [[nodiscard]] Matrix operator^(int n) const {
            if (n == 'T') return this->transpose();
            if (n == 'H') {
                // Hermitian transpose: conj(A^T)
                Matrix result = this->transpose();
                for (size_type i = 0; i < result.rows(); ++i)
                    for (size_type j = 0; j < result.cols(); ++j)
                        result(i, j) = numeric_traits<T>::conj(result(i, j));
                return result;
            }
            if (n == -1) return this->inverse();
            if (n < 0)
                return sangi::power(this->inverse(), static_cast<unsigned int>(-n));
            return sangi::power(*this, static_cast<unsigned int>(n));
        }

        /// noalias() hint: declares that the assignment target does not overlap with the right-hand-side matrices.
        ///
        /// Today operator*(Matrix, Matrix) returns a new Matrix, so A = A * B is always safe.
        /// This is a reserved API for future lazy-evaluation optimizations.
        /// Example: C.noalias() = A * B;
        [[nodiscard]] NoAliasProxy<Matrix> noalias() { return {*this}; }

        /// Solve Ax = b and return x (forward/back substitution safe version). Returns std::nullopt if singular or on size mismatch.
        ///
        /// The internal implementation wraps linalg::solve(Matrix, Vector) from <math/linalg/solvers.hpp>
        /// in a try-catch. Callers therefore need to include
        ///     #include <math/linalg/solvers.hpp>
        /// (it is undefined when only this header is included).
        ///
        /// Primary use case: added so that sangi::Matrix satisfies the root_finding_nd concept requirement
        ///     { m.solve(v) } -> std::convertible_to<std::optional<V>>
        /// directly.
        [[nodiscard]] std::optional<Vector<T>> solve(const Vector<T>& b) const;
   };

   // Fixed-size matrix (shares BaseMatrix<T> as a common parent, allowing direct upcast to LinAlg functions)
   template<typename T, std::size_t Rows, std::size_t Cols>
   class StaticMatrix : public BaseMatrix<T>,
                        public detail::MatrixBase<T, detail::StaticMatrixStorage<T, Rows, Cols>>,
                        public MatExpr<StaticMatrix<T, Rows, Cols>> {
   public:
       // Compile-time dimension traits. Shadow BaseMatrix::static_rows / static_cols
       // (-1) so that LinAlg algorithms can static_assert dimension agreement at build time.
       static constexpr std::ptrdiff_t static_rows = static_cast<std::ptrdiff_t>(Rows);
       static constexpr std::ptrdiff_t static_cols = static_cast<std::ptrdiff_t>(Cols);

       using Base = detail::MatrixBase<T, detail::StaticMatrixStorage<T, Rows, Cols>>;
       using BMBase = BaseMatrix<T>;
       using value_type = typename Base::value_type;
       using size_type = typename Base::size_type;
       using reference = typename Base::reference;
       using const_reference = typename Base::const_reference;
       using iterator = typename Base::iterator;
       using const_iterator = typename Base::const_iterator;

       // Resolve name collisions between BaseMatrix and MatrixBase: prefer the existing MatrixBase API (backward compatibility)
       using Base::operator();
       using Base::at;
       using Base::data;
       using Base::begin;
       using Base::end;
       using Base::cbegin;
       using Base::cend;
       using Base::rows;
       using Base::cols;
       using Base::size;
       using Base::empty;
       using Base::block;
       using Base::is_square;

       // Constructor
       StaticMatrix() : BMBase(), Base() {
           sync_basematrix();
       }

       // Constructor specifying an initial value
       explicit StaticMatrix(const T& value) : BMBase(), Base() {
           this->Base::fill(value);
           sync_basematrix();
       }

       // Constructor using an initializer list (2D)
       StaticMatrix(std::initializer_list<std::initializer_list<T>> init) : BMBase(), Base() {
           size_type i = 0;
           for (const auto& row : init) {
               if (i >= Rows) break;

               size_type j = 0;
               for (const auto& val : row) {
                   if (j >= Cols) break;
                   (*this)(i, j) = val;
                   ++j;
               }
               ++i;
           }
           sync_basematrix();
       }

       // Copy/move constructors and assignments (also sync the BaseMatrix view)
       StaticMatrix(const StaticMatrix& other) : BMBase(), Base(static_cast<const Base&>(other)),
                                                  MatExpr<StaticMatrix<T, Rows, Cols>>() {
           sync_basematrix();
       }
       StaticMatrix(StaticMatrix&& other) noexcept : BMBase(), Base(static_cast<Base&&>(other)),
                                                      MatExpr<StaticMatrix<T, Rows, Cols>>() {
           sync_basematrix();
       }
       StaticMatrix& operator=(const StaticMatrix& other) {
           if (this != &other) {
               Base::operator=(static_cast<const Base&>(other));
               // No sync needed — although the data inside storage_ is swapped, the value of data() remains
               // the address of our own storage_, and rows/cols/strides are unchanged.
               // However, if MatrixBase's copy-assign is designed to swap the entire storage_, then sync
               // is required (for fixed-size storage only the array elements are copied; the ptr is unchanged).
           }
           return *this;
       }
       StaticMatrix& operator=(StaticMatrix&&) noexcept = default;

       // SIMD packet load (linear index)
       auto packet(std::size_t idx) const {
           return PacketTraits<T>::load(this->data() + idx);
       }

       // Construction from an expression template
       template<typename E>
       StaticMatrix(const MatExpr<E>& expr) : BMBase(), Base() {
           assignFromMatExpr(expr.derived());
           sync_basematrix();
       }

       // Assignment from an expression template
       template<typename E>
       StaticMatrix& operator=(const MatExpr<E>& expr) {
           assignFromMatExpr(expr.derived());
           return *this;
       }

   private:
       // Sync the BaseMatrix view from MatrixBase's storage_ (StaticMatrixStorage = std::array).
       // Fixed-size storage lives on the stack and does not have shared_ptr ownership, so use set_view_nonowning.
       void sync_basematrix() noexcept {
           BMBase::set_view_nonowning(this->Base::data(), Rows, Cols, Cols, 1);
       }
   public:

    private:
       template<typename E>
       void assignFromMatExpr(const E& e) {
           constexpr size_type total = Rows * Cols;
           if constexpr (has_simd_packet_v<T> && total >= PacketTraits<T>::size) {
               constexpr size_type PS = PacketTraits<T>::size;
               constexpr size_type vec_end = total - (total % PS);
               T* __restrict dst = this->data();
               for (size_type i = 0; i < vec_end; i += PS)
                   PacketTraits<T>::store(dst + i, e.packet(i));
               for (size_type i = vec_end; i < total; ++i)
                   dst[i] = static_cast<T>(e(i / Cols, i % Cols));
           } else {
               for (size_type i = 0; i < Rows; ++i)
                   for (size_type j = 0; j < Cols; ++j)
                       (*this)(i, j) = static_cast<T>(e(i, j));
           }
       }

    public:

       // Copy construction from another Matrix
       template<typename OtherT, typename OtherStorage>
       StaticMatrix(const detail::MatrixBase<OtherT, OtherStorage>& other) {
           if (other.rows() != Rows || other.cols() != Cols) {
               throw DimensionError("StaticMatrix: dimension mismatch in copy constructor");
           }

           for (size_type i = 0; i < Rows; ++i) {
               for (size_type j = 0; j < Cols; ++j) {
                   (*this)(i, j) = static_cast<T>(other(i, j));
               }
           }
       }

       // Copy construction from a MatrixView
       template<typename MatrixType>
       StaticMatrix(const detail::MatrixView<MatrixType>& view) {
           if (view.rows() != Rows || view.cols() != Cols) {
               throw DimensionError("StaticMatrix: dimension mismatch in copy constructor from view");
           }

           for (size_type i = 0; i < Rows; ++i) {
               for (size_type j = 0; j < Cols; ++j) {
                   (*this)(i, j) = static_cast<T>(view(i, j));
               }
           }
       }

       // Copy construction from a ConstMatrixView
       template<typename MatrixType>
       StaticMatrix(const detail::ConstMatrixView<MatrixType>& view) {
           if (view.rows() != Rows || view.cols() != Cols) {
               throw DimensionError("StaticMatrix: dimension mismatch in copy constructor from const view");
           }

           for (size_type i = 0; i < Rows; ++i) {
               for (size_type j = 0; j < Cols; ++j) {
                   (*this)(i, j) = static_cast<T>(view(i, j));
               }
           }
       }

       // Destructor
       ~StaticMatrix() = default;

       // Add operator
       StaticMatrix& operator+=(const StaticMatrix& rhs) {
           for (size_type i = 0; i < Rows; ++i) {
               for (size_type j = 0; j < Cols; ++j) {
                   (*this)(i, j) += rhs(i, j);
               }
           }

           return *this;
       }

       // Subtract operator
       StaticMatrix& operator-=(const StaticMatrix& rhs) {
           for (size_type i = 0; i < Rows; ++i) {
               for (size_type j = 0; j < Cols; ++j) {
                   (*this)(i, j) -= rhs(i, j);
               }
           }

           return *this;
       }

       // Scalar multiplication operator
       StaticMatrix& operator*=(const T& scalar) {
           for (size_type i = 0; i < Rows; ++i) {
               for (size_type j = 0; j < Cols; ++j) {
                   (*this)(i, j) *= scalar;
               }
           }

           return *this;
       }

       // Scalar division operator
       StaticMatrix& operator/=(const T& scalar) {
           // Division-by-zero check
           if (scalar == T{ 0 }) {
               throw std::invalid_argument("Division by zero");
           }

           for (size_type i = 0; i < Rows; ++i) {
               for (size_type j = 0; j < Cols; ++j) {
                   (*this)(i, j) /= scalar;
               }
           }

           return *this;
       }

       // Compound assignment from an expression template
       template<typename E>
       StaticMatrix& operator+=(const MatExpr<E>& rhs) {
           const E& r = rhs.derived();
           for (size_type i = 0; i < Rows; ++i)
               for (size_type j = 0; j < Cols; ++j)
                   (*this)(i, j) += static_cast<T>(r(i, j));
           return *this;
       }

       template<typename E>
       StaticMatrix& operator-=(const MatExpr<E>& rhs) {
           const E& r = rhs.derived();
           for (size_type i = 0; i < Rows; ++i)
               for (size_type j = 0; j < Cols; ++j)
                   (*this)(i, j) -= static_cast<T>(r(i, j));
           return *this;
       }

       // Obtain the transposed matrix
       StaticMatrix<T, Cols, Rows> transpose() const {
           StaticMatrix<T, Cols, Rows> result;

           for (size_type i = 0; i < Rows; ++i) {
               for (size_type j = 0; j < Cols; ++j) {
                   result(j, i) = (*this)(i, j);
               }
           }

           return result;
       }

       // Specialized methods for square matrices
       T determinant() const requires (Rows == Cols) {
           return sangi::determinant(*this);
       }

       T trace() const requires (Rows == Cols) {
           return sangi::trace(*this);
       }

       static StaticMatrix identity() requires (Rows == Cols) {
           StaticMatrix result;
           result.zero();  // First zero-initialize
           for (size_type i = 0; i < Rows; ++i) {
               result(i, i) = numeric_traits<T>::one();  // Set the diagonal entries to 1
           }
           return result;
       }

       // Frobenius norm computation
       T norm() const {
           return sangi::frobenius_norm(*this);
       }

       // L1 matrix norm (maximum column sum)
       T norm_l1() const {
           return sangi::matrix_norm_l1(*this);
       }

       // L-infinity matrix norm (maximum row sum)
       T norm_linf() const {
           return sangi::matrix_norm_linf(*this);
       }

   };

   //-----------------------------------------------------------------------------
   // Global operators
   //-----------------------------------------------------------------------------

   // Element-wise operations (+, -, *scalar, /scalar, unary -) fully delegate to
   // the MatExpr operators in expr_templates.hpp (expression-template chain enabled),
   // so a compound expression like 2.0*A + B - C/3.0 is evaluated with SIMD in one pass.

   // Matrix-matrix multiplication
   template<typename T>
   Matrix<T> operator*(const Matrix<T>& lhs, const Matrix<T>& rhs) {
       Matrix<T> result(lhs.rows(), rhs.cols());
       multiply(lhs, rhs, result);
       return result;
   }

   // Matrix-matrix multiplication (MatExpr overload: evaluates the expression first, then multiplies)
   template<typename T, typename E>
   Matrix<T> operator*(const Matrix<T>& lhs, const MatExpr<E>& rhs) {
       Matrix<T> rhs_eval = rhs;
       return lhs * rhs_eval;
   }

   template<typename T, typename E>
   Matrix<T> operator*(const MatExpr<E>& lhs, const Matrix<T>& rhs) {
       Matrix<T> lhs_eval = lhs;
       return lhs_eval * rhs;
   }

   // Matrix-vector multiplication
   template<typename T>
   Vector<T> operator*(const Matrix<T>& mat, const Vector<T>& vec) {
       Vector<T> result(mat.rows());
       multiply_vector(mat, vec, result);
       return result;
   }

   // StaticMatrix element-wise operations also fully delegate to ET (via the MatExpr base)

   template<typename T, std::size_t M, std::size_t N, std::size_t P>
   StaticMatrix<T, M, P> operator*(const StaticMatrix<T, M, N>& lhs,
       const StaticMatrix<T, N, P>& rhs) {
       StaticMatrix<T, M, P> result;
       multiply(lhs, rhs, result);
       return result;
   }

   // Matrix-vector multiplication (StaticMatrix * dynamic Vector)
   template<typename T, std::size_t Rows, std::size_t Cols>
   Vector<T> operator*(const StaticMatrix<T, Rows, Cols>& mat, const Vector<T>& vec) {
       Vector<T> result(Rows);
       multiply_vector(mat, vec, result);
       return result;
   }

   // Matrix-vector multiplication (StaticMatrix * StaticVector)
   template<typename T, std::size_t Rows, std::size_t Cols>
   StaticVector<T, Rows> operator*(const StaticMatrix<T, Rows, Cols>& mat,
       const StaticVector<T, Cols>& vec) {
       StaticVector<T, Rows> result;
       for (std::size_t i = 0; i < Rows; ++i) {
           T sum = T(0);
           for (std::size_t j = 0; j < Cols; ++j) {
               sum += mat(i, j) * vec[j];
           }
           result[i] = sum;
       }
       return result;
   }

   // StaticMatrix inverse (2x2 direct formula)
   template<typename T>
   [[nodiscard]] StaticMatrix<T, 2, 2> inverse(const StaticMatrix<T, 2, 2>& m) {
       T det = m(0, 0) * m(1, 1) - m(0, 1) * m(1, 0);
       if (std::abs(det) < std::numeric_limits<T>::epsilon()) {
           throw LinearAlgebraError("inverse: singular 2x2 matrix");
       }
       T inv_det = T(1) / det;
       StaticMatrix<T, 2, 2> result;
       result(0, 0) =  m(1, 1) * inv_det;
       result(0, 1) = -m(0, 1) * inv_det;
       result(1, 0) = -m(1, 0) * inv_det;
       result(1, 1) =  m(0, 0) * inv_det;
       return result;
   }

   // StaticMatrix inverse (3x3 via cofactor matrix)
   template<typename T>
   [[nodiscard]] StaticMatrix<T, 3, 3> inverse(const StaticMatrix<T, 3, 3>& m) {
       T det = m(0, 0) * (m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1))
             - m(0, 1) * (m(1, 0) * m(2, 2) - m(1, 2) * m(2, 0))
             + m(0, 2) * (m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0));
       if (std::abs(det) < std::numeric_limits<T>::epsilon()) {
           throw LinearAlgebraError("inverse: singular 3x3 matrix");
       }
       T inv_det = T(1) / det;
       StaticMatrix<T, 3, 3> result;
       result(0, 0) = (m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1)) * inv_det;
       result(0, 1) = (m(0, 2) * m(2, 1) - m(0, 1) * m(2, 2)) * inv_det;
       result(0, 2) = (m(0, 1) * m(1, 2) - m(0, 2) * m(1, 1)) * inv_det;
       result(1, 0) = (m(1, 2) * m(2, 0) - m(1, 0) * m(2, 2)) * inv_det;
       result(1, 1) = (m(0, 0) * m(2, 2) - m(0, 2) * m(2, 0)) * inv_det;
       result(1, 2) = (m(0, 2) * m(1, 0) - m(0, 0) * m(1, 2)) * inv_det;
       result(2, 0) = (m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0)) * inv_det;
       result(2, 1) = (m(0, 1) * m(2, 0) - m(0, 0) * m(2, 1)) * inv_det;
       result(2, 2) = (m(0, 0) * m(1, 1) - m(0, 1) * m(1, 0)) * inv_det;
       return result;
   }

   //-----------------------------------------------------------------------------
   // Utility functions
   //-----------------------------------------------------------------------------

   // Create a zero matrix
   template<typename T>
   Matrix<T> zeros(std::size_t rows, std::size_t cols) {
       Matrix<T> result(rows, cols);
       result.zero();
       return result;
   }

   // Create an identity matrix
   template<typename T>
   Matrix<T> eye(std::size_t size) {
       return Matrix<T>::identity(size);
   }

   // Create an identity matrix (identity<T>(n) form)
   template<typename T>
   Matrix<T> identity(std::size_t size) {
       return Matrix<T>::identity(size);
   }

   // Create a matrix whose elements are all 1
   template<typename T>
   Matrix<T> ones(std::size_t rows, std::size_t cols) {
       return Matrix<T>(rows, cols, T{ 1 });
   }

   // Special 2x2 matrix — rotation matrix.
   // Returns StaticMatrix<T, 2, 2> for stack allocation and compile-time dimensions.
   // Implicit conversion to Matrix<T> keeps existing callers compatible.
   template<typename T>
   StaticMatrix<T, 2, 2> rotation2D(T angle) {
       const T cos_val = std::cos(angle);
       const T sin_val = std::sin(angle);
       StaticMatrix<T, 2, 2> result;
       result(0, 0) = cos_val;
       result(0, 1) = -sin_val;
       result(1, 0) = sin_val;
       result(1, 1) = cos_val;
       return result;
   }

   // Special 3x3 matrices — rotation matrices around the 3D axes.
   // Returns StaticMatrix<T, 3, 3> for stack allocation and compile-time dimensions.
   template<typename T>
   StaticMatrix<T, 3, 3> rotationX(T angle) {
       const T cos_val = std::cos(angle);
       const T sin_val = std::sin(angle);
       StaticMatrix<T, 3, 3> result;
       result.zero();
       result(0, 0) = T{ 1 };
       result(1, 1) = cos_val;
       result(1, 2) = -sin_val;
       result(2, 1) = sin_val;
       result(2, 2) = cos_val;
       return result;
   }

   template<typename T>
   StaticMatrix<T, 3, 3> rotationY(T angle) {
       const T cos_val = std::cos(angle);
       const T sin_val = std::sin(angle);
       StaticMatrix<T, 3, 3> result;
       result.zero();
       result(0, 0) = cos_val;
       result(0, 2) = sin_val;
       result(1, 1) = T{ 1 };
       result(2, 0) = -sin_val;
       result(2, 2) = cos_val;
       return result;
   }

   template<typename T>
   StaticMatrix<T, 3, 3> rotationZ(T angle) {
       const T cos_val = std::cos(angle);
       const T sin_val = std::sin(angle);
       StaticMatrix<T, 3, 3> result;
       result.zero();
       result(0, 0) = cos_val;
       result(0, 1) = -sin_val;
       result(1, 0) = sin_val;
       result(1, 1) = cos_val;
       result(2, 2) = T{ 1 };
       return result;
   }

   // Matrix concatenation — horizontal
   template<typename T>
   Matrix<T> hconcat(const Matrix<T>& a, const Matrix<T>& b) {
       if (a.rows() != b.rows()) {
           throw DimensionError("hconcat: matrices must have the same number of rows");
       }

       Matrix<T> result(a.rows(), a.cols() + b.cols());

       for (std::size_t i = 0; i < a.rows(); ++i) {
           for (std::size_t j = 0; j < a.cols(); ++j) {
               result(i, j) = a(i, j);
           }
           for (std::size_t j = 0; j < b.cols(); ++j) {
               result(i, a.cols() + j) = b(i, j);
           }
       }

       return result;
   }

   // Matrix concatenation — vertical
   template<typename T>
   Matrix<T> vconcat(const Matrix<T>& a, const Matrix<T>& b) {
       if (a.cols() != b.cols()) {
           throw DimensionError("vconcat: matrices must have the same number of columns");
       }

       Matrix<T> result(a.rows() + b.rows(), a.cols());

       for (std::size_t j = 0; j < a.cols(); ++j) {
           for (std::size_t i = 0; i < a.rows(); ++i) {
               result(i, j) = a(i, j);
           }
           for (std::size_t i = 0; i < b.rows(); ++i) {
               result(a.rows() + i, j) = b(i, j);
           }
       }

       return result;
   }

   // Matrix reshape
   template<typename T>
   Matrix<T> reshape(const Matrix<T>& mat, std::size_t new_rows, std::size_t new_cols) {
       if (mat.cols() == 0 || mat.rows() == 0) {
           if (new_rows == 0 || new_cols == 0) {
               return Matrix<T>(new_rows, new_cols);
           }
           throw DimensionError("reshape: total element count must be preserved");
       }
       if (mat.rows() * mat.cols() != new_rows * new_cols) {
           throw DimensionError("reshape: total element count must be preserved");
       }

       Matrix<T> result(new_rows, new_cols);

       for (std::size_t i = 0; i < mat.rows() * mat.cols(); ++i) {
           const std::size_t old_row = i / mat.cols();
           const std::size_t old_col = i % mat.cols();
           const std::size_t new_row = i / new_cols;
           const std::size_t new_col = i % new_cols;

           result(new_row, new_col) = mat(old_row, old_col);
       }

       return result;
   }

   //-----------------------------------------------------------------------------
   // Matrix-style I/O
   //-----------------------------------------------------------------------------

   // Print a matrix to standard output
   template<typename T>
   void print_matrix(const Matrix<T>& mat, std::ostream& os = std::cout,
       std::size_t width = 12, std::size_t precision = 6) {
       // Save the previous output settings
       std::ios::fmtflags old_flags = os.flags();
       std::streamsize old_precision = os.precision();

       // Set the output format
       os << std::fixed << std::setprecision(precision);

       const std::size_t rows = mat.rows();
       const std::size_t cols = mat.cols();

       for (std::size_t i = 0; i < rows; ++i) {
           for (std::size_t j = 0; j < cols; ++j) {
               os << std::setw(width) << mat(i, j);
               if (j < cols - 1) os << " ";
           }
           os << "\n";
       }

       // Restore the previous output settings
       os.flags(old_flags);
       os.precision(old_precision);
   }

   // Simple formatted output
   template<typename T>
   std::string to_string(const Matrix<T>& mat, std::size_t precision = 4) {
       std::ostringstream ss;
       print_matrix(mat, ss, 8, precision);
       return ss.str();
   }

   // operator<< overload
   template<typename T>
   std::ostream& operator<<(std::ostream& os, const Matrix<T>& mat) {
       print_matrix(mat, os);
       return os;
   }

   // LaTeX output of a matrix
   // style: "pmatrix" (parentheses), "bmatrix" (brackets), "vmatrix" (vertical bars), "matrix" (no delimiters)
   template<typename T>
   [[nodiscard]] std::string toLatex(const Matrix<T>& mat,
       const std::string& style = "pmatrix", std::size_t precision = 6) {
       std::ostringstream oss;
       std::ios::fmtflags old_flags = oss.flags();
       std::streamsize old_prec = oss.precision();
       oss << std::setprecision(precision);
       oss << "\\begin{" << style << "}\n";
       for (std::size_t i = 0; i < mat.rows(); ++i) {
           for (std::size_t j = 0; j < mat.cols(); ++j) {
               if (j > 0) oss << " & ";
               oss << mat(i, j);
           }
           if (i < mat.rows() - 1) oss << " \\\\";
           oss << "\n";
       }
       oss << "\\end{" << style << "}";
       oss.flags(old_flags);
       oss.precision(old_prec);
       return oss.str();
   }

   // Matrix power (left-to-right binary method)
   // Square matrices only. Specializes the generic power() from power.hpp for Matrix.
   template<typename T>
   [[nodiscard]] Matrix<T> power(const Matrix<T>& base, unsigned int exponent) {
       if (base.rows() != base.cols()) {
           throw DimensionError("power(): matrix must be square");
       }
       if (exponent == 0) return Matrix<T>::identity(base.rows());
       if (exponent == 1) return base;

       Matrix<T> y = Matrix<T>::identity(base.rows());
       unsigned int bit = std::bit_width(exponent) - 1;
       while (true) {
           y = y * y;
           if ((exponent >> bit) & 1u) {
               y = y * base;
           }
           if (bit == 0) break;
           --bit;
       }
       return y;
   }

    // =========================================================================
    // MatrixMap — zero-copy view onto external memory (equivalent to Eigen::Map)
    // =========================================================================

    /**
     * @brief Matrix view onto external memory (row-major)
     *
     * @code
     * double data[6] = {1,2,3,4,5,6};
     * auto M = sangi::MatrixMap<double>(data, 2, 3);
     * M(0, 1) = 10.0;  // modifies data[1]
     * @endcode
     */
    template<typename T>
    class MatrixMap {
    public:
        using value_type = T;
        using size_type = std::size_t;

        MatrixMap(T* data, size_type rows, size_type cols)
            : data_(data), rows_(rows), cols_(cols) {}

        size_type rows() const { return rows_; }
        size_type cols() const { return cols_; }
        T* data() { return data_; }
        const T* data() const { return data_; }

        T& operator()(size_type i, size_type j) { return data_[i * cols_ + j]; }
        const T& operator()(size_type i, size_type j) const { return data_[i * cols_ + j]; }

        /// Copy into a Matrix
        operator Matrix<T>() const {
            Matrix<T> m(rows_, cols_);
            for (size_type i = 0; i < rows_; ++i)
                for (size_type j = 0; j < cols_; ++j)
                    m(i, j) = data_[i * cols_ + j];
            return m;
        }

        /// Assignment from a Matrix
        MatrixMap& operator=(const Matrix<T>& m) {
            for (size_type i = 0; i < std::min(rows_, m.rows()); ++i)
                for (size_type j = 0; j < std::min(cols_, m.cols()); ++j)
                    data_[i * cols_ + j] = m(i, j);
            return *this;
        }

    private:
        T* data_;
        size_type rows_, cols_;
    };

    template<typename T>
    class ConstMatrixMap {
    public:
        using value_type = T;
        using size_type = std::size_t;

        ConstMatrixMap(const T* data, size_type rows, size_type cols)
            : data_(data), rows_(rows), cols_(cols) {}

        size_type rows() const { return rows_; }
        size_type cols() const { return cols_; }
        const T* data() const { return data_; }

        const T& operator()(size_type i, size_type j) const { return data_[i * cols_ + j]; }

        operator Matrix<T>() const {
            Matrix<T> m(rows_, cols_);
            for (size_type i = 0; i < rows_; ++i)
                for (size_type j = 0; j < cols_; ++j)
                    m(i, j) = data_[i * cols_ + j];
            return m;
        }

    private:
        const T* data_;
        size_type rows_, cols_;
    };

} // namespace sangi

#ifdef __GNUC__
#pragma GCC diagnostic pop
#endif

// ---------------------------------------------------------------------------
// extern template — suppress re-instantiation for float/double.
// The definitions are provided in src/math/core/matrix.cpp (sangi_matrix.lib).
// Other types T are instantiated from the header as usual.
// ---------------------------------------------------------------------------
#ifndef SANGI_MATRIX_EXPLICIT_INSTANTIATION
namespace sangi {
    extern template class detail::DynamicMatrixStorage<float>;
    extern template class detail::DynamicMatrixStorage<double>;
    extern template class detail::MatrixBase<float, detail::DynamicMatrixStorage<float>>;
    extern template class detail::MatrixBase<double, detail::DynamicMatrixStorage<double>>;
    extern template class Matrix<float>;
    extern template class Matrix<double>;
    extern template Matrix<float> power(const Matrix<float>&, unsigned int);
    extern template Matrix<double> power(const Matrix<double>&, unsigned int);
    extern template std::ostream& operator<<(std::ostream&, const Matrix<float>&);
    extern template std::ostream& operator<<(std::ostream&, const Matrix<double>&);
    extern template std::string toLatex(const Matrix<float>&, const std::string&, std::size_t);
    extern template std::string toLatex(const Matrix<double>&, const std::string&, std::size_t);
    extern template Matrix<float> zeros(std::size_t, std::size_t);
    extern template Matrix<double> zeros(std::size_t, std::size_t);
    extern template Matrix<float> ones(std::size_t, std::size_t);
    extern template Matrix<double> ones(std::size_t, std::size_t);
    extern template Matrix<float> identity(std::size_t);
    extern template Matrix<double> identity(std::size_t);
} // namespace sangi
#endif

#endif // SANGI_MATRIX_HPP