matrix_algorithms Namespace Reference

Classes
struct	EigenResult
	Result of eigenvalue decomposition. More...
struct	GivensQRResult
	Result of Givens QR decomposition. More...
struct	LUResult
	Result of LU decomposition with partial pivoting. More...
struct	QRResult
	Result of Householder QR decomposition. More...

Functions
template<typename MatrixType >
Expected< MatrixType, MatrixStatus >	cholesky (const MatrixType &A, typename MatrixType::value_type tol=typename MatrixType::value_type(PRECISION_TOLERANCE))
	Compute the Cholesky decomposition of a symmetric positive-definite matrix.
template<typename MatrixType >
MatrixType	cholesky_or_die (const MatrixType &A)
	Cholesky decomposition — abort on failure.
template<typename T , my_size_t N>
FusedMatrix< T, N, N >	cholesky_rank1_update (const FusedMatrix< T, N, N > &L, const FusedVector< T, N > &v)
	Rank-1 Cholesky update: compute L' where L'L'ᵀ = LLᵀ + vvᵀ.
template<typename T , my_size_t N>
Expected< EigenResult< T, N >, MatrixStatus >	eigen_jacobi (const FusedMatrix< T, N, N > &A, my_size_t max_iters=100, T tol=T(PRECISION_TOLERANCE))
	Compute eigenvalues and eigenvectors of a symmetric matrix via Jacobi.
template<typename T , my_size_t N>
Expected< LUResult< T, N >, MatrixStatus >	lu (const FusedMatrix< T, N, N > &A, T tol=T(PRECISION_TOLERANCE))
	Compute the LU decomposition of a square matrix with partial pivoting.
template<typename T , my_size_t N>
LUResult< T, N >	lu_or_die (const FusedMatrix< T, N, N > &A)
	LU decomposition — abort on failure.
template<typename T , my_size_t M, my_size_t N>
QRResult< T, M, N >	qr_householder (const FusedMatrix< T, M, N > &A)
	Compute the Householder QR decomposition of a rectangular matrix.
template<typename T , my_size_t M, my_size_t N>
GivensQRResult< T, M, N >	qr_givens (const FusedMatrix< T, M, N > &A)
	Compute the QR decomposition using Givens rotations.
template<typename T , my_size_t N, my_size_t M>
Expected< FusedMatrix< T, N, M >, MatrixStatus >	kalman_gain (const FusedMatrix< T, N, N > &P, const FusedMatrix< T, M, N > &H, const FusedMatrix< T, M, M > &R)
	Compute the Kalman gain K = P·Hᵀ·(H·P·Hᵀ + R)⁻¹.
template<typename T , my_size_t N, my_size_t M>
FusedMatrix< T, N, N >	joseph_update (const FusedMatrix< T, N, M > &K, const FusedMatrix< T, M, N > &H, const FusedMatrix< T, N, N > &P, const FusedMatrix< T, M, M > &R)
	Joseph form covariance update: P' = (I-K·H)·P·(I-K·H)ᵀ + K·R·Kᵀ.
template<typename T , my_size_t N>
Expected< T, MatrixStatus >	condition (const FusedMatrix< T, N, N > &A)
	Compute the condition number of A in the 1-norm: cond₁(A) = ‖A‖₁ · ‖A⁻¹‖₁.
template<typename T , my_size_t N>
FusedVector< T, N >	cross (const FusedVector< T, N > &a, const FusedVector< T, N > &b)
	Compute the cross product of two 3-vectors.
template<typename T , my_size_t N>
T	determinant (const FusedMatrix< T, N, N > &A)
	Compute the determinant of a square matrix.
template<typename T >
FusedMatrix< T, 4, 4 >	homogeneous_inverse (const FusedMatrix< T, 4, 4 > &H)
	Compute the inverse of a 4×4 homogeneous transformation matrix.
template<typename T , my_size_t N>
Expected< FusedMatrix< T, N, N >, MatrixStatus >	inverse (const FusedMatrix< T, N, N > &A)
	Compute the inverse of a square matrix.
template<typename T , my_size_t N>
FusedMatrix< T, N, N >	inverse_or_die (const FusedMatrix< T, N, N > &A)
	Matrix inverse — abort on failure.
template<typename T , my_size_t N>
T	norm1 (const FusedMatrix< T, N, N > &A)
	Compute the 1-norm of a square matrix: max absolute column sum.
template<typename T , my_size_t N>
T	norm2 (const FusedVector< T, N > &v)
	Compute the Euclidean (2-norm) of a vector: √(Σ vᵢ²).
template<typename T , my_size_t N>
FusedMatrix< T, N, N >	symmetric_rank_k_update (const FusedMatrix< T, N, N > &F, const FusedMatrix< T, N, N > &P, const FusedMatrix< T, N, N > &Q)
	Compute the symmetric rank-k update P' = F·P·Fᵀ + Q.
template<typename T , my_size_t N>
FusedMatrix< T, N, N >	symmetric_rank_k_update (const FusedMatrix< T, N, N > &F, const FusedMatrix< T, N, N > &P)
	Compute the symmetric rank-k update P' = F·P·Fᵀ (no noise term).
template<typename T >
FusedMatrix< T, 3, 3 >	skew_symmetric (const FusedVector< T, 3 > &v)
	Construct the 3×3 skew-symmetric matrix [v]× from a 3-vector.
template<typename T >
FusedMatrix< T, 3, 3 >	rodrigues (const FusedVector< T, 3 > &omega, T t=T(1))
	Compute the rotation matrix R = exp(t · [ω]×) via Rodrigues formula.
template<typename T , my_size_t N>
T	trace (const FusedMatrix< T, N, N > &A)
	Compute the trace of a square matrix.
template<typename T , my_size_t N>
Expected< FusedVector< T, N >, MatrixStatus >	cholesky_solve (const FusedMatrix< T, N, N > &A, const FusedVector< T, N > &b)
	Solve Ax = b for symmetric positive-definite A via Cholesky decomposition.
template<typename T , my_size_t M, my_size_t N>
Expected< FusedVector< T, N >, MatrixStatus >	least_squares (const FusedMatrix< T, M, N > &A, const FusedVector< T, M > &b)
	Solve the least squares problem min ‖Ax - b‖² via QR.
template<typename T , my_size_t N>
Expected< FusedVector< T, N >, MatrixStatus >	lu_solve (const FusedMatrix< T, N, N > &A, const FusedVector< T, N > &b)
	Solve Ax = b via LU decomposition with partial pivoting.
template<typename T , my_size_t N>
Expected< FusedVector< T, N >, MatrixStatus >	solve (const FusedMatrix< T, N, N > &A, const FusedVector< T, N > &b)
	Solve Ax = b with automatic algorithm selection.
template<typename T , my_size_t N>
Expected< FusedVector< T, N >, MatrixStatus >	levinson_durbin (const FusedVector< T, N > &r, const FusedVector< T, N > &b)
	Solve a symmetric Toeplitz system Tx = b using Levinson-Durbin.
template<bool UnitDiag = false, typename T , my_size_t N>
Expected< FusedVector< T, N >, MatrixStatus >	forward_substitute (const FusedMatrix< T, N, N > &L, const FusedVector< T, N > &b)
	Solve the lower-triangular system Lx = b by forward substitution.
template<bool UnitDiag = false, typename T , my_size_t N>
Expected< FusedVector< T, N >, MatrixStatus >	back_substitute (const FusedMatrix< T, N, N > &U, const FusedVector< T, N > &b)
	Solve the upper-triangular system Ux = b by back substitution.
template<bool UnitDiag = false, typename T , my_size_t N, my_size_t Ncols>
Expected< FusedMatrix< T, N, Ncols >, MatrixStatus >	forward_substitute (const FusedMatrix< T, N, N > &L, const FusedMatrix< T, N, Ncols > &B)
	Solve the lower-triangular system LX = B for multiple right-hand sides.
template<bool UnitDiag = false, typename T , my_size_t N, my_size_t Ncols>
Expected< FusedMatrix< T, N, Ncols >, MatrixStatus >	back_substitute (const FusedMatrix< T, N, N > &U, const FusedMatrix< T, N, Ncols > &B)
	Solve the upper-triangular system UX = B for multiple right-hand sides.
template<typename T , my_size_t N>
Expected< FusedVector< T, N >, MatrixStatus >	thomas_solve (const FusedVector< T, N > &a, const FusedVector< T, N > &d, const FusedVector< T, N > &c, const FusedVector< T, N > &b)
	Solve a tridiagonal system Ax = b using the Thomas algorithm.

Function Documentation

back_substitute() [1/2]

template<bool UnitDiag = false, typename T , my_size_t N, my_size_t Ncols>

Expected< FusedMatrix< T, N, Ncols >, MatrixStatus > matrix_algorithms::back_substitute	(	const FusedMatrix< T, N, N > &	U,
		const FusedMatrix< T, N, Ncols > &	B
	)

Solve the upper-triangular system UX = B for multiple right-hand sides.

Each column of B is solved independently via back substitution. Uses the generic scalar path (no fixed-size unrolling).

Template Parameters

UnitDiag	If true, the diagonal of U is treated as all ones.
T	Scalar type (deduced).
N	System dimension (deduced).
Ncols	Number of right-hand side columns (deduced).

Parameters

U	Upper-triangular NxN matrix.
B	Right-hand side matrix of size N × Ncols.

Returns

Expected containing solution matrix X, or MatrixStatus::Singular on zero diagonal.

Here is the call graph for this function:

back_substitute() [2/2]

template<bool UnitDiag = false, typename T , my_size_t N>

Expected< FusedVector< T, N >, MatrixStatus > matrix_algorithms::back_substitute	(	const FusedMatrix< T, N, N > &	U,
		const FusedVector< T, N > &	b
	)

Solve the upper-triangular system Ux = b by back substitution.

Compile-time dispatch selects a fully unrolled path for N ∈ {3, 4, 6}; all other sizes use a generic scalar loop.

Template Parameters

UnitDiag	If true, the diagonal of U is treated as all ones. If false, the actual diagonal entries are used.
T	Scalar type (deduced).
N	Matrix/vector dimension (deduced).

Parameters

U	Upper-triangular NxN matrix.
b	Right-hand side vector of length N.

Returns

Expected containing solution x, or MatrixStatus::Singular on zero diagonal.

Here is the call graph for this function:

cholesky()

template<typename MatrixType >

Expected< MatrixType, MatrixStatus > matrix_algorithms::cholesky	(	const MatrixType &	A,
		typename MatrixType::value_type	tol = typename MatrixType::value_type(PRECISION_TOLERANCE)
	)

Compute the Cholesky decomposition of a symmetric positive-definite matrix.

Decomposes A into a lower-triangular matrix L such that A = L · Lᵀ.

Template Parameters

MatrixType

A square FusedMatrix type exposing: isSymmetric() — runtime symmetry check getDim(i) — dimension size along axis i operator()(i, j) — element access value_type — scalar type (float, double)

Parameters

A	Symmetric positive-definite input matrix.
tol	Diagonal tolerance. Elements ≤ tol are rejected as non-positive-definite. Defaults to PRECISION_TOLERANCE.

Returns

Expected containing the lower-triangular factor L on success, or MatrixStatus::NotSymmetric / MatrixStatus::NotPositiveDefinite on failure.

Example:

FusedMatrix<double, 3, 3> A;
// ... fill A as SPD ...
auto result = matrix_algorithms::cholesky(A);
if (!result.has_value()) {
    // handle result.error()
    return;
}
auto& L = result.value();
// L * L^T ≈ A

Here is the call graph for this function:

cholesky_or_die()

template<typename MatrixType >

MatrixType matrix_algorithms::cholesky_or_die

(

const MatrixType &

A

)

Cholesky decomposition — abort on failure.

Convenience wrapper for contexts where failure is unrecoverable (offline computation, test code). Calls MyErrorHandler::error() if the decomposition fails.

Template Parameters

MatrixType

Same requirements as cholesky().

Parameters

A	Symmetric positive-definite input matrix.

Returns

The lower-triangular factor L.

Here is the call graph for this function:

cholesky_rank1_update()

template<typename T , my_size_t N>

FusedMatrix< T, N, N > matrix_algorithms::cholesky_rank1_update	(	const FusedMatrix< T, N, N > &	L,
		const FusedVector< T, N > &	v
	)

Rank-1 Cholesky update: compute L' where L'L'ᵀ = LLᵀ + vvᵀ.

Given a lower-triangular Cholesky factor L and a vector v, computes the updated factor L' in O(N²) via Givens rotations. The input L, v is not modified.

Template Parameters

T	Scalar type (deduced).
N	Matrix/vector dimension (deduced).

Parameters

L	Lower-triangular Cholesky factor (N×N).
v	Update vector (N).

Returns

Updated Cholesky factor L' such that L'L'ᵀ = LLᵀ + vvᵀ.

Here is the call graph for this function:

cholesky_solve()

template<typename T , my_size_t N>

Expected< FusedVector< T, N >, MatrixStatus > matrix_algorithms::cholesky_solve	(	const FusedMatrix< T, N, N > &	A,
		const FusedVector< T, N > &	b
	)

Solve Ax = b for symmetric positive-definite A via Cholesky decomposition.

Decomposes A = LLᵀ, then solves Ly = b (forward substitution) followed by Lᵀx = y (back substitution using L transposed).

The internal back substitution step accesses L(k,i) directly (i.e. Lᵀ(i,k)) to avoid creating a separate transpose. Singular checks are skipped for this step since cholesky() already guarantees a valid L with positive diagonal.

Template Parameters

T	Scalar type (deduced).
N	Matrix/vector dimension (deduced).

Parameters

A	Symmetric positive-definite matrix (N×N).
b	Right-hand side vector (N).

Returns

Expected containing solution x on success, or MatrixStatus error forwarded from cholesky/substitution on failure.

Example:

FusedMatrix<double, 3, 3> A;
FusedVector<double, 3> b;
// ... fill A (SPD) and b ...
auto result = matrix_algorithms::cholesky_solve(A, b);
if (!result.has_value()) {
    // handle result.error()
    return;
}
auto& x = result.value();
// A * x ≈ b

Here is the call graph for this function:

condition()

template<typename T , my_size_t N>

Expected< T, MatrixStatus > matrix_algorithms::condition

(

const FusedMatrix< T, N, N > &

A

)

Compute the condition number of A in the 1-norm: cond₁(A) = ‖A‖₁ · ‖A⁻¹‖₁.

Template Parameters

T	Scalar type (deduced).
N	Matrix dimension (deduced).

Parameters

A	Square input matrix (N×N).

Returns

Expected containing the condition number on success, or MatrixStatus::Singular if A is not invertible.

Here is the call graph for this function:

cross()

template<typename T , my_size_t N>

FusedVector< T, N > matrix_algorithms::cross	(	const FusedVector< T, N > &	a,
		const FusedVector< T, N > &	b
	)

Compute the cross product of two 3-vectors.

Template Parameters

T	Scalar type (deduced). Works on float, double, and integers.
N	Vector dimension (deduced). Must be 3.

Parameters

a	First input vector (3×1).
b	Second input vector (3×1).

Returns

a × b as a 3-vector.

determinant()

template<typename T , my_size_t N>

T matrix_algorithms::determinant

(

const FusedMatrix< T, N, N > &

A

)

Compute the determinant of a square matrix.

Template Parameters

T	Scalar type (deduced).
N	Matrix dimension (deduced).

Parameters

A	Square input matrix (N×N).

Returns

Determinant of A. Zero if A is singular (generic path).

Here is the call graph for this function:

eigen_jacobi()

template<typename T , my_size_t N>

Expected< EigenResult< T, N >, MatrixStatus > matrix_algorithms::eigen_jacobi	(	const FusedMatrix< T, N, N > &	A,
		my_size_t	max_iters = 100,
		T	tol = T(PRECISION_TOLERANCE)
	)

Compute eigenvalues and eigenvectors of a symmetric matrix via Jacobi.

Template Parameters

T	Scalar type (deduced).
N	Matrix dimension (deduced).

Parameters

A	Symmetric input matrix (N×N).
max_iters	Maximum number of sweeps (default 100).
tol	Convergence tolerance for off-diagonal norm (default PRECISION_TOLERANCE).

Returns

Expected containing EigenResult on success, or MatrixStatus error on failure.

Example:

FusedMatrix<double, 3, 3> A;
// ... fill A (symmetric) ...
auto result = matrix_algorithms::eigen_jacobi(A);
if (result.has_value()) {
    auto& eig = result.value();
    // eig.eigenvalues(i) — the i-th eigenvalue
    // eig.eigenvectors column i — the i-th eigenvector
    // V * diag(λ) * Vᵀ ≈ A
}

Here is the call graph for this function:

forward_substitute() [1/2]

template<bool UnitDiag = false, typename T , my_size_t N, my_size_t Ncols>

Expected< FusedMatrix< T, N, Ncols >, MatrixStatus > matrix_algorithms::forward_substitute	(	const FusedMatrix< T, N, N > &	L,
		const FusedMatrix< T, N, Ncols > &	B
	)

Solve the lower-triangular system LX = B for multiple right-hand sides.

Each column of B is solved independently via forward substitution. Uses the generic scalar path (no fixed-size unrolling).

Template Parameters

UnitDiag	If true, the diagonal of L is treated as all ones.
T	Scalar type (deduced).
N	System dimension (deduced).
Ncols	Number of right-hand side columns (deduced).

Parameters

L	Lower-triangular NxN matrix.
B	Right-hand side matrix of size N × Ncols.

Returns

Expected containing solution matrix X, or MatrixStatus::Singular on zero diagonal.

Here is the call graph for this function:

forward_substitute() [2/2]

template<bool UnitDiag = false, typename T , my_size_t N>

Expected< FusedVector< T, N >, MatrixStatus > matrix_algorithms::forward_substitute	(	const FusedMatrix< T, N, N > &	L,
		const FusedVector< T, N > &	b
	)

Solve the lower-triangular system Lx = b by forward substitution.

Compile-time dispatch selects a fully unrolled path for N ∈ {3, 4, 6}; all other sizes use a generic scalar loop.

Template Parameters

UnitDiag	If true, the diagonal of L is treated as all ones (implicit unit diagonal, as produced by LU with partial pivoting). If false, the actual diagonal entries are used.
T	Scalar type (deduced).
N	Matrix/vector dimension (deduced).

Parameters

L	Lower-triangular NxN matrix.
b	Right-hand side vector of length N.

Returns

Expected containing solution x, or MatrixStatus::Singular on zero diagonal.

Here is the call graph for this function:

homogeneous_inverse()

template<typename T >

FusedMatrix< T, 4, 4 > matrix_algorithms::homogeneous_inverse

(

const FusedMatrix< T, 4, 4 > &

H

)

Compute the inverse of a 4×4 homogeneous transformation matrix.

Template Parameters

T	Scalar type (deduced). Must be floating-point.

Parameters

H	4×4 homogeneous transform [R|t; 0 0 0 1].

Returns

H⁻¹ = [Rᵀ | -Rᵀt; 0 0 0 1].

inverse()

template<typename T , my_size_t N>

Expected< FusedMatrix< T, N, N >, MatrixStatus > matrix_algorithms::inverse

(

const FusedMatrix< T, N, N > &

A

)

Compute the inverse of a square matrix.

Template Parameters

T	Scalar type (deduced).
N	Matrix dimension (deduced).

Parameters

A	Square input matrix (N×N).

Returns

Expected containing A⁻¹ on success, or MatrixStatus::Singular if A is not invertible.

Example:

FusedMatrix<double, 3, 3> A;
// ... fill A ...
auto result = matrix_algorithms::inverse(A);
if (!result.has_value()) {
    // handle singular matrix
    return;
}
auto& Ainv = result.value();
// A * Ainv ≈ I

Here is the call graph for this function:

inverse_or_die()

template<typename T , my_size_t N>

FusedMatrix< T, N, N > matrix_algorithms::inverse_or_die

(

const FusedMatrix< T, N, N > &

A

)

Matrix inverse — abort on failure.

Convenience wrapper for contexts where failure is unrecoverable.

Template Parameters

T	Scalar type (deduced).
N	Matrix dimension (deduced).

Parameters

A	Square input matrix (N×N).

Returns

A⁻¹.

Here is the call graph for this function:

joseph_update()

template<typename T , my_size_t N, my_size_t M>

FusedMatrix< T, N, N > matrix_algorithms::joseph_update	(	const FusedMatrix< T, N, M > &	K,
		const FusedMatrix< T, M, N > &	H,
		const FusedMatrix< T, N, N > &	P,
		const FusedMatrix< T, M, M > &	R
	)

Joseph form covariance update: P' = (I-K·H)·P·(I-K·H)ᵀ + K·R·Kᵀ.

Numerically stable covariance update that guarantees symmetry and positive semi-definiteness of the result.

Template Parameters

T	Scalar type (deduced).
N	State dimension (deduced).
M	Measurement dimension (deduced).

Parameters

K	Kalman gain (N×M).
H	Observation matrix (M×N).
P	Prior state covariance (N×N).
R	Measurement noise covariance (M×M).

Returns

Updated covariance P' (N×N).

Here is the call graph for this function:

kalman_gain()

template<typename T , my_size_t N, my_size_t M>

Expected< FusedMatrix< T, N, M >, MatrixStatus > matrix_algorithms::kalman_gain	(	const FusedMatrix< T, N, N > &	P,
		const FusedMatrix< T, M, N > &	H,
		const FusedMatrix< T, M, M > &	R
	)

Compute the Kalman gain K = P·Hᵀ·(H·P·Hᵀ + R)⁻¹.

Template Parameters

T	Scalar type (deduced).
N	State dimension (deduced).
M	Measurement dimension (deduced).

Parameters

P	State covariance (N×N), symmetric positive definite.
H	Observation matrix (M×N).
R	Measurement noise covariance (M×M), symmetric positive definite.

Returns

Expected containing the Kalman gain K (N×M) on success, or MatrixStatus::Singular if the innovation covariance is not invertible.

Here is the call graph for this function:

least_squares()

template<typename T , my_size_t M, my_size_t N>

Expected< FusedVector< T, N >, MatrixStatus > matrix_algorithms::least_squares	(	const FusedMatrix< T, M, N > &	A,
		const FusedVector< T, M > &	b
	)

Solve the least squares problem min ‖Ax - b‖² via QR.

Template Parameters

T	Scalar type (deduced).
M	Number of rows in A (deduced, M ≥ N).
N	Number of columns in A (deduced).

Parameters

A	Input matrix (M×N).
b	Right-hand side vector (M).

Returns

Expected containing the least squares solution x (N) on success, or MatrixStatus::Singular if A is rank-deficient.

Example:

// Fit y = a + b*x to 4 data points
FusedMatrix<double, 4, 2> A;  // [1 x0; 1 x1; 1 x2; 1 x3]
FusedVector<double, 4> b;     // [y0; y1; y2; y3]
auto result = matrix_algorithms::least_squares(A, b);
if (result.has_value()) {
    auto& x = result.value();  // x(0)=intercept, x(1)=slope
}

Here is the call graph for this function:

levinson_durbin()

template<typename T , my_size_t N>

Expected< FusedVector< T, N >, MatrixStatus > matrix_algorithms::levinson_durbin	(	const FusedVector< T, N > &	r,
		const FusedVector< T, N > &	b
	)

Solve a symmetric Toeplitz system Tx = b using Levinson-Durbin.

Template Parameters

T	Scalar type (deduced).
N	System size (deduced).

Parameters

r	First row of the Toeplitz matrix (N elements). r(0) is the diagonal.
b	Right-hand side vector (N).

Returns

Expected containing solution x on success, or MatrixStatus error on failure.

Example:

// Solve [4 2 1; 2 4 2; 1 2 4] x = b
FusedVector<double, 3> r, b;
r(0) = 4; r(1) = 2; r(2) = 1;  // first row defines entire matrix
// ... fill b ...
auto result = matrix_algorithms::levinson_durbin(r, b);

Here is the call graph for this function:

lu()

template<typename T , my_size_t N>

Expected< LUResult< T, N >, MatrixStatus > matrix_algorithms::lu	(	const FusedMatrix< T, N, N > &	A,
		T	tol = T(PRECISION_TOLERANCE)
	)

Compute the LU decomposition of a square matrix with partial pivoting.

Decomposes A into P·A = L·U where L has unit diagonal.

Template Parameters

T	Scalar type (deduced).
N	Matrix dimension (deduced).

Parameters

A	Square input matrix (N×N).
tol	Pivot tolerance. Pivots with |value| ≤ tol are rejected as singular. Defaults to PRECISION_TOLERANCE.

Returns

Expected containing LUResult on success, or MatrixStatus::Singular on zero pivot.

Example:

FusedMatrix<double, 3, 3> A;
// ... fill A ...
auto result = matrix_algorithms::lu(A);
if (!result.has_value()) {
    // handle result.error()
    return;
}
auto& lu = result.value();
// lu.LU contains compact factorization
// lu.perm contains row permutation
// lu.sign contains permutation sign
auto L = lu.L();  // extract L if needed
auto U = lu.U();  // extract U if needed

Here is the call graph for this function:

lu_or_die()

template<typename T , my_size_t N>

LUResult< T, N > matrix_algorithms::lu_or_die

(

const FusedMatrix< T, N, N > &

A

)

LU decomposition — abort on failure.

Convenience wrapper for contexts where failure is unrecoverable.

Template Parameters

T	Scalar type (deduced).
N	Matrix dimension (deduced).

Parameters

A	Square input matrix (N×N).

Returns

LUResult containing factorization.

Here is the call graph for this function:

lu_solve()

template<typename T , my_size_t N>

Expected< FusedVector< T, N >, MatrixStatus > matrix_algorithms::lu_solve	(	const FusedMatrix< T, N, N > &	A,
		const FusedVector< T, N > &	b
	)

Solve Ax = b via LU decomposition with partial pivoting.

Works for any non-singular square matrix.

Template Parameters

T	Scalar type (deduced).
N	Matrix/vector dimension (deduced).

Parameters

A	Square input matrix (N×N).
b	Right-hand side vector (N).

Returns

Expected containing solution x on success, or MatrixStatus error on failure.

Here is the call graph for this function:

norm1()

template<typename T , my_size_t N>

T matrix_algorithms::norm1

(

const FusedMatrix< T, N, N > &

A

)

Compute the 1-norm of a square matrix: max absolute column sum.

Works for any numeric type (not restricted to floating point).

Template Parameters

T	Scalar type (deduced).
N	Matrix dimension (deduced).

Parameters

A	Square input matrix (N×N).

Returns

‖A‖₁ = max_j Σ_i |A(i,j)|.

Here is the call graph for this function:

norm2()

template<typename T , my_size_t N>

T matrix_algorithms::norm2

(

const FusedVector< T, N > &

v

)

Compute the Euclidean (2-norm) of a vector: √(Σ vᵢ²).

Template Parameters

T	Scalar type (deduced).
N	Vector dimension (deduced).

Parameters

v	Input vector (N).

Returns

‖v‖₂ = √(v(0)² + v(1)² + … + v(N−1)²).

Note

Future optimization: with a strided diagonal view, this could use a SIMD dot product kernel (dot(v, v) then sqrt).

Here is the call graph for this function:

qr_givens()

template<typename T , my_size_t M, my_size_t N>

GivensQRResult< T, M, N > matrix_algorithms::qr_givens

(

const FusedMatrix< T, M, N > &

A

)

Compute the QR decomposition using Givens rotations.

Template Parameters

T	Scalar type (deduced).
M	Number of rows (deduced, M ≥ N).
N	Number of columns (deduced).

Parameters

A	Input matrix (M×N).

Returns

GivensQRResult containing Q (M×M) and R (M×N).

Example:

FusedMatrix<double, 4, 3> A;
// ... fill A ...
auto qr = matrix_algorithms::qr_givens(A);
// qr.Q * qr.R ≈ A
// qr.Q is orthogonal: Qᵀ·Q = I

Here is the call graph for this function:

qr_householder()

template<typename T , my_size_t M, my_size_t N>

QRResult< T, M, N > matrix_algorithms::qr_householder

(

const FusedMatrix< T, M, N > &

A

)

Compute the Householder QR decomposition of a rectangular matrix.

Template Parameters

T	Scalar type (deduced).
M	Number of rows (deduced, M ≥ N).
N	Number of columns (deduced).

Parameters

A	Input matrix (M×N).

Returns

QRResult containing compact factorization with Q(), R(), and apply_Qt().

Example:

FusedMatrix<double, 4, 3> A;
// ... fill A ...
auto qr = matrix_algorithms::qr_householder(A);
auto Q = qr.Q();     // 4×4 orthogonal
auto R = qr.R();     // 4×3 upper triangular
// Q * R ≈ A
 
// For least squares (without forming Q):
auto Qtb = qr.apply_Qt(b);  // apply Qᵀ to b
// then back-substitute top N rows of Qtb with R(0:N-1, 0:N-1)

Here is the call graph for this function:

rodrigues()

template<typename T >

FusedMatrix< T, 3, 3 > matrix_algorithms::rodrigues	(	const FusedVector< T, 3 > &	omega,
		T	t = T(1)
	)

Compute the rotation matrix R = exp(t · [ω]×) via Rodrigues formula.

Computes the SO(3) rotation matrix for angular velocity ω over time step t.

Template Parameters

T	Scalar type (deduced).

Parameters

omega	Angular velocity 3-vector (rad/s).
t	Time step (seconds). Default 1.0 (omega is then the rotation vector).

Returns

3×3 rotation matrix R ∈ SO(3).

Here is the call graph for this function:

skew_symmetric()

template<typename T >

FusedMatrix< T, 3, 3 > matrix_algorithms::skew_symmetric

(

const FusedVector< T, 3 > &

v

)

Construct the 3×3 skew-symmetric matrix [v]× from a 3-vector.

Template Parameters

T	Scalar type (deduced).

Parameters

v	3-vector [v₁, v₂, v₃].

Returns

3×3 skew-symmetric matrix such that [v]× · u = v × u.

solve()

template<typename T , my_size_t N>

Expected< FusedVector< T, N >, MatrixStatus > matrix_algorithms::solve	(	const FusedMatrix< T, N, N > &	A,
		const FusedVector< T, N > &	b
	)

Solve Ax = b with automatic algorithm selection.

Attempts Cholesky path first (cheaper for SPD matrices), falls back to LU path for general non-singular matrices.

Template Parameters

T	Scalar type (deduced).
N	Matrix/vector dimension (deduced).

Parameters

A	Square input matrix (N×N).
b	Right-hand side vector (N).

Returns

Expected containing solution x on success, or MatrixStatus error on failure.

Example:

FusedMatrix<double, 3, 3> A;
FusedVector<double, 3> b;
// ... fill A and b ...
auto result = matrix_algorithms::solve(A, b);
if (!result.has_value()) {
    // handle result.error()
    return;
}
auto& x = result.value();
// A * x ≈ b

Here is the call graph for this function:

symmetric_rank_k_update() [1/2]

template<typename T , my_size_t N>

FusedMatrix< T, N, N > matrix_algorithms::symmetric_rank_k_update	(	const FusedMatrix< T, N, N > &	F,
		const FusedMatrix< T, N, N > &	P
	)

Compute the symmetric rank-k update P' = F·P·Fᵀ (no noise term).

Useful when process noise is added separately or is zero.

Template Parameters

T	Scalar type (deduced).
N	Matrix dimension (deduced).

Parameters

F	State transition matrix (N×N).
P	Covariance matrix (N×N), symmetric.

Returns

P' = F·P·Fᵀ.

Here is the call graph for this function:

symmetric_rank_k_update() [2/2]

template<typename T , my_size_t N>

FusedMatrix< T, N, N > matrix_algorithms::symmetric_rank_k_update	(	const FusedMatrix< T, N, N > &	F,
		const FusedMatrix< T, N, N > &	P,
		const FusedMatrix< T, N, N > &	Q
	)

Compute the symmetric rank-k update P' = F·P·Fᵀ + Q.

Template Parameters

T	Scalar type (deduced).
N	Matrix dimension (deduced).

Parameters

F	State transition matrix (N×N).
P	Covariance matrix (N×N), symmetric positive (semi-)definite.
Q	Process noise matrix (N×N), symmetric positive (semi-)definite.

Returns

P' = F·P·Fᵀ + Q.

Here is the call graph for this function:

thomas_solve()

template<typename T , my_size_t N>

Expected< FusedVector< T, N >, MatrixStatus > matrix_algorithms::thomas_solve	(	const FusedVector< T, N > &	a,
		const FusedVector< T, N > &	d,
		const FusedVector< T, N > &	c,
		const FusedVector< T, N > &	b
	)

Solve a tridiagonal system Ax = b using the Thomas algorithm.

Template Parameters

T	Scalar type (deduced).
N	System size (deduced).

Parameters

a	Sub-diagonal vector (N). a(0) is unused.
d	Main diagonal vector (N).
c	Super-diagonal vector (N). c(N−1) is unused.
b	Right-hand side vector (N).

Returns

Expected containing solution x on success, or MatrixStatus::Singular on zero pivot.

Example:

// Solve [-2 1 0; 1 -2 1; 0 1 -2] x = b
FusedVector<double, 3> a, d, c, b;
a(0) = 0; a(1) = 1; a(2) = 1;   // sub-diagonal
d(0) = -2; d(1) = -2; d(2) = -2; // main diagonal
c(0) = 1; c(1) = 1; c(2) = 0;   // super-diagonal
// ... fill b ...
auto result = matrix_algorithms::thomas_solve(a, d, c, b);

Here is the call graph for this function:

trace()

template<typename T , my_size_t N>

T matrix_algorithms::trace

(

const FusedMatrix< T, N, N > &

A

)

Compute the trace of a square matrix.

Template Parameters

T	Scalar type (deduced).
N	Matrix dimension (deduced).

Parameters

A	Square input matrix (N×N).

Returns

Sum of diagonal elements: Σ A(i,i) for i = 0 … N−1.

Here is the call graph for this function: