Merge branch 'ssr'

Author: termoshtt

src/lib.rs

@@ -43,6 +43,7 @@ pub trait SquareMatrix: Matrix {
     /// inverse of matrix
     fn inv(self) -> Result<Self, LinalgError>;
     fn trace(&self) -> Result<Self::Scalar, LinalgError>;
+    fn ssqrt(self) -> Result<Self, LinalgError>;
     fn check_square(&self) -> Result<(), NotSquareError> {
         let (rows, cols) = self.size();
         if rows == cols {
@@ -86,7 +87,7 @@ impl<A: LapackScalar> Matrix for Array<A, (Ix, Ix)> {
     }
 }
 
-impl<A: LapackScalar> SquareMatrix for Array<A, (Ix, Ix)> {
+impl<A: LapackScalar + Float> SquareMatrix for Array<A, (Ix, Ix)> {
     fn eigh(self) -> Result<(Self::Vector, Self), LinalgError> {
         try!(self.check_square());
         let (rows, cols) = self.size();
@@ -107,6 +108,17 @@ impl<A: LapackScalar> SquareMatrix for Array<A, (Ix, Ix)> {
             Ok(m.reversed_axes())
         }
     }
+    fn ssqrt(self) -> Result<Self, LinalgError> {
+        let (n, _) = self.size();
+        let (e, v) = try!(self.eigh());
+        let mut res = Array::zeros((n, n));
+        for i in 0..n {
+            for j in 0..n {
+                res[(i, j)] = e[i].sqrt() * v[(j, i)];
+            }
+        }
+        Ok(v.dot(&res))
+    }
     fn trace(&self) -> Result<Self::Scalar, LinalgError> {
         try!(self.check_square());
         let (n, _) = self.size();

target/doc

@@ -0,0 +1,36 @@
+
+extern crate rand;
+extern crate ndarray;
+extern crate ndarray_rand;
+extern crate ndarray_linalg as linalg;
+
+use ndarray::prelude::*;
+use linalg::SquareMatrix;
+use rand::distributions::*;
+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 ssqrt_symmetric_random() {
+    let r_dist = Range::new(0., 1.);
+    let mut a = Array::<f64, _>::random((3, 3), r_dist);
+    a = a.dot(&a.t());
+    let ar = a.clone().ssqrt().unwrap();
+    all_close(&ar.clone().reversed_axes(), &ar);
+}
+
+#[test]
+fn ssqrt_sqrt_random() {
+    let r_dist = Range::new(0., 1.);
+    let mut a = Array::<f64, _>::random((3, 3), r_dist);
+    a = a.dot(&a.t());
+    let ar = a.clone().ssqrt().unwrap();
+    all_close(&ar.dot(&ar), &a);
+}