Merge pull request #297 from jturner314/fix-factorize_into_inv

Fix Inverse for LUFactorized

Author: termoshtt

lax/src/solve.rs

@@ -42,6 +42,10 @@ macro_rules! impl_solve {
 
             fn inv(l: MatrixLayout, a: &mut [Self], ipiv: &Pivot) -> Result<()> {
                 let (n, _) = l.size();
+                if n == 0 {
+                    // Do nothing for empty matrices.
+                    return Ok(());
+                }
 
                 // calc work size
                 let mut info = 0;

ndarray-linalg/src/solve.rs

@@ -199,9 +199,9 @@ pub trait Solve<A: Scalar> {
 pub struct LUFactorized<S: Data + RawDataClone> {
     /// The factors `L` and `U`; the unit diagonal elements of `L` are not
     /// stored.
-    pub a: ArrayBase<S, Ix2>,
+    a: ArrayBase<S, Ix2>,
     /// The pivot indices that define the permutation matrix `P`.
-    pub ipiv: Pivot,
+    ipiv: Pivot,
 }
 
 impl<A, S> Solve<A> for LUFactorized<S>
@@ -387,8 +387,15 @@ where
     type Output = Array2<A>;
 
     fn inv(&self) -> Result<Array2<A>> {
+        // Preserve the existing layout. This is required to obtain the correct
+        // result, because the result of `A::inv` is layout-dependent.
+        let a = if self.a.is_standard_layout() {
+            replicate(&self.a)
+        } else {
+            replicate(&self.a.t()).reversed_axes()
+        };
         let f = LUFactorized {
-            a: replicate(&self.a),
+            a,
             ipiv: self.ipiv.clone(),
         };
         f.inv_into()

ndarray-linalg/tests/inv.rs

@@ -1,20 +1,103 @@
 use ndarray::*;
 use ndarray_linalg::*;
 
+fn test_inv_random<A>(n: usize, set_f: bool, rtol: A::Real)
+where
+    A: Scalar + Lapack,
+{
+    let a: Array2<A> = random([n; 2].set_f(set_f));
+    let identity = Array2::eye(n);
+    assert_close_l2!(&a.inv().unwrap().dot(&a), &identity, rtol);
+    assert_close_l2!(
+        &a.factorize().unwrap().inv().unwrap().dot(&a),
+        &identity,
+        rtol
+    );
+    assert_close_l2!(
+        &a.clone().factorize_into().unwrap().inv().unwrap().dot(&a),
+        &identity,
+        rtol
+    );
+}
+
+fn test_inv_into_random<A>(n: usize, set_f: bool, rtol: A::Real)
+where
+    A: Scalar + Lapack,
+{
+    let a: Array2<A> = random([n; 2].set_f(set_f));
+    let identity = Array2::eye(n);
+    assert_close_l2!(&a.clone().inv_into().unwrap().dot(&a), &identity, rtol);
+    assert_close_l2!(
+        &a.factorize().unwrap().inv_into().unwrap().dot(&a),
+        &identity,
+        rtol
+    );
+    assert_close_l2!(
+        &a.clone()
+            .factorize_into()
+            .unwrap()
+            .inv_into()
+            .unwrap()
+            .dot(&a),
+        &identity,
+        rtol
+    );
+}
+
+#[test]
+fn inv_empty() {
+    test_inv_random::<f32>(0, false, 0.);
+    test_inv_random::<f64>(0, false, 0.);
+    test_inv_random::<c32>(0, false, 0.);
+    test_inv_random::<c64>(0, false, 0.);
+}
+
+#[test]
+fn inv_random_float() {
+    for n in 1..=8 {
+        for &set_f in &[false, true] {
+            test_inv_random::<f32>(n, set_f, 1e-3);
+            test_inv_random::<f64>(n, set_f, 1e-9);
+        }
+    }
+}
+
+#[test]
+fn inv_random_complex() {
+    for n in 1..=8 {
+        for &set_f in &[false, true] {
+            test_inv_random::<c32>(n, set_f, 1e-3);
+            test_inv_random::<c64>(n, set_f, 1e-9);
+        }
+    }
+}
+
+#[test]
+fn inv_into_empty() {
+    test_inv_into_random::<f32>(0, false, 0.);
+    test_inv_into_random::<f64>(0, false, 0.);
+    test_inv_into_random::<c32>(0, false, 0.);
+    test_inv_into_random::<c64>(0, false, 0.);
+}
+
 #[test]
-fn inv_random() {
-    let a: Array2<f64> = random((3, 3));
-    let ai: Array2<_> = (&a).inv().unwrap();
-    let id = Array::eye(3);
-    assert_close_l2!(&ai.dot(&a), &id, 1e-7);
+fn inv_into_random_float() {
+    for n in 1..=8 {
+        for &set_f in &[false, true] {
+            test_inv_into_random::<f32>(n, set_f, 1e-3);
+            test_inv_into_random::<f64>(n, set_f, 1e-9);
+        }
+    }
 }
 
 #[test]
-fn inv_random_t() {
-    let a: Array2<f64> = random((3, 3).f());
-    let ai: Array2<_> = (&a).inv().unwrap();
-    let id = Array::eye(3);
-    assert_close_l2!(&ai.dot(&a), &id, 1e-7);
+fn inv_into_random_complex() {
+    for n in 1..=8 {
+        for &set_f in &[false, true] {
+            test_inv_into_random::<c32>(n, set_f, 1e-3);
+            test_inv_into_random::<c64>(n, set_f, 1e-9);
+        }
+    }
 }
 
 #[test]