Skip to content

Commit cb4f764

Browse files
termoshtttermoshtt
authored
Merge pull request #207 from rust-ndarray/lax
ndarray_linalg::lapack module as "lax" crate
2 parents 5cce4ac + 44c2819 commit cb4f764

32 files changed

+604
-481
lines changed

‎.gitignore

+1-1
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
# Generated by Cargo
22
# will have compiled files and executables
3-
/target/
3+
target/
44
*.rustfmt
55
rusty-tags.*
66

‎Cargo.toml

+4-1
Original file line numberDiff line numberDiff line change
@@ -1,2 +1,5 @@
11
[workspace]
2-
members = ["ndarray-linalg"]
2+
members = [
3+
"ndarray-linalg",
4+
"lax",
5+
]

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

File renamed without changes.

‎lax/Cargo.toml

+31
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,31 @@
1+
[package]
2+
name = "lax"
3+
version = "0.1.0"
4+
authors = ["Toshiki Teramura <[email protected]>"]
5+
edition = "2018"
6+
7+
[features]
8+
default = []
9+
intel-mkl = ["lapack-src/intel-mkl", "blas-src/intel-mkl"]
10+
netlib = ["lapack-src/netlib", "blas-src/netlib"]
11+
openblas = ["lapack-src/openblas", "blas-src/openblas"]
12+
13+
[dependencies]
14+
thiserror = "1"
15+
cauchy = "0.2"
16+
lapacke = "0.2.0"
17+
num-traits = "0.2"
18+
19+
[dependencies.blas-src]
20+
version = "0.6.1"
21+
default-features = false
22+
23+
[dependencies.lapack-src]
24+
version = "0.6.0"
25+
default-features = false
26+
27+
[dependencies.openblas-src]
28+
version = "0.9.0"
29+
default-features = false
30+
features = ["static"]
31+
optional = true

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

+2-1
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,8 @@
11
//! Cholesky decomposition
22
33
use super::*;
4-
use crate::{error::*, layout::MatrixLayout, types::*};
4+
use crate::{error::*, layout::MatrixLayout};
5+
use cauchy::*;
56

67
pub trait Cholesky_: Sized {
78
/// Cholesky: wrapper of `*potrf`

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

+2-2
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
//! Eigenvalue decomposition for general matrices
22
3-
use super::*;
4-
use crate::{error::*, layout::MatrixLayout, types::*};
3+
use crate::{error::*, layout::MatrixLayout};
4+
use cauchy::*;
55
use num_traits::Zero;
66

77
/// Wraps `*geev` for real/complex

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

+2-1
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,8 @@
11
//! Eigenvalue decomposition for Hermite matrices
22
33
use super::*;
4-
use crate::{error::*, layout::MatrixLayout, types::*};
4+
use crate::{error::*, layout::MatrixLayout};
5+
use cauchy::*;
56
use num_traits::Zero;
67

78
/// Wraps `*syev` for real and `*heev` for complex

‎lax/src/error.rs

+38
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,38 @@
1+
use thiserror::Error;
2+
3+
pub type Result<T> = ::std::result::Result<T, Error>;
4+
5+
#[derive(Error, Debug)]
6+
pub enum Error {
7+
#[error(
8+
"Invalid value for LAPACK subroutine {}-th argument",
9+
-return_code
10+
)]
11+
LapackInvalidValue { return_code: i32 },
12+
13+
#[error(
14+
"Comutational failure in LAPACK subroutine: return_code = {}",
15+
return_code
16+
)]
17+
LapackComputationalFailure { return_code: i32 },
18+
19+
/// Strides of the array is not supported
20+
#[error("Invalid shape")]
21+
InvalidShape,
22+
}
23+
24+
pub trait AsLapackResult {
25+
fn as_lapack_result(self) -> Result<()>;
26+
}
27+
28+
impl AsLapackResult for i32 {
29+
fn as_lapack_result(self) -> Result<()> {
30+
if self > 0 {
31+
return Err(Error::LapackComputationalFailure { return_code: self });
32+
}
33+
if self < 0 {
34+
return Err(Error::LapackInvalidValue { return_code: self });
35+
}
36+
Ok(())
37+
}
38+
}

‎lax/src/layout.rs

+66
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,66 @@
1+
//! Memory layout of matrices
2+
3+
pub type LDA = i32;
4+
pub type LEN = i32;
5+
pub type Col = i32;
6+
pub type Row = i32;
7+
8+
#[derive(Debug, Clone, Copy, PartialEq)]
9+
pub enum MatrixLayout {
10+
C((Row, LDA)),
11+
F((Col, LDA)),
12+
}
13+
14+
impl MatrixLayout {
15+
pub fn size(&self) -> (Row, Col) {
16+
match *self {
17+
MatrixLayout::C((row, lda)) => (row, lda),
18+
MatrixLayout::F((col, lda)) => (lda, col),
19+
}
20+
}
21+
22+
pub fn resized(&self, row: Row, col: Col) -> MatrixLayout {
23+
match *self {
24+
MatrixLayout::C(_) => MatrixLayout::C((row, col)),
25+
MatrixLayout::F(_) => MatrixLayout::F((col, row)),
26+
}
27+
}
28+
29+
pub fn lda(&self) -> LDA {
30+
std::cmp::max(
31+
1,
32+
match *self {
33+
MatrixLayout::C((_, lda)) | MatrixLayout::F((_, lda)) => lda,
34+
},
35+
)
36+
}
37+
38+
pub fn len(&self) -> LEN {
39+
match *self {
40+
MatrixLayout::C((row, _)) => row,
41+
MatrixLayout::F((col, _)) => col,
42+
}
43+
}
44+
45+
pub fn is_empty(&self) -> bool {
46+
self.len() == 0
47+
}
48+
49+
pub fn lapacke_layout(&self) -> lapacke::Layout {
50+
match *self {
51+
MatrixLayout::C(_) => lapacke::Layout::RowMajor,
52+
MatrixLayout::F(_) => lapacke::Layout::ColumnMajor,
53+
}
54+
}
55+
56+
pub fn same_order(&self, other: &MatrixLayout) -> bool {
57+
self.lapacke_layout() == other.lapacke_layout()
58+
}
59+
60+
pub fn toggle_order(&self) -> Self {
61+
match *self {
62+
MatrixLayout::C((row, col)) => MatrixLayout::F((col, row)),
63+
MatrixLayout::F((col, row)) => MatrixLayout::C((row, col)),
64+
}
65+
}
66+
}

‎ndarray-linalg/src/lapack/least_squares.rs renamed to ‎lax/src/least_squares.rs

+4-9
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,7 @@
11
//! Least squares
22
3-
use super::*;
4-
use crate::{error::*, layout::MatrixLayout, types::*};
5-
use ndarray::{ErrorKind, ShapeError};
3+
use crate::{error::*, layout::MatrixLayout};
4+
use cauchy::*;
65
use num_traits::Zero;
76

87
/// Result of LeastSquares
@@ -39,9 +38,7 @@ macro_rules! impl_least_squares {
3938
) -> Result<LeastSquaresOutput<Self>> {
4039
let (m, n) = a_layout.size();
4140
if (m as usize) > b.len() || (n as usize) > b.len() {
42-
return Err(LinalgError::Shape(ShapeError::from_kind(
43-
ErrorKind::IncompatibleShape,
44-
)));
41+
return Err(Error::InvalidShape);
4542
}
4643
let k = ::std::cmp::min(m, n);
4744
let nrhs = 1;
@@ -83,9 +80,7 @@ macro_rules! impl_least_squares {
8380
|| (n as usize) > b.len()
8481
|| a_layout.lapacke_layout() != b_layout.lapacke_layout()
8582
{
86-
return Err(LinalgError::Shape(ShapeError::from_kind(
87-
ErrorKind::IncompatibleShape,
88-
)));
83+
return Err(Error::InvalidShape);
8984
}
9085
let k = ::std::cmp::min(m, n);
9186
let nrhs = b_layout.size().1;

‎lax/src/lib.rs

+150
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,150 @@
1+
//! Linear Algebra eXtension (LAX)
2+
//! ===============================
3+
//!
4+
//! ndarray-free safe Rust wrapper for LAPACK FFI
5+
//!
6+
//! Linear equation, Inverse matrix, Condition number
7+
//! --------------------------------------------------
8+
//!
9+
//! As the property of $A$, several types of triangular factorization are used:
10+
//!
11+
//! - LU-decomposition for general matrix
12+
//! - $PA = LU$, where $L$ is lower matrix, $U$ is upper matrix, and $P$ is permutation matrix
13+
//! - Bunch-Kaufman diagonal pivoting method for nonpositive-definite Hermitian matrix
14+
//! - $A = U D U^\dagger$, where $U$ is upper matrix,
15+
//! $D$ is Hermitian and block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
16+
//!
17+
//! | matrix type | Triangler factorization (TRF) | Solve (TRS) | Inverse matrix (TRI) | Reciprocal condition number (CON) |
18+
//! |:--------------------------------|:------------------------------|:------------|:---------------------|:----------------------------------|
19+
//! | General (GE) | [lu] | [solve] | [inv] | [rcond] |
20+
//! | Symmetric (SY) / Hermitian (HE) | [bk] | [solveh] | [invh] | - |
21+
//!
22+
//! [lu]: solve/trait.Solve_.html#tymethod.lu
23+
//! [solve]: solve/trait.Solve_.html#tymethod.solve
24+
//! [inv]: solve/trait.Solve_.html#tymethod.inv
25+
//! [rcond]: solve/trait.Solve_.html#tymethod.rcond
26+
//!
27+
//! [bk]: solveh/trait.Solveh_.html#tymethod.bk
28+
//! [solveh]: solveh/trait.Solveh_.html#tymethod.solveh
29+
//! [invh]: solveh/trait.Solveh_.html#tymethod.invh
30+
//!
31+
//! Eigenvalue Problem
32+
//! -------------------
33+
//!
34+
//! Solve eigenvalue problem for a matrix $A$
35+
//!
36+
//! $$ Av_i = \lambda_i v_i $$
37+
//!
38+
//! or generalized eigenvalue problem
39+
//!
40+
//! $$ Av_i = \lambda_i B v_i $$
41+
//!
42+
//! | matrix type | Eigenvalue (EV) | Generalized Eigenvalue Problem (EG) |
43+
//! |:--------------------------------|:----------------|:------------------------------------|
44+
//! | General (GE) |[eig] | - |
45+
//! | Symmetric (SY) / Hermitian (HE) |[eigh] |[eigh_generalized] |
46+
//!
47+
//! [eig]: eig/trait.Eig_.html#tymethod.eig
48+
//! [eigh]: eigh/trait.Eigh_.html#tymethod.eigh
49+
//! [eigh_generalized]: eigh/trait.Eigh_.html#tymethod.eigh_generalized
50+
//!
51+
//! Singular Value Decomposition (SVD), Least square problem
52+
//! ----------------------------------------------------------
53+
//!
54+
//! | matrix type | Singular Value Decomposition (SVD) | SVD with divided-and-conquer (SDD) | Least square problem (LSD) |
55+
//! |:-------------|:-----------------------------------|:-----------------------------------|:---------------------------|
56+
//! | General (GE) | [svd] | [svddc] | [least_squares] |
57+
//!
58+
//! [svd]: svd/trait.SVD_.html#tymethod.svd
59+
//! [svddc]: svddck/trait.SVDDC_.html#tymethod.svddc
60+
//! [least_squares]: least_squares/trait.LeastSquaresSvdDivideConquer_.html#tymethod.least_squares
61+
62+
extern crate blas_src;
63+
extern crate lapack_src;
64+
65+
pub mod cholesky;
66+
pub mod eig;
67+
pub mod eigh;
68+
pub mod error;
69+
pub mod layout;
70+
pub mod least_squares;
71+
pub mod opnorm;
72+
pub mod qr;
73+
pub mod solve;
74+
pub mod solveh;
75+
pub mod svd;
76+
pub mod svddc;
77+
pub mod triangular;
78+
pub mod tridiagonal;
79+
80+
pub use self::cholesky::*;
81+
pub use self::eig::*;
82+
pub use self::eigh::*;
83+
pub use self::least_squares::*;
84+
pub use self::opnorm::*;
85+
pub use self::qr::*;
86+
pub use self::solve::*;
87+
pub use self::solveh::*;
88+
pub use self::svd::*;
89+
pub use self::svddc::*;
90+
pub use self::triangular::*;
91+
pub use self::tridiagonal::*;
92+
93+
use cauchy::*;
94+
95+
pub type Pivot = Vec<i32>;
96+
97+
/// Trait for primitive types which implements LAPACK subroutines
98+
pub trait Lapack:
99+
OperatorNorm_
100+
+ QR_
101+
+ SVD_
102+
+ SVDDC_
103+
+ Solve_
104+
+ Solveh_
105+
+ Cholesky_
106+
+ Eig_
107+
+ Eigh_
108+
+ Triangular_
109+
+ Tridiagonal_
110+
{
111+
}
112+
113+
impl Lapack for f32 {}
114+
impl Lapack for f64 {}
115+
impl Lapack for c32 {}
116+
impl Lapack for c64 {}
117+
118+
/// Upper/Lower specification for seveal usages
119+
#[derive(Debug, Clone, Copy)]
120+
#[repr(u8)]
121+
pub enum UPLO {
122+
Upper = b'U',
123+
Lower = b'L',
124+
}
125+
126+
#[derive(Debug, Clone, Copy)]
127+
#[repr(u8)]
128+
pub enum Transpose {
129+
No = b'N',
130+
Transpose = b'T',
131+
Hermite = b'C',
132+
}
133+
134+
#[derive(Debug, Clone, Copy)]
135+
#[repr(u8)]
136+
pub enum NormType {
137+
One = b'O',
138+
Infinity = b'I',
139+
Frobenius = b'F',
140+
}
141+
142+
impl NormType {
143+
pub(crate) fn transpose(self) -> Self {
144+
match self {
145+
NormType::One => NormType::Infinity,
146+
NormType::Infinity => NormType::One,
147+
NormType::Frobenius => NormType::Frobenius,
148+
}
149+
}
150+
}

‎ndarray-linalg/src/lapack/opnorm.rs renamed to ‎lax/src/opnorm.rs

+2-3
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,8 @@
11
//! Operator norms of matrices
22
3-
use lapacke::Layout::ColumnMajor as cm;
4-
53
use crate::layout::MatrixLayout;
6-
use crate::types::*;
4+
use cauchy::*;
5+
use lapacke::Layout::ColumnMajor as cm;
76

87
pub use super::NormType;
98

‎ndarray-linalg/src/lapack/qr.rs renamed to ‎lax/src/qr.rs

+2-2
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
//! QR decomposition
22
3-
use super::*;
4-
use crate::{error::*, layout::MatrixLayout, types::*};
3+
use crate::{error::*, layout::MatrixLayout};
4+
use cauchy::*;
55
use num_traits::Zero;
66
use std::cmp::min;
77

‎ndarray-linalg/src/lapack/solve.rs renamed to ‎lax/src/solve.rs

+2-1
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,8 @@
11
//! Solve linear problem using LU decomposition
22
33
use super::*;
4-
use crate::{error::*, layout::MatrixLayout, types::*};
4+
use crate::{error::*, layout::MatrixLayout};
5+
use cauchy::*;
56
use num_traits::Zero;
67

78
/// Wraps `*getrf`, `*getri`, and `*getrs`

0 commit comments

Comments
 (0)