Add LU tests

nalgebra-sparse/src/csc.rs

@@ -658,6 +658,7 @@ impl<T> CscMatrix<T> {
     /// Assumes that the diagonal of the sparse matrix is all 1.
     pub fn sparse_lower_triangular_solve_sorted(
         &self,
+        // input vector idxs & values
         b_idxs: &[usize],
         b: &[T],
         // idx -> row
@@ -682,24 +683,41 @@ impl<T> CscMatrix<T> {
             let out_pos = out_sparsity_pattern.iter().position(|&p| p == bi).unwrap();
             out[out_pos] = bv;
         }
+        // end init
 
+        // assuming that the output sparsity pattern is sorted
+        // iterate thru
         for (i, &row) in out_sparsity_pattern.iter().enumerate() {
             let col = self.col(row);
             let mut iter = col.row_indices().iter().zip(col.values().iter()).peekable();
             if !assume_unit {
-                match iter.find(|v| *v.0 >= row) {
+                while let Some(n) = iter.peek() {
+                    if *n.0 < row {
+                        iter.next();
+                    } else {
+                        break;
+                    }
+                }
+                match iter.peek() {
                     Some((&r, &l_val)) if r == row => out[i] /= l_val,
+                    // here it now becomes implicitly 0,
+                    // likely this should introduce NaN or some other behavior.
                     _ => {}
                 }
             }
             let mul = out[i];
-            for (offset, &nrow) in out_sparsity_pattern[i..].iter().enumerate() {
-                if nrow <= row {
-                    continue;
+            for (offset, &nrow) in out_sparsity_pattern[i..].iter().enumerate().skip(1) {
+                assert!(nrow > row);
+                while let Some(n) = iter.peek() {
+                    if *n.0 < nrow {
+                        iter.next();
+                    } else {
+                        break;
+                    }
                 }
-                let l_val = match iter.find(|v| *v.0 >= nrow) {
+                let l_val = match iter.peek() {
                     Some((&r, &l_val)) if r == nrow => l_val,
-                    _ => break,
+                    _ => continue,
                 };
                 out[i + offset] -= l_val * mul;
             }

nalgebra-sparse/src/factorization/lu.rs

@@ -1,3 +1,4 @@
+use crate::SparseEntryMut;
 use crate::{CscBuilder, CscMatrix};
 use nalgebra::RealField;
 
@@ -19,14 +20,19 @@ impl<T: RealField + Copy> LeftLookingLUFactorization<T> {
         &self.l_u
     }
 
-    /*
-    /// Returns the upper triangular part of this matrix.
-    pub fn u(&self) -> CscMatrix<T> {
-      self.l_u.lower_triangle()
-      // TODO here need to change the diagonal entries to be 1.
-      todo!();
+    /// Returns the lower triangular part of this matrix.
+    pub fn l(&self) -> CscMatrix<T> {
+        let mut l = self.l_u.lower_triangle();
+        let n = self.l_u.nrows();
+        for i in 0..n {
+            if let SparseEntryMut::NonZero(v) = l.index_entry_mut(i, i) {
+                *v = T::one();
+            } else {
+                unreachable!();
+            }
+        }
+        l
     }
-    */
 
     /// Construct a new sparse LU factorization
     /// from a given CSC matrix.

nalgebra-sparse/tests/unit_tests/cholesky.rs

@@ -9,7 +9,7 @@ use nalgebra::proptest::matrix;
 use proptest::prelude::*;
 use matrixcompare::{assert_matrix_eq, prop_assert_matrix_eq};
 
-fn positive_definite() -> impl Strategy<Value=CscMatrix<f64>> {
+pub fn positive_definite() -> impl Strategy<Value=CscMatrix<f64>> {
     let csc_f64 = csc(value_strategy::<f64>(),
                       PROPTEST_MATRIX_DIM,
                       PROPTEST_MATRIX_DIM,
@@ -114,4 +114,4 @@ fn test_cholesky(a: Matrix5<f64>) {
     let l = DMatrix::from_iterator(l.nrows(), l.ncols(), l.iter().cloned());
     let cs_l_mat = DMatrix::from(&cs_l);
     assert_matrix_eq!(l, cs_l_mat, comp = abs, tol = 1e-12);
-}
+}

nalgebra-sparse/tests/unit_tests/left_looking_lu.proptest-regressions

@@ -0,0 +1,7 @@
+# Seeds for failure cases proptest has generated in the past. It is
+# automatically read and these particular cases re-run before any
+# novel cases are generated.
+#
+# It is recommended to check this file in to source control so that
+# everyone who runs the test benefits from these saved cases.
+cc b8e49f3fa66a3568e1c29a47de5f3d2f67866d59b39329e2a0da1537a9c6d584 # shrinks to matrix = CscMatrix { cs: CsMatrix { sparsity_pattern: SparsityPattern { major_offsets: [0, 4, 8, 11, 16, 21], minor_indices: [0, 1, 3, 4, 0, 1, 3, 4, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4], minor_dim: 5 }, values: [1.0, 0.0, 0.0, 0.0, 0.0, 8.445226477445164, 0.0, 5.072108001361104, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 5.072108001361104, 0.0, 0.0, 4.455405910808415] } }

nalgebra-sparse/tests/unit_tests/left_looking_lu.rs

@@ -0,0 +1,41 @@
+use crate::unit_tests::cholesky::positive_definite;
+
+use nalgebra_sparse::factorization::LeftLookingLUFactorization;
+
+use matrixcompare::{assert_matrix_eq, prop_assert_matrix_eq};
+use proptest::prelude::*;
+
+use nalgebra_sparse::CscMatrix;
+
+proptest! {
+  #[test]
+  fn lu_factorization_correct_for_positive_def_matrices(
+    matrix in positive_definite()
+  ) {
+    let lu = LeftLookingLUFactorization::new(&matrix);
+
+    let l = lu.l();
+    prop_assert!(l.triplet_iter().all(|(i, j, _)| j <= i));
+    let u = lu.u();
+    prop_assert!(u.triplet_iter().all(|(i, j, _)| j >= i));
+    prop_assert_matrix_eq!(l * u, matrix, comp = abs, tol = 1e-7);
+  }
+}
+
+#[test]
+fn minimized_lu() {
+    let major_offsets = vec![0, 1, 3, 4, 5, 8];
+    let minor_indices = vec![0, 1, 4, 2, 3, 1, 2, 4];
+    let values = vec![1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 3.0, 1.0];
+    assert_eq!(minor_indices.len(), values.len());
+    let mat = CscMatrix::try_from_unsorted_csc_data(5, 5, major_offsets, minor_indices, values);
+    let mat = mat.unwrap();
+
+    let lu = LeftLookingLUFactorization::new(&mat);
+
+    let l = lu.l();
+    assert!(l.triplet_iter().all(|(i, j, _)| j <= i));
+    let u = lu.u();
+    assert!(u.triplet_iter().all(|(i, j, _)| j >= i));
+    assert_matrix_eq!(l * u, mat, comp = abs, tol = 1e-7);
+}

nalgebra-sparse/tests/unit_tests/mod.rs

@@ -1,4 +1,7 @@
+// factorization tests
 mod cholesky;
+mod left_looking_lu;
+
 mod convert_serial;
 mod coo;
 mod csc;