How does Matrix store elements internally?

sangi::Matrix uses row-major (C-style) contiguous storage. Elements of row i are stored consecutively in memory, making row traversal cache-friendly. The raw pointer is accessible via data().

What does the operator^ do for matrices?

operator^ provides shorthand for transpose (A ^ 'T'), Hermitian transpose (A ^ 'H'), inverse (A ^ (-1)), and matrix power (A ^ n) using binary exponentiation. Note that ^ has low precedence in C++, so parentheses are required: (A ^ 'T') * B.

What is StaticMatrix and when should I use it?

StaticMatrix is a compile-time fixed-size matrix stored on the stack. It provides the same API as Matrix but avoids heap allocation and enables compile-time size checking for multiplication: StaticMatrix * StaticMatrix yields StaticMatrix .

Can Matrix be used with arbitrary-precision types?

Yes. Matrix is a template that works with any numeric type, including sangi::Rational for exact arithmetic. For example, Matrix can compute the exact inverse of a Hilbert matrix with no floating-point error.

Matrix — Dense Matrix

Overview

Matrix<T> is the dynamic-size dense matrix template class in sangi.

Dynamic size, row-major storage
Full operator overloading — including matrix product and matrix-vector product
Determinant, inverse, transpose, norm, trace
Views — zero-copy references to sub-matrices, rows, and columns
Factory functions — zeros, identity, rotation matrices
Static-size variant: StaticMatrix<T, Rows, Cols>

Header

#include <math/core/matrix.hpp>   // Matrix<T>, StaticMatrix<T, Rows, Cols>

Header-only. No library linking required.

Constructors

Matrix();
Matrix(size_type rows, size_type cols);
Matrix(size_type rows, size_type cols, const T& value);
Matrix(std::initializer_list<std::initializer_list<T>> init);

Signature	Description
`Matrix()`	Empty matrix ($0 \times 0$)
`Matrix(rows, cols)`	Size-specified (zero-initialized)
`Matrix(rows, cols, value)`	Size-specified, all elements initialized to `value`
`Matrix({{...}, {...}})`	2D initializer list

Parameter	Type	Description
`rows`	`size_type`	Number of rows
`cols`	`size_type`	Number of columns
`value`	`const T&`	Initial value

Matrix<double> A(3, 3, 0.0);         // 3x3 zero matrix
Matrix<double> B = {{1, 2}, {3, 4}}; // 2x2

Element Access

reference operator()(size_type row, size_type col);
const_reference operator()(size_type row, size_type col) const;
reference at(size_type row, size_type col);
const_reference at(size_type row, size_type col) const;
pointer data() noexcept;
const_pointer data() const noexcept;

Method	Return Type	Description
`A(i, j)`	`T&` / `const T&`	Assert-checked in debug builds
`A.at(i, j)`	`T&` / `const T&`	Bounds-checked (throws on out-of-range)
`A.data()`	`T` / `const T`	Raw pointer (row-major contiguous layout)

Size & Capacity

size_type rows() const noexcept;
size_type cols() const noexcept;
size_type size() const noexcept;
bool empty() const noexcept;
bool is_square() const noexcept;
void resize(size_type rows, size_type cols);
void resize(size_type rows, size_type cols, const value_type& value);
void reshape(size_type rows, size_type cols);
void clear() noexcept;

Method	Return Type	Description
`rows()`	`size_type`	Number of rows
`cols()`	`size_type`	Number of columns
`size()`	`size_type`	Total number of elements ($= \text{rows} \times \text{cols}$)
`empty()`	`bool`	Whether the matrix is empty
`is_square()`	`bool`	Whether the matrix is square
`resize(r, c)`	`void`	Resize the matrix
`reshape(r, c)`	`void`	Reshape (total element count must match)
`clear()`	`void`	Clear to empty

Arithmetic Operators

// Matrix + Matrix
Matrix<T> operator+(const Matrix<T>& lhs, const Matrix<T>& rhs);
Matrix<T> operator-(const Matrix<T>& lhs, const Matrix<T>& rhs);
Matrix<T> operator*(const Matrix<T>& lhs, const Matrix<T>& rhs);

// Matrix + Scalar
Matrix<T> operator*(const Matrix<T>& mat, const Scalar& scalar);
Matrix<T> operator*(const Scalar& scalar, const Matrix<T>& mat);
Matrix<T> operator/(const Matrix<T>& mat, const Scalar& scalar);
Matrix<T> operator-(const Matrix<T>& mat);  // unary minus

// Matrix + Vector
Vector<T> operator*(const Matrix<T>& mat, const Vector<T>& vec);

// compound assignment
Matrix& operator+=(const Matrix& rhs);
Matrix& operator-=(const Matrix& rhs);
Matrix& operator*=(const T& scalar);
Matrix& operator/=(const T& scalar);

// operator^ (transpose / Hermitian transpose / inverse / power)
Matrix operator^(int n) const;

Operation	Return Type	Description
`A + B`	`Matrix<T>`	Matrix addition (same size). Expression templates eliminate temporaries
`A - B`	`Matrix<T>`	Matrix subtraction
`A * B`	`Matrix<T>`	Matrix product ($A \in \mathbb{R}^{m \times n}$, $B \in \mathbb{R}^{n \times p}$ $\to$ $C \in \mathbb{R}^{m \times p}$)
`A * s`, `s * A`	`Matrix<T>`	Scalar multiplication
`A * v`	`Vector<T>`	Matrix-vector product $y = Ax$
`-A`	`Matrix<T>`	Negate all elements
`A ^ 'T'`	`Matrix<T>`	Transpose $A^\top$
`A ^ 'H'`	`Matrix<T>`	Hermitian transpose $A^*$ (conjugate transpose)
`A ^ (-1)`	`Matrix<T>`	Inverse $A^{-1}$
`A ^ n`	`Matrix<T>`	Power $A^n$ (binary exponentiation). $n < 0$ computes $(A^{-1})^{\|n\|}$

Note: The ^ operator has low precedence in C++. Use parentheses in expressions: (A ^ 'T') * B

auto and expression templates: When you capture the result of A + B or s * A with auto, the actual type is not Matrix<T> but an intermediate expression template node. Passing this intermediate type to functions like logm() or sqrtm() that expect Matrix<T> will cause a compilation error.

// NG: A is MatBinOp<...>, logm() cannot deduce T
auto A = matmul(B.transpose(), B) + Matrix<double>::identity(n);
auto L = logm(A);  // compile error

// OK: Matrix<double> evaluates the expression
Matrix<double> A = matmul(B.transpose(), B) + Matrix<double>::identity(n);
auto L = logm(A);  // OK

Use Matrix<T> explicitly when passing expression results to other functions.

Matrix<double> A = {{1, 2}, {3, 4}};
Matrix<double> B = {{5, 6}, {7, 8}};
auto C = A * B;  // {{19, 22}, {43, 50}}

Vector<double> v = {1, 2};
auto w = A * v;  // {5, 11}

auto At   = A ^ 'T';   // transpose
auto Ainv = A ^ (-1);   // inverse
auto A3   = A ^ 3;      // A*A*A

Matrix Operations

Matrix transpose() const;
T determinant() const;
[[nodiscard]] Matrix inverse() const;
T trace() const;
T norm() const;
void zero();
void fill(const T& value);
void identity();
static Matrix identity(size_type size);

Method	Return Type	Description
`transpose()`	`Matrix`	Transpose $A^\top$
`determinant()`	`T`	Determinant $\det(A)$ (LU decomposition, $O(n^3)$)
`inverse()`	`Matrix`	Inverse $A^{-1}$ (Gauss-Jordan, $O(n^3)$)
`trace()`	`T`	Trace $\operatorname{tr}(A) = \sum_i a_{ii}$
`norm()`	`T`	Frobenius norm $\\|A\\|_F = \sqrt{\sum_{ij} a_{ij}^2}$
`zero()`	`void`	Set all elements to zero
`fill(value)`	`void`	Set all elements to `value`
`identity()`	`void`	Set to identity matrix (square matrices only)
`identity(n)`	`Matrix`	$n \times n$ identity matrix (static)

Matrix<double> A = {{1, 2}, {3, 4}};
std::cout << A.determinant() << std::endl;  // -2
auto Ainv = A.inverse();
auto I = A * Ainv;  // approximately identity matrix

Low-Level Free Functions

Output-parameter form. Customizable via the ComputePolicy template parameter.

template<...> void add(const MatrixA& a, const MatrixB& b, MatrixC& c);
template<...> void subtract(const MatrixA& a, const MatrixB& b, MatrixC& c);
template<...> void scale(const MatrixA& a, const Scalar& alpha, MatrixB& b);
template<...> void multiply(const MatrixA& a, const MatrixB& b, MatrixC& c);
template<...> void multiply_vector(const MatrixA& a, const VectorX& x, VectorY& y);
template<...> void transpose(const MatrixA& a, MatrixB& b);
template<typename MatrixA> auto trace(const MatrixA& a) -> typename MatrixA::value_type;
template<typename MatrixA> auto frobenius_norm(const MatrixA& a) -> typename MatrixA::value_type;
template<typename MatrixA> auto determinant(const MatrixA& a) -> typename MatrixA::value_type;
template<typename T> Matrix<T> inverse(const Matrix<T>& a);

Row/Column Operations & Views

template<typename VectorType> void get_row(size_type row_idx, VectorType& dest) const;
template<typename VectorType> void get_column(size_type col_idx, VectorType& dest) const;
template<typename VectorType> void set_row(size_type row_idx, const VectorType& src);
template<typename VectorType> void set_column(size_type col_idx, const VectorType& src);
MatrixView<MatrixBase> submatrix(size_type r, size_type c, size_type nr, size_type nc);
ConstMatrixView<MatrixBase> submatrix(size_type r, size_type c, size_type nr, size_type nc) const;
MatrixView<MatrixBase> row(size_type row_idx);
MatrixView<MatrixBase> column(size_type col_idx);

Method	Return Type	Description
`get_row(i, dest)`	`void`	Copy row $i$ into a vector
`get_column(j, dest)`	`void`	Copy column $j$ into a vector
`set_row(i, src)`	`void`	Set row $i$
`set_column(j, src)`	`void`	Set column $j$
`submatrix(r, c, nr, nc)`	`MatrixView`	Sub-matrix view (zero-copy)
`row(i)`	`MatrixView`	View of row $i$
`column(j)`	`MatrixView`	View of column $j$

Factory Functions

template<typename T> Matrix<T> zeros(std::size_t rows, std::size_t cols);
template<typename T> Matrix<T> ones(std::size_t rows, std::size_t cols);
template<typename T> Matrix<T> eye(std::size_t size);
template<typename T> Matrix<T> rotation2D(T angle);
template<typename T> Matrix<T> rotationX(T angle);
template<typename T> Matrix<T> rotationY(T angle);
template<typename T> Matrix<T> rotationZ(T angle);
template<typename T> Matrix<T> hconcat(const Matrix<T>& a, const Matrix<T>& b);
template<typename T> Matrix<T> vconcat(const Matrix<T>& a, const Matrix<T>& b);
template<typename T> [[nodiscard]] Matrix<T> power(const Matrix<T>& base, unsigned int exponent);

Function	Return Type	Description
`zeros<T>(r, c)`	`Matrix<T>`	Zero matrix
`ones<T>(r, c)`	`Matrix<T>`	All elements set to $1$
`eye<T>(n)`	`Matrix<T>`	$n \times n$ identity matrix
`rotation2D(θ)`	`Matrix<T>`	2D rotation matrix ($2 \times 2$)
`rotationX(θ)`	`Matrix<T>`	3D rotation about X-axis ($3 \times 3$)
`rotationY(θ)`	`Matrix<T>`	3D rotation about Y-axis
`rotationZ(θ)`	`Matrix<T>`	3D rotation about Z-axis
`hconcat(A, B)`	`Matrix<T>`	Horizontal concatenation $[A \mid B]$
`vconcat(A, B)`	`Matrix<T>`	Vertical concatenation $\begin{bmatrix}A \\ B\end{bmatrix}$
`power(A, n)`	`Matrix<T>`	Matrix power $A^n$ (binary exponentiation, $O(n^3 \log k)$)

auto I = eye<double>(3);
auto R = rotationZ<double>(M_PI / 4);  // Z-axis 45° rotation

Output

template<typename T> std::ostream& operator<<(std::ostream& os, const Matrix<T>& mat);
template<typename T> std::string to_string(const Matrix<T>& mat, std::size_t precision = 4);
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);

Function	Return Type	Description
`operator<<`	`std::ostream&`	Stream output
`to_string(A, prec)`	`std::string`	String conversion
`print_matrix(A, os, w, p)`	`void`	Formatted output (width and precision)

LaTeX Output

Free functions in the sangi namespace that convert a matrix or vector into a LaTeX-formatted string. These are not member functions.

Matrix LaTeX Output

// style: "pmatrix" (parentheses), "bmatrix" (square brackets), "vmatrix" (vertical bars), "matrix" (no delimiters)
template<typename T>
std::string toLatex(const Matrix<T>& mat,
    const std::string& style = "pmatrix", std::size_t precision = 6);

Vector LaTeX Output (Column Vector)

template<typename T>
std::string toLatex(const Vector<T>& v,
    const std::string& style = "pmatrix", std::size_t precision = 6);

Parameters

Parameter	Description
`style`	LaTeX environment name. `"pmatrix"` = parentheses, `"bmatrix"` = square brackets, `"vmatrix"` = vertical bars, `"matrix"` = no delimiters. Default: `"pmatrix"`
`precision`	Number of decimal places. Default: 6

Examples

Matrix<double> A({{1, 2}, {3, 4}});
std::cout << toLatex(A) << std::endl;
// \begin{pmatrix}
// 1 & 2 \\
// 3 & 4
// \end{pmatrix}

Vector<double> v({1, 2, 3});
std::cout << toLatex(v, "bmatrix") << std::endl;
// \begin{bmatrix}
// 1 \\
// 2 \\
// 3
// \end{bmatrix}

StaticMatrix<T, Rows, Cols>

Compile-time fixed-size variant. Allocated on the stack with a fixed size. Provides the same API as the dynamic Matrix<T>.

StaticMatrix<double, 3, 3> A;
constexpr auto r = A.static_rows();  // 3
constexpr auto c = A.static_cols();  // 3

Specific Methods

static constexpr size_type static_rows();
static constexpr size_type static_cols();
StaticMatrix<T, Cols, Rows> transpose() const;
T determinant() const requires (Rows == Cols);
T trace() const requires (Rows == Cols);
static StaticMatrix identity() requires (Rows == Cols);

Multiplication between StaticMatrix instances is size-checked at compile time: StaticMatrix<T, M, P> = StaticMatrix<T, M, N> * StaticMatrix<T, N, P>.

Examples

Basic Matrix Creation and Operations

#include <math/core/matrix.hpp>
#include <iostream>
using namespace sangi;

int main() {
    // initializer_list construction
    Matrix<double> A = {{1, 2, 3},
                         {4, 5, 6},
                         {7, 8, 10}};

    // transpose, trace, determinant
    auto At = A.transpose();
    std::cout << "trace(A) = " << A.trace() << std::endl;        // 16
    std::cout << "det(A)   = " << A.determinant() << std::endl;  // -3

    // matrix multiplication
    Matrix<double> B = {{1, 0}, {0, 1}, {1, 1}};
    auto C = A * B;  // 3x2 result
    print_matrix(C);

    // scalar operations
    auto D = 2.0 * A - A;  // == A

    // factory functions
    auto I = eye<double>(3);
    auto R = rotationZ<double>(M_PI / 6);  // Z-axis 30 deg rotation
}

Solving Linear Systems with solve()

#include <math/core/matrix.hpp>
#include <math/linalg/solve.hpp>
#include <iostream>
using namespace sangi;

int main() {
    // Solve Ax = b
    Matrix<double> A = {{2, 1, -1},
                         {-3, -1, 2},
                         {-2, 1, 2}};
    Vector<double> b = {8, -11, -3};

    auto x = solve(A, b);
    std::cout << "x = " << x << std::endl;  // {2, 3, -1}

    // verify: A * x == b
    auto residual = A * x - b;
    std::cout << "residual norm = " << norm(residual) << std::endl;  // ~0
}

Exact Inverse of a Hilbert Matrix with Matrix<Rational>

#include <math/core/matrix.hpp>
#include <math/core/mp/Rational.hpp>
#include <iostream>
using namespace sangi;

int main() {
    const int n = 4;
    Matrix<Rational> H(n, n);
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            H(i, j) = Rational(1, i + j + 1);

    // exact inverse -- no floating-point error
    auto Hinv = H.inverse();
    auto I = H * Hinv;

    // verify: every diagonal element is exactly 1
    for (int i = 0; i < n; ++i)
        std::cout << "I(" << i << "," << i << ") = " << I(i, i) << std::endl;
    // Output: 1/1 for all

    print_matrix(Hinv);
}

Related Mathematical Background

The following articles explain the mathematical concepts underlying matrix operations.

Matrix Norms — Definitions and properties of various norms
Gaussian Elimination — Fundamental method for solving linear systems
Gauss-Jordan Elimination — Matrix inversion