Added is_finite check for 4x4 matrix inversion

benches/lib.rs

@@ -6,8 +6,9 @@ extern crate rand_package as rand;
 #[macro_use]
 extern crate criterion;
 
-use na::DMatrix;
+use na::{DMatrix, Matrix, Scalar};
 use rand::Rng;
+use rand_distr::{Distribution, Standard};
 use rand_isaac::IsaacRng;
 
 pub mod core;
@@ -20,6 +21,19 @@ fn reproducible_dmatrix(nrows: usize, ncols: usize) -> DMatrix<f64> {
     DMatrix::<f64>::from_fn(nrows, ncols, |_, _| rng.gen())
 }
 
+fn reproducible_matrix<T: Scalar + Default, const R: usize, const C: usize>(
+) -> Matrix<T, na::Const<R>, na::Const<C>, na::ArrayStorage<T, R, C>>
+where
+    na::ArrayStorage<T, R, C>: Default,
+    Standard: Distribution<T>,
+{
+    use rand::SeedableRng;
+    let mut rng = IsaacRng::seed_from_u64(0);
+    let mut m: Matrix<T, na::Const<R>, na::Const<C>, na::ArrayStorage<T, R, C>> = Matrix::default();
+    m.iter_mut().for_each(|x| *x = rng.gen());
+    return m;
+}
+
 criterion_main!(
     core::matrix,
     core::vector,
@@ -34,4 +48,5 @@ criterion_main!(
     linalg::solve,
     linalg::svd,
     linalg::symmetric_eigen,
+    linalg::matrix_inversion
 );

benches/linalg/inverse.rs

@@ -0,0 +1,19 @@
+use na::LU;
+
+use crate::reproducible_matrix;
+
+fn inverse_4(bh: &mut criterion::Criterion) {
+    let m = reproducible_matrix::<f64, 4, 4>();
+    bh.bench_function("4x4_matrix_inverse_det", move |bh| {
+        bh.iter(|| std::hint::black_box(m.try_inverse().unwrap()))
+    });
+}
+
+fn inverse_lu(bh: &mut criterion::Criterion) {
+    let m = reproducible_matrix::<f64, 4, 4>();
+    bh.bench_function("4x4_matrix_inverse_lu", move |bh| {
+        bh.iter(|| std::hint::black_box(LU::new(m).try_inverse().unwrap()))
+    });
+}
+
+criterion_group!(matrix_inversion, inverse_4, inverse_lu);

benches/linalg/mod.rs

@@ -2,6 +2,7 @@ pub use self::bidiagonal::bidiagonal;
 pub use self::cholesky::cholesky;
 pub use self::full_piv_lu::full_piv_lu;
 pub use self::hessenberg::hessenberg;
+pub use self::inverse::matrix_inversion;
 pub use self::lu::lu;
 pub use self::qr::qr;
 pub use self::schur::schur;
@@ -13,6 +14,7 @@ mod bidiagonal;
 mod cholesky;
 mod full_piv_lu;
 mod hessenberg;
+mod inverse;
 mod lu;
 mod qr;
 mod schur;

src/linalg/inverse.rs

@@ -10,6 +10,11 @@ use crate::linalg::lu;
 impl<T: ComplexField, D: Dim, S: Storage<T, D, D>> SquareMatrix<T, D, S> {
     /// Attempts to invert this square matrix.
     ///
+    /// # Returns
+    ///
+    /// - `Some<SquareMatrix>`: The inverted matrix, if the inversion succeeds and all entries are finite.
+    /// - `None`: If the inversion fails or produces `Inf` or `NaN` values.
+    ///
     /// # Panics
     ///
     /// Panics if `self` isn’t a square matrix.
@@ -20,7 +25,7 @@ impl<T: ComplexField, D: Dim, S: Storage<T, D, D>> SquareMatrix<T, D, S> {
         DefaultAllocator: Allocator<D, D>,
     {
         let mut me = self.into_owned();
-        if me.try_inverse_mut() {
+        if me.try_inverse_mut() && me.iter().all(|x| x.is_finite()) {
             Some(me)
         } else {
             None
@@ -32,6 +37,8 @@ impl<T: ComplexField, D: Dim, S: StorageMut<T, D, D>> SquareMatrix<T, D, S> {
     /// Attempts to invert this square matrix in-place. Returns `false` and leaves `self` untouched if
     /// inversion fails.
     ///
+    /// Calling this function on near singular matrices may lead to`Inf` or `NaN` values!
+    ///
     /// # Panics
     ///
     /// Panics if `self` isn’t a square matrix.

tests/linalg/inverse.rs

@@ -140,3 +140,70 @@ fn matrix5_try_inverse_scaled_identity() {
 
     assert_relative_eq!(a_inv, expected_inverse);
 }
+
+#[test]
+#[rustfmt::skip]
+fn matrix2_try_inverse_singular(){
+    let m1 = Matrix2::new(
+        1.0, 0.0,
+         0.0,0.0
+    );
+    assert!(m1.try_inverse().is_none());
+}
+
+#[test]
+#[rustfmt::skip]
+fn matrix2_try_inverse_near_singular(){
+    let m1 = Matrix2::new(
+        1.0, 1e1, 
+        0.0,  f64::MIN_POSITIVE);
+    if let Some(inv) = m1.try_inverse(){
+        assert!(inv.iter().all(|x|x.is_finite()));
+    }
+}
+
+#[test]
+#[rustfmt::skip]
+fn matrix3_try_inverse_singular(){
+    let m1 = Matrix3::new(
+        1.0, 0.0, 0.0, 
+        0.0, 1.0, 0.0,
+        0.0, 0.0, 0.0);
+    assert!(m1.try_inverse().is_none());
+}
+
+#[test]
+#[rustfmt::skip]
+fn matrix3_try_inverse_near_singular(){
+    let m1 = Matrix3::new(
+        1.0, 1e1, 1e1, 
+        0.0, 1.0, 1e1,
+        0.0, 0.0, f64::MIN_POSITIVE);
+    if let Some(inv) = m1.try_inverse(){
+        assert!(inv.iter().all(|x|x.is_finite()));
+    }
+}
+
+#[test]
+#[rustfmt::skip]
+fn matrix4_try_inverse_singular(){
+    let m1 = Matrix4::new(
+        1.0, 0.0, 0.0, 0.0, 
+        0.0, 1.0, 0.0, 0.0,
+        0.0, 0.0, 1.0, 0.0,
+        0.0, 0.0, 0.0, 0.0 );
+    assert!(m1.try_inverse().is_none());
+}
+
+#[test]
+#[rustfmt::skip]
+fn matrix4_try_inverse_near_singular(){
+    let m1 = Matrix4::new(
+        1.0, 1e1, 1e1, 1e1, 
+        0.0, 1.0, 1e1, 1e1,
+        0.0, 0.0, 1.0, 1e1,
+        0.0, 0.0, 0.0, f64::MIN_POSITIVE );
+    if let Some(inv) = m1.try_inverse(){
+        assert!(inv.iter().all(|x|x.is_finite()));
+    }
+}