Implement SparseLU factorization

Add `solve_upper_triangular` to `CsrMatrix`

This allows a sparse matrix to be used for efficient solving with a dense LU decomposition.

nalgebra-sparse/src/csc.rs

@@ -12,7 +12,7 @@ use crate::csr::CsrMatrix;
 use crate::pattern::{SparsityPattern, SparsityPatternFormatError, SparsityPatternIter};
 use crate::{SparseEntry, SparseEntryMut, SparseFormatError, SparseFormatErrorKind};
 
-use nalgebra::Scalar;
+use nalgebra::{RealField, Scalar};
 use num_traits::One;
 use std::slice::{Iter, IterMut};
 
@@ -553,6 +553,87 @@ impl<T> CscMatrix<T> {
         self.filter(|i, j, _| i >= j)
     }
 
+    /// Solves a lower triangular system, `self` is a matrix of NxN, and `b` is a column vector of size N
+    /// Assuming that b is dense.
+    // TODO add an option here for assuming diagonal is one.
+    pub fn dense_lower_triangular_solve(&self, b: &[T], out: &mut [T], unit_diagonal: bool)
+    where
+        T: RealField + Copy,
+    {
+        assert_eq!(self.nrows(), self.ncols());
+        assert_eq!(self.ncols(), b.len());
+        assert_eq!(out.len(), b.len());
+        out.copy_from_slice(b);
+        let n = b.len();
+
+        for i in 0..n {
+            let mul = out[i];
+            let col = self.col(i);
+            for (&ri, &v) in col.row_indices().iter().zip(col.values().iter()) {
+                // ensure that only using the lower part
+                if ri == i {
+                    if !unit_diagonal {
+                        out[ri] /= v;
+                    }
+                } else if ri > i {
+                    out[ri] -= v * mul;
+                }
+            }
+        }
+    }
+
+    /// Solves a sparse lower triangular system `Ax = b`, with both the matrix and vector
+    /// sparse.
+    /// sparsity_idxs should be precomputed using the sparse_lower_triangle.
+    /// Assumes that the diagonal of the sparse matrix is all 1.
+    pub fn sparse_lower_triangular_solve(
+        &self,
+        b_idxs: &[usize],
+        b: &[T],
+        // idx -> row
+        // for now, is permitted to be unsorted
+        // TODO maybe would be better to enforce sorted, but would have to sort internally.
+        out_sparsity_pattern: &[usize],
+        out: &mut [T],
+        assume_unit: bool,
+    ) where
+        T: RealField + Copy,
+    {
+        assert_eq!(self.nrows(), self.ncols());
+        assert_eq!(b.len(), b_idxs.len());
+        assert!(b_idxs.iter().all(|&bi| bi < self.ncols()));
+
+        assert_eq!(out_sparsity_pattern.len(), out.len());
+        assert!(out_sparsity_pattern.iter().all(|&i| i < self.ncols()));
+
+        // initialize out with b
+        for (&bv, &bi) in b.iter().zip(b_idxs.iter()) {
+            let out_pos = out_sparsity_pattern.iter().position(|&p| p == bi).unwrap();
+            out[out_pos] = bv;
+        }
+
+        for (i, &row) in out_sparsity_pattern.iter().enumerate() {
+            let mul = out[i];
+            let col = self.col(row);
+            for (ni, &nrow) in out_sparsity_pattern.iter().enumerate() {
+                if nrow < row {
+                    continue;
+                }
+                // TODO in a sorted version may be able to iterate without
+                // having the cost of binary search at each iteration
+                let l_val = col
+                    .get_entry(nrow)
+                    .map(|v| v.into_value())
+                    .unwrap_or_else(T::zero);
+                if !assume_unit && (ni == i || nrow == row) {
+                    out[i] /= l_val;
+                } else if nrow > row {
+                    out[ni] -= l_val * mul;
+                }
+            }
+        }
+    }
+
     /// Returns the diagonal of the matrix as a sparse matrix.
     #[must_use]
     pub fn diagonal_as_csc(&self) -> Self

nalgebra-sparse/src/csr.rs

@@ -12,7 +12,7 @@ use crate::csc::CscMatrix;
 use crate::pattern::{SparsityPattern, SparsityPatternFormatError, SparsityPatternIter};
 use crate::{SparseEntry, SparseEntryMut, SparseFormatError, SparseFormatErrorKind};
 
-use nalgebra::Scalar;
+use nalgebra::{DMatrix, DMatrixView, RealField, Scalar};
 use num_traits::One;
 
 use std::slice::{Iter, IterMut};
@@ -573,6 +573,55 @@ impl<T> CsrMatrix<T> {
     {
         CscMatrix::from(self).transpose_as_csr()
     }
+
+    /// Solves the equation `Ax = b`, treating `self` as an upper triangular matrix.
+    /// If `A` is not upper triangular, elements in the lower triangle below the diagonal
+    /// will be ignored.
+    ///
+    /// If `m` has zeros along the diagonal, returns `None`.
+    /// Panics:
+    /// Panics if `A` and `b` have incompatible shapes, specifically if `b`
+    /// has a different number of rows and than `a`.
+    pub fn solve_upper_triangular<'a>(&self, b: impl Into<DMatrixView<'a, T>>) -> Option<DMatrix<T>>
+    where
+        T: RealField + Scalar,
+    {
+        // https://www.nicolasboumal.net/papers/MAT321_Lecture_notes_Boumal_2019.pdf
+        // page 48
+        let b: DMatrixView<'a, T> = b.into();
+        assert_eq!(b.nrows(), self.nrows());
+
+        let out_cols = b.ncols();
+        let out_rows = self.nrows();
+
+        let mut out = DMatrix::zeros(out_rows, out_cols);
+        for r in (0..out_rows).rev() {
+            let row = self.row(r);
+            // only take upper triangle elements
+            let mut row_iter = row
+                .col_indices()
+                .iter()
+                .copied()
+                .zip(row.values().iter())
+                .filter(|&(c, _)| c >= r);
+
+            let (c, div) = row_iter.next()?;
+            // This implies there is a 0 on the diagonal
+            if c != r || div.is_zero() {
+                return None;
+            }
+            for c in 0..out_cols {
+                let numer = b.index((r, c)).clone();
+                let numer = numer
+                    - row_iter
+                        .clone()
+                        .map(|(a_col, val)| val.clone() * b.index((a_col, c)).clone())
+                        .fold(T::zero(), |acc, n| acc + n);
+                *out.index_mut((r, c)) = numer / div.clone();
+            }
+        }
+        Some(out)
+    }
 }
 
 /// Convert pattern format errors into more meaningful CSR-specific errors.

nalgebra-sparse/src/factorization/lu.rs

@@ -0,0 +1,16 @@
+use crate::CscMatrix;
+use nalgebra::RealField;
+
+pub struct LeftLookingLUFactorization<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!();
+  }
+}

nalgebra-sparse/src/factorization/mod.rs

@@ -4,3 +4,7 @@
 mod cholesky;
 
 pub use cholesky::*;
+
+//mod lu;
+
+//pub use lu::*;

nalgebra-sparse/src/lib.rs

@@ -279,4 +279,11 @@ impl<'a, T: Clone + Zero> SparseEntryMut<'a, T> {
             SparseEntryMut::Zero => T::zero(),
         }
     }
+    /// If the entry is nonzero, returns `Some(&mut value)`, otherwise returns `None`.
+    pub fn nonzero(self) -> Option<&'a mut T> {
+        match self {
+            SparseEntryMut::NonZero(v) => Some(v),
+            SparseEntryMut::Zero => None,
+        }
+    }
 }

nalgebra-sparse/src/pattern.rs

@@ -61,6 +61,14 @@ impl SparsityPattern {
             minor_dim,
         }
     }
+    /// Creates the sparsity pattern of an identity matrix of size `n`.
+    pub fn identity(n: usize) -> Self {
+        Self {
+            major_offsets: (0..=n).collect(),
+            minor_indices: (0..n).collect(),
+            minor_dim: n,
+        }
+    }
 
     /// The offsets for the major dimension.
     #[inline]
@@ -109,6 +117,13 @@ impl SparsityPattern {
     pub fn lane(&self, major_index: usize) -> &[usize] {
         self.get_lane(major_index).unwrap()
     }
+    /// Returns an iterator over all lanes in this sparsity pattern, in order.
+    /// Does not omit empty lanes.
+    #[inline]
+    #[must_use]
+    pub fn lanes(&self) -> impl Iterator<Item = &[usize]> {
+        (0..self.major_offsets.len() - 1).map(move |r| self.lane(r))
+    }
 
     /// Get the lane at the given index, or `None` if out of bounds.
     #[inline]
@@ -289,6 +304,133 @@ impl SparsityPattern {
         )
         .expect("Internal error: Transpose should never fail.")
     }
+
+    /// Computes the output sparsity pattern of `x` in `Ax = b`.
+    /// where A's nonzero pattern is given by `self` and the non-zero indices
+    /// of vector `b` are specified as a slice.
+    /// 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> {
+        assert!(b.iter().all(|&i| i < self.major_dim()));
+        let mut out = vec![];
+
+        // From a given starting column, traverses and finds all reachable indices.
+        fn reach(sp: &SparsityPattern, j: usize, out: &mut Vec<usize>) {
+            // already traversed
+            if out.contains(&j) {
+                return;
+            }
+
+            out.push(j);
+            for &i in sp.lane(j) {
+                reach(sp, i, out);
+            }
+        }
+
+        for &i in b {
+            reach(&self, i, &mut out);
+        }
+
+        out
+    }
+
+    /// Left-looking Sparse LU decomposition from Gilbert/Peierls.
+    /// returns the sparsity pattern of the output
+    pub fn left_looking_lu_decomposition(&self) -> SparsityPattern {
+        assert_eq!(self.major_dim(), self.minor_dim());
+        let n = self.minor_dim();
+        let mut sp = SparsityPatternBuilder::new(n, n);
+        for col in 0..n {
+            let mut x = sp
+                .valid_partial()
+                .sparse_lower_triangular_solve(self.lane(col));
+            x.sort_unstable();
+            for row in x {
+                assert!(sp.insert(col, row).is_ok());
+            }
+        }
+        sp.build()
+    }
+}
+
+/// A builder that allows for constructing a sparsity pattern.
+/// It requires elements to be added in sorted order. Specifically,
+/// For each element the major must be >= the previous element's major.
+/// If the major is the same, the minor must be in ascending order.
+#[derive(Debug, Clone, PartialEq)]
+pub struct SparsityPatternBuilder {
+    buf: SparsityPattern,
+    major_dim: usize,
+}
+
+///
+#[derive(Debug, Copy, Clone, PartialEq, Eq)]
+pub enum BuilderInsertError {
+    ///
+    MajorTooLow,
+    ///
+    MinorTooLow,
+}
+
+impl SparsityPatternBuilder {
+    /// Constructs a new empty builder.
+    pub fn new(major_dim: usize, minor_dim: usize) -> Self {
+        Self {
+            buf: SparsityPattern {
+                major_offsets: vec![0],
+                minor_indices: vec![],
+                minor_dim,
+            },
+            major_dim,
+        }
+    }
+    /// Allows for general assignment of indices
+    pub fn insert(&mut self, maj: usize, min: usize) -> Result<(), BuilderInsertError> {
+        assert!(maj < self.major_dim);
+        assert!(min < self.buf.minor_dim);
+
+        let curr_major = self.buf.major_dim();
+
+        // cannot go backwards in major
+        if maj < curr_major {
+            return Err(BuilderInsertError::MajorTooLow);
+        }
+        // cannot go backwards in minor
+        if maj == curr_major
+            && !self.buf.minor_indices.is_empty()
+            && min <= *self.buf.minor_indices.last().unwrap()
+        {
+            return Err(BuilderInsertError::MinorTooLow);
+        }
+        // add any advances in row.
+        for _ in curr_major..maj {
+            self.buf.major_offsets.push(self.buf.minor_indices.len());
+        }
+        self.buf.minor_indices.push(min);
+        Ok(())
+    }
+    /// Returns a valid partial sparsity pattern.
+    /// All the major lanes up to the current insertion will be completed.
+    pub(crate) fn valid_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());
+        }
+        assert_eq!(self.buf.major_dim(), self.major_dim);
+        self.buf
+    }
 }
 
 /// Error type for `SparsityPattern` format errors.