Merge pull request #20 from termoshtt/memory_layout

termoshtt · web-flow · commit 86f38a74483f · 2016-11-23T01:40:14.000+09:00
unify memory layout management
diff --git a/.gitignore b/.gitignore
@@ -5,3 +5,6 @@
 # Remove Cargo.lock from gitignore if creating an executable, leave it for libraries
 # More information here http://doc.crates.io/guide.html#cargotoml-vs-cargolock
 Cargo.lock
+
+# cargo fmt
+*.bk
diff --git a/rustfmt.toml b/rustfmt.toml
@@ -1 +1,2 @@
 max_width = 120
+use_try_shorthand = true
diff --git a/src/eigh.rs b/src/eigh.rs
@@ -1,30 +1,20 @@
 //! Implement eigenvalue decomposition of Hermite matrix
 
-use lapack::fortran::*;
+use lapack::c::*;
 use num_traits::Zero;
 
 use error::LapackError;
 
 pub trait ImplEigh: Sized {
-    fn eigh(n: usize, mut a: Vec<Self>) -> Result<(Vec<Self>, Vec<Self>), LapackError>;
+    fn eigh(layout: Layout, n: usize, mut a: Vec<Self>) -> Result<(Vec<Self>, Vec<Self>), LapackError>;
 }
 
 macro_rules! impl_eigh {
     ($scalar:ty, $syev:path) => {
 impl ImplEigh for $scalar {
-    fn eigh(n: usize, mut a: Vec<Self>) -> Result<(Vec<Self>, Vec<Self>), LapackError> {
+    fn eigh(layout: Layout, n: usize, mut a: Vec<Self>) -> Result<(Vec<Self>, Vec<Self>), LapackError> {
         let mut w = vec![Self::zero(); n];
-        let mut work = vec![Self::zero(); 4 * n];
-        let mut info = 0;
-        $syev(b'V',
-              b'U',
-              n as i32,
-              &mut a,
-              n as i32,
-              &mut w,
-              &mut work,
-              4 * n as i32,
-              &mut info);
+        let info = $syev(layout, b'V', b'U', n as i32, &mut a, n as i32, &mut w);
         if info == 0 {
             Ok((w, a))
         } else {
diff --git a/src/error.rs b/src/error.rs
@@ -2,6 +2,7 @@
 
 use std::error;
 use std::fmt;
+use ndarray::Ixs;
 
 #[derive(Debug)]
 pub struct LapackError {
@@ -44,17 +45,37 @@ impl error::Error for NotSquareError {
     }
 }
 
+#[derive(Debug)]
+pub struct StrideError {
+    pub s0: Ixs,
+    pub s1: Ixs,
+}
+
+impl fmt::Display for StrideError {
+    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
+        write!(f, "invalid stride: s0={}, s1={}", self.s0, self.s1)
+    }
+}
+
+impl error::Error for StrideError {
+    fn description(&self) -> &str {
+        "invalid stride"
+    }
+}
+
 #[derive(Debug)]
 pub enum LinalgError {
     NotSquare(NotSquareError),
     Lapack(LapackError),
+    Stride(StrideError),
 }
 
 impl fmt::Display for LinalgError {
     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),
         }
     }
 }
@@ -64,6 +85,7 @@ impl error::Error for LinalgError {
         match *self {
             LinalgError::NotSquare(ref err) => err.description(),
             LinalgError::Lapack(ref err) => err.description(),
+            LinalgError::Stride(ref err) => err.description(),
         }
     }
 }
@@ -79,3 +101,9 @@ impl From<LapackError> for LinalgError {
         LinalgError::Lapack(err)
     }
 }
+
+impl From<StrideError> for LinalgError {
+    fn from(err: StrideError) -> LinalgError {
+        LinalgError::Stride(err)
+    }
+}
diff --git a/src/hermite.rs b/src/hermite.rs
@@ -29,16 +29,20 @@ impl<A> HermiteMatrix for Array<A, Ix2>
     where A: ImplQR + ImplSVD + ImplNorm + ImplSolve + ImplEigh + ImplCholesky + LinalgScalar + Float + Debug
 {
     fn eigh(self) -> Result<(Self::Vector, Self), LinalgError> {
-        try!(self.check_square());
+        self.check_square()?;
+        let layout = self.layout()?;
         let (rows, cols) = self.size();
-        let (w, a) = try!(ImplEigh::eigh(rows, self.into_raw_vec()));
+        let (w, a) = ImplEigh::eigh(layout, rows, self.into_raw_vec())?;
         let ea = Array::from_vec(w);
-        let va = Array::from_vec(a).into_shape((rows, cols)).unwrap().reversed_axes();
+        let va = match layout {
+            Layout::ColumnMajor => Array::from_vec(a).into_shape((rows, cols)).unwrap().reversed_axes(),
+            Layout::RowMajor => Array::from_vec(a).into_shape((rows, cols)).unwrap(),
+        };
         Ok((ea, va))
     }
     fn ssqrt(self) -> Result<Self, LinalgError> {
         let (n, _) = self.size();
-        let (e, v) = try!(self.eigh());
+        let (e, v) = self.eigh()?;
         let mut res = Array::zeros((n, n));
         for i in 0..n {
             for j in 0..n {
@@ -48,11 +52,10 @@ impl<A> HermiteMatrix for Array<A, Ix2>
         Ok(v.dot(&res))
     }
     fn cholesky(self) -> Result<Self, LinalgError> {
-        try!(self.check_square());
-        println!("layout = {:?}", self.layout());
+        self.check_square()?;
         let (n, _) = self.size();
-        let layout = self.layout();
-        let a = try!(ImplCholesky::cholesky(layout, n, self.into_raw_vec()));
+        let layout = self.layout()?;
+        let a = ImplCholesky::cholesky(layout, n, self.into_raw_vec())?;
         let mut c = match layout {
             Layout::RowMajor => Array::from_vec(a).into_shape((n, n)).unwrap(),
             Layout::ColumnMajor => Array::from_vec(a).into_shape((n, n)).unwrap().reversed_axes(),
diff --git a/src/matrix.rs b/src/matrix.rs
@@ -6,7 +6,7 @@ use ndarray::prelude::*;
 use ndarray::LinalgScalar;
 use lapack::c::Layout;
 
-use error::LapackError;
+use error::{LinalgError, StrideError};
 use qr::ImplQR;
 use svd::ImplSVD;
 use norm::ImplNorm;
@@ -20,19 +20,19 @@ pub trait Matrix: Sized {
     /// number of (rows, columns)
     fn size(&self) -> (usize, usize);
     /// Layout (C/Fortran) of matrix
-    fn layout(&self) -> Layout;
+    fn layout(&self) -> Result<Layout, StrideError>;
     /// Operator norm for L-1 norm
     fn norm_1(&self) -> Self::Scalar;
     /// Operator norm for L-inf norm
     fn norm_i(&self) -> Self::Scalar;
     /// Frobenius norm
     fn norm_f(&self) -> Self::Scalar;
     /// singular-value decomposition (SVD)
-    fn svd(self) -> Result<(Self, Self::Vector, Self), LapackError>;
+    fn svd(self) -> Result<(Self, Self::Vector, Self), LinalgError>;
     /// QR decomposition
-    fn qr(self) -> Result<(Self, Self), LapackError>;
+    fn qr(self) -> Result<(Self, Self), LinalgError>;
     /// LU decomposition
-    fn lu(self) -> Result<(Self::Permutator, Self, Self), LapackError>;
+    fn lu(self) -> Result<(Self::Permutator, Self, Self), LinalgError>;
     /// permutate matrix (inplace)
     fn permutate(&mut self, p: &Self::Permutator);
     /// permutate matrix (outplace)
@@ -52,12 +52,18 @@ impl<A> Matrix for Array<A, Ix2>
     fn size(&self) -> (usize, usize) {
         (self.rows(), self.cols())
     }
-    fn layout(&self) -> Layout {
+    fn layout(&self) -> Result<Layout, StrideError> {
         let strides = self.strides();
+        if min(strides[0], strides[1]) != 1 {
+            return Err(StrideError {
+                s0: strides[0],
+                s1: strides[1],
+            });;
+        }
         if strides[0] < strides[1] {
-            Layout::ColumnMajor
+            Ok(Layout::ColumnMajor)
         } else {
-            Layout::RowMajor
+            Ok(Layout::RowMajor)
         }
     }
     fn norm_1(&self) -> Self::Scalar {
@@ -82,35 +88,24 @@ impl<A> Matrix for Array<A, Ix2>
         let (m, n) = self.size();
         ImplNorm::norm_f(m, n, self.clone().into_raw_vec())
     }
-    fn svd(self) -> Result<(Self, Self::Vector, Self), LapackError> {
-        let strides = self.strides();
-        let (m, n) = if strides[0] > strides[1] {
-            self.size()
-        } else {
-            let (n, m) = self.size();
-            (m, n)
-        };
-        let (u, s, vt) = try!(ImplSVD::svd(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()?;
+        let (u, s, vt) = ImplSVD::svd(layout, m, n, self.clone().into_raw_vec())?;
         let sv = Array::from_vec(s);
-        if strides[0] > strides[1] {
-            let ua = Array::from_vec(u).into_shape((n, n)).unwrap();
-            let va = Array::from_vec(vt).into_shape((m, m)).unwrap();
-            Ok((va, sv, ua))
-        } else {
-            let ua = Array::from_vec(u).into_shape((n, n)).unwrap().reversed_axes();
-            let va = Array::from_vec(vt).into_shape((m, m)).unwrap().reversed_axes();
-            Ok((ua, sv, va))
+        let ua = Array::from_vec(u).into_shape((n, n)).unwrap();
+        let va = Array::from_vec(vt).into_shape((m, m)).unwrap();
+        match layout {
+            Layout::RowMajor => Ok((ua, sv, va)),
+            Layout::ColumnMajor => Ok((ua.reversed_axes(), sv, va.reversed_axes())),
         }
     }
-    fn qr(self) -> Result<(Self, Self), LapackError> {
+    fn qr(self) -> Result<(Self, Self), LinalgError> {
         let (n, m) = self.size();
         let strides = self.strides();
         let k = min(n, m);
-        let (q, r) = if strides[0] < strides[1] {
-            try!(ImplQR::qr(m, n, self.clone().into_raw_vec()))
-        } else {
-            try!(ImplQR::lq(n, m, self.clone().into_raw_vec()))
-        };
+        let layout = self.layout()?;
+        let (q, r) = ImplQR::qr(layout, m, n, self.clone().into_raw_vec())?;
         let (qa, ra) = if strides[0] < strides[1] {
             (Array::from_vec(q).into_shape((m, n)).unwrap().reversed_axes(),
              Array::from_vec(r).into_shape((m, n)).unwrap().reversed_axes())
@@ -136,23 +131,14 @@ impl<A> Matrix for Array<A, Ix2>
         }
         Ok((qm, rm))
     }
-    fn lu(self) -> Result<(Self::Permutator, Self, Self), LapackError> {
+    fn lu(self) -> Result<(Self::Permutator, Self, Self), LinalgError> {
         let (n, m) = self.size();
-        println!("n={}, m={}", n, m);
         let k = min(n, m);
-        let (p, mut a) = match self.layout() {
-            Layout::ColumnMajor => {
-                println!("ColumnMajor");
-                let (p, l) = ImplSolve::lu(self.layout(), n, m, self.clone().into_raw_vec())?;
-                (p, Array::from_vec(l).into_shape((m, n)).unwrap().reversed_axes())
-            }
-            Layout::RowMajor => {
-                println!("RowMajor");
-                let (p, l) = ImplSolve::lu(self.layout(), n, m, self.clone().into_raw_vec())?;
-                (p, Array::from_vec(l).into_shape((n, m)).unwrap())
-            }
+        let (p, l) = ImplSolve::lu(self.layout()?, n, m, self.clone().into_raw_vec())?;
+        let mut a = match self.layout()? {
+            Layout::ColumnMajor => Array::from_vec(l).into_shape((m, n)).unwrap().reversed_axes(),
+            Layout::RowMajor => Array::from_vec(l).into_shape((n, m)).unwrap(),
         };
-        println!("a (after LU) = \n{:?}", &a);
         let mut lm = Array::zeros((n, k));
         for ((i, j), val) in lm.indexed_iter_mut() {
             if i > j {
@@ -171,7 +157,6 @@ impl<A> Matrix for Array<A, Ix2>
         } else {
             a
         };
-        println!("am = \n{:?}", am);
         Ok((p, lm, am))
     }
     fn permutate(&mut self, ipiv: &Self::Permutator) {
diff --git a/src/norm.rs b/src/norm.rs
@@ -1,7 +1,6 @@
 //! Implement Norms for matrices
 
-use lapack::fortran::*;
-use num_traits::Zero;
+use lapack::c::*;
 
 pub trait ImplNorm: Sized {
     fn norm_1(m: usize, n: usize, mut a: Vec<Self>) -> Self;
@@ -13,16 +12,13 @@ macro_rules! impl_norm {
     ($scalar:ty, $lange:path) => {
 impl ImplNorm for $scalar {
     fn norm_1(m: usize, n: usize, mut a: Vec<Self>) -> Self {
-        let mut work = Vec::<Self>::new();
-        $lange(b'o', m as i32, n as i32, &mut a, m as i32, &mut work)
+        $lange(Layout::ColumnMajor, b'o', m as i32, n as i32, &mut a, m as i32)
     }
     fn norm_i(m: usize, n: usize, mut a: Vec<Self>) -> Self {
-        let mut work = vec![Self::zero(); m];
-        $lange(b'i', m as i32, n as i32, &mut a, m as i32, &mut work)
+        $lange(Layout::ColumnMajor, b'i', m as i32, n as i32, &mut a, m as i32)
     }
     fn norm_f(m: usize, n: usize, mut a: Vec<Self>) -> Self {
-        let mut work = Vec::<Self>::new();
-        $lange(b'f', m as i32, n as i32, &mut a, m as i32, &mut work)
+        $lange(Layout::ColumnMajor, b'f', m as i32, n as i32, &mut a, m as i32)
     }
 }
 }} // end macro_rules
diff --git a/src/qr.rs b/src/qr.rs
diff --git a/src/square.rs b/src/square.rs
diff --git a/src/svd.rs b/src/svd.rs

Original file line number	Diff line number	Diff line change
`@@ -1 +1,2 @@`
`1`	`1`	`max_width = 120`
	`2`	`+use_try_shorthand = true`
Original file line number	Diff line number	Diff line change
`@@ -2,6 +2,7 @@`
`2`	`2`
`3`	`3`	`use std::error;`
`4`	`4`	`use std::fmt;`
	`5`	`+use ndarray::Ixs;`
`5`	`6`
`6`	`7`	`#[derive(Debug)]`
`7`	`8`	`pub struct LapackError {`
`@@ -44,17 +45,37 @@ impl error::Error for NotSquareError {`
`44`	`45`	`}`
`45`	`46`	`}`
`46`	`47`
	`48`	`+#[derive(Debug)]`
	`49`	`+pub struct StrideError {`
	`50`	`+ pub s0: Ixs,`
	`51`	`+ pub s1: Ixs,`
	`52`	`+}`
	`53`	`+`
	`54`	`+impl fmt::Display for StrideError {`
	`55`	`+ fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {`
	`56`	`+ write!(f, "invalid stride: s0={}, s1={}", self.s0, self.s1)`
	`57`	`+ }`
	`58`	`+}`
	`59`	`+`
	`60`	`+impl error::Error for StrideError {`
	`61`	`+ fn description(&self) -> &str {`
	`62`	`+ "invalid stride"`
	`63`	`+ }`
	`64`	`+}`
	`65`	`+`
`47`	`66`	`#[derive(Debug)]`
`48`	`67`	`pub enum LinalgError {`
`49`	`68`	`NotSquare(NotSquareError),`
`50`	`69`	`Lapack(LapackError),`
	`70`	`+ Stride(StrideError),`
`51`	`71`	`}`
`52`	`72`
`53`	`73`	`impl fmt::Display for LinalgError {`
`54`	`74`	`fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {`
`55`	`75`	`match *self {`
`56`	`76`	`LinalgError::NotSquare(ref err) => err.fmt(f),`
`57`	`77`	`LinalgError::Lapack(ref err) => err.fmt(f),`
	`78`	`+ LinalgError::Stride(ref err) => err.fmt(f),`
`58`	`79`	`}`
`59`	`80`	`}`
`60`	`81`	`}`
`@@ -64,6 +85,7 @@ impl error::Error for LinalgError {`
`64`	`85`	`match *self {`
`65`	`86`	`LinalgError::NotSquare(ref err) => err.description(),`
`66`	`87`	`LinalgError::Lapack(ref err) => err.description(),`
	`88`	`+ LinalgError::Stride(ref err) => err.description(),`
`67`	`89`	`}`
`68`	`90`	`}`
`69`	`91`	`}`
`@@ -79,3 +101,9 @@ impl From<LapackError> for LinalgError {`
`79`	`101`	`LinalgError::Lapack(err)`
`80`	`102`	`}`
`81`	`103`	`}`
	`104`	`+`
	`105`	`+impl From<StrideError> for LinalgError {`
	`106`	`+ fn from(err: StrideError) -> LinalgError {`
	`107`	`+ LinalgError::Stride(err)`
	`108`	`+ }`
	`109`	`+}`