Merge pull request #10 from termoshtt/norm

Implement Matrix and Vector norms

Author: termoshtt

Cargo.toml

@@ -6,3 +6,4 @@ authors = ["Toshiki Teramura <[email protected]>"]
 [dependencies]
 ndarray = "*"
 lapack = "*"
+num-traits = "*"

src/lapack_binding.rs

@@ -3,14 +3,18 @@ extern crate lapack;
 
 use self::lapack::fortran::*;
 use error::LapackError;
+use ndarray::LinalgScalar;
 
 /// Eigenvalue decomposition for Hermite matrix
-pub trait Eigh: Sized {
+pub trait LapackScalar: LinalgScalar {
     /// execute *syev subroutine
     fn eigh(row_size: usize, matrix: Vec<Self>) -> Result<(Vec<Self>, Vec<Self>), LapackError>;
+    fn norm_1(rows: usize, cols: usize, matrix: Vec<Self>) -> Self;
+    fn norm_i(rows: usize, cols: usize, matrix: Vec<Self>) -> Self;
+    fn norm_f(rows: usize, cols: usize, matrix: Vec<Self>) -> Self;
 }
 
-impl Eigh for f64 {
+impl LapackScalar for f64 {
     fn eigh(n: usize, mut a: Vec<Self>) -> Result<(Vec<Self>, Vec<Self>), LapackError> {
         let mut w = vec![0.0; n ];
         let mut work = vec![0.0; 4 * n ];
@@ -30,9 +34,22 @@ impl Eigh for f64 {
             Err(From::from(info))
         }
     }
+
+    fn norm_1(m: usize, n: usize, mut a: Vec<Self>) -> Self {
+        let mut work = Vec::<Self>::new();
+        dlange(b'o', m as i32, n as i32, &mut a, m as i32, &mut work)
+    }
+    fn norm_i(m: usize, n: usize, mut a: Vec<Self>) -> Self {
+        let mut work = vec![0.0; m];
+        dlange(b'i', m as i32, n as i32, &mut a, m as i32, &mut work)
+    }
+    fn norm_f(m: usize, n: usize, mut a: Vec<Self>) -> Self {
+        let mut work = Vec::<Self>::new();
+        dlange(b'f', m as i32, n as i32, &mut a, m as i32, &mut work)
+    }
 }
 
-impl Eigh for f32 {
+impl LapackScalar for f32 {
     fn eigh(n: usize, mut a: Vec<Self>) -> Result<(Vec<Self>, Vec<Self>), LapackError> {
         let mut w = vec![0.0; n];
         let mut work = vec![0.0; 4 * n];
@@ -52,4 +69,16 @@ impl Eigh for f32 {
             Err(From::from(info))
         }
     }
+    fn norm_1(m: usize, n: usize, mut a: Vec<Self>) -> Self {
+        let mut work = Vec::<Self>::new();
+        slange(b'o', m as i32, n as i32, &mut a, m as i32, &mut work)
+    }
+    fn norm_i(m: usize, n: usize, mut a: Vec<Self>) -> Self {
+        let mut work = vec![0.0; m];
+        slange(b'i', m as i32, n as i32, &mut a, m as i32, &mut work)
+    }
+    fn norm_f(m: usize, n: usize, mut a: Vec<Self>) -> Self {
+        let mut work = Vec::<Self>::new();
+        slange(b'f', m as i32, n as i32, &mut a, m as i32, &mut work)
+    }
 }

src/lib.rs

@@ -1,18 +1,44 @@
 
 extern crate ndarray;
+extern crate num_traits;
 
 pub mod lapack_binding;
 pub mod error;
 
 use ndarray::prelude::*;
-use ndarray::LinalgScalar;
-use lapack_binding::Eigh;
+use ndarray::{LinalgScalar, DataOwned};
+use num_traits::float::Float;
+use lapack_binding::LapackScalar;
 use error::{LinalgError, NotSquareError};
 
+pub trait Vector {
+    type Scalar;
+    fn norm(&self) -> Self::Scalar;
+}
+
+impl<A: Float + LinalgScalar> Vector for Array<A, Ix> {
+    type Scalar = A;
+    fn norm(&self) -> Self::Scalar {
+        self.dot(&self).sqrt()
+    }
+}
+
+pub fn norm<S, A>(v: &ArrayBase<S, Ix>) -> A
+    where S: DataOwned<Elem = A>,
+          A: Float + LinalgScalar
+{
+    let n2 = &v.dot(&v);
+    n2.sqrt()
+}
+
 pub trait Matrix: Sized {
+    type Scalar;
     type Vector;
     /// number of rows and cols
     fn size(&self) -> (usize, usize);
+    fn norm_1(&self) -> Self::Scalar;
+    fn norm_i(&self) -> Self::Scalar;
+    fn norm_f(&self) -> Self::Scalar;
     // fn svd(self) -> (Self, Self::Vector, Self);
 }
 
@@ -35,20 +61,41 @@ pub trait SquareMatrix: Matrix {
     }
 }
 
-impl<A> Matrix for Array<A, (Ix, Ix)> {
+impl<A: LapackScalar> Matrix for Array<A, (Ix, Ix)> {
+    type Scalar = A;
     type Vector = Array<A, Ix>;
     fn size(&self) -> (usize, usize) {
         (self.rows(), self.cols())
     }
+    fn norm_1(&self) -> Self::Scalar {
+        let (m, n) = self.size();
+        let strides = self.strides();
+        if strides[0] > strides[1] {
+            LapackScalar::norm_i(n, m, self.clone().into_raw_vec())
+        } else {
+            LapackScalar::norm_1(m, n, self.clone().into_raw_vec())
+        }
+    }
+    fn norm_i(&self) -> Self::Scalar {
+        let (m, n) = self.size();
+        let strides = self.strides();
+        if strides[0] > strides[1] {
+            LapackScalar::norm_1(n, m, self.clone().into_raw_vec())
+        } else {
+            LapackScalar::norm_i(m, n, self.clone().into_raw_vec())
+        }
+    }
+    fn norm_f(&self) -> Self::Scalar {
+        let (m, n) = self.size();
+        LapackScalar::norm_f(m, n, self.clone().into_raw_vec())
+    }
 }
 
-impl<A> SquareMatrix for Array<A, (Ix, Ix)>
-    where A: Eigh + LinalgScalar
-{
+impl<A: LapackScalar> SquareMatrix for Array<A, (Ix, Ix)> {
     fn eigh(self) -> Result<(Self::Vector, Self), LinalgError> {
         try!(self.check_square());
         let (rows, cols) = self.size();
-        let (w, a) = try!(Eigh::eigh(rows, self.into_raw_vec()));
+        let (w, a) = try!(LapackScalar::eigh(rows, self.into_raw_vec()));
         let ea = Array::from_vec(w);
         let va = Array::from_vec(a).into_shape((rows, cols)).unwrap().reversed_axes();
         Ok((ea, va))

tests/norm.rs

@@ -0,0 +1,75 @@
+
+extern crate ndarray;
+extern crate ndarray_linalg as linalg;
+extern crate num_traits;
+
+use ndarray::prelude::*;
+use linalg::{Matrix, Vector};
+use num_traits::float::Float;
+
+fn assert_almost_eq(a: f64, b: f64) {
+    let rel_dev = (a - b).abs() / (a.abs() + b.abs());
+    if rel_dev > 1.0e-7 {
+        panic!("a={:?}, b={:?} are not almost equal", a, b);
+    }
+}
+
+#[test]
+fn vector_norm() {
+    let a = Array::range(1., 10., 1.);
+    assert_almost_eq(a.norm(), 285.0.sqrt());
+}
+
+#[test]
+fn matrix_norm_square() {
+    let a = Array::range(1., 10., 1.).into_shape((3, 3)).unwrap();
+    assert_almost_eq(a.norm_1(), 18.0);
+    assert_almost_eq(a.norm_i(), 24.0);
+    assert_almost_eq(a.norm_f(), 285.0.sqrt());
+}
+
+#[test]
+fn matrix_norm_square_t() {
+    let a = Array::range(1., 10., 1.).into_shape((3, 3)).unwrap().reversed_axes();
+    assert_almost_eq(a.norm_1(), 24.0);
+    assert_almost_eq(a.norm_i(), 18.0);
+    assert_almost_eq(a.norm_f(), 285.0.sqrt());
+}
+
+#[test]
+fn matrix_norm_3x4() {
+    let a = Array::range(1., 13., 1.).into_shape((3, 4)).unwrap();
+    assert_almost_eq(a.norm_1(), 24.0);
+    assert_almost_eq(a.norm_i(), 42.0);
+    assert_almost_eq(a.norm_f(), 650.0.sqrt());
+}
+
+#[test]
+fn matrix_norm_3x4_t() {
+    let a = Array::range(1., 13., 1.)
+        .into_shape((3, 4))
+        .unwrap()
+        .reversed_axes();
+    assert_almost_eq(a.norm_1(), 42.0);
+    assert_almost_eq(a.norm_i(), 24.0);
+    assert_almost_eq(a.norm_f(), 650.0.sqrt());
+}
+
+#[test]
+fn matrix_norm_4x3() {
+    let a = Array::range(1., 13., 1.).into_shape((4, 3)).unwrap();
+    assert_almost_eq(a.norm_1(), 30.0);
+    assert_almost_eq(a.norm_i(), 33.0);
+    assert_almost_eq(a.norm_f(), 650.0.sqrt());
+}
+
+#[test]
+fn matrix_norm_4x3_t() {
+    let a = Array::range(1., 13., 1.)
+        .into_shape((4, 3))
+        .unwrap()
+        .reversed_axes();
+    assert_almost_eq(a.norm_1(), 33.0);
+    assert_almost_eq(a.norm_i(), 30.0);
+    assert_almost_eq(a.norm_f(), 650.0.sqrt());
+}