Struct Matrix 

pub struct Matrix<const N: usize, const M: usize> {
    pub rows: [[f64; M]; N],
}

A matrix with N rows and M columns, allocated on the stack.

Example

use hoomd_linear_algebra::matrix::Matrix;

let a = Matrix {
    rows: [[-1.0, 4.0, -6.0], [2.0, -3.0, 1.0]],
};

Fields

rows: [[f64; M]; N]

The elements of the matrix

Implementations

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

pub fn transpose(&self) -> Matrix<M, N>

Interchange the rows and columns of matrix.

\mathbf{A}^\top_{ji} = \mathbf{A}_{ij}

Examples

use hoomd_linear_algebra::matrix::Matrix;

let m = Matrix {
    rows: [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]],
};
let m_t = m.transpose();
assert_eq!(m_t.rows, [[1.0, 4.0], [2.0, 5.0], [3.0, 6.0]]);

pub fn map_rows<F>(self, f: F) -> Self

where F: FnMut([f64; M]) -> [f64; M],

Apply a function to an Matrix by rows.

Examples

use hoomd_linear_algebra::{Full, GeneralMatrix, matrix::Matrix22};

let m = Matrix22::full(3.0);
assert_eq!(
    m.map_rows(|v| [v[0] + 2.0, v[1]]),
    Matrix22 {
        rows: [[5.0, 3.0], [5.0, 3.0]]
    }
);

pub fn map_columns<F>(self, f: F) -> Self

where F: FnMut([f64; N]) -> [f64; N],

Apply a function to an Matrix by columns.

Examples

use hoomd_linear_algebra::{Full, GeneralMatrix, matrix::Matrix22};

let m = Matrix22::full(3.0);
assert_eq!(
    m.map_columns(|v| [v[0] + 2.0, v[1]]),
    Matrix22 {
        rows: [[5.0, 5.0], [3.0, 3.0]]
    }
);

pub fn map_elements<F>(self, f: F) -> Self

where F: Fn(f64) -> f64,

Apply a function to Matrix elements.

Examples

use hoomd_linear_algebra::{Full, GeneralMatrix, matrix::Matrix33};

let m = Matrix33::full(3.0);
assert_eq!(m.map_elements(|x| x + 2.0), m + Matrix33::full(2.0));

pub fn iter_elements(&self) -> impl Iterator<Item = f64>

Iterate over every element in the Matrix.

The iterator yields all matrix elements in row major order.

Examples

use hoomd_linear_algebra::{SquareMatrix, matrix::Matrix22};

let x = Matrix22 {
    rows: [[1.0, 2.0], [3.0, 4.0]],
};

let mut iterator = x.iter_elements();
assert_eq!(iterator.next(), Some(1.0));
assert_eq!(iterator.next(), Some(2.0));
assert_eq!(iterator.next(), Some(3.0));
assert_eq!(iterator.next(), Some(4.0));
assert_eq!(iterator.next(), None);

pub fn iter_elements_mut(&mut self) -> impl Iterator<Item = &mut f64>

Iterate over every element in the Matrix.

The iterator yields all matrix elements in row major order.

Examples

use hoomd_linear_algebra::{SquareMatrix, matrix::Matrix22};

let mut x = Matrix22 {
    rows: [[1.0, 2.0], [3.0, 4.0]],
};

let x_copy = x.clone();
let mut iterator = x.iter_elements_mut();
iterator.for_each(|x| *x *= 2.0);
assert_eq!(x, x_copy * 2.0);

pub fn iter_rows(&self) -> impl Iterator<Item = [f64; M]>

Iterate over matrix rows.

Examples

use hoomd_linear_algebra::{SquareMatrix, matrix::Matrix22};

let x = Matrix22 {
    rows: [[1.0, 2.0], [3.0, 4.0]],
};
let mut iterator = x.iter_rows();
assert_eq!(iterator.next(), Some([1.0, 2.0]));
assert_eq!(iterator.next(), Some([3.0, 4.0]));
assert_eq!(iterator.next(), None);

pub fn fold_elements<B, F>(self, init: B, f: F) -> B

where F: FnMut(B, f64) -> B,

Folds every element into an accumulator by applying an operation.

fold_elements takes two arguments: an initial value, and a closure with two arguments: an ‘accumulator’, and an element. The closure returns the value that the accumulator should have for the next iteration.

The initial value is the value the accumulator will have on the first call. After applying this closure to every element of the flattened iterator, fold_elements returns the accumulator.

Examples

use hoomd_linear_algebra::{Full, GeneralMatrix, matrix::Matrix22};

let m = Matrix22::full(3.0);
// Sum the elements of a matrix
assert_eq!(m.fold_elements(0.0, |acc, x| acc + x), 3.0 * 4.0);

pub fn fold_rows<B, F>(self, init: B, f: F) -> B

where F: FnMut(B, [f64; M]) -> B,

Folds every row into an accumulator by applying an operation.

fold_rows takes two arguments: an initial value, and a closure with two arguments: an ‘accumulator’, and an element. The closure returns the value that the accumulator should have for the next iteration.

The initial value is the value the accumulator will have on the first call. After applying this closure to every element of the flattened iterator, fold_rows returns the accumulator.

Examples

use hoomd_linear_algebra::{GeneralMatrix, matrix::Matrix22};

let m = Matrix22 {
    rows: [[1.0, 2.0], [3.0, 4.0]],
};
// Average the columns of a matrix
let n_rows = m.n_rows() as f64;
assert_eq!(
    m.fold_rows([0.0; 2], |acc, x| [acc[0] + x[0], acc[1] + x[1]])
        .map(|x| x / n_rows),
    [(1.0 + 3.0) / 2.0, (2.0 + 4.0) / 2.0]
);

pub const fn n_rows(&self) -> usize

Get the number of rows in the Matrix.

Examples

use hoomd_linear_algebra::matrix::Matrix;

let m = Matrix {
    rows: [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]],
};
assert_eq!(m.n_rows(), 2);

pub const fn n_columns(&self) -> usize

Get the number of columns in the Matrix.

Examples

use hoomd_linear_algebra::matrix::Matrix;

let m = Matrix {
    rows: [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]],
};
assert_eq!(m.n_columns(), 3);

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

pub fn with_diagonal(diagonal: [f64; N]) -> Self

Construct a square matrix with the given diagonal.

diagonal[i] is the matrix element $A_{ii}$. All off-diagonal elements are 0.

Example

use hoomd_linear_algebra::matrix::Matrix;

let a = Matrix::with_diagonal([2.0, -3.0]);

assert_eq!(a.rows, [[2.0, 0.0], [0.0, -3.0]]);

pub fn determinant(&self) -> f64

Compute the signed hypervolume of the hyperparallelepiped defined by a matrix.

This implementation uses the Laplace expansion, which is optimal for small matrices but will be extremely slow for large matrices due to its $O(N!)$ complexity.

Example

use hoomd_linear_algebra::{SquareMatrix, matrix::Matrix22};

let identity = Matrix22::identity();
assert_eq!(identity.determinant(), 1.0);

let scaled = identity * 2.0;
assert_eq!(scaled.determinant(), 2.0 * 2.0);

pub fn trace(&self) -> f64

Compute the sum of diagonal elements of a square matrix.

Examples

use hoomd_linear_algebra::{SquareMatrix, matrix::Matrix22};

let identity = Matrix22::identity();
assert_eq!(identity.trace(), 2.0);

let scaled = identity * 3.0;
assert_eq!(scaled.trace(), 3.0 + 3.0);

pub fn powi(&self, n: u32) -> Self

Compute a matrix to an integer power

\mathbf{A}^n = \prod_{i=1}^n \mathbf{A}

Examples

use hoomd_linear_algebra::{Full, GeneralMatrix, MatMul, matrix::Matrix22};

let matrix = Matrix22::full(2.0);

assert_eq!(matrix.powi(2), matrix.matmul(&matrix));

assert_eq!(matrix.powi(2).powi(2), matrix.powi(4));

pub fn diagonal(&self) -> DiagonalMatrix<N>

Extract the diagonal elements from a square matrix.

This method returns a DiagonalMatrix<N> containing the diagonal elements of the input matrix, where the element at position (i, i) is taken from the input matrix. All off-diagonal elements are ignored.

Examples

use hoomd_linear_algebra::matrix::Matrix33;
let a = Matrix33 {
    rows: [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]],
};
let b = a.diagonal();
assert_eq!(b.elements, [1.0, 5.0, 9.0]);

impl Matrix<2, 2>

pub fn svd(&self) -> (Self, DiagonalMatrix<2>, Self)

Decompose a Matrix22 into a rotation, a scaling, and a second rotation.

\mathbf{A} = \mathbf{U} \boldsymbol{\Sigma} \mathbf{V}^\intercal

This implementation is based on the math in Blinn 1996, and ensures good (but not optimal) numerical stability. For certain pathological inputs, preconditioning the matrix could provide a benefit in numerical stability.

svd sets all singular values to be positive.

Examples

use hoomd_linear_algebra::{
    MatMul,
    matrix::{DiagonalMatrix, Matrix22},
};
let m = Matrix22 {
    rows: [[1.0, 2.0], [3.0, 4.0]],
};
let (u, s, vt) = m.svd();
let m_recon = u.matmul(&s.to_dense()).matmul(&vt);
for i in 0..2 {
    for j in 0..2 {
        assert!((m.rows[i][j] - m_recon.rows[i][j]).abs() < 1e-9);
    }
}

impl Matrix<3, 3>

pub fn svd(&self) -> (Self, DiagonalMatrix<3>, Self)

Decompose a Matrix33 into a rotation, a scaling, and a second rotation.

\mathbf{A} = \mathbf{U} \boldsymbol{\Sigma} \mathbf{V}^\top

This implementation is based on the method described by McAdams 2011, which is a fast variant of the Jacobi iteration method.

The method ensures that U and V are pure rotation matrices (determinant = 1). As a result, the third singular value may be negative. For a conventional SVD with non-negative singular values, the sign can be absorbed into U.

Examples

use hoomd_linear_algebra::{
    MatMul, SquareMatrix,
    matrix::{DiagonalMatrix, Matrix33},
};
let m = Matrix33 {
    rows: [[1.0, 2.0, 3.0], [0.0, 1.0, 4.0], [5.0, 6.0, 0.0]],
};
let (u, s, vt) = m.svd();
let m_recon = u.matmul(&s.to_dense()).matmul(&vt);
for i in 0..3 {
    for j in 0..3 {
        assert!((m.rows[i][j] - m_recon.rows[i][j]).abs() < 1e-9);
    }
}

Trait Implementations

impl<const N: usize, const M: usize> Add for Matrix<N, M>

type Output = Matrix<N, M>

The resulting type after applying the + operator.

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

Performs the + operation.

impl<const N: usize, const M: usize> AddAssign for Matrix<N, M>

fn add_assign(&mut self, rhs: Self)

Performs the += operation.

impl<const N: usize, const M: usize> Clone for Matrix<N, M>

fn clone(&self) -> Matrix<N, M>

Returns a duplicate of the value.

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

Performs copy-assignment from source.

impl<const N: usize, const M: usize> Debug for Matrix<N, M>

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

Formats the value using the given formatter.

impl<'de, const N: usize, const M: usize> Deserialize<'de> for Matrix<N, M>

fn deserialize<__D>(__deserializer: __D) -> Result<Self, __D::Error>

where __D: Deserializer<'de>,

Deserialize this value from the given Serde deserializer.

impl<const N: usize, const M: usize> Display for Matrix<N, M>

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

Formats the value using the given formatter.

impl<const N: usize, const M: usize> Full for Matrix<N, M>

fn full(value: f64) -> Self

Construct a matrix with the same value in every element.

Examples

use hoomd_linear_algebra::{Full, matrix::Matrix22};

let m = Matrix22::full(5.0);
assert_eq!(m.rows, [[5.0, 5.0], [5.0, 5.0]]);

impl<const N: usize, const M: usize> GeneralMatrix for Matrix<N, M>

fn zeros() -> Self

Fill a matrix with zeros.

Examples

use hoomd_linear_algebra::{GeneralMatrix, matrix::Matrix22};

let m = Matrix22::zeros();
assert_eq!(m.rows, [[0.0, 0.0], [0.0, 0.0]]);

fn shape(&self) -> (usize, usize)

Get the shape of a matrix (n_rows,n_columns).

impl<const N: usize, const M: usize> Index<(usize, usize)> for Matrix<N, M>

fn index(&self, index: (usize, usize)) -> &f64

Access matrix elements..

Elements are indexed by (row, column).

Examples

use hoomd_linear_algebra::matrix::Matrix;

let rows = [[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]];
let a = Matrix { rows };
assert_eq!(a[(0, 1)], rows[0][1]);
assert_eq!(a[(2, 1)], 6.0);
assert_eq!(a[(1, 1)], 4.0);

type Output = f64

The returned type after indexing.

impl<const N: usize, const M: usize> IndexMut<(usize, usize)> for Matrix<N, M>

fn index_mut(&mut self, index: (usize, usize)) -> &mut f64

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

impl Invertible for Matrix<2, 2>

fn inverse(&self) -> Option<Self>

Compute the inverse of a matrix. Will be None if the matrix is not invertible.

This implementation uses a closed form solution for the matrix inverse.

Examples

use hoomd_linear_algebra::{Invertible, matrix::Matrix22};

let m = Matrix22 {
    rows: [[1.0, 2.0], [3.0, 4.0]],
};
let m_inv = m.inverse().unwrap();
assert_eq!(m_inv.rows, [[-2.0, 1.0], [1.5, -0.5]]);

impl Invertible for Matrix<3, 3>

fn inverse(&self) -> Option<Self>

Compute the inverse of a matrix. Will be None if the matrix is not invertible.

This implementation uses a closed form solution for the matrix inverse based on the cross product of rows.

Example

use hoomd_linear_algebra::{Invertible, SquareMatrix, matrix::Matrix};

let m = Matrix::identity() * 5.0;
let m_inv = m.inverse();

assert_eq!(m_inv, Some(Matrix::with_diagonal([1.0 / 5.0; 3])));

impl Invertible for Matrix<4, 4>

fn inverse(&self) -> Option<Self>

Compute the inverse of a matrix. Will be None if the matrix is not invertible.

This implementation uses a closed form solution for the matrix inverse based on the Cayley–Hamilton method.

Examples

use hoomd_linear_algebra::{
    Invertible, MatMul, SquareMatrix, matrix::Matrix44,
};

let m = Matrix44::identity();
let m_inv = m.inverse().unwrap();
assert_eq!(m_inv.rows, m.rows);

impl<const N: usize, const M: usize> MatMul<DiagonalMatrix<M>> for Matrix<N, M>

fn matmul(&self, rhs: &DiagonalMatrix<M>) -> Self::Output

Matrix-diagonal matrix multiplication.

This is equivalent to scaling each column of a Matrix by the corresponding element in a DiagonalMatrix.

Examples

use hoomd_linear_algebra::{
    Full, GeneralMatrix, MatMul,
    matrix::{DiagonalMatrix, Matrix22},
};

let diag = DiagonalMatrix {
    elements: [3.0, 4.0],
};
let a = Matrix22::full(1.0).matmul(&diag);
assert_eq!(a[(0, 1)], 4.0);
assert_eq!(a[(1, 0)], 3.0);

type Output = Matrix<N, M>

The type of the output matrix.

impl<const N: usize, const M: usize, const K: usize> MatMul<Matrix<M, K>> for Matrix<N, M>

fn matmul(&self, rhs: &Matrix<M, K>) -> Self::Output

Matrix-matrix multiplication.

Examples

use hoomd_linear_algebra::{
    MatMul,
    matrix::{Matrix, Matrix22},
};

let a = Matrix22 {
    rows: [[1.0, 2.0], [3.0, 4.0]],
};
let b = Matrix22 {
    rows: [[5.0, 6.0], [7.0, 8.0]],
};
let c = a.matmul(&b);
assert_eq!(c.rows, [[19.0, 22.0], [43.0, 50.0]]);

type Output = Matrix<N, K>

The type of the output matrix.

impl<const N: usize, const M: usize> Mul<Matrix<N, M>> for f64

fn mul(self, rhs: Self::Output) -> Self::Output

Matrix-scalar multiplication.

Examples

use hoomd_linear_algebra::{Full, GeneralMatrix, matrix::Matrix22};

let matrix = Matrix22::full(2.0);
let scalar = 3.0;
assert_eq!(scalar * matrix, matrix * scalar);

type Output = Matrix<N, M>

The resulting type after applying the * operator.

impl<const N: usize, const M: usize> Mul<f64> for Matrix<N, M>

fn mul(self, rhs: f64) -> Self

Matrix-scalar multiplication.

Examples

use hoomd_linear_algebra::{Full, GeneralMatrix, matrix::Matrix22};

let matrix = Matrix22::full(2.0);
let scalar = 2.0;
assert_eq!(matrix * scalar, matrix + matrix);

type Output = Matrix<N, M>

The resulting type after applying the * operator.

impl<const N: usize, const M: usize> MulAssign<f64> for Matrix<N, M>

fn mul_assign(&mut self, rhs: f64)

Matrix-scalar multiplication assignment.

Examples

use hoomd_linear_algebra::{Full, GeneralMatrix, matrix::Matrix22};

let mut matrix = Matrix22::full(2.0);
let matrix_copy = matrix.clone();
matrix *= 3.0;
assert_eq!(matrix, matrix_copy * 3.0);

impl<const N: usize, const M: usize> Neg for Matrix<N, M>

fn neg(self) -> Self

Matrix negation.

Examples

use hoomd_linear_algebra::{Full, GeneralMatrix, matrix::Matrix22};

let matrix = Matrix22::full(5.0);
assert_eq!(-matrix, Matrix22::zeros() - matrix);

type Output = Matrix<N, M>

The resulting type after applying the - operator.

impl<const N: usize, const M: usize> PartialEq for Matrix<N, M>

fn eq(&self, other: &Matrix<N, M>) -> 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 N: usize> QuadraticForm<N> for Matrix<N, N>

fn compute_quadratic_form(&self, x: &[f64; N]) -> f64

Evaluate the quadratic form.

impl<const N: usize, const M: usize> Serialize for Matrix<N, M>

fn serialize<__S>(&self, __serializer: __S) -> Result<__S::Ok, __S::Error>

where __S: Serializer,

Serialize this value into the given Serde serializer.

impl<const N: usize> SquareMatrix for Matrix<N, N>

fn identity() -> Self

Examples

use hoomd_linear_algebra::{SquareMatrix, matrix::Matrix22};

let m = Matrix22::identity();
assert_eq!(m.rows, [[1.0, 0.0], [0.0, 1.0]]);

impl<const N: usize, const M: usize> Sub for Matrix<N, M>

type Output = Matrix<N, M>

The resulting type after applying the - operator.

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

Performs the - operation.

impl<const N: usize, const M: usize> SubAssign for Matrix<N, M>

fn sub_assign(&mut self, rhs: Self)

Performs the -= operation.

impl Copy for Matrix<1, 1>

impl Copy for Matrix<1, 2>

impl Copy for Matrix<1, 3>

impl Copy for Matrix<1, 4>

impl Copy for Matrix<2, 1>

impl Copy for Matrix<2, 2>

impl Copy for Matrix<2, 3>

impl Copy for Matrix<2, 4>

impl Copy for Matrix<3, 1>

impl Copy for Matrix<3, 2>

impl Copy for Matrix<3, 3>

impl Copy for Matrix<3, 4>

impl Copy for Matrix<4, 1>

impl Copy for Matrix<4, 2>

impl Copy for Matrix<4, 3>

impl Copy for Matrix<4, 4>

impl<const N: usize, const M: usize> StructuralPartialEq for Matrix<N, M>