Introduce Field, RealField

Author: termoshtt

src/assert.rs

@@ -1,32 +1,30 @@
 //! Assertions for array
 
-use std::iter::Sum;
-use num_traits::Float;
 use ndarray::*;
 
 use super::types::*;
 use super::vector::*;
 
 pub fn rclose<A, Tol>(test: A, truth: A, rtol: Tol) -> Result<Tol, Tol>
-    where A: LinalgScalar + Absolute<Output = Tol>,
-          Tol: Float
+    where A: Field + Absolute<Output = Tol>,
+          Tol: RealField
 {
     let dev = (test - truth).abs() / truth.abs();
     if dev < rtol { Ok(dev) } else { Err(dev) }
 }
 
 pub fn aclose<A, Tol>(test: A, truth: A, atol: Tol) -> Result<Tol, Tol>
-    where A: LinalgScalar + Absolute<Output = Tol>,
-          Tol: Float
+    where A: Field + Absolute<Output = Tol>,
+          Tol: RealField
 {
     let dev = (test - truth).abs();
     if dev < atol { Ok(dev) } else { Err(dev) }
 }
 
 /// check two arrays are close in maximum norm
 pub fn close_max<A, Tol, S1, S2, D>(test: &ArrayBase<S1, D>, truth: &ArrayBase<S2, D>, atol: Tol) -> Result<Tol, Tol>
-    where A: LinalgScalar + Absolute<Output = Tol>,
-          Tol: Float + Sum,
+    where A: Field + Absolute<Output = Tol>,
+          Tol: RealField,
           S1: Data<Elem = A>,
           S2: Data<Elem = A>,
           D: Dimension
@@ -37,8 +35,8 @@ pub fn close_max<A, Tol, S1, S2, D>(test: &ArrayBase<S1, D>, truth: &ArrayBase<S
 
 /// check two arrays are close in L1 norm
 pub fn close_l1<A, Tol, S1, S2, D>(test: &ArrayBase<S1, D>, truth: &ArrayBase<S2, D>, rtol: Tol) -> Result<Tol, Tol>
-    where A: LinalgScalar + Absolute<Output = Tol>,
-          Tol: Float + Sum,
+    where A: Field + Absolute<Output = Tol>,
+          Tol: RealField,
           S1: Data<Elem = A>,
           S2: Data<Elem = A>,
           D: Dimension
@@ -49,8 +47,8 @@ pub fn close_l1<A, Tol, S1, S2, D>(test: &ArrayBase<S1, D>, truth: &ArrayBase<S2
 
 /// check two arrays are close in L2 norm
 pub fn close_l2<A, Tol, S1, S2, D>(test: &ArrayBase<S1, D>, truth: &ArrayBase<S2, D>, rtol: Tol) -> Result<Tol, Tol>
-    where A: LinalgScalar + Absolute<Output = Tol>,
-          Tol: Float + Sum,
+    where A: Field + Absolute<Output = Tol>,
+          Tol: RealField,
           S1: Data<Elem = A>,
           S2: Data<Elem = A>,
           D: Dimension

src/impl2/triangular.rs

@@ -3,6 +3,7 @@
 use lapack::c;
 
 use error::*;
+use types::*;
 use layout::Layout;
 use super::{UPLO, Transpose, into_result};
 
@@ -18,11 +19,14 @@ pub trait Triangular_: Sized {
     fn solve_triangular(al: Layout, bl: Layout, UPLO, Diag, a: &[Self], b: &mut [Self]) -> Result<()>;
 }
 
-impl Triangular_ for f64 {
+macro_rules! impl_triangular {
+    ($scalar:ty, $trtri:path, $trtrs:path) => {
+
+impl Triangular_ for $scalar {
     fn inv_triangular(l: Layout, uplo: UPLO, diag: Diag, a: &mut [Self]) -> Result<()> {
         let (n, _) = l.size();
         let lda = l.lda();
-        let info = c::dtrtri(l.lapacke_layout(), uplo as u8, diag as u8, n, a, lda);
+        let info = $trtri(l.lapacke_layout(), uplo as u8, diag as u8, n, a, lda);
         into_result(info, ())
     }
 
@@ -31,16 +35,14 @@ impl Triangular_ for f64 {
         let lda = al.lda();
         let nrhs = bl.len();
         let ldb = bl.lda();
-        let info = c::dtrtrs(al.lapacke_layout(),
-                             uplo as u8,
-                             Transpose::No as u8,
-                             diag as u8,
-                             n,
-                             nrhs,
-                             a,
-                             lda,
-                             &mut b,
-                             ldb);
+        let info = $trtrs(al.lapacke_layout(), uplo as u8, Transpose::No as u8, diag as u8, n, nrhs, a, lda, &mut b, ldb);
         into_result(info, ())
     }
 }
+
+}} // impl_triangular!
+
+impl_triangular!(f64, c::dtrtri, c::dtrtrs);
+impl_triangular!(f32, c::strtri, c::strtrs);
+impl_triangular!(c64, c::ztrtri, c::ztrtrs);
+impl_triangular!(c32, c::ctrtri, c::ctrtrs);

src/prelude.rs

@@ -13,3 +13,4 @@ pub use opnorm::*;
 pub use solve::*;
 pub use eigh::*;
 pub use cholesky::*;
+pub use impl2::LapackScalar;

src/triangular.rs

@@ -7,6 +7,8 @@ use super::layout::*;
 use super::error::*;
 use super::impl2::*;
 
+pub use super::impl2::Diag;
+
 /// solve a triangular system with upper triangular matrix
 pub trait SolveTriangular<Rhs> {
     type Output;

src/types.rs

@@ -1,11 +1,39 @@
 
-pub use num_complex::Complex32 as c32;
-pub use num_complex::Complex64 as c64;
-use num_complex::Complex;
-use num_traits::Float;
 use std::ops::*;
+use std::fmt::Debug;
+use std::iter::Sum;
+use num_complex::Complex;
+use num_traits::*;
 use rand::Rng;
 use rand::distributions::*;
+use ndarray::LinalgScalar;
+
+use super::impl2::LapackScalar;
+
+pub use num_complex::Complex32 as c32;
+pub use num_complex::Complex64 as c64;
+
+macro_rules! trait_alias {
+    ($name:ident: $($t:ident),*) => {
+
+pub trait $name : $($t +)* {}
+
+impl<T> $name for T where T: $($t +)* {}
+
+}} // trait_alias!
+
+trait_alias!(Field: LapackScalar,
+             LinalgScalar,
+             AssociatedReal,
+             AssociatedComplex,
+             Absolute,
+             SquareRoot,
+             Conjugate,
+             RandNormal,
+             Sum,
+             Debug);
+
+trait_alias!(RealField: Field, Float);
 
 pub trait AssociatedReal: Sized {
     type Real: Float + Mul<Self, Output = Self>;
@@ -16,13 +44,17 @@ pub trait AssociatedComplex: Sized {
 
 /// Field with norm
 pub trait Absolute {
-    type Output: Float;
+    type Output: RealField;
     fn squared(&self) -> Self::Output;
     fn abs(&self) -> Self::Output {
         self.squared().sqrt()
     }
 }
 
+pub trait SquareRoot {
+    fn sqrt(&self) -> Self;
+}
+
 pub trait Conjugate: Copy {
     fn conj(self) -> Self;
 }
@@ -70,6 +102,18 @@ impl Absolute for $complex {
     }
 }
 
+impl SquareRoot for $real {
+    fn sqrt(&self) -> Self {
+        Float::sqrt(*self)
+    }
+}
+
+impl SquareRoot for $complex {
+    fn sqrt(&self) -> Self {
+        Complex::sqrt(self)
+    }
+}
+
 impl Conjugate for $real {
     fn conj(self) -> Self {
         self

src/util.rs

@@ -1,6 +1,5 @@
 //! misc utilities
 
-use std::iter::Sum;
 use ndarray::*;
 use num_traits::Float;
 use std::ops::Div;
@@ -15,9 +14,9 @@ pub enum NormalizeAxis {
 
 /// normalize in L2 norm
 pub fn normalize<A, S, T>(mut m: ArrayBase<S, Ix2>, axis: NormalizeAxis) -> (ArrayBase<S, Ix2>, Vec<T>)
-    where A: LinalgScalar + Absolute<Output = T> + Div<T, Output = A>,
+    where A: Field + Absolute<Output = T> + Div<T, Output = A>,
           S: DataMut<Elem = A>,
-          T: Float + Sum
+          T: Field + Float
 {
     let mut ms = Vec::new();
     for mut v in m.axis_iter_mut(Axis(axis as usize)) {

src/vector.rs

@@ -1,6 +1,5 @@
 //! Define trait for vectors
 
-use std::iter::Sum;
 use ndarray::*;
 use num_traits::Float;
 
@@ -24,8 +23,8 @@ pub trait Norm {
 }
 
 impl<A, S, D, T> Norm for ArrayBase<S, D>
-    where A: LinalgScalar + Absolute<Output = T>,
-          T: Float + Sum,
+    where A: Field + Absolute<Output = T>,
+          T: Field + Float,
           S: Data<Elem = A>,
           D: Dimension
 {