Cholesky decomposition rewrite (#150)

* Add minimal internal_utils module

* Basic implementation of Cholesky

* Rewrite lu unpack to use nullify from internal_utils

* Add Cholesky benchmarks

* Take raw_slice_mut() in nullify

* Rewrite Cholesky in terms of dot()

Benchmarks show significant performance gains.

* Cholesky: validate diagonal entries

* Implement det() for Cholesky

* Implement transpose_back_substitution

* Implement Cholesky::solve

* Add benchmarks for Cholesky::solve

* Rewrite transpose_back_sub.. for better data access pattern

* Enable previously disabled test

* Implement Cholesky::inverse()

* Documentation for Cholesky

* Fix errors in Cholesky doctests

* Remove stray ///

* Remove confusing comment

* Add Cholesky to decomposition module overview

* Deprecate .cholesky()

* Add warning to Cholesky::decompose

* Rewrite Cholesky solve/inverse to return Result

* Explain determinant of empty matrix to Cholesky/PartialPivLu

* Reorder Cholesky and FullPivLu in decomposition docs

benches/lib.rs

@@ -11,6 +11,7 @@ pub mod linalg {
 	mod matrix;
 	mod svd;
 	mod lu;
+	mod cholesky;
 	mod norm;
 	mod triangular;
 	mod permutation;

benches/linalg/cholesky.rs

@@ -0,0 +1,69 @@
+use rulinalg::matrix::Matrix;
+use rulinalg::matrix::decomposition::{Cholesky, Decomposition};
+use test::Bencher;
+
+#[bench]
+fn cholesky_decompose_unpack_100x100(b: &mut Bencher) {
+    let n = 100;
+    let x = Matrix::<f64>::identity(n);
+    b.iter(|| {
+        // Assume that the cost of cloning x is roughly
+        // negligible in comparison with the cost of LU
+        Cholesky::decompose(x.clone()).expect("Matrix is invertible")
+                                      .unpack()
+    })
+}
+
+#[bench]
+fn cholesky_decompose_unpack_500x500(b: &mut Bencher) {
+    let n = 500;
+    let x = Matrix::<f64>::identity(n);
+    b.iter(|| {
+        // Assume that the cost of cloning x is roughly
+        // negligible in comparison with the cost of LU
+        Cholesky::decompose(x.clone()).expect("Matrix is invertible")
+                                      .unpack()
+    })
+}
+
+#[bench]
+fn cholesky_100x100(b: &mut Bencher) {
+    // Benchmark for legacy cholesky(). Remove when
+    // cholesky() has been removed.
+    let n = 100;
+    let x = Matrix::<f64>::identity(n);
+    b.iter(|| {
+        x.cholesky().expect("Matrix is invertible")
+    })
+}
+
+#[bench]
+fn cholesky_500x500(b: &mut Bencher) {
+    // Benchmark for legacy cholesky(). Remove when
+    // cholesky() has been removed.
+    let n = 500;
+    let x = Matrix::<f64>::identity(n);
+    b.iter(|| {
+        x.cholesky().expect("Matrix is invertible")
+    })
+}
+
+#[bench]
+fn cholesky_solve_1000x1000(b: &mut Bencher) {
+    let n = 1000;
+    let x = Matrix::identity(n);
+    let cholesky = Cholesky::decompose(x).unwrap();
+    b.iter(|| {
+        cholesky.solve(vector![0.0; n])
+    });
+}
+
+#[bench]
+fn cholesky_solve_100x100(b: &mut Bencher) {
+    let n = 100;
+    let x = Matrix::identity(n);
+    let cholesky = Cholesky::decompose(x).unwrap();
+    b.iter(|| {
+        cholesky.solve(vector![0.0; n])
+    });
+}

src/internal_utils.rs

@@ -0,0 +1,52 @@
+use matrix::BaseMatrixMut;
+use libnum::Zero;
+
+pub fn nullify_lower_triangular_part<T, M>(matrix: &mut M)
+    where T: Zero, M: BaseMatrixMut<T> {
+    for (i, mut row) in matrix.row_iter_mut().enumerate() {
+        for element in row.raw_slice_mut().iter_mut().take(i) {
+            *element = T::zero();
+        }
+    }
+}
+
+pub fn nullify_upper_triangular_part<T, M>(matrix: &mut M)
+    where T: Zero, M: BaseMatrixMut<T> {
+    for (i, mut row) in matrix.row_iter_mut().enumerate() {
+        for element in row.raw_slice_mut().iter_mut().skip(i + 1) {
+            *element = T::zero();
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::nullify_lower_triangular_part;
+    use super::nullify_upper_triangular_part;
+
+    #[test]
+    fn nullify_lower_triangular_part_examples() {
+        let mut x = matrix![1.0, 2.0, 3.0;
+                            4.0, 5.0, 6.0;
+                            7.0, 8.0, 9.0];
+        nullify_lower_triangular_part(&mut x);
+        assert_matrix_eq!(x, matrix![
+            1.0, 2.0, 3.0;
+            0.0, 5.0, 6.0;
+            0.0, 0.0, 9.0
+        ]);
+    }
+
+    #[test]
+    fn nullify_upper_triangular_part_examples() {
+        let mut x = matrix![1.0, 2.0, 3.0;
+                            4.0, 5.0, 6.0;
+                            7.0, 8.0, 9.0];
+        nullify_upper_triangular_part(&mut x);
+        assert_matrix_eq!(x, matrix![
+            1.0, 0.0, 0.0;
+            4.0, 5.0, 0.0;
+            7.0, 8.0, 9.0
+        ]);
+    }
+}

src/lib.rs

@@ -100,6 +100,8 @@ pub mod vector;
 pub mod ulp;
 pub mod norm;
 
+mod internal_utils;
+
 #[cfg(test)]
 mod testsupport;