Add implementation of the left-looking cholesky decomposition.

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<Item = (usize, &'a mut N)>;
+
+    fn column_entries_mut(&'a mut self, j: usize) -> Self::ColumnEntriesMut;
+}
+
 pub trait CsStorage<N, R, C = U1>: 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<N, R, C = U1>: for<'a> CsStorageIter<'a, N, R, C> {
     fn len(&self) -> usize;
 }
 
-pub trait CsStorageMut<N, R, C = U1>: CsStorage<N, R, C> {
+pub trait CsStorageMut<N, R, C = U1>:
-    /*
+    CsStorage<N, R, C> + for<'a> CsStorageIterMut<'a, N, R, C>
-    /// 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);
-    */
 }
 
 #[derive(Clone, Debug)]
@@ -133,6 +133,27 @@ where
     }
 }
 
+impl<'a, N: Scalar, R: Dim, C: Dim> CsStorageIterMut<'a, N, R, C> for CsVecStorage<N, R, C>
+where
+    DefaultAllocator: Allocator<usize, C>,
+{
+    type ColumnEntriesMut = iter::Zip<iter::Cloned<slice::Iter<'a, usize>>, 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<N: Scalar, R: Dim, C: Dim> CsStorageMut<N, R, C> for CsVecStorage<N, R, C> where
+    DefaultAllocator: Allocator<usize, C>
+{
+}
+
 /*
 pub struct CsSliceStorage<'a, N: Scalar, R: Dim, C: DimAdd<U1>> {
     shape: (R, C),

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<N, D, D>) -> 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!(

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;
 

tests/sparse/cs_cholesky.rs

@@ -51,5 +51,5 @@ fn test_cholesky(a: Matrix5<f32>) {
     println!("{}", l);
     println!("{}", cs_l);
 
-    assert_eq!(l, cs_l);
+    assert_relative_eq!(l, cs_l);
 }