Merge branch 'cholesky'

Author: termoshtt

src/cholesky.rs

@@ -0,0 +1,24 @@
+
+use lapack::c::*;
+use error::LapackError;
+
+pub trait ImplCholesky: Sized {
+    fn cholesky(layout: Layout, n: usize, a: Vec<Self>) -> Result<Vec<Self>, LapackError>;
+}
+
+macro_rules! impl_cholesky {
+    ($scalar:ty, $potrf:path) => {
+impl ImplCholesky for $scalar {
+    fn cholesky(layout: Layout, n: usize, mut a: Vec<Self>) -> Result<Vec<Self>, LapackError> {
+        let info = $potrf(layout, b'U', n as i32, &mut a, n as i32);
+        if info == 0 {
+            Ok(a)
+        } else {
+            Err(From::from(info))
+        }
+    }
+}
+}} // end macro_rules
+
+impl_cholesky!(f64, dpotrf);
+impl_cholesky!(f32, spotrf);

src/hermite.rs

@@ -4,6 +4,7 @@ use std::fmt::Debug;
 use ndarray::prelude::*;
 use ndarray::LinalgScalar;
 use num_traits::float::Float;
+use lapack::c::Layout;
 
 use matrix::Matrix;
 use square::SquareMatrix;
@@ -13,17 +14,20 @@ use qr::ImplQR;
 use svd::ImplSVD;
 use norm::ImplNorm;
 use solve::ImplSolve;
+use cholesky::ImplCholesky;
 
 /// Methods for Hermite matrix
 pub trait HermiteMatrix: SquareMatrix + Matrix {
     /// eigenvalue decomposition
     fn eigh(self) -> Result<(Self::Vector, Self), LinalgError>;
     /// symmetric square root of Hermite matrix
     fn ssqrt(self) -> Result<Self, LinalgError>;
+    /// Cholesky factorization
+    fn cholesky(self) -> Result<Self, LinalgError>;
 }
 
 impl<A> HermiteMatrix for Array<A, (Ix, Ix)>
-    where A: ImplQR + ImplSVD + ImplNorm + ImplSolve + ImplEigh + LinalgScalar + Float + Debug
+    where A: ImplQR + ImplSVD + ImplNorm + ImplSolve + ImplEigh + ImplCholesky + LinalgScalar + Float + Debug
 {
     fn eigh(self) -> Result<(Self::Vector, Self), LinalgError> {
         try!(self.check_square());
@@ -44,4 +48,21 @@ impl<A> HermiteMatrix for Array<A, (Ix, Ix)>
         }
         Ok(v.dot(&res))
     }
+    fn cholesky(self) -> Result<Self, LinalgError> {
+        try!(self.check_square());
+        println!("layout = {:?}", self.layout());
+        let (n, _) = self.size();
+        let layout = self.layout();
+        let a = try!(ImplCholesky::cholesky(layout, n, self.into_raw_vec()));
+        let mut c = match layout {
+          Layout::RowMajor => Array::from_vec(a).into_shape((n, n)).unwrap(),
+          Layout::ColumnMajor => Array::from_vec(a).into_shape((n, n)).unwrap().reversed_axes(),
+        };
+        for ((i, j), val) in c.indexed_iter_mut() {
+            if i > j {
+                *val = A::zero();
+            }
+        }
+        Ok(c)
+    }
 }

src/lib.rs

@@ -41,3 +41,4 @@ pub mod svd;
 pub mod eigh;
 pub mod norm;
 pub mod solve;
+pub mod cholesky;

tests/cholesky.rs

@@ -0,0 +1,42 @@
+
+extern crate rand;
+extern crate ndarray;
+extern crate ndarray_rand;
+extern crate ndarray_linalg;
+
+use rand::distributions::*;
+use ndarray::prelude::*;
+use ndarray_linalg::prelude::*;
+use ndarray_rand::RandomExt;
+
+fn all_close(a: Array<f64, (Ix, Ix)>, b: Array<f64, (Ix, Ix)>) {
+    if !a.all_close(&b, 1.0e-7) {
+        panic!("\nTwo matrices are not equal:\na = \n{:?}\nb = \n{:?}\n",
+               a,
+               b);
+    }
+}
+
+#[test]
+fn cholesky() {
+    let r_dist = Range::new(0., 1.);
+    let mut a = Array::<f64, _>::random((3, 3), r_dist);
+    a = a.dot(&a.t());
+    println!("a = \n{:?}", a);
+    let c = a.clone().cholesky().unwrap();
+    println!("c = \n{:?}", c);
+    println!("cc = \n{:?}", c.t().dot(&c));
+    all_close(c.t().dot(&c), a);
+}
+
+#[test]
+fn cholesky_t() {
+    let r_dist = Range::new(0., 1.);
+    let mut a = Array::<f64, _>::random((3, 3), r_dist);
+    a = a.dot(&a.t()).reversed_axes();
+    println!("a = \n{:?}", a);
+    let c = a.clone().cholesky().unwrap();
+    println!("c = \n{:?}", c);
+    println!("cc = \n{:?}", c.t().dot(&c));
+    all_close(c.t().dot(&c), a);
+}