Struct Matrix 

pub struct Matrix<const ROWS: usize, const COLS: usize, T = f64> { /* private fields */ }

A ROWS×COLS matrix stored inline on the stack in row-major order.

use multicalc::linear_algebra::{Matrix, Vector};
let a = Matrix::new([[1.0, 2.0], [3.0, 4.0]]);
let b = Matrix::new([[5.0, 6.0], [7.0, 8.0]]);

assert_eq!(a[(0, 1)], 2.0);
assert_eq!((a + b).into_array(), [[6.0, 8.0], [10.0, 12.0]]);
assert_eq!((b - a).into_array(), [[4.0, 4.0], [4.0, 4.0]]);
assert_eq!((-a).into_array(), [[-1.0, -2.0], [-3.0, -4.0]]);
assert_eq!((a * 2.0).into_array(), [[2.0, 4.0], [6.0, 8.0]]);
assert_eq!((a * b).into_array(), [[19.0, 22.0], [43.0, 50.0]]);
assert_eq!(a * Vector::new([1.0, 1.0]), Vector::new([3.0, 7.0]));

Implementations

impl<const N: usize, T: Numeric> Matrix<N, N, T>

pub fn cholesky(self) -> Result<Cholesky<N, T>, CalcError>

Factorizes self as L·Lᵀ by the Cholesky–Banachiewicz algorithm.

Only the lower triangle is read; self is assumed symmetric. Returns CalcError::NotPositiveDefinite if a diagonal radicand is not strictly positive — the matrix is not positive definite — rather than taking a root of it.

use multicalc::linear_algebra::Matrix;
let a = Matrix::<3, 3>::new([[4.0, 12.0, -16.0], [12.0, 37.0, -43.0], [-16.0, -43.0, 98.0]]);
let l = a.cholesky().unwrap().l();
// L·Lᵀ == A.
let prod = l * l.transpose();
for r in 0..3 {
    for c in 0..3 {
        assert!((prod[(r, c)] - a[(r, c)]).abs() < 1e-12);
    }
}

impl<const N: usize, T: Numeric> Matrix<N, N, T>

pub fn lu(self) -> Result<Lu<N, T>, CalcError>

Factorizes self by Doolittle LU with partial pivoting.

Returns CalcError::SingularMatrix if a pivot column is entirely zero — the largest available pivot is zero — rather than dividing by it.

use multicalc::linear_algebra::Matrix;
let a = Matrix::<3, 3>::new([[2.0, 1.0, 1.0], [4.0, 3.0, 3.0], [8.0, 7.0, 9.0]]);
let f = a.lu().unwrap();
// P·A == L·U.
let (l, u, perm) = (f.l(), f.u(), f.permutation());
let pa = Matrix::<3, 3>::from_fn(|i, c| a[(perm[i], c)]);
let prod = l * u;
for i in 0..3 {
    for c in 0..3 {
        assert!((pa[(i, c)] - prod[(i, c)]).abs() < 1e-12);
    }
}

pub fn solve(self, b: Vector<N, T>) -> Result<Vector<N, T>, CalcError>

Solves A·x = b for x, factorizing self by LU.

A one-call convenience over Matrix::lu followed by Lu::solve. Returns CalcError::SingularMatrix if self is singular. To solve several right-hand sides, factor once with Matrix::lu and reuse the result. For a symmetric positive-definite matrix, Matrix::cholesky is faster.

use multicalc::linear_algebra::{Matrix, Vector};
let a = Matrix::<3, 3>::new([[2.0, 1.0, 1.0], [4.0, 3.0, 3.0], [8.0, 7.0, 9.0]]);
let x = a.solve(Vector::new([7.0, 19.0, 49.0])).unwrap();
assert!((x[0] - 1.0).abs() < 1e-12);
assert!((x[1] - 2.0).abs() < 1e-12);
assert!((x[2] - 3.0).abs() < 1e-12);

impl<const ROWS: usize, const COLS: usize, T> Matrix<ROWS, COLS, T>

pub const fn new(data: [[T; COLS]; ROWS]) -> Self

Wraps a row-major array of rows into a matrix.

pub fn from_fn(f: impl FnMut(usize, usize) -> T) -> Self

Builds a matrix by calling f with each (row, column) index.

use multicalc::linear_algebra::Matrix;
let m = Matrix::<2, 2>::from_fn(|r, c| (r * 2 + c) as f64);
assert_eq!(m.into_array(), [[0.0, 1.0], [2.0, 3.0]]);

pub const fn as_slice_rows(&self) -> &[[T; COLS]; ROWS]

Borrows the rows.

use multicalc::linear_algebra::Matrix;
assert_eq!(Matrix::new([[1.0, 2.0]]).as_slice_rows(), &[[1.0, 2.0]]);

pub fn into_array(self) -> [[T; COLS]; ROWS]

Consumes the matrix, returning its rows.

impl<const ROWS: usize, const COLS: usize, T: Copy> Matrix<ROWS, COLS, T>

pub fn try_from_row_slice(slice: &[T]) -> Option<Self>

Builds a matrix from a row-major slice, or None if slice.len() is not ROWS * COLS.

use multicalc::linear_algebra::Matrix;
assert!(Matrix::<2, 2>::try_from_row_slice(&[1.0, 2.0, 3.0, 4.0]).is_some());
assert!(Matrix::<2, 2>::try_from_row_slice(&[1.0, 2.0, 3.0]).is_none());

pub fn row(&self, r: usize) -> Vector<COLS, T>

Copies row r into a vector. Panics if r >= ROWS.

use multicalc::linear_algebra::{Matrix, Vector};
let m = Matrix::new([[1.0, 2.0], [3.0, 4.0]]);
assert_eq!(m.row(1), Vector::new([3.0, 4.0]));

pub fn column(&self, c: usize) -> Vector<ROWS, T>

Copies column c into a vector. Panics if c >= COLS.

use multicalc::linear_algebra::{Matrix, Vector};
let m = Matrix::new([[1.0, 2.0], [3.0, 4.0]]);
assert_eq!(m.column(0), Vector::new([1.0, 3.0]));

impl<const ROWS: usize, const COLS: usize, T: Numeric> Matrix<ROWS, COLS, T>

pub fn zeros() -> Self

The zero matrix.

use multicalc::linear_algebra::Matrix;
let m: Matrix<2, 3> = Matrix::zeros();
assert_eq!(m.into_array(), [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]);

pub fn scale(self, scalar: T) -> Self

Multiplies every element by scalar.

use multicalc::linear_algebra::Matrix;
assert_eq!(Matrix::new([[1.0, 2.0]]).scale(3.0).into_array(), [[3.0, 6.0]]);

pub fn transpose(self) -> Matrix<COLS, ROWS, T>

The transpose, with rows and columns swapped.

use multicalc::linear_algebra::Matrix;
let m = Matrix::new([[1.0, 2.0, 3.0]]);
assert_eq!(m.transpose().into_array(), [[1.0], [2.0], [3.0]]);

impl<const N: usize, T: Numeric> Matrix<N, N, T>

pub fn identity() -> Self

The N×N identity matrix.

use multicalc::linear_algebra::Matrix;
let i: Matrix<3, 3> = Matrix::identity();
assert_eq!(i.into_array(), [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]);

impl<T: Numeric> Matrix<2, 2, T>

pub fn determinant(self) -> T

The determinant m[(0,0)]*m[(1,1)] - m[(0,1)]*m[(1,0)].

use multicalc::linear_algebra::Matrix;
assert_eq!(Matrix::new([[1.0, 2.0], [3.0, 4.0]]).determinant(), -2.0);

pub fn inverse(self) -> Result<Self, CalcError>

The inverse, or CalcError::SingularMatrix if the matrix is singular (determinant() == 0).

use multicalc::linear_algebra::Matrix;
let m: Matrix<2, 2> = Matrix::new([[4.0, 7.0], [2.0, 6.0]]);
let p = m * m.inverse().unwrap();
assert!((p[(0, 0)] - 1.0).abs() < 1e-12 && (p[(1, 1)] - 1.0).abs() < 1e-12);
assert!(Matrix::<2, 2>::new([[1.0, 2.0], [2.0, 4.0]]).inverse().is_err());

impl<T: Numeric> Matrix<3, 3, T>

pub fn determinant(self) -> T

The determinant, by cofactor expansion along the first row.

use multicalc::linear_algebra::Matrix;
let i: Matrix<3, 3> = Matrix::identity();
assert_eq!(i.determinant(), 1.0);

pub fn inverse(self) -> Result<Self, CalcError>

The inverse, or CalcError::SingularMatrix if the matrix is singular (determinant() == 0).

use multicalc::linear_algebra::Matrix;
let i: Matrix<3, 3> = Matrix::identity();
assert_eq!(i.inverse().unwrap().into_array(), i.into_array());

impl<T: Numeric> Matrix<4, 4, T>

pub fn determinant(self) -> T

The determinant, as the Laplace expansion along the first two rows.

use multicalc::linear_algebra::Matrix;
let m = Matrix::<4, 4>::new([
    [2.0, 1.0, 1.0, 1.0],
    [0.0, 3.0, 1.0, 1.0],
    [0.0, 0.0, 4.0, 1.0],
    [0.0, 0.0, 0.0, 5.0],
]);
assert_eq!(m.determinant(), 120.0);

pub fn inverse(self) -> Result<Self, CalcError>

The inverse, or CalcError::SingularMatrix if the matrix is singular (determinant() == 0).

Built from the adjugate, the transpose of the cofactor matrix, scaled by 1/det.

use multicalc::linear_algebra::Matrix;
let m = Matrix::<4, 4>::new([
    [2.0, 1.0, 1.0, 1.0],
    [0.0, 3.0, 1.0, 1.0],
    [0.0, 0.0, 4.0, 1.0],
    [0.0, 0.0, 0.0, 5.0],
]);
let p = m * m.inverse().unwrap();
for r in 0..4 {
    for c in 0..4 {
        let expected = if r == c { 1.0 } else { 0.0 };
        assert!((p[(r, c)] - expected).abs() < 1e-12);
    }
}

impl<const M: usize, const N: usize, T: Numeric> Matrix<M, N, T>

pub fn svd(self) -> Result<Svd<M, N, T>, CalcError>

Decomposes self as U · diag(σ) · Vᵀ by one-sided Jacobi (thin form, M ≥ N).

U has orthonormal columns, the σ are non-negative and descending, and V has orthonormal columns. Returns CalcError::Underdetermined for a wide matrix (M < N) — transpose it first — or CalcError::NonFiniteValue if any entry is not finite.

use multicalc::linear_algebra::Matrix;
let a = Matrix::<3, 2>::new([[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]]);
let svd = a.svd().unwrap();
let (u, s, v) = (svd.u(), svd.singular_values(), svd.v());
// U · diag(σ) · Vᵀ == A.
for r in 0..3 {
    for c in 0..2 {
        let mut acc = 0.0;
        for k in 0..2 {
            acc += u[(r, k)] * s[k] * v[(c, k)];
        }
        assert!((acc - a[(r, c)]).abs() < 1e-12);
    }
}

A wide matrix (M < N) has no thin form here; take the SVD of its transpose, whose singular values are the same. Its pseudo-inverse then follows from A⁺ = ((Aᵀ)⁺)ᵀ, which Matrix::pseudo_inverse applies for any shape.

use multicalc::linear_algebra::Matrix;
// For a wide matrix, decompose its transpose.
let a = Matrix::<2, 3>::new([[1.0, 0.0, 2.0], [0.0, 1.0, 1.0]]);
let at = a.transpose();
let svd = at.svd().unwrap();
let (u, s, v) = (svd.u(), svd.singular_values(), svd.v());
for r in 0..3 {
    for c in 0..2 {
        let mut acc = 0.0;
        for k in 0..2 {
            acc += u[(r, k)] * s[k] * v[(c, k)];
        }
        assert!((acc - at[(r, c)]).abs() < 1e-12);
    }
}

pub fn pseudo_inverse(self) -> Result<Matrix<N, M, T>, CalcError>

The Moore–Penrose pseudo-inverse of self, for any shape.

Tall or square inputs go straight through Matrix::svd; a wide input (M < N) is handled as ((Aᵀ)⁺)ᵀ. Returns CalcError::NonFiniteValue if any entry is not finite.

use multicalc::linear_algebra::Matrix;
// A wide matrix, handled through the transpose route.
let a = Matrix::<2, 3>::new([[1.0, 0.0, 2.0], [0.0, 1.0, 1.0]]);
let pinv = a.pseudo_inverse().unwrap();
let recon = a * pinv * a; // A·A⁺·A == A
for r in 0..2 {
    for c in 0..3 {
        assert!((recon[(r, c)] - a[(r, c)]).abs() < 1e-12);
    }
}

Trait Implementations

impl<const ROWS: usize, const COLS: usize, T: Numeric> Add for Matrix<ROWS, COLS, T>

type Output = Matrix<ROWS, COLS, T>

The resulting type after applying the + operator.

fn add(self, rhs: Self) -> Self

Performs the + operation.

impl<const ROWS: usize, const COLS: usize, T: Numeric> AddAssign for Matrix<ROWS, COLS, T>

fn add_assign(&mut self, rhs: Self)

Performs the += operation.

impl<const ROWS: usize, const COLS: usize, T: Clone> Clone for Matrix<ROWS, COLS, T>

fn clone(&self) -> Matrix<ROWS, COLS, T>

Returns a duplicate of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl<const ROWS: usize, const COLS: usize, T: Copy> Copy for Matrix<ROWS, COLS, T>

impl<const ROWS: usize, const COLS: usize, T: Debug> Debug for Matrix<ROWS, COLS, T>

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.

impl<const ROWS: usize, const COLS: usize, T> From<[[T; COLS]; ROWS]> for Matrix<ROWS, COLS, T>

fn from(data: [[T; COLS]; ROWS]) -> Self

Converts to this type from the input type.

impl<const ROWS: usize, const COLS: usize, T> Index<(usize, usize)> for Matrix<ROWS, COLS, T>

type Output = T

The returned type after indexing.

fn index(&self, (row, col): (usize, usize)) -> &T

Performs the indexing (container[index]) operation.

impl<const ROWS: usize, const COLS: usize, T> IndexMut<(usize, usize)> for Matrix<ROWS, COLS, T>

fn index_mut(&mut self, (row, col): (usize, usize)) -> &mut T

Performs the mutable indexing (container[index]) operation.

impl<const ROWS: usize, const COLS: usize, const C2: usize, T: Numeric> Mul<Matrix<COLS, C2, T>> for Matrix<ROWS, COLS, T>

type Output = Matrix<ROWS, C2, T>

The resulting type after applying the * operator.

fn mul(self, rhs: Matrix<COLS, C2, T>) -> Matrix<ROWS, C2, T>

Performs the * operation.

impl<const ROWS: usize, const COLS: usize, T: Numeric> Mul<T> for Matrix<ROWS, COLS, T>

type Output = Matrix<ROWS, COLS, T>

The resulting type after applying the * operator.

fn mul(self, scalar: T) -> Self

Performs the * operation.

impl<const ROWS: usize, const COLS: usize, T: Numeric> Mul<Vector<COLS, T>> for Matrix<ROWS, COLS, T>

type Output = Vector<ROWS, T>

The resulting type after applying the * operator.

fn mul(self, rhs: Vector<COLS, T>) -> Vector<ROWS, T>

Performs the * operation.

impl<const ROWS: usize, const COLS: usize, T: Numeric> Neg for Matrix<ROWS, COLS, T>

type Output = Matrix<ROWS, COLS, T>

The resulting type after applying the - operator.

fn neg(self) -> Self

Performs the unary - operation.

impl<const ROWS: usize, const COLS: usize, T: PartialEq> PartialEq for Matrix<ROWS, COLS, T>

fn eq(&self, other: &Matrix<ROWS, COLS, T>) -> bool

Tests for self and other values to be equal, and is used by ==.

fn ne(&self, other: &Rhs) -> bool

Tests for !=. The default implementation is almost always sufficient, and should not be overridden without very good reason.

impl<const ROWS: usize, const COLS: usize, T: PartialEq> StructuralPartialEq for Matrix<ROWS, COLS, T>

impl<const ROWS: usize, const COLS: usize, T: Numeric> Sub for Matrix<ROWS, COLS, T>

type Output = Matrix<ROWS, COLS, T>

The resulting type after applying the - operator.

fn sub(self, rhs: Self) -> Self

Performs the - operation.

impl<const ROWS: usize, const COLS: usize, T: Numeric> SubAssign for Matrix<ROWS, COLS, T>

fn sub_assign(&mut self, rhs: Self)

Performs the -= operation.