Start working on actual LU factorization

nalgebra-sparse/src/cs.rs

@@ -700,15 +700,48 @@ pub struct CsBuilder<T> {
 }
 
 impl<T> CsBuilder<T> {
+    /// Constructs a new CsBuilder of the given size.
     pub fn new(major_dim: usize, minor_dim: usize) -> Self {
         Self {
             sparsity_builder: SparsityPatternBuilder::new(major_dim, minor_dim),
             values: vec![],
         }
     }
+    /// Given an existing CsMatrix, allows for modification by converting it into a builder.
+    pub fn from_mat(mat: CsMatrix<T>) -> Self {
+        let CsMatrix {
+            sparsity_pattern,
+            values,
+        } = mat;
+
+        CsBuilder {
+            sparsity_builder: SparsityPatternBuilder::from(sparsity_pattern),
+            values,
+        }
+    }
+    /// Backtracks to a given major index
+    pub fn revert_to_major(&mut self, maj: usize) -> bool {
+        if !self.sparsity_builder.revert_to_major(maj) {
+            return false;
+        }
+
+        self.values.truncate(self.sparsity_builder.num_entries());
+        true
+    }
     pub fn insert(&mut self, maj: usize, min: usize, val: T) -> Result<(), BuilderInsertError> {
         self.sparsity_builder.insert(maj, min)?;
         self.values.push(val);
         Ok(())
     }
+    pub fn build(self) -> CsMatrix<T> {
+        let CsBuilder {
+            sparsity_builder,
+            values,
+        } = self;
+        let sparsity_pattern = sparsity_builder.build();
+        CsMatrix {
+            sparsity_pattern,
+            values,
+        }
+    }
 }

nalgebra-sparse/src/csc.rs

@@ -937,8 +937,20 @@ impl<T> CscBuilder<T> {
     pub fn new(rows: usize, cols: usize) -> Self {
         Self(CsBuilder::new(cols, rows))
     }
+    /// Convert back from a matrix to a CscBuilder.
+    pub fn from_mat(mat: CscMatrix<T>) -> Self {
+        Self(CsBuilder::from_mat(mat.cs))
+    }
+    /// Backtracks back to column `col`, deleting all entries ahead of it.
+    pub fn revert_to_col(&mut self, col: usize) -> bool {
+        self.0.revert_to_major(col)
+    }
     /// Inserts a value into the builder. Must be called in ascending col, row order.
     pub fn insert(&mut self, row: usize, col: usize, val: T) -> Result<(), BuilderInsertError> {
         self.0.insert(col, row, val)
     }
+    /// Converts this builder into a valid CscMatrix.
+    pub fn build(self) -> CscMatrix<T> {
+        CscMatrix { cs: self.0.build() }
+    }
 }

nalgebra-sparse/src/factorization/lu.rs

@@ -1,16 +1,67 @@
-use crate::CscMatrix;
+use crate::{CscBuilder, CscMatrix};
 use nalgebra::RealField;
 
+/// Constructs an LU Factorization using a left-looking approach.
+/// This means it will construct each column, starting from the leftmost one.
 pub struct LeftLookingLUFactorization<T> {
-  /// A single matrix stores both the lower and upper triangular components
-  l_u: CscMatrix<T>
+    /// A single matrix stores both the lower and upper triangular components
+    l_u: CscMatrix<T>,
 }
 
-impl<T: RealField> LeftLookingLUFactorization<T> {
-  /// Construct a new sparse LU factorization
-  /// from a given CSC matrix.
-  pub fn new(a: &CscMatrix<T>) -> Self {
-    assert_eq!(a.nrows(), a.ncols());
-    todo!();
-  }
+impl<T: RealField + Copy> LeftLookingLUFactorization<T> {
+    /// Returns the upper triangular part of this matrix.
+    pub fn u(&self) -> CscMatrix<T> {
+        self.l_u.upper_triangle()
+    }
+
+    /*
+    /// 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!();
+    }
+    */
+
+    /// Construct a new sparse LU factorization
+    /// from a given CSC matrix.
+    pub fn new(a: &CscMatrix<T>) -> Self {
+        assert_eq!(a.nrows(), a.ncols());
+        let n = a.nrows();
+
+        // this initially starts as an identity  matrix.
+        // but the ones are all implicit.
+        let mut csc_builder = CscBuilder::new(n, n);
+
+        let mut val_buf = vec![];
+        let mut pat_buf = vec![];
+
+        for (i, col) in a.col_iter().enumerate() {
+            let curr_mat = csc_builder.build();
+
+            curr_mat
+                .pattern()
+                .sparse_lower_triangular_solve(col.row_indices(), &mut pat_buf);
+            pat_buf.sort_unstable();
+            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(
+                col.row_indices(),
+                col.values(),
+                &pat_buf,
+                &mut val_buf,
+                true,
+            );
+
+            // convert builder back to matrix
+            csc_builder = CscBuilder::from_mat(curr_mat);
+            assert!(csc_builder.revert_to_col(i));
+            let _ = x;
+            todo!();
+        }
+
+        //csc_builder.build()
+        todo!();
+    }
 }

nalgebra-sparse/src/factorization/mod.rs

@@ -2,9 +2,7 @@
 //!
 //! Currently, the only factorization provided here is the [`CscCholesky`] factorization.
 mod cholesky;
+mod lu;
 
 pub use cholesky::*;
-
-//mod lu;
-
-//pub use lu::*;
+pub use lu::*;

nalgebra-sparse/src/pattern.rs

@@ -311,9 +311,9 @@ impl SparsityPattern {
     /// The output is not necessarily in sorted order, but is topological sort order.
     /// Treats `self` as lower triangular, even if there are elements in the upper triangle.
     /// Acts as if b is one major lane (i.e. CSC matrix and one column)
-    pub fn sparse_lower_triangular_solve(&self, b: &[usize]) -> Vec<usize> {
+    pub fn sparse_lower_triangular_solve(&self, b: &[usize], out: &mut Vec<usize>) {
         assert!(b.iter().all(|&i| i < self.major_dim()));
-        let mut out = vec![];
+        out.clear();
 
         // From a given starting column, traverses and finds all reachable indices.
         fn reach(sp: &SparsityPattern, j: usize, out: &mut Vec<usize>) {
@@ -329,10 +329,8 @@ impl SparsityPattern {
         }
 
         for &i in b {
-            reach(&self, i, &mut out);
+            reach(&self, i, out);
         }
-
-        out
     }
 
     /// Left-looking Sparse LU decomposition from Gilbert/Peierls.
@@ -341,12 +339,12 @@ impl SparsityPattern {
         assert_eq!(self.major_dim(), self.minor_dim());
         let n = self.minor_dim();
         let mut sp = SparsityPatternBuilder::new(n, n);
+        let mut x = vec![];
         for col in 0..n {
-            let mut x = sp
-                .valid_partial()
-                .sparse_lower_triangular_solve(self.lane(col));
+            sp.valid_partial()
+                .sparse_lower_triangular_solve(self.lane(col), &mut x);
             x.sort_unstable();
-            for row in x {
+            for &row in &x {
                 assert!(sp.insert(col, row).is_ok());
             }
         }
@@ -385,6 +383,25 @@ impl SparsityPatternBuilder {
             major_dim,
         }
     }
+    /// The number of entries into this sparsity pattern
+    pub fn num_entries(&self) -> usize {
+        self.buf.minor_indices.last().copied().unwrap_or(0)
+    }
+
+    /// 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);
@@ -418,19 +435,58 @@ 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
     /// output the values.
     #[inline]
     pub fn build(mut self) -> SparsityPattern {
-        let curr_major = self.buf.major_dim();
-        for _ in curr_major..self.major_dim {
-            self.buf.major_offsets.push(self.buf.minor_indices.len());
-        }
+        self.buf
+            .major_offsets
+            .resize(self.major_dim + 1, self.buf.minor_indices.len());
         assert_eq!(self.buf.major_dim(), self.major_dim);
         self.buf
     }
+
+    /// 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() {
+            return false;
+        }
+        let last = self.buf.major_offsets[maj + 1];
+        self.buf.major_offsets.truncate(maj + 1);
+        self.buf.minor_indices.truncate(last);
+        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.
+    pub fn from(sp: SparsityPattern) -> Self {
+        SparsityPatternBuilder {
+            major_dim: sp.major_dim(),
+            buf: sp,
+        }
+    }
 }
 
 /// Error type for `SparsityPattern` format errors.

nalgebra-sparse/tests/sparsity.rs

@@ -16,7 +16,9 @@ fn lower_sparse_solve() {
     // just a smaller number so it's easier to debug
     let n = 8;
     let speye = SparsityPattern::identity(n);
-    assert_eq!(speye.sparse_lower_triangular_solve(&[0, 5]), vec![0, 5]);
+    let mut buf = vec![];
+    speye.sparse_lower_triangular_solve(&[0, 5], &mut buf);
+    assert_eq!(buf, vec![0, 5]);
 
     // test case from
     // https://www.youtube.com/watch?v=1dGRTOwBkQs&ab_channel=TimDavis
@@ -49,10 +51,8 @@ fn lower_sparse_solve() {
     for ij in sp.entries() {
         assert!(indices.contains(&ij));
     }
-    assert_eq!(
-        sp.sparse_lower_triangular_solve(&[3, 5]),
-        vec![3, 8, 11, 12, 13, 5, 9, 10]
-    );
+    sp.sparse_lower_triangular_solve(&[3, 5], &mut buf);
+    assert_eq!(buf, vec![3, 8, 11, 12, 13, 5, 9, 10]);
 }
 
 #[test]
@@ -65,3 +65,28 @@ fn test_builder() {
     assert!(builder.insert(1, 0).is_ok());
     assert!(builder.insert(1, 0).is_err());
 }
+
+#[test]
+fn test_builder_reset() {
+    let mut builder = SparsityPatternBuilder::new(4, 4);
+    for i in 0..3 {
+        assert!(builder.insert(i, i + 1).is_ok());
+    }
+    let out = builder.build();
+    assert_eq!(out.major_dim(), 4);
+    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);
+    }
+    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);
+}