ndarray_linalg::lapack module as "lax" crate

.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",
+]

lax/Cargo.toml

@@ -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

ndarray-linalg/src/lapack/cholesky.rs → lax/src/cholesky.rs

@@ -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/eig.rs → lax/src/eig.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/eigh.rs → lax/src/eigh.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

lax/src/error.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/layout.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)),
+        }
+    }
+}

ndarray-linalg/src/lapack/least_squares.rs → lax/src/least_squares.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;

lax/src/lib.rs

@@ -0,0 +1,150 @@
+//! Linear Algebra eXtension (LAX)
+//! ===============================
+//!
+//! ndarray-free safe Rust wrapper for LAPACK FFI
+//!
+//! Linear equation, Inverse matrix, Condition number
+//! --------------------------------------------------
+//!
+//! As the property of $A$, several types of triangular factorization are used:
+//!
+//! - LU-decomposition for general matrix
+//!   - $PA = LU$, where $L$ is lower matrix, $U$ is upper matrix, and $P$ is permutation matrix
+//! - Bunch-Kaufman diagonal pivoting method for nonpositive-definite Hermitian matrix
+//!   - $A = U D U^\dagger$, where $U$ is upper matrix,
+//!     $D$ is Hermitian and block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
+//!
+//! | matrix type                     | Triangler factorization (TRF) | Solve (TRS) | Inverse matrix (TRI) | Reciprocal condition number (CON) |
+//! |:--------------------------------|:------------------------------|:------------|:---------------------|:----------------------------------|
+//! | General (GE)                    | [lu]                          | [solve]     | [inv]                | [rcond]                           |
+//! | Symmetric (SY) / Hermitian (HE) | [bk]                          | [solveh]    | [invh]               | -                                 |
+//!
+//! [lu]:    solve/trait.Solve_.html#tymethod.lu
+//! [solve]: solve/trait.Solve_.html#tymethod.solve
+//! [inv]:   solve/trait.Solve_.html#tymethod.inv
+//! [rcond]: solve/trait.Solve_.html#tymethod.rcond
+//!
+//! [bk]:     solveh/trait.Solveh_.html#tymethod.bk
+//! [solveh]: solveh/trait.Solveh_.html#tymethod.solveh
+//! [invh]:   solveh/trait.Solveh_.html#tymethod.invh
+//!
+//! Eigenvalue Problem
+//! -------------------
+//!
+//! Solve eigenvalue problem for a matrix $A$
+//!
+//! $$ Av_i = \lambda_i v_i $$
+//!
+//! or generalized eigenvalue problem
+//!
+//! $$ Av_i = \lambda_i B v_i $$
+//!
+//! | matrix type                     | Eigenvalue (EV) | Generalized Eigenvalue Problem (EG) |
+//! |:--------------------------------|:----------------|:------------------------------------|
+//! | General (GE)                    |[eig]            | -                                   |
+//! | Symmetric (SY) / Hermitian (HE) |[eigh]           |[eigh_generalized]                   |
+//!
+//! [eig]:              eig/trait.Eig_.html#tymethod.eig
+//! [eigh]:             eigh/trait.Eigh_.html#tymethod.eigh
+//! [eigh_generalized]: eigh/trait.Eigh_.html#tymethod.eigh_generalized
+//!
+//! Singular Value Decomposition (SVD), Least square problem
+//! ----------------------------------------------------------
+//!
+//! | matrix type  | Singular Value Decomposition (SVD) | SVD with divided-and-conquer (SDD) | Least square problem (LSD) |
+//! |:-------------|:-----------------------------------|:-----------------------------------|:---------------------------|
+//! | General (GE) | [svd]                              | [svddc]                            | [least_squares]            |
+//!
+//! [svd]:   svd/trait.SVD_.html#tymethod.svd
+//! [svddc]: svddck/trait.SVDDC_.html#tymethod.svddc
+//! [least_squares]: least_squares/trait.LeastSquaresSvdDivideConquer_.html#tymethod.least_squares
+
+extern crate blas_src;
+extern crate lapack_src;
+
+pub mod cholesky;
+pub mod eig;
+pub mod eigh;
+pub mod error;
+pub mod layout;
+pub mod least_squares;
+pub mod opnorm;
+pub mod qr;
+pub mod solve;
+pub mod solveh;
+pub mod svd;
+pub mod svddc;
+pub mod triangular;
+pub mod tridiagonal;
+
+pub use self::cholesky::*;
+pub use self::eig::*;
+pub use self::eigh::*;
+pub use self::least_squares::*;
+pub use self::opnorm::*;
+pub use self::qr::*;
+pub use self::solve::*;
+pub use self::solveh::*;
+pub use self::svd::*;
+pub use self::svddc::*;
+pub use self::triangular::*;
+pub use self::tridiagonal::*;
+
+use cauchy::*;
+
+pub type Pivot = Vec<i32>;
+
+/// Trait for primitive types which implements LAPACK subroutines
+pub trait Lapack:
+    OperatorNorm_
+    + QR_
+    + SVD_
+    + SVDDC_
+    + Solve_
+    + Solveh_
+    + Cholesky_
+    + Eig_
+    + Eigh_
+    + Triangular_
+    + Tridiagonal_
+{
+}
+
+impl Lapack for f32 {}
+impl Lapack for f64 {}
+impl Lapack for c32 {}
+impl Lapack for c64 {}
+
+/// Upper/Lower specification for seveal usages
+#[derive(Debug, Clone, Copy)]
+#[repr(u8)]
+pub enum UPLO {
+    Upper = b'U',
+    Lower = b'L',
+}
+
+#[derive(Debug, Clone, Copy)]
+#[repr(u8)]
+pub enum Transpose {
+    No = b'N',
+    Transpose = b'T',
+    Hermite = b'C',
+}
+
+#[derive(Debug, Clone, Copy)]
+#[repr(u8)]
+pub enum NormType {
+    One = b'O',
+    Infinity = b'I',
+    Frobenius = b'F',
+}
+
+impl NormType {
+    pub(crate) fn transpose(self) -> Self {
+        match self {
+            NormType::One => NormType::Infinity,
+            NormType::Infinity => NormType::One,
+            NormType::Frobenius => NormType::Frobenius,
+        }
+    }
+}

ndarray-linalg/src/lapack/opnorm.rs → lax/src/opnorm.rs

@@ -1,9 +1,8 @@
 //! Operator norms of matrices
 
-use lapacke::Layout::ColumnMajor as cm;
-
 use crate::layout::MatrixLayout;
-use crate::types::*;
+use cauchy::*;
+use lapacke::Layout::ColumnMajor as cm;
 
 pub use super::NormType;
 

ndarray-linalg/src/lapack/qr.rs → lax/src/qr.rs

@@ -1,7 +1,7 @@
 //! QR decomposition
 
-use super::*;
-use crate::{error::*, layout::MatrixLayout, types::*};
+use crate::{error::*, layout::MatrixLayout};
+use cauchy::*;
 use num_traits::Zero;
 use std::cmp::min;
 

ndarray-linalg/src/lapack/solve.rs → lax/src/solve.rs

@@ -1,7 +1,8 @@
 //! Solve linear problem using LU decomposition
 
 use super::*;
-use crate::{error::*, layout::MatrixLayout, types::*};
+use crate::{error::*, layout::MatrixLayout};
+use cauchy::*;
 use num_traits::Zero;
 
 /// Wraps `*getrf`, `*getri`, and `*getrs`