From cd9c3baeada0e7442041a3c04e53b741c6d6516b Mon Sep 17 00:00:00 2001 From: Andreas Longva Date: Tue, 19 Jan 2021 15:20:01 +0100 Subject: [PATCH] Clean up CscCholesky --- nalgebra-sparse/src/factorization/cholesky.rs | 103 ++++++------------ 1 file changed, 36 insertions(+), 67 deletions(-) diff --git a/nalgebra-sparse/src/factorization/cholesky.rs b/nalgebra-sparse/src/factorization/cholesky.rs index 3e9c032e..9dfd20ff 100644 --- a/nalgebra-sparse/src/factorization/cholesky.rs +++ b/nalgebra-sparse/src/factorization/cholesky.rs @@ -4,10 +4,9 @@ use crate::pattern::SparsityPattern; use crate::csc::CscMatrix; use core::{mem, iter}; -use nalgebra::{U1, VectorN, Dynamic, Scalar, RealField}; -use num_traits::Zero; +use nalgebra::{Scalar, RealField}; use std::sync::Arc; -use std::ops::Add; +use std::fmt::{Display, Formatter}; pub struct CscSymbolicCholesky { // Pattern of the original matrix that was decomposed @@ -21,37 +20,11 @@ impl CscSymbolicCholesky { pub fn factor(pattern: &Arc) -> Self { assert_eq!(pattern.major_dim(), pattern.minor_dim(), "Major and minor dimensions must be the same (square matrix)."); - - // TODO: Temporary stopgap solution to make things work until we can refactor - #[derive(Copy, Clone, PartialEq, Eq, Debug)] - struct DummyVal; - impl Zero for DummyVal { - fn zero() -> Self { - DummyVal - } - - fn is_zero(&self) -> bool { - true - } - } - - impl Add for DummyVal { - type Output = Self; - - fn add(self, rhs: DummyVal) -> Self::Output { - rhs - } - } - - let dummy_vals = vec![DummyVal; pattern.nnz()]; - let dummy_csc = CscMatrix::try_from_pattern_and_values(Arc::clone(pattern), dummy_vals) - .unwrap(); - let (l, u) = nonzero_pattern(&dummy_csc); - // TODO: Don't clone unnecessarily + let (l_pattern, u_pattern) = nonzero_pattern(&*pattern); Self { m_pattern: Arc::clone(pattern), - l_pattern: l.pattern().as_ref().clone(), - u_pattern: u.pattern().as_ref().clone() + l_pattern, + u_pattern, } } @@ -70,10 +43,20 @@ pub struct CscCholesky { } #[derive(Debug, PartialEq, Eq, Clone)] +#[non_exhaustive] pub enum CholeskyError { - + /// The matrix is not positive definite. + NotPositiveDefinite, } +impl Display for CholeskyError { + fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { + write!(f, "Matrix is not positive definite") + } +} + +impl std::error::Error for CholeskyError {} + impl CscCholesky { pub fn factor(matrix: &CscMatrix) -> Result { @@ -181,9 +164,7 @@ impl CscCholesky { *self.work_x.get_unchecked_mut(p) = T::zero(); } } else { - // self.ok = false; - // TODO: Return indefinite error (i.e. encountered non-positive diagonal - unimplemented!() + return Err(CholeskyError::NotPositiveDefinite); } } } @@ -196,8 +177,8 @@ impl CscCholesky { -fn reach( - m: &CscMatrix, +fn reach( + pattern: &SparsityPattern, j: usize, max_j: usize, tree: &[usize], @@ -211,7 +192,7 @@ fn reach( let mut tmp = Vec::new(); let mut res = Vec::new(); - for &irow in m.col(j).row_indices() { + for &irow in pattern.lane(j) { let mut curr = irow; while curr != usize::max_value() && curr <= max_j && !marks[curr] { marks[curr] = true; @@ -223,57 +204,45 @@ fn reach( mem::swap(&mut tmp, &mut res); } - // TODO: Is this right? res.sort_unstable(); out.append(&mut res); } -fn nonzero_pattern( - m: &CscMatrix -) -> (CscMatrix, CscMatrix) { - // TODO: In order to stay as faithful as possible to the original implementation, - // we here return full matrices, whereas we actually only need to construct sparsity patterns - +fn nonzero_pattern(m: &SparsityPattern) -> (SparsityPattern, SparsityPattern) { let etree = elimination_tree(m); - let (nrows, ncols) = (m.nrows(), m.ncols()); + // Note: We assume CSC, therefore rows == minor and cols == major + let (nrows, ncols) = (m.minor_dim(), m.major_dim()); let mut rows = Vec::with_capacity(m.nnz()); - // TODO: Use a Vec here instead - let mut cols = unsafe { VectorN::new_uninitialized_generic(Dynamic::new(nrows), U1) }; + let mut col_offsets = Vec::with_capacity(ncols + 1); let mut marks = Vec::new(); // NOTE: the following will actually compute the non-zero pattern of // the transpose of l. + col_offsets.push(0); for i in 0..nrows { - cols[i] = rows.len(); reach(m, i, i, &etree, &mut marks, &mut rows); + col_offsets.push(rows.len()); } - // TODO: Get rid of this in particular - let mut vals = Vec::with_capacity(rows.len()); - unsafe { - vals.set_len(rows.len()); - } - vals.shrink_to_fit(); + let u_pattern = SparsityPattern::try_from_offsets_and_indices(nrows, ncols, col_offsets, rows) + .unwrap(); - // TODO: Remove this unnecessary conversion by using Vec throughout - let mut cols: Vec<_> = cols.iter().cloned().collect(); - cols.push(rows.len()); + // TODO: Avoid this transpose? + let l_pattern = u_pattern.transpose(); - let u = CscMatrix::try_from_csc_data(nrows, ncols, cols, rows, vals).unwrap(); - // TODO: Avoid this transpose - let l = u.transpose(); - - (l, u) + (l_pattern, u_pattern) } -fn elimination_tree(m: &CscMatrix) -> Vec { - let nrows = m.nrows(); +fn elimination_tree(pattern: &SparsityPattern) -> Vec { + // Note: The pattern is assumed to of a CSC matrix, so the number of rows is + // given by the minor dimension + let nrows = pattern.minor_dim(); let mut forest: Vec<_> = iter::repeat(usize::max_value()).take(nrows).collect(); let mut ancestor: Vec<_> = iter::repeat(usize::max_value()).take(nrows).collect(); for k in 0..nrows { - for &irow in m.col(k).row_indices() { + for &irow in pattern.lane(k) { let mut i = irow; while i < k {