Merge pull request #32 from termoshtt/arraybase

Layout module

Author: termoshtt

Cargo.toml

@@ -15,6 +15,8 @@ openblas = ["blas/openblas", "lapack/openblas"]
 netlib   = ["blas/netlib", "lapack/netlib"]
 
 [dependencies]
+derive-new = "0.4"
+enum-error-derive = "0.1"
 num-traits  = "0.1"
 num-complex = "0.1"
 ndarray = { version = "0.9",  default-features = false, features = ["blas"] }

src/error.rs

@@ -4,7 +4,18 @@ use std::error;
 use std::fmt;
 use ndarray::{Ixs, ShapeError};
 
-#[derive(Debug)]
+pub type Result<T> = ::std::result::Result<T, LinalgError>;
+
+#[derive(Debug, EnumError)]
+pub enum LinalgError {
+    NotSquare(NotSquareError),
+    Lapack(LapackError),
+    Stride(StrideError),
+    MemoryCont(MemoryContError),
+    Shape(ShapeError),
+}
+
+#[derive(Debug, new)]
 pub struct LapackError {
     pub return_code: i32,
 }
@@ -27,10 +38,10 @@ impl From<i32> for LapackError {
     }
 }
 
-#[derive(Debug)]
+#[derive(Debug, new)]
 pub struct NotSquareError {
-    pub rows: usize,
-    pub cols: usize,
+    pub rows: i32,
+    pub cols: i32,
 }
 
 impl fmt::Display for NotSquareError {
@@ -45,7 +56,7 @@ impl error::Error for NotSquareError {
     }
 }
 
-#[derive(Debug)]
+#[derive(Debug, new)]
 pub struct StrideError {
     pub s0: Ixs,
     pub s1: Ixs,
@@ -63,55 +74,17 @@ impl error::Error for StrideError {
     }
 }
 
-#[derive(Debug)]
-pub enum LinalgError {
-    NotSquare(NotSquareError),
-    Lapack(LapackError),
-    Stride(StrideError),
-    Shape(ShapeError),
-}
+#[derive(Debug, new)]
+pub struct MemoryContError {}
 
-impl fmt::Display for LinalgError {
+impl fmt::Display for MemoryContError {
     fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
-        match *self {
-            LinalgError::NotSquare(ref err) => err.fmt(f),
-            LinalgError::Lapack(ref err) => err.fmt(f),
-            LinalgError::Stride(ref err) => err.fmt(f),
-            LinalgError::Shape(ref err) => err.fmt(f),
-        }
+        write!(f, "Memory is not contiguous")
     }
 }
 
-impl error::Error for LinalgError {
+impl error::Error for MemoryContError {
     fn description(&self) -> &str {
-        match *self {
-            LinalgError::NotSquare(ref err) => err.description(),
-            LinalgError::Lapack(ref err) => err.description(),
-            LinalgError::Stride(ref err) => err.description(),
-            LinalgError::Shape(ref err) => err.description(),
-        }
-    }
-}
-
-impl From<NotSquareError> for LinalgError {
-    fn from(err: NotSquareError) -> LinalgError {
-        LinalgError::NotSquare(err)
-    }
-}
-
-impl From<LapackError> for LinalgError {
-    fn from(err: LapackError) -> LinalgError {
-        LinalgError::Lapack(err)
-    }
-}
-
-impl From<StrideError> for LinalgError {
-    fn from(err: StrideError) -> LinalgError {
-        LinalgError::Stride(err)
-    }
-}
-impl From<ShapeError> for LinalgError {
-    fn from(err: ShapeError) -> LinalgError {
-        LinalgError::Shape(err)
+        "Memory is not contiguous"
     }
 }

src/impl2/mod.rs

@@ -0,0 +1,6 @@
+
+pub mod opnorm;
+pub use self::opnorm::*;
+
+pub trait LapackScalar: OperatorNorm_ {}
+impl<A> LapackScalar for A where A: OperatorNorm_ {}

src/impl2/opnorm.rs

@@ -0,0 +1,45 @@
+//! Implement Operator norms for matrices
+
+use lapack::c;
+use lapack::c::Layout::ColumnMajor as cm;
+
+use types::*;
+use layout::*;
+
+#[repr(u8)]
+pub enum NormType {
+    One = b'o',
+    Infinity = b'i',
+    Frobenius = b'f',
+}
+
+impl NormType {
+    fn transpose(self) -> Self {
+        match self {
+            NormType::One => NormType::Infinity,
+            NormType::Infinity => NormType::One,
+            NormType::Frobenius => NormType::Frobenius,
+        }
+    }
+}
+
+pub trait OperatorNorm_: AssociatedReal {
+    fn opnorm(NormType, Layout, &[Self]) -> Self::Real;
+}
+
+macro_rules! impl_opnorm {
+    ($scalar:ty, $lange:path) => {
+impl OperatorNorm_ for $scalar {
+    fn opnorm(t: NormType, l: Layout, a: &[Self]) -> Self::Real {
+        match l {
+            Layout::F((col, lda)) => $lange(cm, t as u8, lda, col, a, lda),
+            Layout::C((row, lda)) => $lange(cm, t.transpose() as u8, lda, row, a, lda),
+        }
+    }
+}
+}} // impl_opnorm!
+
+impl_opnorm!(f64, c::dlange);
+impl_opnorm!(f32, c::slange);
+impl_opnorm!(c64, c::zlange);
+impl_opnorm!(c32, c::clange);

src/impls/mod.rs

@@ -3,6 +3,5 @@ pub mod outer;
 pub mod qr;
 pub mod svd;
 pub mod eigh;
-pub mod opnorm;
 pub mod solve;
 pub mod cholesky;

src/impls/opnorm.rs

@@ -0,0 +1,62 @@
+
+use ndarray::*;
+
+use super::error::*;
+
+pub type LDA = i32;
+pub type Col = i32;
+pub type Row = i32;
+
+pub enum Layout {
+    C((Row, LDA)),
+    F((Col, LDA)),
+}
+
+impl Layout {
+    pub fn size(&self) -> (Row, Col) {
+        match *self {
+            Layout::C((row, lda)) => (row, lda),
+            Layout::F((col, lda)) => (lda, col),
+        }
+    }
+}
+
+pub trait AllocatedArray {
+    type Scalar;
+    fn layout(&self) -> Result<Layout>;
+    fn square_layout(&self) -> Result<Layout>;
+    fn as_allocated(&self) -> Result<&[Self::Scalar]>;
+}
+
+impl<A, S> AllocatedArray for ArrayBase<S, Ix2>
+    where S: Data<Elem = A>
+{
+    type Scalar = A;
+
+    fn layout(&self) -> Result<Layout> {
+        let strides = self.strides();
+        if ::std::cmp::min(strides[0], strides[1]) != 1 {
+            return Err(StrideError::new(strides[0], strides[1]).into());
+        }
+        if strides[0] > strides[1] {
+            Ok(Layout::C((self.rows() as i32, self.cols() as i32)))
+        } else {
+            Ok(Layout::F((self.cols() as i32, self.rows() as i32)))
+        }
+    }
+
+    fn square_layout(&self) -> Result<Layout> {
+        let l = self.layout()?;
+        let (n, m) = l.size();
+        if n == m {
+            Ok(l)
+        } else {
+            Err(NotSquareError::new(n, m).into())
+        }
+    }
+
+    fn as_allocated(&self) -> Result<&[A]> {
+        let slice = self.as_slice_memory_order().ok_or(MemoryContError::new())?;
+        Ok(slice)
+    }
+}

src/layout.rs

@@ -37,9 +37,18 @@ extern crate num_traits;
 extern crate num_complex;
 #[macro_use(s)]
 extern crate ndarray;
+#[macro_use]
+extern crate enum_error_derive;
+#[macro_use]
+extern crate derive_new;
 
-pub mod impls;
+pub mod types;
 pub mod error;
+pub mod layout;
+pub mod impls;
+pub mod impl2;
+
+pub mod traits;
 
 pub mod vector;
 pub mod matrix;

src/lib.rs

@@ -8,11 +8,10 @@ use lapack::c::Layout;
 use super::error::{LinalgError, StrideError};
 use super::impls::qr::ImplQR;
 use super::impls::svd::ImplSVD;
-use super::impls::opnorm::ImplOpNorm;
 use super::impls::solve::ImplSolve;
 
-pub trait MFloat: ImplQR + ImplSVD + ImplOpNorm + ImplSolve + NdFloat {}
-impl<A: ImplQR + ImplSVD + ImplOpNorm + ImplSolve + NdFloat> MFloat for A {}
+pub trait MFloat: ImplQR + ImplSVD + ImplSolve + NdFloat {}
+impl<A: ImplQR + ImplSVD + ImplSolve + NdFloat> MFloat for A {}
 
 /// Methods for general matrices
 pub trait Matrix: Sized {
@@ -23,12 +22,6 @@ pub trait Matrix: Sized {
     fn size(&self) -> (usize, usize);
     /// Layout (C/Fortran) of matrix
     fn layout(&self) -> Result<Layout, StrideError>;
-    /// Operator norm for L-1 norm
-    fn opnorm_1(&self) -> Self::Scalar;
-    /// Operator norm for L-inf norm
-    fn opnorm_i(&self) -> Self::Scalar;
-    /// Frobenius norm
-    fn opnorm_f(&self) -> Self::Scalar;
     /// singular-value decomposition (SVD)
     fn svd(self) -> Result<(Self, Self::Vector, Self), LinalgError>;
     /// QR decomposition
@@ -84,28 +77,6 @@ impl<A: MFloat> Matrix for Array<A, Ix2> {
     fn layout(&self) -> Result<Layout, StrideError> {
         check_layout(self.strides())
     }
-    fn opnorm_1(&self) -> Self::Scalar {
-        let (m, n) = self.size();
-        let strides = self.strides();
-        if strides[0] > strides[1] {
-            ImplOpNorm::opnorm_i(n, m, self.clone().into_raw_vec())
-        } else {
-            ImplOpNorm::opnorm_1(m, n, self.clone().into_raw_vec())
-        }
-    }
-    fn opnorm_i(&self) -> Self::Scalar {
-        let (m, n) = self.size();
-        let strides = self.strides();
-        if strides[0] > strides[1] {
-            ImplOpNorm::opnorm_1(n, m, self.clone().into_raw_vec())
-        } else {
-            ImplOpNorm::opnorm_i(m, n, self.clone().into_raw_vec())
-        }
-    }
-    fn opnorm_f(&self) -> Self::Scalar {
-        let (m, n) = self.size();
-        ImplOpNorm::opnorm_f(m, n, self.clone().into_raw_vec())
-    }
     fn svd(self) -> Result<(Self, Self::Vector, Self), LinalgError> {
         let (n, m) = self.size();
         let layout = self.layout()?;
@@ -192,18 +163,6 @@ impl<A: MFloat> Matrix for RcArray<A, Ix2> {
     fn layout(&self) -> Result<Layout, StrideError> {
         check_layout(self.strides())
     }
-    fn opnorm_1(&self) -> Self::Scalar {
-        // XXX unnecessary clone
-        self.to_owned().opnorm_1()
-    }
-    fn opnorm_i(&self) -> Self::Scalar {
-        // XXX unnecessary clone
-        self.to_owned().opnorm_i()
-    }
-    fn opnorm_f(&self) -> Self::Scalar {
-        // XXX unnecessary clone
-        self.to_owned().opnorm_f()
-    }
     fn svd(self) -> Result<(Self, Self::Vector, Self), LinalgError> {
         let (u, s, v) = self.into_owned().svd()?;
         Ok((u.into_shared(), s.into_shared(), v.into_shared()))

src/matrix.rs

@@ -5,3 +5,4 @@ pub use hermite::HermiteMatrix;
 pub use triangular::*;
 pub use util::*;
 pub use assert::*;
+pub use traits::*;