Huge refactor to flatten everything

but no API change I think

src/decompose.rs

@@ -69,7 +69,7 @@ pub struct LUDecomposition<T: Copy, const N: usize> {
     /// 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

src/identities.rs

@@ -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())
+    }
+}

src/lib.rs

@@ -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
+    }
+}

src/math.rs

@@ -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()
+    }
+}

src/matrix/mod.rs

@@ -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)
-    }
-}

src/mod.rs

@@ -0,0 +1 @@
+

src/ops.rs

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

tests/decompose.rs

@@ -4,9 +4,8 @@ mod common;
 use common::Approx;
 use generic_parameterize::parameterize;
 use num_traits::real::Real;
-use num_traits::{Float, One, Signed, Zero};
+use num_traits::Zero;
 use std::fmt::Debug;
-use std::iter::{Product, Sum};
 use vector_victor::decompose::{LUDecompose, LUDecomposition, Parity};
 use vector_victor::{Matrix, Vector};