drewcassidy/vector-victor
1
0
Fork 0
mirror of https://github.com/drewcassidy/vector-victor.git synced 2024-09-01 14:58:35 +00:00
Code Issues Packages Projects Releases Wiki Activity

Compare commits

base: drewcassidy:bc1b3f199dd7e120e0671142a3533c1d904935aa
Branches Tags
drewcassidy:develop
drewcassidy:main
drewcassidy:0.1.2
drewcassidy:0.1.1
drewcassidy:0.1.0
..
compare: drewcassidy:e2a2bc7529a5b3d0e22075bdd2f9238af5e83d00
Branches Tags
drewcassidy:develop
drewcassidy:main
drewcassidy:0.1.2
drewcassidy:0.1.1
drewcassidy:0.1.0

2 Commits
bc1b3f199d ... e2a2bc7529

Author SHA1 Message Date
Andrew Cassidy
e2a2bc7529 Huge refactor to flatten everything 
but no API change I think
 2023-05-21 22:52:47 -07:00
Andrew Cassidy
2b303892f7 Refactor LUDecompose slightly 2023-05-21 20:55:37 -07:00
9 changed files with 737 additions and 604 deletions
Show Stats Expand all files Collapse all files

96
src/decompose.rs
View File

@ -1,6 +1,7 @@
use crate::util::checked_inv;
use crate::{Matrix, Vector};
use num_traits::real::Real;
use num_traits::Signed;
use std::iter::{Product, Sum};
use std::ops::{Mul, Neg, Not};
@ -37,7 +38,7 @@ impl Not for Parity {
}
}
/// The result of the [LU decomposition](LUDecomposable::lu) of a matrix.
/// The result of the [LU decomposition](LUDecompose::lu) of a matrix.
///
/// This struct provides a convenient way to reuse one LU decomposition to solve multiple
/// matrix equations. You likely do not need to worry about its contents.
@ -46,29 +47,29 @@ impl Not for Parity {
/// on wikipedia for more information
#[derive(Copy, Clone, Debug, PartialEq)]
pub struct LUDecomposition<T: Copy, const N: usize> {
/// The $L$ and $U$ matrices combined into one
/// The $bbL$ and $bbU$ matrices combined into one
///
/// for example if
///
/// $ U = [[u_{11}, u_{12}, cdots, u_{1n} ],
/// [0, u_{22}, cdots, u_{2n} ],
/// [vdots, vdots, ddots, vdots ],
/// [0, 0, cdots, u_{mn} ]] $
/// $ bbU = [[u_{11}, u_{12}, cdots, u_{1n} ],
/// [0, u_{22}, cdots, u_{2n} ],
/// [vdots, vdots, ddots, vdots ],
/// [0, 0, cdots, u_{mn} ]] $
/// and
/// $ L = [[1, 0, cdots, 0 ],
/// [l_{21}, 1, cdots, 0 ],
/// [vdots, vdots, ddots, vdots ],
/// [l_{m1}, l_{m2}, cdots, 1 ]] $,
/// $ bbL = [[1, 0, cdots, 0 ],
/// [l_{21}, 1, cdots, 0 ],
/// [vdots, vdots, ddots, vdots ],
/// [l_{m1}, l_{m2}, cdots, 1 ]] $,
/// then
/// $ LU = [[u_{11}, u_{12}, cdots, u_{1n} ],
/// [l_{21}, u_{22}, cdots, u_{2n} ],
/// [vdots, vdots, ddots, vdots ],
/// [l_{m1}, l_{m2}, cdots, u_{mn} ]] $
/// $ bb{LU} = [[u_{11}, u_{12}, cdots, u_{1n} ],
/// [l_{21}, u_{22}, cdots, u_{2n} ],
/// [vdots, vdots, ddots, vdots ],
/// [l_{m1}, l_{m2}, cdots, u_{mn} ]] $
///
/// note that the diagonals of the $L$ matrix are always 1, so no information is lost
/// note that the diagonals of the $bbL$ matrix are always 1, so no information is lost
pub lu: Matrix<T, N, N>,
/// The indices of the permutation matrix $P$, such that $PxxA$ = $LxxU$
/// The indices of the permutation matrix $bbP$, such that $bbP xx bbA$ = $bbL xx bbU$
///
/// The permutation matrix rearranges the rows of the original matrix in order to produce
/// the LU decomposition. This makes calculation simpler, but makes the result
@ -79,13 +80,10 @@ pub struct LUDecomposition<T: Copy, const N: usize> {
pub parity: Parity,
}
impl<T: Copy + Default, const N: usize> LUDecomposition<T, N>
where
T: Real + Default + Sum + Product,
{
/// Solve for $x$ in $M xx x = b$, where $M$ is the original matrix this is a decomposition of.
impl<T: Copy + Default + Real, const N: usize> LUDecomposition<T, N> {
/// Solve for $x$ in $bbM xx x = b$, where $bbM$ is the original matrix this is a decomposition of.
///
/// This is equivalent to [`LUDecomposable::solve`] while allowing the LU decomposition
/// This is equivalent to [`LUDecompose::solve`] while allowing the LU decomposition
/// to be reused
#[must_use]
pub fn solve<const M: usize>(&self, b: &Matrix<T, N, M>) -> Matrix<T, N, M> {
@ -123,17 +121,17 @@ where
/// Calculate the determinant $|M|$ of the matrix $M$.
/// If the matrix is singular, the determinant is 0.
///
/// This is equivalent to [`LUDecomposable::det`] while allowing the LU decomposition
/// This is equivalent to [`LUDecompose::det`] while allowing the LU decomposition
/// to be reused
pub fn det(&self) -> T {
self.parity * self.lu.diagonals().product()
self.parity * self.lu.diagonals().fold(T::one(), T::mul)
}
/// Calculate the inverse of the original matrix, such that $MxxM^{-1} = I$
/// Calculate the inverse of the original matrix, such that $bbM xx bbM^{-1} = bbI$
///
/// This is equivalent to [`Matrix::inverse`] while allowing the LU decomposition to be reused
/// This is equivalent to [`Matrix::inv`] while allowing the LU decomposition to be reused
#[must_use]
pub fn inverse(&self) -> Matrix<T, N, N> {
pub fn inv(&self) -> Matrix<T, N, N> {
return self.solve(&Matrix::<T, N, N>::identity());
}
@ -160,17 +158,14 @@ where
///
/// See [LU decomposition](https://en.wikipedia.org/wiki/LU_decomposition)
/// on wikipedia for more information
pub trait LUDecomposable<T, const N: usize>
where
T: Copy + Default + Real + Product + Sum,
{
pub trait LUDecompose<T: Copy, const N: usize> {
/// return this matrix's [`LUDecomposition`], or [`None`] if the matrix is singular.
/// This can be used to solve for multiple results
///
/// ```
/// # use vector_victor::decompose::LUDecomposable;
/// # use vector_victor::decompose::LUDecompose;
/// # use vector_victor::{Matrix, Vector};
/// let m = Matrix::new([[1.0,3.0],[2.0,4.0]]);
/// let m = Matrix::mat([[1.0,3.0],[2.0,4.0]]);
/// let lu = m.lu().expect("Cannot decompose a signular matrix");
///
/// let b = Vector::vec([7.0,10.0]);
@ -183,34 +178,35 @@ where
#[must_use]
fn lu(&self) -> Option<LUDecomposition<T, N>>;
/// Calculate the inverse of the matrix, such that $MxxM^{-1} = I$, or [`None`] if the matrix is singular.
/// Calculate the inverse of the matrix, such that $bbMxxbbM^{-1} = bbI$,
/// or [`None`] if the matrix is singular.
///
/// ```
/// # use vector_victor::decompose::LUDecomposable;
/// # use vector_victor::decompose::LUDecompose;
/// # use vector_victor::Matrix;
/// let m = Matrix::new([[1.0,3.0],[2.0,4.0]]);
/// let mi = m.inverse().expect("Cannot invert a singular matrix");
/// let m = Matrix::mat([[1.0,3.0],[2.0,4.0]]);
/// let mi = m.inv().expect("Cannot invert a singular matrix");
///
/// assert_eq!(mi, Matrix::new([[-2.0, 1.5],[1.0, -0.5]]), "unexpected inverse matrix");
/// assert_eq!(mi, Matrix::mat([[-2.0, 1.5],[1.0, -0.5]]), "unexpected inverse matrix");
///
/// // multiplying a matrix by its inverse yields the identity matrix
/// assert_eq!(m.mmul(&mi), Matrix::identity())
/// ```
#[must_use]
fn inverse(&self) -> Option<Matrix<T, N, N>>;
fn inv(&self) -> Option<Matrix<T, N, N>>;
/// Calculate the determinant $|M|$ of the matrix $M$.
/// If the matrix is singular, the determinant is 0
#[must_use]
fn det(&self) -> T;
/// Solve for $x$ in $M xx x = b$
/// Solve for $x$ in $bbM xx x = b$
///
/// ```
/// # use vector_victor::decompose::LUDecomposable;
/// # use vector_victor::decompose::LUDecompose;
/// # use vector_victor::{Matrix, Vector};
///
/// let m = Matrix::new([[1.0,3.0],[2.0,4.0]]);
/// let m = Matrix::mat([[1.0,3.0],[2.0,4.0]]);
/// let b = Vector::vec([7.0,10.0]);
/// let x = m.solve(&b).expect("Cannot solve a singular matrix");
///
@ -219,26 +215,26 @@ where
/// ```
///
/// $x$ does not need to be a column-vector, it can also be a 2D matrix. For example,
/// the following is another way to calculate the [inverse](LUDecomposable::inverse()) by solving for the identity matrix $I$.
/// the following is another way to calculate the [inverse](LUDecompose::inv()) by solving for the identity matrix $I$.
///
/// ```
/// # use vector_victor::decompose::LUDecomposable;
/// # use vector_victor::decompose::LUDecompose;
/// # use vector_victor::{Matrix, Vector};
///
/// let m = Matrix::new([[1.0,3.0],[2.0,4.0]]);
/// let m = Matrix::mat([[1.0,3.0],[2.0,4.0]]);
/// let i = Matrix::<f64,2,2>::identity();
/// let mi = m.solve(&i).expect("Cannot solve a singular matrix");
///
/// assert_eq!(mi, Matrix::new([[-2.0, 1.5],[1.0, -0.5]]));
/// assert_eq!(mi, Matrix::mat([[-2.0, 1.5],[1.0, -0.5]]));
/// assert_eq!(m.mmul(&mi), i, "M x M^-1 = I");
/// ```
#[must_use]
fn solve<const M: usize>(&self, b: &Matrix<T, N, M>) -> Option<Matrix<T, N, M>>;
}
impl<T, const N: usize> LUDecomposable<T, N> for Matrix<T, N, N>
impl<T, const N: usize> LUDecompose<T, N> for Matrix<T, N, N>
where
T: Copy + Default + Real + Sum + Product,
T: Copy + Default + Real + Sum + Product + Signed,
{
fn lu(&self) -> Option<LUDecomposition<T, N>> {
// Implementation from Numerical Recipes §2.3
@ -300,7 +296,7 @@ where
return Some(LUDecomposition { lu, idx, parity });
}
fn inverse(&self) -> Option<Matrix<T, N, N>> {
fn inv(&self) -> Option<Matrix<T, N, N>> {
match N {
1 => Some(Self::fill(checked_inv(self[0])?)),
2 => {
@ -311,7 +307,7 @@ where
result[(0, 1)] = -self[(0, 1)];
Some(result * checked_inv(self.det())?)
}
_ => Some(self.lu()?.inverse()),
_ => Some(self.lu()?.inv()),
}
}

66
src/identities.rs Normal file
View File

@ -0,0 +1,66 @@
use crate::Matrix;
use num_traits::{Bounded, One, Zero};
// Identity
impl<T: Copy + Zero + One, const N: usize> Matrix<T, N, N> {
/// Create an identity matrix, a square matrix where the diagonals are 1 and all other elements
/// are 0.
/// for example,
///
/// $bbI = [[1,0,0],[0,1,0],[0,0,1]]$
///
/// Matrix multiplication between a matrix and the identity matrix always results in itself
///
/// $bbA xx bbI = bbA$
///
/// # Examples
/// ```
/// # use vector_victor::Matrix;
/// let i = Matrix::<i32,3,3>::identity();
/// assert_eq!(i, Matrix::mat([[1,0,0],[0,1,0],[0,0,1]]))
/// ```
///
/// Note that the identity only exists for matrices that are square, so this doesnt work:
/// ```compile_fail
/// # use vector_victor::Matrix;
/// let i = Matrix::<i32,4,2>::identity();
/// ```
#[must_use]
pub fn identity() -> Self {
let mut result = Self::zero();
for i in 0..N {
result[(i, i)] = T::one();
}
return result;
}
}
// Zero
impl<T: Copy + Zero, const M: usize, const N: usize> Zero for Matrix<T, M, N> {
fn zero() -> Self {
Matrix::fill(T::zero())
}
fn is_zero(&self) -> bool {
self.elements().all(|e| e.is_zero())
}
}
// One
impl<T: Copy + One, const M: usize, const N: usize> One for Matrix<T, M, N> {
fn one() -> Self {
Matrix::fill(T::one())
}
}
// min_value and max_value
// LowerBounded and UpperBounded are automatically implemented from this
impl<T: Copy + Bounded, const N: usize, const M: usize> Bounded for Matrix<T, N, M> {
fn min_value() -> Self {
Self::fill(T::min_value())
}
fn max_value() -> Self {
Self::fill(T::max_value())
}
}

0
src/matrix/index.rs → src/index.rs
View File

419
src/lib.rs
View File

@ -1,7 +1,422 @@
extern crate core;
use index::Index2D;
use std::cmp::min;
use std::fmt::Debug;
use std::iter::{zip, Flatten};
use std::ops::{Index, IndexMut};
pub mod decompose;
mod matrix;
mod identities;
pub mod index;
mod math;
mod ops;
mod util;
pub use matrix::{Matrix, Vector};
/// A 2D array of values which can be operated upon.
///
/// Matrices have a fixed size known at compile time
#[derive(Debug, Copy, Clone, PartialEq)]
pub struct Matrix<T, const M: usize, const N: usize>
where
T: Copy,
{
data: [[T; N]; M], // Row-Major order
}
/// An alias for a [Matrix] with a single column
pub type Vector<T, const N: usize> = Matrix<T, N, 1>;
// CONSTRUCTORS
// Default
impl<T: Copy + Default, const M: usize, const N: usize> Default for Matrix<T, M, N> {
fn default() -> Self {
Matrix::fill(T::default())
}
}
// Matrix constructors
impl<T: Copy, const M: usize, const N: usize> Matrix<T, M, N> {
/// Generate a new matrix from a 2D Array
///
/// # Arguments
///
/// * `data`: A 2D array of elements to copy into the new matrix
///
/// returns: Matrix<T, M, N>
///
/// # Examples
///
/// ```
/// # use vector_victor::Matrix;
/// let a = Matrix::mat([[1,2,3,4];4]);
/// ```
#[must_use]
pub fn mat(data: [[T; N]; M]) -> Self {
assert!(M > 0, "Matrix must have at least 1 row");
assert!(N > 0, "Matrix must have at least 1 column");
Matrix::<T, M, N> { data }
}
/// Generate a new matrix from a single scalar
///
/// # Arguments
///
/// * `scalar`: Scalar value to copy into the new matrix.
///
/// returns: Matrix<T, M, N>
///
/// # Examples
///
/// ```
/// # use vector_victor::Matrix;
/// let my_matrix = Matrix::<i32,4,4>::fill(5);
/// // is equivalent to
/// assert_eq!(my_matrix, Matrix::mat([[5;4];4]))
/// ```
#[must_use]
pub fn fill(scalar: T) -> Matrix<T, M, N> {
assert!(M > 0, "Matrix must have at least 1 row");
assert!(N > 0, "Matrix must have at least 1 column");
Matrix::<T, M, N> {
data: [[scalar; N]; M],
}
}
/// Create a matrix from an iterator of vectors
///
/// # Arguments
///
/// * `iter`: iterator of vectors to copy into rows
///
/// returns: Matrix<T, M, N>
///
/// # Examples
///
/// ```
/// # use vector_victor::Matrix;
/// let my_matrix = Matrix::mat([[1,2,3],[4,5,6]]);
/// let transpose : Matrix<_,3,2>= Matrix::from_rows(my_matrix.cols());
/// assert_eq!(transpose, Matrix::mat([[1,4],[2,5],[3,6]]))
/// ```
#[must_use]
pub fn from_rows<I>(iter: I) -> Self
where
I: IntoIterator<Item = Vector<T, N>>,
Self: Default,
{
let mut result = Self::default();
for (m, row) in iter.into_iter().enumerate().take(M) {
result.set_row(m, &row)
}
result
}
/// Create a matrix from an iterator of vectors
///
/// # Arguments
///
/// * `iter`: iterator of vectors to copy into columns
///
/// returns: Matrix<T, M, N>
///
/// # Examples
///
/// ```
/// # use vector_victor::Matrix;
/// let my_matrix = Matrix::mat([[1,2,3],[4,5,6]]);
/// let transpose : Matrix<_,3,2>= Matrix::from_cols(my_matrix.rows());
/// assert_eq!(transpose, Matrix::mat([[1,4],[2,5],[3,6]]))
/// ```
#[must_use]
pub fn from_cols<I>(iter: I) -> Self
where
I: IntoIterator<Item = Vector<T, M>>,
Self: Default,
{
let mut result = Self::default();
for (n, col) in iter.into_iter().enumerate().take(N) {
result.set_col(n, &col)
}
result
}
}
// Vector constructor
impl<T: Copy, const N: usize> Vector<T, N> {
/// Create a vector from a 1D array.
/// Note that vectors are always column vectors unless explicitly instantiated as row vectors
///
/// # Examples
/// ```
/// # use vector_victor::{Matrix, Vector};
/// let my_vector = Vector::vec([1,2,3,4]);
/// // is equivalent to
/// assert_eq!(my_vector, Matrix::mat([[1],[2],[3],[4]]));
/// ```
pub fn vec(data: [T; N]) -> Self {
assert!(N > 0, "Vector must have at least 1 element");
return Vector::<T, N> {
data: data.map(|e| [e]),
};
}
}
// ACCESSORS AND MUTATORS
impl<T: Copy, const M: usize, const N: usize> Matrix<T, M, N> {
/// Returns an iterator over the elements of the matrix in row-major order.
///
/// # Examples
/// ```
/// # use vector_victor::Matrix;
/// let my_matrix = Matrix::mat([[1,2],[3,4]]);
/// assert!(vec![1,2,3,4].iter().eq(my_matrix.elements()))
/// ```
#[must_use]
pub fn elements<'a>(&'a self) -> impl Iterator<Item = &'a T> + 'a {
self.data.iter().flatten()
}
/// Returns a mutable iterator over the elements of the matrix in row-major order.
#[must_use]
pub fn elements_mut<'a>(&'a mut self) -> impl Iterator<Item = &'a mut T> + 'a {
self.data.iter_mut().flatten()
}
/// returns an iterator over the elements along the diagonal of a matrix
#[must_use]
pub fn diagonals<'s>(&'s self) -> impl Iterator<Item = T> + 's {
(0..min(N, M)).map(|n| self[(n, n)])
}
/// Returns an iterator over the elements directly below the diagonal of a matrix
#[must_use]
pub fn subdiagonals<'s>(&'s self) -> impl Iterator<Item = T> + 's {
(0..min(N, M) - 1).map(|n| self[(n, n + 1)])
}
/// Returns a reference to the element at that position in the matrix, or `None` if out of bounds.
///
/// # Examples
///
/// ```
/// # use vector_victor::Matrix;
/// let my_matrix = Matrix::mat([[1,2],[3,4]]);
///
/// // element at index 2 is the same as the element at (row 1, column 0).
/// assert_eq!(my_matrix.get(2), my_matrix.get((1,0)));
///
/// // my_matrix.get() is equivalent to my_matrix[],
/// // but returns an Option instead of panicking
/// assert_eq!(my_matrix.get(2), Some(&my_matrix[2]));
///
/// // index 4 is out of range, so get(4) returns None.
/// assert_eq!(my_matrix.get(4), None);
/// ```
#[inline]
#[must_use]
pub fn get(&self, index: impl Index2D) -> Option<&T> {
let (m, n) = index.to_2d(M, N)?;
Some(&self.data[m][n])
}
/// Returns a mutable reference to the element at that position in the matrix, or `None` if out of bounds.
#[inline]
#[must_use]
pub fn get_mut(&mut self, index: impl Index2D) -> Option<&mut T> {
let (m, n) = index.to_2d(M, N)?;
Some(&mut self.data[m][n])
}
/// Returns a row of the matrix. or [None] if index is out of bounds
///
/// # Examples
///
/// ```
/// # use vector_victor::{Matrix, Vector};
/// let my_matrix = Matrix::mat([[1,2],[3,4]]);
///
/// // row at index 1
/// assert_eq!(my_matrix.row(1), Vector::vec([3,4]));
/// ```
#[inline]
#[must_use]
pub fn row(&self, m: usize) -> Vector<T, N> {
assert!(
m < M,
"Row index {} out of bounds for {}x{} matrix",
m,
M,
N
);
Vector::<T, N>::vec(self.data[m])
}
#[inline]
pub fn set_row(&mut self, m: usize, val: &Vector<T, N>) {
assert!(
m < M,
"Row index {} out of bounds for {}x{} matrix",
m,
M,
N
);
for n in 0..N {
self.data[m][n] = val.data[n][0];
}
}
pub fn pivot_row(&mut self, m1: usize, m2: usize) {
let tmp = self.row(m2);
self.set_row(m2, &self.row(m1));
self.set_row(m1, &tmp);
}
#[inline]
#[must_use]
pub fn col(&self, n: usize) -> Vector<T, M> {
assert!(
n < N,
"Column index {} out of bounds for {}x{} matrix",
n,
M,
N
);
Vector::<T, M>::vec(self.data.map(|r| r[n]))
}
#[inline]
pub fn set_col(&mut self, n: usize, val: &Vector<T, M>) {
assert!(
n < N,
"Column index {} out of bounds for {}x{} matrix",
n,
M,
N
);
for m in 0..M {
self.data[m][n] = val.data[m][0];
}
}
pub fn pivot_col(&mut self, n1: usize, n2: usize) {
let tmp = self.col(n2);
self.set_col(n2, &self.col(n1));
self.set_col(n1, &tmp);
}
#[must_use]
pub fn rows<'a>(&'a self) -> impl Iterator<Item = Vector<T, N>> + 'a {
(0..M).map(|m| self.row(m))
}
#[must_use]
pub fn cols<'a>(&'a self) -> impl Iterator<Item = Vector<T, M>> + 'a {
(0..N).map(|n| self.col(n))
}
#[must_use]
pub fn permute_rows(&self, ms: &Vector<usize, M>) -> Self
where
T: Default,
{
Self::from_rows(ms.elements().map(|&m| self.row(m)))
}
#[must_use]
pub fn permute_cols(&self, ns: &Vector<usize, N>) -> Self
where
T: Default,
{
Self::from_cols(ns.elements().map(|&n| self.col(n)))
}
pub fn transpose(&self) -> Matrix<T, N, M>
where
Matrix<T, N, M>: Default,
{
Matrix::<T, N, M>::from_rows(self.cols())
}
}
// Index
impl<I, T, const M: usize, const N: usize> Index<I> for Matrix<T, M, N>
where
I: Index2D,
T: Copy,
{
type Output = T;
#[inline(always)]
fn index(&self, index: I) -> &Self::Output {
self.get(index).expect(&*format!(
"index {:?} out of range for {}x{} Matrix",
index, M, N
))
}
}
// IndexMut
impl<I, T, const M: usize, const N: usize> IndexMut<I> for Matrix<T, M, N>
where
I: Index2D,
T: Copy,
{
#[inline(always)]
fn index_mut(&mut self, index: I) -> &mut Self::Output {
self.get_mut(index).expect(&*format!(
"index {:?} out of range for {}x{} Matrix",
index, M, N
))
}
}
// CONVERSIONS
// Convert from 2D Array (equivalent to new)
impl<T: Copy, const M: usize, const N: usize> From<[[T; N]; M]> for Matrix<T, M, N> {
fn from(data: [[T; N]; M]) -> Self {
Self::mat(data)
}
}
// Convert from 1D Array (equivalent to vec)
impl<T: Copy, const M: usize> From<[T; M]> for Vector<T, M> {
fn from(data: [T; M]) -> Self {
Self::vec(data)
}
}
// Convert from scalar (equivalent to fill)
impl<T: Copy, const M: usize, const N: usize> From<T> for Matrix<T, M, N> {
fn from(scalar: T) -> Self {
Self::fill(scalar)
}
}
// IntoIter
impl<T: Copy, const M: usize, const N: usize> IntoIterator for Matrix<T, M, N> {
type Item = T;
type IntoIter = Flatten<std::array::IntoIter<[T; N], M>>;
fn into_iter(self) -> Self::IntoIter {
self.data.into_iter().flatten()
}
}
// FromIterator
impl<T: Copy, const M: usize, const N: usize> FromIterator<T> for Matrix<T, M, N>
where
Self: Default,
{
fn from_iter<I: IntoIterator<Item = T>>(iter: I) -> Self {
let mut result: Self = Default::default();
for (l, r) in zip(result.elements_mut(), iter) {
*l = r;
}
result
}
}

182
src/math.rs Normal file
View File

@ -0,0 +1,182 @@
use crate::{Matrix, Vector};
use num_traits::real::Real;
use num_traits::{Inv, Num, NumOps, One, Pow, Signed, Zero};
use std::iter::{zip, Product, Sum};
use std::ops::{Add, Mul};
/// Operations for column vectors
impl<T: Copy, const N: usize> Vector<T, N> {
/// Compute the dot product of two vectors, otherwise known as the scalar product.
/// This is the sum of the elementwise product, or in math terms
///
/// $vec(a) * vec(b) = sum_(i=1)^n a_i b_i = a_1 b_1 + a_2 b_2 + ... + a_n b_n$
///
/// for example, $[[1],[2],[3]] * [[4],[5],[6]] = (1 * 4) + (2 * 5) + (3 * 6) = 32$
///
/// For vectors in euclidean space, this has the property that it is equal to the magnitudes of
/// the vectors times the cosine of the angle between them.
///
/// $vec(a) * vec(b) = |vec(a)| |vec(b)| cos(theta)$
///
/// this also gives it the special property that the dot product of a vector and itself is the
/// square of its magnitude. You may recognize the 2D version as the
/// [pythagorean theorem](https://en.wikipedia.org/wiki/Pythagorean_theorem).
///
/// see [dot product](https://en.wikipedia.org/wiki/Dot_product) on Wikipedia for more
/// information.
pub fn dot<R>(&self, rhs: &R) -> T
where
for<'s> &'s Self: Mul<&'s R, Output = Self>,
T: Sum<T>,
{
(self * rhs).elements().cloned().sum()
}
pub fn sqrmag(&self) -> T
where
for<'s> &'s Self: Mul<&'s Self, Output = Self>,
T: Sum<T>,
{
self.dot(self)
}
pub fn mag(&self) -> T
where
T: Sum<T> + Mul<T> + Real,
{
self.sqrmag().sqrt()
}
pub fn normalized(&self) -> Option<Self>
where
T: Sum<T> + Mul<T> + Real,
{
match self.mag() {
mag if mag.abs() < T::epsilon() => None,
mag => Some(self / mag),
}
}
}
/// Cross product operations for column vectors in $RR^3$
impl<T: Copy> Vector<T, 3> {
pub fn cross_r<R: Copy>(&self, rhs: &Vector<R, 3>) -> Self
where
T: NumOps<R> + NumOps,
{
Self::vec([
(self[1] * rhs[2]) - (self[2] * rhs[1]),
(self[2] * rhs[0]) - (self[0] * rhs[2]),
(self[0] * rhs[1]) - (self[1] * rhs[0]),
])
}
pub fn cross_l<R: Copy>(&self, rhs: &Vector<R, 3>) -> Vector<R, 3>
where
R: NumOps<T> + NumOps,
{
rhs.cross_r(self)
}
}
/// Operations for Matrices
impl<T: Copy, const M: usize, const N: usize> Matrix<T, M, N> {
pub fn mmul<R: Copy, const P: usize>(&self, rhs: &Matrix<R, N, P>) -> Matrix<T, M, P>
where
T: Default + NumOps<R> + Sum,
{
let mut result: Matrix<T, M, P> = Default::default();
for (m, a) in self.rows().enumerate() {
for (n, b) in rhs.cols().enumerate() {
result[(m, n)] = a.dot(&b)
}
}
return result;
}
/// Computes the absolute value of each element of the matrix
pub fn abs(&self) -> Self
where
T: Signed + Default,
{
self.elements().map(|&x| x.abs()).collect()
}
/// Computes the sign of each element of the matrix
pub fn signum(&self) -> Self
where
T: Signed + Default,
{
self.elements().map(|&x| x.signum()).collect()
}
/// Raises every element to the power of rhs, where rhs is either a scalar or a matrix of exponents
pub fn pow<R, O>(self, rhs: R) -> O
where
Self: Pow<R, Output = O>,
{
Pow::pow(self, rhs)
}
}
// Sum up matrices
impl<T: Copy, const M: usize, const N: usize> Sum for Matrix<T, M, N>
where
Self: Zero + Add<Output = Self>,
{
fn sum<I: Iterator<Item = Self>>(iter: I) -> Self {
iter.fold(Self::zero(), Self::add)
}
}
// Product of matrices
impl<T: Copy, const M: usize, const N: usize> Product for Matrix<T, M, N>
where
Self: One + Mul<Output = Self>,
{
fn product<I: Iterator<Item = Self>>(iter: I) -> Self {
iter.fold(Self::one(), Self::mul)
}
}
/// Inverse trait. Note that this is the elementwise inverse, not the matrix inverse.
/// For the inverse matrix see [`LUDecomposable::inv()`](crate::decompose::LUDecompose::inv())
impl<T: Copy + Inv<Output = T> + Default, const M: usize, const N: usize> Inv for Matrix<T, M, N> {
type Output = Self;
fn inv(self) -> Self::Output {
self.elements().map(|t| t.inv()).collect()
}
}
/// Pow for $Matrix^{scalar}$
impl<T, R, O, const M: usize, const N: usize> Pow<R> for Matrix<T, M, N>
where
T: Copy + Pow<R, Output = O>,
R: Copy + Num,
O: Copy + Default,
{
type Output = Matrix<O, M, N>;
fn pow(self, rhs: R) -> Self::Output {
self.elements().map(|&x| x.pow(rhs)).collect()
}
}
/// Pow for $Matrix^{Matrix}$
impl<T, R, O, const M: usize, const N: usize> Pow<Matrix<R, M, N>> for Matrix<T, M, N>
where
T: Copy + Pow<R, Output = O>,
R: Copy,
O: Copy + Default,
{
type Output = Matrix<O, M, N>;
fn pow(self, rhs: Matrix<R, M, N>) -> Self::Output {
zip(self.elements(), rhs.elements())
.map(|(x, &r)| x.pow(r))
.collect()
}
}

533
src/matrix/mod.rs
View File

@ -1,533 +0,0 @@
use index::Index2D;
use num_traits::{NumOps, One, Zero};
use std::fmt::Debug;
use std::iter::{zip, Flatten, Product, Sum};
use std::ops::{Add, AddAssign, Deref, DerefMut, Index, IndexMut, Mul, MulAssign, Neg};
pub mod ops;
mod index;
/// A 2D array of values which can be operated upon.
///
/// Matrices have a fixed size known at compile time
#[derive(Debug, Copy, Clone, PartialEq)]
pub struct Matrix<T, const M: usize, const N: usize>
where
T: Copy,
{
data: [[T; N]; M], // Row-Major order
}
/// An alias for a [Matrix] with a single column
pub type Vector<T, const N: usize> = Matrix<T, N, 1>;
// Simple access functions that only require T be copyable
impl<T: Copy, const M: usize, const N: usize> Matrix<T, M, N> {
/// Generate a new matrix from a 2D Array
///
/// # Arguments
///
/// * `data`: A 2D array of elements to copy into the new matrix
///
/// returns: Matrix<T, M, N>
///
/// # Examples
///
/// ```
/// # use vector_victor::Matrix;
/// let a = Matrix::new([[1,2,3,4];4]);
/// ```
#[must_use]
pub fn new(data: [[T; N]; M]) -> Self {
assert!(M > 0, "Matrix must have at least 1 row");
assert!(N > 0, "Matrix must have at least 1 column");
Matrix::<T, M, N> { data }
}
/// Generate a new matrix from a single scalar
///
/// # Arguments
///
/// * `scalar`: Scalar value to copy into the new matrix.
///
/// returns: Matrix<T, M, N>
///
/// # Examples
///
/// ```
/// # use vector_victor::Matrix;
/// let my_matrix = Matrix::<i32,4,4>::fill(5);
/// // is equivalent to
/// assert_eq!(my_matrix, Matrix::new([[5;4];4]))
/// ```
#[must_use]
pub fn fill(scalar: T) -> Matrix<T, M, N> {
assert!(M > 0, "Matrix must have at least 1 row");
assert!(N > 0, "Matrix must have at least 1 column");
Matrix::<T, M, N> {
data: [[scalar; N]; M],
}
}
/// Create a matrix from an iterator of vectors
///
/// # Arguments
///
/// * `iter`: iterator of vectors to copy into rows
///
/// returns: Matrix<T, M, N>
///
/// # Examples
///
/// ```
/// # use vector_victor::Matrix;
/// let my_matrix = Matrix::new([[1,2,3],[4,5,6]]);
/// let transpose : Matrix<_,3,2>= Matrix::from_rows(my_matrix.cols());
/// assert_eq!(transpose, Matrix::new([[1,4],[2,5],[3,6]]))
/// ```
#[must_use]
pub fn from_rows<I>(iter: I) -> Self
where
I: IntoIterator<Item = Vector<T, N>>,
Self: Default,
{
let mut result = Self::default();
for (m, row) in iter.into_iter().enumerate().take(M) {
result.set_row(m, &row)
}
result
}
/// Create a matrix from an iterator of vectors
///
/// # Arguments
///
/// * `iter`: iterator of vectors to copy into columns
///
/// returns: Matrix<T, M, N>
///
/// # Examples
///
/// ```
/// # use vector_victor::Matrix;
/// let my_matrix = Matrix::new([[1,2,3],[4,5,6]]);
/// let transpose : Matrix<_,3,2>= Matrix::from_cols(my_matrix.rows());
/// assert_eq!(transpose, Matrix::new([[1,4],[2,5],[3,6]]))
/// ```
#[must_use]
pub fn from_cols<I>(iter: I) -> Self
where
I: IntoIterator<Item = Vector<T, M>>,
Self: Default,
{
let mut result = Self::default();
for (n, col) in iter.into_iter().enumerate().take(N) {
result.set_col(n, &col)
}
result
}
/// Returns an iterator over the elements of the matrix in row-major order.
///
/// # Examples
/// ```
/// # use vector_victor::Matrix;
/// let my_matrix = Matrix::new([[1,2],[3,4]]);
/// assert!(vec![1,2,3,4].iter().eq(my_matrix.elements()))
/// ```
#[must_use]
pub fn elements<'a>(&'a self) -> impl Iterator<Item = &'a T> + 'a {
self.data.iter().flatten()
}
/// Returns a mutable iterator over the elements of the matrix in row-major order.
#[must_use]
pub fn elements_mut<'a>(&'a mut self) -> impl Iterator<Item = &'a mut T> + 'a {
self.data.iter_mut().flatten()
}
/// Returns a reference to the element at that position in the matrix, or `None` if out of bounds.
///
/// # Examples
///
/// ```
/// # use vector_victor::Matrix;
/// let my_matrix = Matrix::new([[1,2],[3,4]]);
///
/// // element at index 2 is the same as the element at (row 1, column 0).
/// assert_eq!(my_matrix.get(2), my_matrix.get((1,0)));
///
/// // my_matrix.get() is equivalent to my_matrix[],
/// // but returns an Option instead of panicking
/// assert_eq!(my_matrix.get(2), Some(&my_matrix[2]));
///
/// // index 4 is out of range, so get(4) returns None.
/// assert_eq!(my_matrix.get(4), None);
/// ```
#[inline]
#[must_use]
pub fn get(&self, index: impl Index2D) -> Option<&T> {
let (m, n) = index.to_2d(M, N)?;
Some(&self.data[m][n])
}
/// Returns a mutable reference to the element at that position in the matrix, or `None` if out of bounds.
#[inline]
#[must_use]
pub fn get_mut(&mut self, index: impl Index2D) -> Option<&mut T> {
let (m, n) = index.to_2d(M, N)?;
Some(&mut self.data[m][n])
}
/// Returns a row of the matrix. or [None] if index is out of bounds
///
/// # Examples
///
/// ```
/// # use vector_victor::{Matrix, Vector};
/// let my_matrix = Matrix::new([[1,2],[3,4]]);
///
/// // row at index 1
/// assert_eq!(my_matrix.row(1), Vector::vec([3,4]));
/// ```
#[inline]
#[must_use]
pub fn row(&self, m: usize) -> Vector<T, N> {
assert!(
m < M,
"Row index {} out of bounds for {}x{} matrix",
m,
M,
N
);
Vector::<T, N>::vec(self.data[m])
}
#[inline]
pub fn set_row(&mut self, m: usize, val: &Vector<T, N>) {
assert!(
m < M,
"Row index {} out of bounds for {}x{} matrix",
m,
M,
N
);
for n in 0..N {
self.data[m][n] = val.data[n][0];
}
}
pub fn pivot_row(&mut self, m1: usize, m2: usize) {
let tmp = self.row(m2);
self.set_row(m2, &self.row(m1));
self.set_row(m1, &tmp);
}
#[inline]
#[must_use]
pub fn col(&self, n: usize) -> Vector<T, M> {
assert!(
n < N,
"Column index {} out of bounds for {}x{} matrix",
n,
M,
N
);
Vector::<T, M>::vec(self.data.map(|r| r[n]))
}
#[inline]
pub fn set_col(&mut self, n: usize, val: &Vector<T, M>) {
assert!(
n < N,
"Column index {} out of bounds for {}x{} matrix",
n,
M,
N
);
for m in 0..M {
self.data[m][n] = val.data[m][0];
}
}
pub fn pivot_col(&mut self, n1: usize, n2: usize) {
let tmp = self.col(n2);
self.set_col(n2, &self.col(n1));
self.set_col(n1, &tmp);
}
#[must_use]
pub fn rows<'a>(&'a self) -> impl Iterator<Item = Vector<T, N>> + 'a {
(0..M).map(|m| self.row(m))
}
#[must_use]
pub fn cols<'a>(&'a self) -> impl Iterator<Item = Vector<T, M>> + 'a {
(0..N).map(|n| self.col(n))
}
#[must_use]
pub fn permute_rows(&self, ms: &Vector<usize, M>) -> Self
where
T: Default,
{
Self::from_rows(ms.elements().map(|&m| self.row(m)))
}
#[must_use]
pub fn permute_cols(&self, ns: &Vector<usize, N>) -> Self
where
T: Default,
{
Self::from_cols(ns.elements().map(|&n| self.col(n)))
}
pub fn transpose(&self) -> Matrix<T, N, M>
where
Matrix<T, N, M>: Default,
{
Matrix::<T, N, M>::from_rows(self.cols())
}
pub fn abs(&self) -> Self
where
T: Default + PartialOrd + Zero + Neg<Output = T>,
{
self.elements()
.map(|&x| match x > T::zero() {
true => x,
false => -x,
})
.collect()
}
}
// 1D vector implementations
impl<T: Copy, const N: usize> Vector<T, N> {
/// Create a vector from a 1D array.
/// Note that vectors are always column vectors unless explicitly instantiated as row vectors
///
/// # Examples
/// ```
/// # use vector_victor::{Matrix, Vector};
/// let my_vector = Vector::vec([1,2,3,4]);
/// // is equivalent to
/// assert_eq!(my_vector, Matrix::new([[1],[2],[3],[4]]));
/// ```
pub fn vec(data: [T; N]) -> Self {
assert!(N > 0, "Vector must have at least 1 element");
return Vector::<T, N> {
data: data.map(|e| [e]),
};
}
pub fn dot<R>(&self, rhs: &R) -> T
where
for<'s> &'s Self: Mul<&'s R, Output = Self>,
T: Sum<T>,
{
(self * rhs).elements().cloned().sum()
}
}
// Cross Product
impl<T: Copy> Vector<T, 3> {
pub fn cross_r<R: Copy>(&self, rhs: &Vector<R, 3>) -> Self
where
T: NumOps<R> + NumOps,
{
Self::vec([
(self[1] * rhs[2]) - (self[2] * rhs[1]),
(self[2] * rhs[0]) - (self[0] * rhs[2]),
(self[0] * rhs[1]) - (self[1] * rhs[0]),
])
}
pub fn cross_l<R: Copy>(&self, rhs: &Vector<R, 3>) -> Vector<R, 3>
where
R: NumOps<T> + NumOps,
{
rhs.cross_r(self)
}
}
//Matrix Multiplication
impl<T: Copy, const M: usize, const N: usize> Matrix<T, M, N> {
pub fn mmul<R: Copy, const P: usize>(&self, rhs: &Matrix<R, N, P>) -> Matrix<T, M, P>
where
T: Default + NumOps<R> + Sum,
{
let mut result: Matrix<T, M, P> = Default::default();
for (m, a) in self.rows().enumerate() {
for (n, b) in rhs.cols().enumerate() {
result[(m, n)] = a.dot(&b)
}
}
return result;
}
}
// Square matrix implementations
impl<T: Copy, const N: usize> Matrix<T, N, N> {
/// Create an identity matrix
#[must_use]
pub fn identity() -> Self
where
T: Zero + One,
{
let mut result = Self::zero();
for i in 0..N {
result[(i, i)] = T::one();
}
return result;
}
/// returns an iterator over the elements along the diagonal of a square matrix
#[must_use]
pub fn diagonals<'s>(&'s self) -> impl Iterator<Item = T> + 's {
(0..N).map(|n| self[(n, n)])
}
/// Returns an iterator over the elements directly below the diagonal of a square matrix
#[must_use]
pub fn subdiagonals<'s>(&'s self) -> impl Iterator<Item = T> + 's {
(0..N - 1).map(|n| self[(n, n + 1)])
}
}
// Index
impl<I, T, const M: usize, const N: usize> Index<I> for Matrix<T, M, N>
where
I: Index2D,
T: Copy,
{
type Output = T;
#[inline(always)]
fn index(&self, index: I) -> &Self::Output {
self.get(index).expect(&*format!(
"index {:?} out of range for {}x{} Matrix",
index, M, N
))
}
}
// IndexMut
impl<I, T, const M: usize, const N: usize> IndexMut<I> for Matrix<T, M, N>
where
I: Index2D,
T: Copy,
{
#[inline(always)]
fn index_mut(&mut self, index: I) -> &mut Self::Output {
self.get_mut(index).expect(&*format!(
"index {:?} out of range for {}x{} Matrix",
index, M, N
))
}
}
// Default
impl<T: Copy + Default, const M: usize, const N: usize> Default for Matrix<T, M, N> {
fn default() -> Self {
Matrix::fill(T::default())
}
}
// Zero
impl<T: Copy + Zero, const M: usize, const N: usize> Zero for Matrix<T, M, N> {
fn zero() -> Self {
Matrix::fill(T::zero())
}
fn is_zero(&self) -> bool {
self.elements().all(|e| e.is_zero())
}
}
// One
impl<T: Copy + One, const M: usize, const N: usize> One for Matrix<T, M, N> {
fn one() -> Self {
Matrix::fill(T::one())
}
}
impl<T: Copy, const M: usize, const N: usize> From<[[T; N]; M]> for Matrix<T, M, N> {
fn from(data: [[T; N]; M]) -> Self {
Self::new(data)
}
}
impl<T: Copy, const M: usize> From<[T; M]> for Vector<T, M> {
fn from(data: [T; M]) -> Self {
Self::vec(data)
}
}
impl<T: Copy, const M: usize, const N: usize> From<T> for Matrix<T, M, N> {
fn from(scalar: T) -> Self {
Self::fill(scalar)
}
}
// deref 1x1 matrices to a scalar automatically
impl<T: Copy> Deref for Matrix<T, 1, 1> {
type Target = T;
fn deref(&self) -> &Self::Target {
&self.data[0][0]
}
}
// deref 1x1 matrices to a mutable scalar automatically
impl<T: Copy> DerefMut for Matrix<T, 1, 1> {
fn deref_mut(&mut self) -> &mut Self::Target {
&mut self.data[0][0]
}
}
// IntoIter
impl<T: Copy, const M: usize, const N: usize> IntoIterator for Matrix<T, M, N> {
type Item = T;
type IntoIter = Flatten<std::array::IntoIter<[T; N], M>>;
fn into_iter(self) -> Self::IntoIter {
self.data.into_iter().flatten()
}
}
// FromIterator
impl<T: Copy, const M: usize, const N: usize> FromIterator<T> for Matrix<T, M, N>
where
Self: Default,
{
fn from_iter<I: IntoIterator<Item = T>>(iter: I) -> Self {
let mut result: Self = Default::default();
for (l, r) in zip(result.elements_mut(), iter) {
*l = r;
}
result
}
}
impl<T: Copy + AddAssign, const M: usize, const N: usize> Sum for Matrix<T, M, N>
where
Self: Zero + Add<Output = Self>,
{
fn sum<I: Iterator<Item = Self>>(iter: I) -> Self {
iter.fold(Self::zero(), Self::add)
}
}
impl<T: Copy + MulAssign, const M: usize, const N: usize> Product for Matrix<T, M, N>
where
Self: One + Mul<Output = Self>,
{
fn product<I: Iterator<Item = Self>>(iter: I) -> Self {
iter.fold(Self::one(), Self::mul)
}
}

1
src/mod.rs Normal file
View File

@ -0,0 +1 @@

2
src/matrix/ops.rs → src/ops.rs
View File

@ -1,4 +1,4 @@
use crate::matrix::Matrix;
use crate::Matrix;
use num_traits::Num;
// borrowed from the auto_ops crate

42
tests/decompose.rs
View File

@ -6,16 +6,18 @@ use generic_parameterize::parameterize;
use num_traits::real::Real;
use num_traits::Zero;
use std::fmt::Debug;
use std::iter::{Product, Sum};
use vector_victor::decompose::Parity::Even;
use vector_victor::decompose::{LUDecomposable, LUDecomposition};
use vector_victor::decompose::{LUDecompose, LUDecomposition, Parity};
use vector_victor::{Matrix, Vector};
#[parameterize(S = (f32, f64), M = [1,2,3,4])]
#[test]
/// The LU decomposition of the identity matrix should produce
/// the identity matrix with no permutations and parity 1
fn test_lu_identity<S: Default + Approx + Real + Debug + Product + Sum, const M: usize>() {
fn test_lu_identity<S, const M: usize>()
where
Matrix<S, M, M>: LUDecompose<S, M>,
S: Copy + Real + Debug + Approx + Default,
{
// let a: Matrix<f32, 3, 3> = Matrix::<f32, 3, 3>::identity();
let i = Matrix::<S, M, M>::identity();
let ones = Vector::<S, M>::fill(S::one());
@ -26,7 +28,7 @@ fn test_lu_identity<S: Default + Approx + Real + Debug + Product + Sum, const M:
(0..M).eq(idx.elements().cloned()),
"Incorrect permutation matrix",
);
assert_eq!(parity, Even, "Incorrect permutation parity");
assert_eq!(parity, Parity::Even, "Incorrect permutation parity");
// Check determinant calculation which uses LU decomposition
assert_approx!(
@ -37,7 +39,7 @@ fn test_lu_identity<S: Default + Approx + Real + Debug + Product + Sum, const M:
// Check inverse calculation with uses LU decomposition
assert_eq!(
i.inverse(),
i.inv(),
Some(i),
"Identity matrix should be its own inverse"
);
@ -54,7 +56,11 @@ fn test_lu_identity<S: Default + Approx + Real + Debug + Product + Sum, const M:
#[parameterize(S = (f32, f64), M = [2,3,4])]
#[test]
/// The LU decomposition of any singular matrix should be `None`
fn test_lu_singular<S: Default + Real + Debug + Product + Sum, const M: usize>() {
fn test_lu_singular<S, const M: usize>()
where
Matrix<S, M, M>: LUDecompose<S, M>,
S: Copy + Real + Debug + Approx + Default,
{
// let a: Matrix<f32, 3, 3> = Matrix::<f32, 3, 3>::identity();
let mut a = Matrix::<S, M, M>::zero();
let ones = Vector::<S, M>::fill(S::one());
@ -66,7 +72,7 @@ fn test_lu_singular<S: Default + Real + Debug + Product + Sum, const M: usize>()
S::zero(),
"Singular matrix should have determinant of zero"
);
assert_eq!(a.inverse(), None, "Singular matrix should have no inverse");
assert_eq!(a.inv(), None, "Singular matrix should have no inverse");
assert_eq!(
a.solve(&ones),
None,
@ -76,7 +82,7 @@ fn test_lu_singular<S: Default + Real + Debug + Product + Sum, const M: usize>()
#[test]
fn test_lu_2x2() {
let a = Matrix::new([[1.0, 2.0], [3.0, 0.0]]);
let a = Matrix::mat([[1.0, 2.0], [3.0, 0.0]]);
let decomp = a.lu().expect("Singular matrix encountered");
// the decomposition is non-unique, due to the combination of lu and idx.
// Instead of checking the exact value, we only check the results.
@ -90,16 +96,16 @@ fn test_lu_2x2() {
assert_approx!(a.det(), decomp.det());
assert_approx!(
a.inverse().unwrap(),
Matrix::new([[0.0, 2.0], [3.0, -1.0]]) * (1.0 / 6.0)
a.inv().unwrap(),
Matrix::mat([[0.0, 2.0], [3.0, -1.0]]) * (1.0 / 6.0)
);
assert_approx!(a.inverse().unwrap(), decomp.inverse());
assert_approx!(a.inverse().unwrap().inverse().unwrap(), a)
assert_approx!(a.inv().unwrap(), decomp.inv());
assert_approx!(a.inv().unwrap().inv().unwrap(), a)
}
#[test]
fn test_lu_3x3() {
let a = Matrix::new([[1.0, -5.0, 8.0], [1.0, -2.0, 1.0], [2.0, -1.0, -4.0]]);
let a = Matrix::mat([[1.0, -5.0, 8.0], [1.0, -2.0, 1.0], [2.0, -1.0, -4.0]]);
let decomp = a.lu().expect("Singular matrix encountered");
let (l, u) = decomp.separate();
@ -109,9 +115,9 @@ fn test_lu_3x3() {
assert_approx!(a.det(), decomp.det());
assert_approx!(
a.inverse().unwrap(),
Matrix::new([[9.0, -28.0, 11.0], [6.0, -20.0, 7.0], [3.0, -9.0, 3.0]]) * (1.0 / 3.0)
a.inv().unwrap(),
Matrix::mat([[9.0, -28.0, 11.0], [6.0, -20.0, 7.0], [3.0, -9.0, 3.0]]) * (1.0 / 3.0)
);
assert_approx!(a.inverse().unwrap(), decomp.inverse());
assert_approx!(a.inverse().unwrap().inverse().unwrap(), a)
assert_approx!(a.inv().unwrap(), decomp.inv());
assert_approx!(a.inv().unwrap().inv().unwrap(), a)
}