diff --git a/Cargo.toml b/Cargo.toml index 071f00e9..bc60de63 100644 --- a/Cargo.toml +++ b/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"] } diff --git a/src/error.rs b/src/error.rs index efbb7c27..d98df7b3 100644 --- a/src/error.rs +++ b/src/error.rs @@ -4,7 +4,18 @@ use std::error; use std::fmt; use ndarray::{Ixs, ShapeError}; -#[derive(Debug)] +pub type Result = ::std::result::Result; + +#[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 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 for LinalgError { - fn from(err: NotSquareError) -> LinalgError { - LinalgError::NotSquare(err) - } -} - -impl From for LinalgError { - fn from(err: LapackError) -> LinalgError { - LinalgError::Lapack(err) - } -} - -impl From for LinalgError { - fn from(err: StrideError) -> LinalgError { - LinalgError::Stride(err) - } -} -impl From for LinalgError { - fn from(err: ShapeError) -> LinalgError { - LinalgError::Shape(err) + "Memory is not contiguous" } } diff --git a/src/impl2/mod.rs b/src/impl2/mod.rs new file mode 100644 index 00000000..89a38823 --- /dev/null +++ b/src/impl2/mod.rs @@ -0,0 +1,6 @@ + +pub mod opnorm; +pub use self::opnorm::*; + +pub trait LapackScalar: OperatorNorm_ {} +impl LapackScalar for A where A: OperatorNorm_ {} diff --git a/src/impl2/opnorm.rs b/src/impl2/opnorm.rs new file mode 100644 index 00000000..e2687ebf --- /dev/null +++ b/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); diff --git a/src/impls/mod.rs b/src/impls/mod.rs index 8239548b..7a97a974 100644 --- a/src/impls/mod.rs +++ b/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; diff --git a/src/impls/opnorm.rs b/src/impls/opnorm.rs deleted file mode 100644 index 22a9f11a..00000000 --- a/src/impls/opnorm.rs +++ /dev/null @@ -1,27 +0,0 @@ -//! Implement Norms for matrices - -use lapack::c::*; - -pub trait ImplOpNorm: Sized { - fn opnorm_1(m: usize, n: usize, a: Vec) -> Self; - fn opnorm_i(m: usize, n: usize, a: Vec) -> Self; - fn opnorm_f(m: usize, n: usize, a: Vec) -> Self; -} - -macro_rules! impl_opnorm { - ($scalar:ty, $lange:path) => { -impl ImplOpNorm for $scalar { - fn opnorm_1(m: usize, n: usize, mut a: Vec) -> Self { - $lange(Layout::ColumnMajor, b'o', m as i32, n as i32, &mut a, m as i32) - } - fn opnorm_i(m: usize, n: usize, mut a: Vec) -> Self { - $lange(Layout::ColumnMajor, b'i', m as i32, n as i32, &mut a, m as i32) - } - fn opnorm_f(m: usize, n: usize, mut a: Vec) -> Self { - $lange(Layout::ColumnMajor, b'f', m as i32, n as i32, &mut a, m as i32) - } -} -}} // end macro_rules - -impl_opnorm!(f64, dlange); -impl_opnorm!(f32, slange); diff --git a/src/layout.rs b/src/layout.rs new file mode 100644 index 00000000..66e56c19 --- /dev/null +++ b/src/layout.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; + fn square_layout(&self) -> Result; + fn as_allocated(&self) -> Result<&[Self::Scalar]>; +} + +impl AllocatedArray for ArrayBase + where S: Data +{ + type Scalar = A; + + fn layout(&self) -> Result { + 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 { + 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) + } +} diff --git a/src/lib.rs b/src/lib.rs index 2816df3c..1d4f537e 100644 --- a/src/lib.rs +++ b/src/lib.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; diff --git a/src/matrix.rs b/src/matrix.rs index e90911a0..ea6c5f79 100644 --- a/src/matrix.rs +++ b/src/matrix.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 MFloat for A {} +pub trait MFloat: ImplQR + ImplSVD + ImplSolve + NdFloat {} +impl 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; - /// 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 Matrix for Array { fn layout(&self) -> Result { 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 Matrix for RcArray { fn layout(&self) -> Result { 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())) diff --git a/src/prelude.rs b/src/prelude.rs index c445c0f7..9fb296b2 100644 --- a/src/prelude.rs +++ b/src/prelude.rs @@ -5,3 +5,4 @@ pub use hermite::HermiteMatrix; pub use triangular::*; pub use util::*; pub use assert::*; +pub use traits::*; diff --git a/src/square.rs b/src/square.rs index aedde160..709abd9a 100644 --- a/src/square.rs +++ b/src/square.rs @@ -25,8 +25,8 @@ pub trait SquareMatrix: Matrix { Ok(()) } else { Err(NotSquareError { - rows: rows, - cols: cols, + rows: rows as i32, + cols: cols as i32, }) } } diff --git a/src/traits.rs b/src/traits.rs new file mode 100644 index 00000000..5dada9c6 --- /dev/null +++ b/src/traits.rs @@ -0,0 +1,36 @@ + +pub use impl2::LapackScalar; +pub use impl2::NormType; + +use ndarray::*; + +use super::types::*; +use super::error::*; +use super::layout::*; + +pub trait OperationNorm { + type Output; + fn opnorm(&self, t: NormType) -> Self::Output; + fn opnorm_one(&self) -> Self::Output { + self.opnorm(NormType::One) + } + fn opnorm_inf(&self) -> Self::Output { + self.opnorm(NormType::Infinity) + } + fn opnorm_fro(&self) -> Self::Output { + self.opnorm(NormType::Frobenius) + } +} + +impl OperationNorm for ArrayBase + where A: LapackScalar + AssociatedReal, + S: Data +{ + type Output = Result; + + fn opnorm(&self, t: NormType) -> Self::Output { + let l = self.layout()?; + let a = self.as_allocated()?; + Ok(A::opnorm(t, l, a)) + } +} diff --git a/src/types.rs b/src/types.rs new file mode 100644 index 00000000..2d0393d8 --- /dev/null +++ b/src/types.rs @@ -0,0 +1,31 @@ + +pub use num_complex::Complex32 as c32; +pub use num_complex::Complex64 as c64; +use num_traits::Float; +use std::ops::*; + +pub trait AssociatedReal: Sized { + type Real: Float + Mul; +} +pub trait AssociatedComplex: Sized { + type Complex; +} + +macro_rules! impl_assoc { + ($real:ty, $complex:ty) => { +impl AssociatedReal for $real { + type Real = $real; +} +impl AssociatedReal for $complex { + type Real = $real; +} +impl AssociatedComplex for $real { + type Complex = $complex; +} +impl AssociatedComplex for $complex { + type Complex = $complex; +} +}} // impl_assoc! + +impl_assoc!(f64, c64); +impl_assoc!(f32, c32); diff --git a/tests/opnorm.rs b/tests/opnorm.rs index 0569a1b0..eda418a1 100644 --- a/tests/opnorm.rs +++ b/tests/opnorm.rs @@ -5,9 +5,12 @@ macro_rules! impl_test { #[test] fn $funcname() { let a = $a; - assert_rclose!(a.opnorm_1(), $op1, 1e-7); - assert_rclose!(a.opnorm_i(), $opi, 1e-7); - assert_rclose!(a.opnorm_f(), $opf, 1e-7); + println!("ONE = {:?}", a.opnorm_one()); + println!("INF = {:?}", a.opnorm_inf()); + println!("FRO = {:?}", a.opnorm_fro()); + assert_rclose!(a.opnorm_fro().unwrap(), $opf, 1e-7; "Frobenius norm"); + assert_rclose!(a.opnorm_one().unwrap(), $op1, 1e-7; "One norm"); + assert_rclose!(a.opnorm_inf().unwrap(), $opi, 1e-7; "Infinity norm"); } }} // impl_test