Complete basic version of sparse LU factorization

nalgebra-sparse/src/factorization/lu.rs

@@ -14,6 +14,11 @@ impl<T: RealField + Copy> LeftLookingLUFactorization<T> {
         self.l_u.upper_triangle()
     }
 
+    /// Returns the joint L\U matrix. Here, `L` implicitly has 1 along the diagonal.
+    pub fn lu(&self) -> &CscMatrix<T> {
+        &self.l_u
+    }
+
     /*
     /// Returns the upper triangular part of this matrix.
     pub fn u(&self) -> CscMatrix<T> {
@@ -36,7 +41,7 @@ impl<T: RealField + Copy> LeftLookingLUFactorization<T> {
         let mut val_buf = vec![];
         let mut pat_buf = vec![];
 
-        for (i, col) in a.col_iter().enumerate() {
+        for (ci, col) in a.col_iter().enumerate() {
             let curr_mat = csc_builder.build();
 
             curr_mat
@@ -46,7 +51,7 @@ impl<T: RealField + Copy> LeftLookingLUFactorization<T> {
             val_buf.resize_with(pat_buf.len(), T::zero);
 
             // Solve the current column, assuming that it is lower triangular
-            let x = curr_mat.sparse_lower_triangular_solve(
+            curr_mat.sparse_lower_triangular_solve_sorted(
                 col.row_indices(),
                 col.values(),
                 &pat_buf,
@@ -56,12 +61,26 @@ impl<T: RealField + Copy> LeftLookingLUFactorization<T> {
 
             // convert builder back to matrix
             csc_builder = CscBuilder::from_mat(curr_mat);
-            assert!(csc_builder.revert_to_col(i));
-            let _ = x;
-            todo!();
+            assert!(csc_builder.revert_to_col(ci));
+            let mut ukk = T::zero();
+            for (row, val) in pat_buf.drain(..).zip(val_buf.drain(..)) {
+                use std::cmp::Ordering;
+                let val = match row.cmp(&ci) {
+                    Ordering::Less => val,
+                    Ordering::Equal => {
+                        ukk = val;
+                        val
+                    }
+                    Ordering::Greater => {
+                        assert_ne!(ukk, T::zero());
+                        val / ukk
+                    }
+                };
+                assert_eq!(csc_builder.insert(row, ci, val), Ok(()));
+            }
         }
 
-        //csc_builder.build()
-        todo!();
+        let l_u = csc_builder.build();
+        Self { l_u }
     }
 }

nalgebra-sparse/src/pattern.rs

@@ -383,25 +383,11 @@ impl SparsityPatternBuilder {
             major_dim,
         }
     }
-    /// The number of entries into this sparsity pattern
+    /// The number of non-zero entries inserted into `self`.
     pub fn num_entries(&self) -> usize {
-        self.buf.minor_indices.last().copied().unwrap_or(0)
+        self.buf.minor_indices.len()
     }
 
-    /// The current major dimension of this sparsity pattern builder.
-    pub fn major_dim(&self) -> usize {
-        self.buf.major_dim()
-    }
-    /// The current minor dimension of this sparsity pattern builder.
-    /// Note this doesn't represent the number of minor dimensions,
-    /// But the last major dim's current max minor idx.
-    pub fn current_minor(&self) -> usize {
-        let maj = match self.buf.major_offsets.last() {
-            None => return 0,
-            Some(&maj) => maj,
-        };
-        self.buf.minor_indices.get(maj).copied().unwrap_or(0)
-    }
     /// Allows for general assignment of indices
     pub fn insert(&mut self, maj: usize, min: usize) -> Result<(), BuilderInsertError> {
         assert!(maj < self.major_dim);
@@ -415,6 +401,7 @@ impl SparsityPatternBuilder {
         }
         // cannot go backwards in minor
         if maj == curr_major
+            && *self.buf.major_offsets.last().unwrap() < self.buf.minor_indices.len()
             && !self.buf.minor_indices.is_empty()
             && min <= *self.buf.minor_indices.last().unwrap()
         {
@@ -435,14 +422,6 @@ impl SparsityPatternBuilder {
         }
         &self.buf
     }
-    /// Instead of creating a matrix of the given input size, constructs a matrix
-    /// which is only filled up to the current row.
-    pub fn into_partial(mut self) -> SparsityPattern {
-        if *self.buf.major_offsets.last().unwrap() != self.buf.minor_indices.len() {
-            self.buf.major_offsets.push(self.buf.minor_indices.len());
-        }
-        self.buf
-    }
     /// Consumes self and outputs the constructed `SparsityPattern`.
     /// If elements were added to the last major, but `advance_major`
     /// was not called, will implicitly call `advance_major` then
@@ -455,29 +434,24 @@ impl SparsityPatternBuilder {
         assert_eq!(self.buf.major_dim(), self.major_dim);
         self.buf
     }
+    /// Returns the current major being modified by `self`.
+    pub fn current_major(&self) -> usize {
+        assert!(!self.buf.major_offsets.is_empty());
+        self.buf.major_offsets.len() - 1
+    }
 
     /// Reverts the major index of `self` back to `maj`, deleting any entries ahead of it.
     /// Preserves entries in `maj`.
     pub fn revert_to_major(&mut self, maj: usize) -> bool {
-        if maj >= self.buf.major_dim() {
+        // preserve maj + 1 elements in self
+        if self.buf.major_offsets.len() + 1 <= maj {
             return false;
         }
         let last = self.buf.major_offsets[maj + 1];
         self.buf.major_offsets.truncate(maj + 1);
-        self.buf.minor_indices.truncate(last);
+        self.buf.minor_indices.truncate(last + 1);
         true
     }
-    /// Backtracks along the minor dimension, allowing for adding different elements along the
-    /// current major axis.
-    pub fn revert_current_minor(&mut self, min: usize) -> bool {
-        if min >= self.current_minor() {
-            return false;
-        }
-        while min < self.current_minor() {
-            self.buf.minor_indices.pop();
-        }
-        return true;
-    }
 
     /// Allows for rebuilding part of a sparsity pattern, assuming that
     /// items after maj_start have not been filled in.

nalgebra-sparse/tests/lu_factorization.rs

@@ -0,0 +1,46 @@
+use nalgebra_sparse::factorization::LeftLookingLUFactorization;
+use nalgebra_sparse::CscBuilder;
+
+#[test]
+fn test_basic_lu_factorization() {
+    let n = 5;
+    let mut a = CscBuilder::new(n, n);
+    for i in 0..n {
+        assert!(a.insert(i, i, 1.).is_ok());
+    }
+    // construct an identity matrix as a basic test
+    let a = a.build();
+
+    let lu_fact = LeftLookingLUFactorization::new(&a);
+
+    assert_eq!(lu_fact.u(), a);
+}
+
+#[test]
+fn test_basic_lu_factorization_with_one_more_entry() {
+    let n = 3;
+    let mut a = CscBuilder::new(n, n);
+    for i in 0..n {
+        assert!(a.insert(i, i, if i == 0 { 0.5 } else { 1. }).is_ok());
+        if i == 0 {
+            assert!(a.insert(1, 0, 1.).is_ok());
+        }
+    }
+    // construct an identity matrix as a basic test
+    let a = a.build();
+
+    let lu_fact = LeftLookingLUFactorization::new(&a);
+
+    let mut ground_truth = CscBuilder::new(n, n);
+    for i in 0..n {
+        assert!(ground_truth
+            .insert(i, i, if i == 0 { 0.5 } else { 1. })
+            .is_ok());
+        if i == 0 {
+            assert!(ground_truth.insert(1, 0, 2.).is_ok());
+        }
+    }
+    let gt = ground_truth.build();
+
+    assert_eq!(lu_fact.lu(), &gt);
+}

nalgebra-sparse/tests/sparsity.rs

@@ -77,16 +77,12 @@ fn test_builder_reset() {
     assert_eq!(out.minor_dim(), 4);
     let mut builder = SparsityPatternBuilder::from(out);
     for i in (0..=2).rev() {
-        builder.revert_to_major(i);
-        assert_eq!(builder.major_dim(), i);
-        assert_eq!(builder.current_minor(), i + 1);
+        assert!(builder.revert_to_major(i));
+        assert_eq!(builder.current_major(), i);
     }
     let out = builder.build();
-    assert_eq!(out.major_dim(), 4);
-    assert_eq!(out.minor_dim(), 4);
 
     let mut builder = SparsityPatternBuilder::from(out);
-    builder.revert_to_major(1);
-    assert_eq!(builder.major_dim(), 1);
-    assert_eq!(builder.current_minor(), 0);
+    assert!(builder.revert_to_major(1));
+    assert_eq!(builder.current_major(), 1);
 }