Add rref and use Gauss-Jordan Elimination for Inverse

src/err.rs

@@ -19,8 +19,7 @@ impl Display for MatrixError {
         match self {
             MatrixError::IndexOutOfBounds(idx) => write!(
                 f,
-                "Tried to access a matrix at index `{}`, which is out of bounds.",
-                idx
+                "Tried to access a matrix at index `{idx}`, which is out of bounds.",
             )?,
         }
         Ok(())
@@ -33,12 +32,11 @@ impl Display for DimensionError {
             DimensionError::InvalidDimensions => {
                 write!(f, "Dimensions with a size of less than 1 are invalid.")?
             }
-            DimensionError::NoMatch(dims, bad_dims, op) => {write!(
+            DimensionError::NoMatch(dims, bad_dims, op) => write!(
                 f,
-                "Dimensions of two matrices do not match in the correct way. Cannot {} {} matrix with {} matrix.",
-                op, dims, bad_dims
-            )?}
-            DimensionError::InvalidInputDimensions(input_len, correct_len) => {write!(f, "Invalid input dimensions. Input has length {}, but should have length {}.", input_len, correct_len)?}
+                "Dimensions of two matrices do not match in the correct way. Cannot {op} {dims} matrix with {bad_dims} matrix.",
+            )?,
+            DimensionError::InvalidInputDimensions(input_len, correct_len) => write!(f, "Invalid input dimensions. Input has length {input_len}, but should have length {correct_len}.")?,
             DimensionError::NoSquare => {
                 write!(f, "Not a square matrix. Rows and cols need to be the same.")?
             }

src/mat/_mat/mat_impl.rs

@@ -168,6 +168,7 @@ where
     }
 
     /// Creates a diagonal matrix with dimensions `dim x dim` and initial entries specified in `entries`.
+    #[allow(clippy::manual_memcpy)]
     pub fn diag_with(dim: usize, entries: &[T]) -> Result<Matrix<T>, DimensionError> {
         if entries.len() != dim {
             return Err(DimensionError::InvalidInputDimensions(entries.len(), dim));
@@ -308,6 +309,71 @@ where
         }
         Matrix::<T>::from_vec(self.cols(), self.rows(), vec).unwrap()
     }
+
+    /// Find Reduced Row Echelon Form.
+    ///
+    /// # Example
+    ///
+    /// ```
+    /// # use libmat::mat::Matrix;
+    /// # use libmat::matrix;
+    /// let mat_a = matrix!{1, 2, 3, 4; 5, 6, 7, 8; 9, 10, 11, 12};
+    /// // 1  2  3  4
+    /// // 5  6  7  8
+    /// // 9 10 11 12
+    /// let mat_b = matrix!{1, 0, -1, -2; 0, 1, 2, 3; 0, 0, 0, 0};
+    /// // 1  0 -1 -2
+    /// // 0  1  2  3
+    /// // 0  0  0  0
+    /// assert_eq!(mat_a.rref(), mat_b);
+    pub fn rref(&self) -> Matrix<T>
+    where
+        T: sign::Signed + std::ops::DivAssign + std::ops::SubAssign + Clone + Zero + One,
+    {
+        let mut mat = self.clone();
+        let mut col = 0;
+        let mut row = 0;
+        while row < mat.rows() && col < mat.cols() {
+            if mat[row][col].is_zero() {
+                // find non-zero
+                for r in row..mat.rows() {
+                    if !mat[r][r].is_zero() {
+                        // swap r -> row
+                        for (i, item) in mat[row].to_vec().iter().cloned().enumerate() {
+                            mat[row][i] = mat[r][i].clone();
+                            mat[r][i] = item;
+                        }
+                        break;
+                    }
+                }
+            }
+            if mat[row][col].is_zero() {
+                col += 1;
+                continue;
+            }
+            // ensure first item is 1
+            if !mat[row][col].is_one() {
+                let val = mat[row][col].clone();
+                for c in col..mat.cols() {
+                    mat[row][c] /= val.clone();
+                }
+            }
+            // reduce all other rows
+            for r in 0..mat.rows() {
+                if mat[r][col].is_zero() || r == row {
+                    continue;
+                }
+                let val = mat[r][col].clone();
+                for c in col..mat.cols() {
+                    let x = mat[row][c].clone();
+                    mat[r][c] -= val.clone() * x;
+                }
+            }
+            row += 1;
+            col += 1;
+        }
+        mat
+    }
 }
 
 impl<T> From<Vector<T>> for Matrix<T>
@@ -346,7 +412,7 @@ impl<T> Matrix<T> {
         self.matrix[self.cols() * i.into() + j.into()].clone()
     }
 
-    pub fn entry_mut<'a>(&'a mut self, i: impl Into<usize>, j: impl Into<usize>) -> &'a mut T {
+    pub fn entry_mut(&mut self, i: impl Into<usize>, j: impl Into<usize>) -> &mut T {
         let cols = self.cols();
         &mut self.matrix[cols * i.into() + j.into()]
     }

src/mat/_mat/mat_traits.rs

@@ -16,11 +16,11 @@ where
             for j in 0..self.cols() {
                 let n = &self.matrix[i * self.cols() + j];
                 if j == self.cols() - 1 && i == self.rows() - 1 {
-                    write!(f, "{}", n)?;
+                    write!(f, "{n}")?;
                 } else if j == self.cols() - 1 {
-                    writeln!(f, "{}", n)?;
+                    writeln!(f, "{n}")?;
                 } else {
-                    write!(f, "{}\t", n)?;
+                    write!(f, "{n}\t")?;
                 }
             }
         }
@@ -30,7 +30,7 @@ where
 
 impl<T> Inv for Matrix<T>
 where
-    T: One + Zero + Clone + Signed + PartialOrd + std::iter::Sum + std::ops::DivAssign,
+    T: Signed + PartialOrd + std::ops::DivAssign + std::ops::SubAssign + Clone + Zero + One,
 {
     type Output = Result<Option<Matrix<T>>, DimensionError>;
 
@@ -44,47 +44,47 @@ where
     /// # use num_traits::ops::inv::Inv;
     /// # use libmat::err::DimensionError;
     /// # fn main() -> Result<(), DimensionError> {
-    /// let mat_a: Matrix<f32> = matrix!{{0.0,-1.0,2.0},{1.0,2.0,0.0},{2.0,1.0,0.0}};
+    /// let mat_a: Matrix<f32> = matrix!{{1.0,2.0,3.0},{0.0,1.0,4.0},{-5.0,-6.0,0.0}};
     /// let mat_c: Matrix<i32> = matrix!{{1,0,0},{0,1,0},{0,0,0}}; // not invertible
-    /// let mat_b = matrix!{{0.0, -1.0/3.0, 2.0/3.0}, {0.0, 2.0/3.0, -1.0/3.0}, {1.0/2.0, 1.0/3.0, -1.0/6.0}};
+    /// let mat_b = matrix!{{-24.0, 18.0, -5.0}, {20.0, -15.0, 4.0}, {-5.0, 4.0, -1.0}};
     /// assert_eq!(mat_a.inv()?, Some(mat_b));
     /// assert_eq!(mat_c.inv()?, None);
     /// # Ok(()) }
     /// ```
     fn inv(self) -> Self::Output {
-        if let Some((mat, p)) = self.lupdecompose()? {
-            let dim = mat.rows();
-            let mut mat_inv = Matrix::<T>::zero(dim, dim).unwrap();
-            for j in 0..dim {
-                for i in 0..dim {
-                    mat_inv[i][j] = {
-                        if p[i] == j {
-                            T::one()
-                        } else {
-                            T::zero()
-                        }
-                    };
-
-                    for k in 0..i {
-                        mat_inv[i][j] =
-                            mat_inv[i][j].clone() - mat[i][k].clone() * mat_inv[k][j].clone();
-                    }
-                }
+        // This uses the Gauss-Jordan Elimination method
+        if !self.is_square() {
+            return Err(DimensionError::NoSquare);
+        }
+        let dim = self.rows();
+        let mut mat = Matrix::zero(dim, dim*2)?;
 
-                for i in (0..dim).rev() {
-                    for k in (i + 1)..dim {
-                        mat_inv[i][j] =
-                            mat_inv[i][j].clone() - mat[i][k].clone() * mat_inv[k][j].clone();
-                    }
-                    mat_inv[i][j] /= mat[i][i].clone();
+        for i in 0..dim {
+            mat[i][..dim].clone_from_slice(&self[i][..dim]);
+            mat[i][i+dim] = T::one();
+        }
+        mat = mat.rref();
+        for i in 0..dim {
+            for j in 0..dim {
+                if mat[i][j] != if i == j { T::one() } else { T::zero() } {
+                    return Ok(None);
                 }
             }
-            if (p[dim] - dim) % 2 != 0 {
-                mat_inv.matrix.reverse();
-            }
-            Ok(Some(mat_inv))
-        } else {
-            Ok(None)
         }
+
+        let mut inv = Matrix::zero(dim, dim)?;
+        for i in 0..dim {
+            inv[i][..dim].clone_from_slice(&mat[i][dim..2*dim])
+        }
+        Ok(Some(inv))
+    }
+}
+
+impl<T> IntoIterator for Matrix<T> {
+    type Item = T;
+    type IntoIter = std::vec::IntoIter<Self::Item>;
+
+    fn into_iter(self) -> Self::IntoIter {
+        self.matrix.into_iter()
     }
 }

src/mat/smat/smat_impl.rs

@@ -194,6 +194,7 @@ where
     }
 
     /// Creates a diagonal matrix with initial entries specified in `entries`.
+    #[allow(clippy::manual_memcpy)]
     pub fn diag_with(entries: &[T]) -> SMatrix<T, N, N>
     where
         T: One + Copy + Zero + std::iter::Sum,

src/mat/smat/smat_traits.rs

@@ -15,7 +15,7 @@ where
         while let Some(r) = rs.next() {
             let mut es = r.iter().peekable();
             while let Some(e) = es.next() {
-                write!(f, "{}", e)?;
+                write!(f, "{e}")?;
                 if rs.peek().is_some() {
                     write!(f, "")?;
                 } else if es.peek().is_some() {

src/mat/vec/vec_ops.rs

@@ -86,8 +86,7 @@ where
     T: AddAssign + Clone,
 {
     fn add_assign(&mut self, rhs: T) {
-        self.iter_mut()
-            .for_each(|a| *a += rhs.clone());
+        self.iter_mut().for_each(|a| *a += rhs.clone());
     }
 }
 
@@ -171,8 +170,7 @@ where
     T: SubAssign + Clone,
 {
     fn sub_assign(&mut self, rhs: T) {
-        self.iter_mut()
-            .for_each(|a| *a -= rhs.clone());
+        self.iter_mut().for_each(|a| *a -= rhs.clone());
     }
 }
 
@@ -255,9 +253,7 @@ where
     fn mul(self, mat: Matrix<T>) -> Self::Output {
         let vector: Vector<T> = self;
         let mat_v: Matrix<T> = vector.into();
-        println!("{}\n\n{}", mat_v, mat);
         let res = (mat_v * mat)?;
-        println!("\n{}", res);
         Ok(res.into())
     }
 }

tests/complex_tests.rs

@@ -9,3 +9,12 @@ fn double_inverse() -> Result<(), DimensionError> {
     assert_eq!(mat_b.inv()?, Some(mat_a));
     Ok(())
 }
+
+#[test]
+fn double_inverse_3() -> Result<(), DimensionError> {
+    let mat_a = matrix! {{0.0, 1.0, 2.0}, {1.0, 2.0, 3.0}, {3.0, 1.0, 1.0}};
+    let mat_b = matrix! {{0.5, -0.5, 0.5}, {-4.0, 3.0, -1.0}, {2.5, -1.5, 0.5}};
+    assert_eq!(mat_a.clone().inv()?, Some(mat_b.clone()));
+    assert_eq!(mat_b.inv()?, Some(mat_a));
+    Ok(())
+}

tests/determinant.rs

@@ -34,7 +34,10 @@ fn some_dets() -> Result<(), DimensionError> {
             8, 6, 1, 0, 1, 9, 5, 9, 9, 9, 0, 8, 4, 3, 4, 0, 5, 6, 5, 1, 0, 9, 4, 6, 4, 9, 8, 3, 5,
             1, 10, 6, 3, 10, 7, 4, 9, 2, 0, 1, 2, 1, 6, 8, 7, 3, 2, 9, 1, 7, 1, 4, 4, 9, 0, 0, 7,
             6, 4, 0, 10, 4, 5, 9,
-        ].iter().map(|x| (*x as i16).into()).collect(),
+        ]
+        .iter()
+        .map(|x| (*x as i16).into())
+        .collect(),
     )?;
     assert_eq!(b.det()?, -15546220_f32);
     Ok(())