Merge pull request #207 from rust-ndarray/lax

ndarray_linalg::lapack module as "lax" crate

Author: termoshtt

.gitignore

@@ -1,6 +1,6 @@
 # Generated by Cargo
 # will have compiled files and executables
-/target/
+target/
 *.rustfmt
 rusty-tags.*
 

Cargo.toml

@@ -1,2 +1,5 @@
 [workspace]
-members = ["ndarray-linalg"]
+members = [
+  "ndarray-linalg",
+  "lax",
+]

ndarray-linalg/katex-header.html katex-header.htmlndarray-linalg/katex-header.html renamed to katex-header.html

@@ -0,0 +1,31 @@
+[package]
+name = "lax"
+version = "0.1.0"
+authors = ["Toshiki Teramura <toshiki.teramura@gmail.com>"]
+edition = "2018"
+
+[features]
+default    = []
+intel-mkl  = ["lapack-src/intel-mkl", "blas-src/intel-mkl"]
+netlib     = ["lapack-src/netlib", "blas-src/netlib"]
+openblas   = ["lapack-src/openblas", "blas-src/openblas"]
+
+[dependencies]
+thiserror = "1"
+cauchy = "0.2"
+lapacke = "0.2.0"
+num-traits = "0.2"
+
+[dependencies.blas-src]
+version = "0.6.1"
+default-features = false
+
+[dependencies.lapack-src]
+version = "0.6.0"
+default-features = false
+
+[dependencies.openblas-src]
+version = "0.9.0"
+default-features = false
+features = ["static"]
+optional = true

lax/Cargo.toml

@@ -1,7 +1,8 @@
 //! Cholesky decomposition
 
 use super::*;
-use crate::{error::*, layout::MatrixLayout, types::*};
+use crate::{error::*, layout::MatrixLayout};
+use cauchy::*;
 
 pub trait Cholesky_: Sized {
     /// Cholesky: wrapper of `*potrf`

ndarray-linalg/src/lapack/cholesky.rs lax/src/cholesky.rsndarray-linalg/src/lapack/cholesky.rs renamed to lax/src/cholesky.rs

@@ -1,7 +1,7 @@
 //! Eigenvalue decomposition for general matrices
 
-use super::*;
-use crate::{error::*, layout::MatrixLayout, types::*};
+use crate::{error::*, layout::MatrixLayout};
+use cauchy::*;
 use num_traits::Zero;
 
 /// Wraps `*geev` for real/complex

ndarray-linalg/src/lapack/eig.rs lax/src/eig.rsndarray-linalg/src/lapack/eig.rs renamed to lax/src/eig.rs

@@ -1,7 +1,8 @@
 //! Eigenvalue decomposition for Hermite matrices
 
 use super::*;
-use crate::{error::*, layout::MatrixLayout, types::*};
+use crate::{error::*, layout::MatrixLayout};
+use cauchy::*;
 use num_traits::Zero;
 
 /// Wraps `*syev` for real and `*heev` for complex

ndarray-linalg/src/lapack/eigh.rs lax/src/eigh.rsndarray-linalg/src/lapack/eigh.rs renamed to lax/src/eigh.rs

@@ -0,0 +1,38 @@
+use thiserror::Error;
+
+pub type Result<T> = ::std::result::Result<T, Error>;
+
+#[derive(Error, Debug)]
+pub enum Error {
+    #[error(
+        "Invalid value for LAPACK subroutine {}-th argument",
+        -return_code
+    )]
+    LapackInvalidValue { return_code: i32 },
+
+    #[error(
+        "Comutational failure in LAPACK subroutine: return_code = {}",
+        return_code
+    )]
+    LapackComputationalFailure { return_code: i32 },
+
+    /// Strides of the array is not supported
+    #[error("Invalid shape")]
+    InvalidShape,
+}
+
+pub trait AsLapackResult {
+    fn as_lapack_result(self) -> Result<()>;
+}
+
+impl AsLapackResult for i32 {
+    fn as_lapack_result(self) -> Result<()> {
+        if self > 0 {
+            return Err(Error::LapackComputationalFailure { return_code: self });
+        }
+        if self < 0 {
+            return Err(Error::LapackInvalidValue { return_code: self });
+        }
+        Ok(())
+    }
+}

lax/src/error.rs

@@ -0,0 +1,66 @@
+//! Memory layout of matrices
+
+pub type LDA = i32;
+pub type LEN = i32;
+pub type Col = i32;
+pub type Row = i32;
+
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub enum MatrixLayout {
+    C((Row, LDA)),
+    F((Col, LDA)),
+}
+
+impl MatrixLayout {
+    pub fn size(&self) -> (Row, Col) {
+        match *self {
+            MatrixLayout::C((row, lda)) => (row, lda),
+            MatrixLayout::F((col, lda)) => (lda, col),
+        }
+    }
+
+    pub fn resized(&self, row: Row, col: Col) -> MatrixLayout {
+        match *self {
+            MatrixLayout::C(_) => MatrixLayout::C((row, col)),
+            MatrixLayout::F(_) => MatrixLayout::F((col, row)),
+        }
+    }
+
+    pub fn lda(&self) -> LDA {
+        std::cmp::max(
+            1,
+            match *self {
+                MatrixLayout::C((_, lda)) | MatrixLayout::F((_, lda)) => lda,
+            },
+        )
+    }
+
+    pub fn len(&self) -> LEN {
+        match *self {
+            MatrixLayout::C((row, _)) => row,
+            MatrixLayout::F((col, _)) => col,
+        }
+    }
+
+    pub fn is_empty(&self) -> bool {
+        self.len() == 0
+    }
+
+    pub fn lapacke_layout(&self) -> lapacke::Layout {
+        match *self {
+            MatrixLayout::C(_) => lapacke::Layout::RowMajor,
+            MatrixLayout::F(_) => lapacke::Layout::ColumnMajor,
+        }
+    }
+
+    pub fn same_order(&self, other: &MatrixLayout) -> bool {
+        self.lapacke_layout() == other.lapacke_layout()
+    }
+
+    pub fn toggle_order(&self) -> Self {
+        match *self {
+            MatrixLayout::C((row, col)) => MatrixLayout::F((col, row)),
+            MatrixLayout::F((col, row)) => MatrixLayout::C((row, col)),
+        }
+    }
+}

lax/src/layout.rs

@@ -1,8 +1,7 @@
 //! Least squares
 
-use super::*;
-use crate::{error::*, layout::MatrixLayout, types::*};
-use ndarray::{ErrorKind, ShapeError};
+use crate::{error::*, layout::MatrixLayout};
+use cauchy::*;
 use num_traits::Zero;
 
 /// Result of LeastSquares
@@ -39,9 +38,7 @@ macro_rules! impl_least_squares {
             ) -> Result<LeastSquaresOutput<Self>> {
                 let (m, n) = a_layout.size();
                 if (m as usize) > b.len() || (n as usize) > b.len() {
-                    return Err(LinalgError::Shape(ShapeError::from_kind(
-                        ErrorKind::IncompatibleShape,
-                    )));
+                    return Err(Error::InvalidShape);
                 }
                 let k = ::std::cmp::min(m, n);
                 let nrhs = 1;
@@ -83,9 +80,7 @@ macro_rules! impl_least_squares {
                     || (n as usize) > b.len()
                     || a_layout.lapacke_layout() != b_layout.lapacke_layout()
                 {
-                    return Err(LinalgError::Shape(ShapeError::from_kind(
-                        ErrorKind::IncompatibleShape,
-                    )));
+                    return Err(Error::InvalidShape);
                 }
                 let k = ::std::cmp::min(m, n);
                 let nrhs = b_layout.size().1;