From c3e8112d5ea13f1e135272e65e38a0c1570e0ddd Mon Sep 17 00:00:00 2001 From: sebcrozet Date: Sun, 4 Nov 2018 07:10:43 +0100 Subject: [PATCH] Add implementation of the left-looking cholesky decomposition. --- src/sparse/cs_matrix.rs | 39 +++++++++++++---- src/sparse/cs_matrix_cholesky.rs | 72 +++++++++++++++++++++++++++++++- src/sparse/mod.rs | 2 +- tests/sparse/cs_cholesky.rs | 2 +- 4 files changed, 102 insertions(+), 13 deletions(-) diff --git a/src/sparse/cs_matrix.rs b/src/sparse/cs_matrix.rs index ec4cc5b0..6403f572 100644 --- a/src/sparse/cs_matrix.rs +++ b/src/sparse/cs_matrix.rs @@ -20,6 +20,12 @@ pub trait CsStorageIter<'a, N, R, C = U1> { fn column_entries(&'a self, j: usize) -> Self::ColumnEntries; } +pub trait CsStorageIterMut<'a, N: 'a, R, C = U1> { + type ColumnEntriesMut: Iterator; + + fn column_entries_mut(&'a mut self, j: usize) -> Self::ColumnEntriesMut; +} + pub trait CsStorage: for<'a> CsStorageIter<'a, N, R, C> { fn shape(&self) -> (R, C); unsafe fn row_index_unchecked(&self, i: usize) -> usize; @@ -30,15 +36,9 @@ pub trait CsStorage: for<'a> CsStorageIter<'a, N, R, C> { fn len(&self) -> usize; } -pub trait CsStorageMut: CsStorage { - /* - /// Sets the length of this column without initializing its values and row indices. - /// - /// If the given length is larger than the current one, uninitialized entries are - /// added at the end of the column `i`. This will effectively shift all the matrix entries - /// of the columns at indices `j` with `j > i`. - fn set_column_len(&mut self, i: usize, len: usize); - */ +pub trait CsStorageMut: + CsStorage + for<'a> CsStorageIterMut<'a, N, R, C> +{ } #[derive(Clone, Debug)] @@ -133,6 +133,27 @@ where } } +impl<'a, N: Scalar, R: Dim, C: Dim> CsStorageIterMut<'a, N, R, C> for CsVecStorage +where + DefaultAllocator: Allocator, +{ + type ColumnEntriesMut = iter::Zip>, slice::IterMut<'a, N>>; + + #[inline] + fn column_entries_mut(&'a mut self, j: usize) -> Self::ColumnEntriesMut { + let rng = self.column_range(j); + self.i[rng.clone()] + .iter() + .cloned() + .zip(self.vals[rng].iter_mut()) + } +} + +impl CsStorageMut for CsVecStorage where + DefaultAllocator: Allocator +{ +} + /* pub struct CsSliceStorage<'a, N: Scalar, R: Dim, C: DimAdd> { shape: (R, C), diff --git a/src/sparse/cs_matrix_cholesky.rs b/src/sparse/cs_matrix_cholesky.rs index 3b4185a8..01a6b081 100644 --- a/src/sparse/cs_matrix_cholesky.rs +++ b/src/sparse/cs_matrix_cholesky.rs @@ -8,7 +8,7 @@ use std::slice; use allocator::Allocator; use constraint::{AreMultipliable, DimEq, SameNumberOfRows, ShapeConstraint}; -use sparse::{CsMatrix, CsStorage, CsStorageIter, CsVecStorage, CsVector}; +use sparse::{CsMatrix, CsStorage, CsStorageIter, CsStorageIterMut, CsVecStorage, CsVector}; use storage::{Storage, StorageMut}; use {DefaultAllocator, Dim, Matrix, MatrixMN, Real, Scalar, Vector, VectorN, U1}; @@ -39,7 +39,7 @@ where /// Computes the cholesky decomposition of the sparse matrix `m`. pub fn new(m: &CsMatrix) -> Self { let mut me = Self::new_symbolic(m); - let _ = me.decompose(&m.data.vals); + let _ = me.decompose_left_looking(&m.data.vals); me } /// Perform symbolic analysis for the given matrix. @@ -86,6 +86,74 @@ where } } + pub fn decompose_left_looking(&mut self, values: &[N]) -> bool { + assert!( + values.len() >= self.original_i.len(), + "The set of values is too small." + ); + + let n = self.l.nrows(); + + // Reset `work_c` to the column pointers of `l`. + self.work_c.copy_from(&self.l.data.p); + + unsafe { + for k in 0..n { + // Scatter the k-th column of the original matrix with the values provided. + let range_k = + *self.original_p.get_unchecked(k)..*self.original_p.get_unchecked(k + 1); + + *self.work_x.vget_unchecked_mut(k) = N::zero(); + for p in range_k.clone() { + let irow = *self.original_i.get_unchecked(p); + + if irow >= k { + *self.work_x.vget_unchecked_mut(irow) = *values.get_unchecked(p); + } + } + + for j in self.u.data.column_row_indices(k) { + let factor = -*self + .l + .data + .vals + .get_unchecked(*self.work_c.vget_unchecked(j)); + *self.work_c.vget_unchecked_mut(j) += 1; + + if j < k { + for (z, val) in self.l.data.column_entries(j) { + if z >= k { + *self.work_x.vget_unchecked_mut(z) += val * factor; + } + } + } + } + + let diag = *self.work_x.vget_unchecked(k); + + if diag > N::zero() { + let denom = diag.sqrt(); + *self + .l + .data + .vals + .get_unchecked_mut(*self.l.data.p.vget_unchecked(k)) = denom; + + for (p, val) in self.l.data.column_entries_mut(k) { + *val = *self.work_x.vget_unchecked(p) / denom; + *self.work_x.vget_unchecked_mut(p) = N::zero(); + } + } else { + self.ok = false; + return false; + } + } + } + + self.ok = true; + true + } + // Performs the numerical Cholesky decomposition given the set of numerical values. pub fn decompose(&mut self, values: &[N]) -> bool { assert!( diff --git a/src/sparse/mod.rs b/src/sparse/mod.rs index 6ce898e5..411e133b 100644 --- a/src/sparse/mod.rs +++ b/src/sparse/mod.rs @@ -1,5 +1,5 @@ pub use self::cs_matrix::{ - CsMatrix, CsStorage, CsStorageIter, CsStorageMut, CsVecStorage, CsVector, + CsMatrix, CsStorage, CsStorageIter, CsStorageIterMut, CsStorageMut, CsVecStorage, CsVector, }; pub use self::cs_matrix_cholesky::CsCholesky; diff --git a/tests/sparse/cs_cholesky.rs b/tests/sparse/cs_cholesky.rs index 9c199737..aebefacb 100644 --- a/tests/sparse/cs_cholesky.rs +++ b/tests/sparse/cs_cholesky.rs @@ -51,5 +51,5 @@ fn test_cholesky(a: Matrix5) { println!("{}", l); println!("{}", cs_l); - assert_eq!(l, cs_l); + assert_relative_eq!(l, cs_l); }