Struct Matrix

Source

pub struct Matrix<const D: usize> { /* private fields */ }

Expand description

Fixed-size square matrix D×D, stored inline.

Implementations§

Source §

impl<const D: usize> Matrix<D>

Source

pub fn det_sign_exact(&self) -> Result<i8, LaError>

Exact determinant sign using adaptive-precision arithmetic.

Returns 1 if det > 0, -1 if det < 0, and 0 if det == 0 (singular).

For D ≤ 4, a fast f64 filter is tried first: det_direct() is compared against a conservative error bound derived from the matrix permanent. If the f64 result clearly exceeds the bound, the sign is returned immediately without allocating. Otherwise (and always for D ≥ 5), the Bareiss algorithm runs in exact BigRational arithmetic.

§When to use

Use this when the sign of the determinant must be correct regardless of floating-point conditioning (e.g. geometric predicates on near-degenerate configurations). For well-conditioned matrices the fast filter resolves the sign without touching BigRational, so the overhead is minimal.

§Examples

use la_stack::prelude::*;

let m = Matrix::<3>::from_rows([
    [1.0, 2.0, 3.0],
    [4.0, 5.0, 6.0],
    [7.0, 8.0, 9.0],
]);
// This matrix is singular (row 3 = row 1 + row 2 in exact arithmetic).
assert_eq!(m.det_sign_exact().unwrap(), 0);

assert_eq!(Matrix::<3>::identity().det_sign_exact().unwrap(), 1);

§Errors

Returns LaError::NonFinite if any matrix entry is NaN or infinite.

Source §

impl<const D: usize> Matrix<D>

Source

pub const fn from_rows(rows: [[f64; D]; D]) -> Self

Construct from row-major storage.

§Examples

use la_stack::prelude::*;

let m = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]);
assert_eq!(m.get(0, 1), Some(2.0));

Source

pub const fn zero() -> Self

All-zeros matrix.

§Examples

use la_stack::prelude::*;

let z = Matrix::<2>::zero();
assert_eq!(z.get(1, 1), Some(0.0));

Source

pub const fn identity() -> Self

Identity matrix.

§Examples

use la_stack::prelude::*;

let i = Matrix::<3>::identity();
assert_eq!(i.get(0, 0), Some(1.0));
assert_eq!(i.get(0, 1), Some(0.0));
assert_eq!(i.get(2, 2), Some(1.0));

Source

pub const fn get(&self, r: usize, c: usize) -> Option<f64>

Get an element with bounds checking.

§Examples

use la_stack::prelude::*;

let m = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]);
assert_eq!(m.get(1, 0), Some(3.0));
assert_eq!(m.get(2, 0), None);

Source

pub const fn set(&mut self, r: usize, c: usize, value: f64) -> bool

Set an element with bounds checking.

Returns true if the index was in-bounds.

§Examples

use la_stack::prelude::*;

let mut m = Matrix::<2>::zero();
assert!(m.set(0, 1, 2.5));
assert_eq!(m.get(0, 1), Some(2.5));
assert!(!m.set(10, 0, 1.0));

Source

pub fn inf_norm(&self) -> f64

Infinity norm (maximum absolute row sum).

§Examples

use la_stack::prelude::*;

let m = Matrix::<2>::from_rows([[1.0, -2.0], [3.0, 4.0]]);
assert!((m.inf_norm() - 7.0).abs() <= 1e-12);

Source

pub fn lu(self, tol: f64) -> Result<Lu<D>, LaError>

Compute an LU decomposition with partial pivoting.

§Examples

use la_stack::prelude::*;

let a = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]);
let lu = a.lu(DEFAULT_PIVOT_TOL)?;

let b = Vector::<2>::new([5.0, 11.0]);
let x = lu.solve_vec(b)?.into_array();

assert!((x[0] - 1.0).abs() <= 1e-12);
assert!((x[1] - 2.0).abs() <= 1e-12);

§Errors

Returns LaError::Singular if, for some column k, the largest-magnitude candidate pivot in that column satisfies |pivot| <= tol (so no numerically usable pivot exists). Returns LaError::NonFinite if NaN/∞ is detected during factorization.

Source

pub fn ldlt(self, tol: f64) -> Result<Ldlt<D>, LaError>

Compute an LDLT factorization (A = L D Lᵀ) without pivoting.

This is intended for symmetric positive definite (SPD) and positive semi-definite (PSD) matrices such as Gram matrices.

§Examples

use la_stack::prelude::*;

let a = Matrix::<2>::from_rows([[4.0, 2.0], [2.0, 3.0]]);
let ldlt = a.ldlt(DEFAULT_SINGULAR_TOL)?;

// det(A) = 8
assert!((ldlt.det() - 8.0).abs() <= 1e-12);

// Solve A x = b
let b = Vector::<2>::new([1.0, 2.0]);
let x = ldlt.solve_vec(b)?.into_array();
assert!((x[0] - (-0.125)).abs() <= 1e-12);
assert!((x[1] - 0.75).abs() <= 1e-12);

§Errors

Returns LaError::Singular if, for some step k, the required diagonal entry d = D[k,k] is <= tol (non-positive or too small). This treats PSD degeneracy (and indefinite inputs) as singular/degenerate. Returns LaError::NonFinite if NaN/∞ is detected during factorization.

Source

pub const fn det_direct(&self) -> Option<f64>

Closed-form determinant for dimensions 0–4, bypassing LU factorization.

Returns Some(det) for D ∈ {0, 1, 2, 3, 4}, None for D ≥ 5. D = 0 returns Some(1.0) (empty product). This is a const fn (Rust 1.94+) and uses fused multiply-add (mul_add) for improved accuracy and performance.

For a determinant that works for any dimension (falling back to LU for D ≥ 5), use det.

§Examples

use la_stack::prelude::*;

let m = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]);
assert!((m.det_direct().unwrap() - (-2.0)).abs() <= 1e-12);

// D = 0 is the empty product.
assert_eq!(Matrix::<0>::zero().det_direct(), Some(1.0));

// D ≥ 5 returns None.
assert!(Matrix::<5>::identity().det_direct().is_none());

Source

pub fn det(self, tol: f64) -> Result<f64, LaError>

Determinant, using closed-form formulas for D ≤ 4 and LU decomposition for D ≥ 5.

For D ∈ {1, 2, 3, 4}, this bypasses LU factorization entirely for a significant speedup (see det_direct). The tol parameter is only used by the LU fallback path for D ≥ 5.

§Examples

use la_stack::prelude::*;

let det = Matrix::<3>::identity().det(DEFAULT_PIVOT_TOL)?;
assert!((det - 1.0).abs() <= 1e-12);

§Errors

Returns LaError::NonFinite if the result contains NaN or infinity. For D ≥ 5, propagates LU factorization errors (e.g. LaError::Singular).

Source

pub fn det_errbound(&self) -> Option<f64>

Conservative absolute error bound for det_direct().

Returns Some(bound) such that |det_direct() - det_exact| ≤ bound, or None for D ≥ 5 where no fast bound is available.

For D ≤ 4, the bound is derived from the matrix permanent using Shewchuk-style error analysis. For D = 0 or 1, returns Some(0.0) since the determinant computation is exact (no arithmetic).

This method does NOT require the exact feature — the bounds use pure f64 arithmetic and are useful for custom adaptive-precision logic.

§When to use

Use this to build adaptive-precision logic: if |det_direct()| > bound, the f64 sign is provably correct. Otherwise fall back to exact arithmetic.

§Examples

use la_stack::prelude::*;

let m = Matrix::<3>::from_rows([
    [1.0, 2.0, 3.0],
    [4.0, 5.0, 6.0],
    [7.0, 8.0, 9.0],
]);
let bound = m.det_errbound().unwrap();
let det_approx = m.det_direct().unwrap();
// If |det_approx| > bound, the sign is guaranteed correct.

§Adaptive precision pattern (requires `exact` feature)

use la_stack::prelude::*;

let m = Matrix::<3>::identity();
if let Some(bound) = m.det_errbound() {
    let det = m.det_direct().unwrap();
    if det.abs() > bound {
        // f64 sign is guaranteed correct
        let sign = det.signum() as i8;
    } else {
        // Fall back to exact arithmetic (requires `exact` feature)
        let sign = m.det_sign_exact().unwrap();
    }
} else {
    // D ≥ 5: no fast filter, use exact directly
    let sign = m.det_sign_exact().unwrap();
}