diff --git a/.gitignore b/.gitignore index 686005d3..c90915c2 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,6 @@ # Generated by Cargo # will have compiled files and executables -/target/ +target/ *.rustfmt rusty-tags.* diff --git a/Cargo.toml b/Cargo.toml index 546f4980..ba1ae403 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,2 +1,5 @@ [workspace] -members = ["ndarray-linalg"] +members = [ + "ndarray-linalg", + "lax", +] diff --git a/ndarray-linalg/katex-header.html b/katex-header.html similarity index 100% rename from ndarray-linalg/katex-header.html rename to katex-header.html diff --git a/lax/Cargo.toml b/lax/Cargo.toml new file mode 100644 index 00000000..e42c8d74 --- /dev/null +++ b/lax/Cargo.toml @@ -0,0 +1,31 @@ +[package] +name = "lax" +version = "0.1.0" +authors = ["Toshiki Teramura "] +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 diff --git a/ndarray-linalg/src/lapack/cholesky.rs b/lax/src/cholesky.rs similarity index 97% rename from ndarray-linalg/src/lapack/cholesky.rs rename to lax/src/cholesky.rs index 80a6ab29..ef9473b4 100644 --- a/ndarray-linalg/src/lapack/cholesky.rs +++ b/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` diff --git a/ndarray-linalg/src/lapack/eig.rs b/lax/src/eig.rs similarity index 98% rename from ndarray-linalg/src/lapack/eig.rs rename to lax/src/eig.rs index 4e80cb9e..e8f9a8f8 100644 --- a/ndarray-linalg/src/lapack/eig.rs +++ b/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 diff --git a/ndarray-linalg/src/lapack/eigh.rs b/lax/src/eigh.rs similarity index 97% rename from ndarray-linalg/src/lapack/eigh.rs rename to lax/src/eigh.rs index 509c08a4..a992a67e 100644 --- a/ndarray-linalg/src/lapack/eigh.rs +++ b/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 diff --git a/lax/src/error.rs b/lax/src/error.rs new file mode 100644 index 00000000..fb4b9838 --- /dev/null +++ b/lax/src/error.rs @@ -0,0 +1,38 @@ +use thiserror::Error; + +pub type Result = ::std::result::Result; + +#[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(()) + } +} diff --git a/lax/src/layout.rs b/lax/src/layout.rs new file mode 100644 index 00000000..13caca2c --- /dev/null +++ b/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)), + } + } +} diff --git a/ndarray-linalg/src/lapack/least_squares.rs b/lax/src/least_squares.rs similarity index 89% rename from ndarray-linalg/src/lapack/least_squares.rs rename to lax/src/least_squares.rs index 28e7980e..eaf496f2 100644 --- a/ndarray-linalg/src/lapack/least_squares.rs +++ b/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> { 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; diff --git a/lax/src/lib.rs b/lax/src/lib.rs new file mode 100644 index 00000000..769910db --- /dev/null +++ b/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; + +/// 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, + } + } +} diff --git a/ndarray-linalg/src/lapack/opnorm.rs b/lax/src/opnorm.rs similarity index 97% rename from ndarray-linalg/src/lapack/opnorm.rs rename to lax/src/opnorm.rs index 55b377d5..305629c1 100644 --- a/ndarray-linalg/src/lapack/opnorm.rs +++ b/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; diff --git a/ndarray-linalg/src/lapack/qr.rs b/lax/src/qr.rs similarity index 95% rename from ndarray-linalg/src/lapack/qr.rs rename to lax/src/qr.rs index 4afb8f43..6c26273d 100644 --- a/ndarray-linalg/src/lapack/qr.rs +++ b/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; diff --git a/ndarray-linalg/src/lapack/solve.rs b/lax/src/solve.rs similarity index 98% rename from ndarray-linalg/src/lapack/solve.rs rename to lax/src/solve.rs index 89762cf1..67af6409 100644 --- a/ndarray-linalg/src/lapack/solve.rs +++ b/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` diff --git a/ndarray-linalg/src/lapack/solveh.rs b/lax/src/solveh.rs similarity index 97% rename from ndarray-linalg/src/lapack/solveh.rs rename to lax/src/solveh.rs index b7466455..2f851d42 100644 --- a/ndarray-linalg/src/lapack/solveh.rs +++ b/lax/src/solveh.rs @@ -3,7 +3,8 @@ //! See also [the manual of dsytrf](http://www.netlib.org/lapack/lapack-3.1.1/html/dsytrf.f.html) use super::*; -use crate::{error::*, layout::MatrixLayout, types::*}; +use crate::{error::*, layout::MatrixLayout}; +use cauchy::*; pub trait Solveh_: Sized { /// Bunch-Kaufman: wrapper of `*sytrf` and `*hetrf` diff --git a/ndarray-linalg/src/lapack/svd.rs b/lax/src/svd.rs similarity index 97% rename from ndarray-linalg/src/lapack/svd.rs rename to lax/src/svd.rs index 2bdd145a..47a9a1be 100644 --- a/ndarray-linalg/src/lapack/svd.rs +++ b/lax/src/svd.rs @@ -1,7 +1,7 @@ //! Singular-value decomposition -use super::*; -use crate::{error::*, layout::MatrixLayout, types::*}; +use crate::{error::*, layout::MatrixLayout}; +use cauchy::*; use num_traits::Zero; #[repr(u8)] diff --git a/ndarray-linalg/src/lapack/svddc.rs b/lax/src/svddc.rs similarity index 78% rename from ndarray-linalg/src/lapack/svddc.rs rename to lax/src/svddc.rs index 850d5a15..84f8394b 100644 --- a/ndarray-linalg/src/lapack/svddc.rs +++ b/lax/src/svddc.rs @@ -1,7 +1,22 @@ use super::*; -use crate::{error::*, layout::MatrixLayout, svddc::UVTFlag, types::*}; +use crate::{error::*, layout::MatrixLayout}; +use cauchy::*; use num_traits::Zero; +/// Specifies how many of the columns of *U* and rows of *V*ᵀ are computed and returned. +/// +/// For an input array of shape *m*×*n*, the following are computed: +#[derive(Clone, Copy, Eq, PartialEq)] +#[repr(u8)] +pub enum UVTFlag { + /// All *m* columns of *U* and all *n* rows of *V*ᵀ. + Full = b'A', + /// The first min(*m*,*n*) columns of *U* and the first min(*m*,*n*) rows of *V*ᵀ. + Some = b'S', + /// No columns of *U* or rows of *V*ᵀ. + None = b'N', +} + pub trait SVDDC_: Scalar { unsafe fn svddc(l: MatrixLayout, jobz: UVTFlag, a: &mut [Self]) -> Result>; } diff --git a/ndarray-linalg/src/lapack/triangular.rs b/lax/src/triangular.rs similarity index 97% rename from ndarray-linalg/src/lapack/triangular.rs rename to lax/src/triangular.rs index 98967bf1..e725fa9e 100644 --- a/ndarray-linalg/src/lapack/triangular.rs +++ b/lax/src/triangular.rs @@ -1,7 +1,8 @@ //! Implement linear solver and inverse matrix use super::*; -use crate::{error::*, layout::MatrixLayout, types::*}; +use crate::{error::*, layout::MatrixLayout}; +use cauchy::*; #[derive(Debug, Clone, Copy)] #[repr(u8)] diff --git a/lax/src/tridiagonal.rs b/lax/src/tridiagonal.rs new file mode 100644 index 00000000..4eb8ff13 --- /dev/null +++ b/lax/src/tridiagonal.rs @@ -0,0 +1,218 @@ +//! Implement linear solver using LU decomposition +//! for tridiagonal matrix + +use super::*; +use crate::{error::*, layout::MatrixLayout}; +use cauchy::*; +use num_traits::Zero; +use std::ops::{Index, IndexMut}; + +/// Represents a tridiagonal matrix as 3 one-dimensional vectors. +/// +/// ```text +/// [d0, u1, 0, ..., 0, +/// l1, d1, u2, ..., +/// 0, l2, d2, +/// ... ..., u{n-1}, +/// 0, ..., l{n-1}, d{n-1},] +/// ``` +#[derive(Clone, PartialEq)] +pub struct Tridiagonal { + /// layout of raw matrix + pub l: MatrixLayout, + /// (n-1) sub-diagonal elements of matrix. + pub dl: Vec, + /// (n) diagonal elements of matrix. + pub d: Vec, + /// (n-1) super-diagonal elements of matrix. + pub du: Vec, +} + +impl Tridiagonal { + fn opnorm_one(&self) -> A::Real { + let mut col_sum: Vec = self.d.iter().map(|val| val.abs()).collect(); + for i in 0..col_sum.len() { + if i < self.dl.len() { + col_sum[i] += self.dl[i].abs(); + } + if i > 0 { + col_sum[i] += self.du[i - 1].abs(); + } + } + let mut max = A::Real::zero(); + for &val in &col_sum { + if max < val { + max = val; + } + } + max + } +} + +/// Represents the LU factorization of a tridiagonal matrix `A` as `A = P*L*U`. +#[derive(Clone, PartialEq)] +pub struct LUFactorizedTridiagonal { + /// A tridiagonal matrix which consists of + /// - l : layout of raw matrix + /// - dl: (n-1) multipliers that define the matrix L. + /// - d : (n) diagonal elements of the upper triangular matrix U. + /// - du: (n-1) elements of the first super-diagonal of U. + pub a: Tridiagonal, + /// (n-2) elements of the second super-diagonal of U. + pub du2: Vec, + /// The pivot indices that define the permutation matrix `P`. + pub ipiv: Pivot, + + a_opnorm_one: A::Real, +} + +impl Index<(i32, i32)> for Tridiagonal { + type Output = A; + #[inline] + fn index(&self, (row, col): (i32, i32)) -> &A { + let (n, _) = self.l.size(); + assert!( + std::cmp::max(row, col) < n, + "ndarray: index {:?} is out of bounds for array of shape {}", + [row, col], + n + ); + match row - col { + 0 => &self.d[row as usize], + 1 => &self.dl[col as usize], + -1 => &self.du[row as usize], + _ => panic!( + "ndarray-linalg::tridiagonal: index {:?} is not tridiagonal element", + [row, col] + ), + } + } +} + +impl Index<[i32; 2]> for Tridiagonal { + type Output = A; + #[inline] + fn index(&self, [row, col]: [i32; 2]) -> &A { + &self[(row, col)] + } +} + +impl IndexMut<(i32, i32)> for Tridiagonal { + #[inline] + fn index_mut(&mut self, (row, col): (i32, i32)) -> &mut A { + let (n, _) = self.l.size(); + assert!( + std::cmp::max(row, col) < n, + "ndarray: index {:?} is out of bounds for array of shape {}", + [row, col], + n + ); + match row - col { + 0 => &mut self.d[row as usize], + 1 => &mut self.dl[col as usize], + -1 => &mut self.du[row as usize], + _ => panic!( + "ndarray-linalg::tridiagonal: index {:?} is not tridiagonal element", + [row, col] + ), + } + } +} + +impl IndexMut<[i32; 2]> for Tridiagonal { + #[inline] + fn index_mut(&mut self, [row, col]: [i32; 2]) -> &mut A { + &mut self[(row, col)] + } +} + +/// Wraps `*gttrf`, `*gtcon` and `*gttrs` +pub trait Tridiagonal_: Scalar + Sized { + /// Computes the LU factorization of a tridiagonal `m x n` matrix `a` using + /// partial pivoting with row interchanges. + unsafe fn lu_tridiagonal(a: Tridiagonal) -> Result>; + + unsafe fn rcond_tridiagonal(lu: &LUFactorizedTridiagonal) -> Result; + + unsafe fn solve_tridiagonal( + lu: &LUFactorizedTridiagonal, + bl: MatrixLayout, + t: Transpose, + b: &mut [Self], + ) -> Result<()>; +} + +macro_rules! impl_tridiagonal { + ($scalar:ty, $gttrf:path, $gtcon:path, $gttrs:path) => { + impl Tridiagonal_ for $scalar { + unsafe fn lu_tridiagonal( + mut a: Tridiagonal, + ) -> Result> { + let (n, _) = a.l.size(); + let mut du2 = vec![Zero::zero(); (n - 2) as usize]; + let mut ipiv = vec![0; n as usize]; + // We have to calc one-norm before LU factorization + let a_opnorm_one = a.opnorm_one(); + $gttrf(n, &mut a.dl, &mut a.d, &mut a.du, &mut du2, &mut ipiv) + .as_lapack_result()?; + Ok(LUFactorizedTridiagonal { + a, + du2, + ipiv, + a_opnorm_one, + }) + } + + unsafe fn rcond_tridiagonal(lu: &LUFactorizedTridiagonal) -> Result { + let (n, _) = lu.a.l.size(); + let ipiv = &lu.ipiv; + let mut rcond = Self::Real::zero(); + $gtcon( + NormType::One as u8, + n, + &lu.a.dl, + &lu.a.d, + &lu.a.du, + &lu.du2, + ipiv, + lu.a_opnorm_one, + &mut rcond, + ) + .as_lapack_result()?; + Ok(rcond) + } + + unsafe fn solve_tridiagonal( + lu: &LUFactorizedTridiagonal, + bl: MatrixLayout, + t: Transpose, + b: &mut [Self], + ) -> Result<()> { + let (n, _) = lu.a.l.size(); + let (_, nrhs) = bl.size(); + let ipiv = &lu.ipiv; + let ldb = bl.lda(); + $gttrs( + lu.a.l.lapacke_layout(), + t as u8, + n, + nrhs, + &lu.a.dl, + &lu.a.d, + &lu.a.du, + &lu.du2, + ipiv, + b, + ldb, + ) + .as_lapack_result()?; + Ok(()) + } + } + }; +} // impl_tridiagonal! + +impl_tridiagonal!(f64, lapacke::dgttrf, lapacke::dgtcon, lapacke::dgttrs); +impl_tridiagonal!(f32, lapacke::sgttrf, lapacke::sgtcon, lapacke::sgttrs); +impl_tridiagonal!(c64, lapacke::zgttrf, lapacke::zgtcon, lapacke::zgttrs); +impl_tridiagonal!(c32, lapacke::cgttrf, lapacke::cgtcon, lapacke::cgttrs); diff --git a/ndarray-linalg/Cargo.toml b/ndarray-linalg/Cargo.toml index 6347345c..717a6d09 100644 --- a/ndarray-linalg/Cargo.toml +++ b/ndarray-linalg/Cargo.toml @@ -13,20 +13,15 @@ readme = "README.md" categories = ["algorithms", "science"] [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"] -serde-1 = ["ndarray/serde-1", "num-complex/serde"] - -static = ["openblas-static"] -openblas-static = ["openblas", "openblas-src"] +default = [] +intel-mkl = ["lax/intel-mkl"] +netlib = ["lax/netlib"] +openblas = ["lax/openblas"] [dependencies] -lapacke = "0.2.0" -num-traits = "0.2.11" cauchy = "0.2.2" num-complex = "0.2.4" +num-traits = "0.2.11" rand = "0.5" thiserror = "1.0.20" @@ -35,19 +30,10 @@ version = "0.13.0" features = ["blas", "approx"] default-features = false -[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" +[dependencies.lax] +version = "0.1.0" +path = "../lax" default-features = false -features = ["static"] -optional = true [dev-dependencies] paste = "0.1.9" diff --git a/ndarray-linalg/src/convert.rs b/ndarray-linalg/src/convert.rs index 52ce0dd8..a4f37d59 100644 --- a/ndarray-linalg/src/convert.rs +++ b/ndarray-linalg/src/convert.rs @@ -35,7 +35,15 @@ pub fn into_matrix(l: MatrixLayout, a: Vec) -> Result where S: DataOwned, { - Ok(ArrayBase::from_shape_vec(l.as_shape(), a)?) + match l { + MatrixLayout::C((row, col)) => { + Ok(ArrayBase::from_shape_vec((row as usize, col as usize), a)?) + } + MatrixLayout::F((col, row)) => Ok(ArrayBase::from_shape_vec( + (row as usize, col as usize).f(), + a, + )?), + } } fn uninitialized(l: MatrixLayout) -> ArrayBase @@ -43,7 +51,14 @@ where A: Copy, S: DataOwned, { - unsafe { ArrayBase::uninitialized(l.as_shape()) } + match l { + MatrixLayout::C((row, col)) => unsafe { + ArrayBase::uninitialized((row as usize, col as usize)) + }, + MatrixLayout::F((col, row)) => unsafe { + ArrayBase::uninitialized((row as usize, col as usize).f()) + }, + } } pub fn replicate(a: &ArrayBase) -> ArrayBase diff --git a/ndarray-linalg/src/error.rs b/ndarray-linalg/src/error.rs index 35042b51..893a7a85 100644 --- a/ndarray-linalg/src/error.rs +++ b/ndarray-linalg/src/error.rs @@ -12,17 +12,9 @@ pub enum LinalgError { #[error("Not square: rows({}) != cols({})", rows, cols)] NotSquare { rows: i32, cols: i32 }, - #[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 }, + /// LAPACK subroutine returns non-zero code + #[error(transparent)] + Lapack(#[from] lapack::error::Error), /// Strides of the array is not supported #[error("invalid stride: s0={}, s1={}", s0, s1)] diff --git a/ndarray-linalg/src/lapack/mod.rs b/ndarray-linalg/src/lapack/mod.rs deleted file mode 100644 index 5c5728d0..00000000 --- a/ndarray-linalg/src/lapack/mod.rs +++ /dev/null @@ -1,105 +0,0 @@ -//! Define traits wrapping LAPACK routines - -#![allow(clippy::missing_safety_doc)] - -pub mod cholesky; -pub mod eig; -pub mod eigh; -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 super::error::*; -use super::types::*; - -pub type Pivot = Vec; - -/// 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 {} - -trait AsLapackResult { - fn as_lapack_result(self) -> Result<()>; -} - -impl AsLapackResult for i32 { - fn as_lapack_result(self) -> Result<()> { - if self > 0 { - return Err(LinalgError::LapackComputationalFailure { return_code: self }); - } - if self < 0 { - return Err(LinalgError::LapackInvalidValue { return_code: self }); - } - Ok(()) - } -} - -/// 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, - } - } -} diff --git a/ndarray-linalg/src/lapack/tridiagonal.rs b/ndarray-linalg/src/lapack/tridiagonal.rs deleted file mode 100644 index 417ac0f8..00000000 --- a/ndarray-linalg/src/lapack/tridiagonal.rs +++ /dev/null @@ -1,91 +0,0 @@ -//! Implement linear solver using LU decomposition -//! for tridiagonal matrix - -use super::*; -use crate::{error::*, layout::MatrixLayout, opnorm::*, tridiagonal::*, types::*}; -use num_traits::Zero; - -/// Wraps `*gttrf`, `*gtcon` and `*gttrs` -pub trait Tridiagonal_: Scalar + Sized { - /// Computes the LU factorization of a tridiagonal `m x n` matrix `a` using - /// partial pivoting with row interchanges. - unsafe fn lu_tridiagonal(a: &mut Tridiagonal) -> Result<(Vec, Self::Real, Pivot)>; - /// Estimates the the reciprocal of the condition number of the tridiagonal matrix in 1-norm. - unsafe fn rcond_tridiagonal(lu: &LUFactorizedTridiagonal) -> Result; - unsafe fn solve_tridiagonal( - lu: &LUFactorizedTridiagonal, - bl: MatrixLayout, - t: Transpose, - b: &mut [Self], - ) -> Result<()>; -} - -macro_rules! impl_tridiagonal { - ($scalar:ty, $gttrf:path, $gtcon:path, $gttrs:path) => { - impl Tridiagonal_ for $scalar { - unsafe fn lu_tridiagonal( - a: &mut Tridiagonal, - ) -> Result<(Vec, Self::Real, Pivot)> { - let (n, _) = a.l.size(); - let anom = a.opnorm_one()?; - let mut du2 = vec![Zero::zero(); (n - 2) as usize]; - let mut ipiv = vec![0; n as usize]; - $gttrf(n, &mut a.dl, &mut a.d, &mut a.du, &mut du2, &mut ipiv) - .as_lapack_result()?; - Ok((du2, anom, ipiv)) - } - - unsafe fn rcond_tridiagonal(lu: &LUFactorizedTridiagonal) -> Result { - let (n, _) = lu.a.l.size(); - let ipiv = &lu.ipiv; - let anorm = lu.anom; - let mut rcond = Self::Real::zero(); - $gtcon( - NormType::One as u8, - n, - &lu.a.dl, - &lu.a.d, - &lu.a.du, - &lu.du2, - ipiv, - anorm, - &mut rcond, - ) - .as_lapack_result()?; - Ok(rcond) - } - - unsafe fn solve_tridiagonal( - lu: &LUFactorizedTridiagonal, - bl: MatrixLayout, - t: Transpose, - b: &mut [Self], - ) -> Result<()> { - let (n, _) = lu.a.l.size(); - let (_, nrhs) = bl.size(); - let ipiv = &lu.ipiv; - let ldb = bl.lda(); - $gttrs( - lu.a.l.lapacke_layout(), - t as u8, - n, - nrhs, - &lu.a.dl, - &lu.a.d, - &lu.a.du, - &lu.du2, - ipiv, - b, - ldb, - ) - .as_lapack_result()?; - Ok(()) - } - } - }; -} // impl_tridiagonal! - -impl_tridiagonal!(f64, lapacke::dgttrf, lapacke::dgtcon, lapacke::dgttrs); -impl_tridiagonal!(f32, lapacke::sgttrf, lapacke::sgtcon, lapacke::sgttrs); -impl_tridiagonal!(c64, lapacke::zgttrf, lapacke::zgtcon, lapacke::zgttrs); -impl_tridiagonal!(c32, lapacke::cgttrf, lapacke::cgtcon, lapacke::cgttrs); diff --git a/ndarray-linalg/src/layout.rs b/ndarray-linalg/src/layout.rs index 3710063b..9e43369b 100644 --- a/ndarray-linalg/src/layout.rs +++ b/ndarray-linalg/src/layout.rs @@ -1,80 +1,9 @@ //! Memory layout of matrices -use ndarray::*; - use super::error::*; +use ndarray::*; -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 as_shape(&self) -> Shape { - match *self { - MatrixLayout::C((row, col)) => (row as usize, col as usize).into_shape(), - MatrixLayout::F((col, row)) => (row as usize, col as usize).f().into_shape(), - } - } - - pub fn toggle_order(&self) -> Self { - match *self { - MatrixLayout::C((row, col)) => MatrixLayout::F((col, row)), - MatrixLayout::F((col, row)) => MatrixLayout::C((row, col)), - } - } -} +pub use lapack::layout::MatrixLayout; pub trait AllocatedArray { type Elem; diff --git a/ndarray-linalg/src/least_squares.rs b/ndarray-linalg/src/least_squares.rs index 549c8e08..9fab0062 100644 --- a/ndarray-linalg/src/least_squares.rs +++ b/ndarray-linalg/src/least_squares.rs @@ -719,7 +719,6 @@ mod tests { // Testing error cases // use crate::layout::MatrixLayout; - use ndarray::ErrorKind; #[test] fn test_incompatible_shape_error_on_mismatching_num_rows() { @@ -727,12 +726,7 @@ mod tests { let b: Array1 = array![1., 2.]; let res = a.least_squares(&b); match res { - Err(err) => match err { - LinalgError::Shape(shape_error) => { - assert_eq!(shape_error.kind(), ErrorKind::IncompatibleShape) - } - _ => panic!("Expected ShapeError"), - }, + Err(LinalgError::Lapack(err)) if matches!(err, lapack::error::Error::InvalidShape) => {} _ => panic!("Expected Err()"), } } @@ -745,12 +739,7 @@ mod tests { let res = a.least_squares(&b); match res { - Err(err) => match err { - LinalgError::Shape(shape_error) => { - assert_eq!(shape_error.kind(), ErrorKind::IncompatibleShape) - } - _ => panic!("Expected ShapeError"), - }, + Err(LinalgError::Lapack(err)) if matches!(err, lapack::error::Error::InvalidShape) => {} _ => panic!("Expected Err()"), } } diff --git a/ndarray-linalg/src/lib.rs b/ndarray-linalg/src/lib.rs index b34e21df..0a22591f 100644 --- a/ndarray-linalg/src/lib.rs +++ b/ndarray-linalg/src/lib.rs @@ -48,8 +48,7 @@ #[macro_use] extern crate ndarray; -extern crate blas_src; -extern crate lapack_src; +extern crate lax as lapack; pub mod assert; pub mod cholesky; @@ -61,7 +60,6 @@ pub mod error; pub mod generate; pub mod inner; pub mod krylov; -pub mod lapack; pub mod layout; pub mod least_squares; pub mod lobpcg; diff --git a/ndarray-linalg/src/lobpcg/lobpcg.rs b/ndarray-linalg/src/lobpcg/lobpcg.rs index a86f868e..c942a743 100644 --- a/ndarray-linalg/src/lobpcg/lobpcg.rs +++ b/ndarray-linalg/src/lobpcg/lobpcg.rs @@ -337,7 +337,9 @@ pub fn lobpcg< // if this fails (or the algorithm was restarted), then just use span{R, X} let result = p_ap .as_ref() - .ok_or(LinalgError::LapackComputationalFailure { return_code: 1 }) + .ok_or(LinalgError::Lapack( + lapack::error::Error::LapackComputationalFailure { return_code: 1 }, + )) .and_then(|(active_p, active_ap)| { let xap = x.t().dot(active_ap); let rap = r.t().dot(active_ap); diff --git a/ndarray-linalg/src/solve.rs b/ndarray-linalg/src/solve.rs index 60eb070b..566511f3 100644 --- a/ndarray-linalg/src/solve.rs +++ b/ndarray-linalg/src/solve.rs @@ -476,7 +476,8 @@ where self.ensure_square()?; match self.factorize() { Ok(fac) => fac.sln_det(), - Err(LinalgError::LapackComputationalFailure { .. }) => { + Err(LinalgError::Lapack(e)) if matches!(e, lapack::error::Error::LapackComputationalFailure {..}) => + { // The determinant is zero. Ok((A::zero(), A::Real::neg_infinity())) } @@ -494,7 +495,8 @@ where self.ensure_square()?; match self.factorize_into() { Ok(fac) => fac.sln_det_into(), - Err(LinalgError::LapackComputationalFailure { .. }) => { + Err(LinalgError::Lapack(e)) if matches!(e, lapack::error::Error::LapackComputationalFailure { .. }) => + { // The determinant is zero. Ok((A::zero(), A::Real::neg_infinity())) } @@ -538,11 +540,11 @@ where { fn rcond(&self) -> Result { unsafe { - A::rcond( + Ok(A::rcond( self.a.layout()?, self.a.as_allocated()?, self.a.opnorm_one()?, - ) + )?) } } } diff --git a/ndarray-linalg/src/solveh.rs b/ndarray-linalg/src/solveh.rs index 481757e6..c05e2bde 100644 --- a/ndarray-linalg/src/solveh.rs +++ b/ndarray-linalg/src/solveh.rs @@ -423,7 +423,8 @@ where fn sln_deth(&self) -> Result<(A::Real, A::Real)> { match self.factorizeh() { Ok(fac) => Ok(fac.sln_deth()), - Err(LinalgError::LapackComputationalFailure { .. }) => { + Err(LinalgError::Lapack(e)) if matches!(e, lapack::error::Error::LapackComputationalFailure {..}) => + { // Determinant is zero. Ok((A::Real::zero(), A::Real::neg_infinity())) } @@ -447,7 +448,8 @@ where fn sln_deth_into(self) -> Result<(A::Real, A::Real)> { match self.factorizeh_into() { Ok(fac) => Ok(fac.sln_deth_into()), - Err(LinalgError::LapackComputationalFailure { .. }) => { + Err(LinalgError::Lapack(e)) if matches!(e, lapack::error::Error::LapackComputationalFailure {..}) => + { // Determinant is zero. Ok((A::Real::zero(), A::Real::neg_infinity())) } diff --git a/ndarray-linalg/src/svddc.rs b/ndarray-linalg/src/svddc.rs index f2c102fd..22f3ae0c 100644 --- a/ndarray-linalg/src/svddc.rs +++ b/ndarray-linalg/src/svddc.rs @@ -7,19 +7,7 @@ use super::error::*; use super::layout::*; use super::types::*; -/// Specifies how many of the columns of *U* and rows of *V*ᵀ are computed and returned. -/// -/// For an input array of shape *m*×*n*, the following are computed: -#[derive(Clone, Copy, Eq, PartialEq)] -#[repr(u8)] -pub enum UVTFlag { - /// All *m* columns of *U* and all *n* rows of *V*ᵀ. - Full = b'A', - /// The first min(*m*,*n*) columns of *U* and the first min(*m*,*n*) rows of *V*ᵀ. - Some = b'S', - /// No columns of *U* or rows of *V*ᵀ. - None = b'N', -} +pub use lapack::svddc::UVTFlag; /// Singular-value decomposition of matrix (copying) by divide-and-conquer pub trait SVDDC { diff --git a/ndarray-linalg/src/tridiagonal.rs b/ndarray-linalg/src/tridiagonal.rs index 63dbf3d5..1fb0d558 100644 --- a/ndarray-linalg/src/tridiagonal.rs +++ b/ndarray-linalg/src/tridiagonal.rs @@ -2,94 +2,15 @@ //! & //! Methods for tridiagonal matrices -use std::ops::{Index, IndexMut}; - -use cauchy::Scalar; -use ndarray::*; -use num_traits::One; - use super::convert::*; use super::error::*; use super::lapack::*; use super::layout::*; +use cauchy::Scalar; +use ndarray::*; +use num_traits::One; -/// Represents a tridiagonal matrix as 3 one-dimensional vectors. -/// This struct also holds the layout of the raw matrix. -#[derive(Clone, PartialEq)] -pub struct Tridiagonal { - /// layout of raw matrix - pub l: MatrixLayout, - /// (n-1) sub-diagonal elements of matrix. - pub dl: Vec, - /// (n) diagonal elements of matrix. - pub d: Vec, - /// (n-1) super-diagonal elements of matrix. - pub du: Vec, -} - -pub trait TridiagIndex { - fn to_tuple(&self) -> (i32, i32); -} -impl TridiagIndex for [Ix; 2] { - fn to_tuple(&self) -> (i32, i32) { - (self[0] as i32, self[1] as i32) - } -} - -impl Index for Tridiagonal -where - A: Scalar, - I: TridiagIndex, -{ - type Output = A; - #[inline] - fn index(&self, index: I) -> &A { - let (n, _) = self.l.size(); - let (row, col) = index.to_tuple(); - assert!( - std::cmp::max(row, col) < n, - "ndarray: index {:?} is out of bounds for array of shape {}", - [row, col], - n - ); - match row - col { - 0 => &self.d[row as usize], - 1 => &self.dl[col as usize], - -1 => &self.du[row as usize], - _ => panic!( - "ndarray-linalg::tridiagonal: index {:?} is not tridiagonal element", - [row, col] - ), - } - } -} - -impl IndexMut for Tridiagonal -where - A: Scalar, - I: TridiagIndex, -{ - #[inline] - fn index_mut(&mut self, index: I) -> &mut A { - let (n, _) = self.l.size(); - let (row, col) = index.to_tuple(); - assert!( - std::cmp::max(row, col) < n, - "ndarray: index {:?} is out of bounds for array of shape {}", - [row, col], - n - ); - match row - col { - 0 => &mut self.d[row as usize], - 1 => &mut self.dl[col as usize], - -1 => &mut self.du[row as usize], - _ => panic!( - "ndarray-linalg::tridiagonal: index {:?} is not tridiagonal element", - [row, col] - ), - } - } -} +pub use lapack::tridiagonal::{LUFactorizedTridiagonal, Tridiagonal}; /// An interface for making a Tridiagonal struct. pub trait ExtractTridiagonal { @@ -188,23 +109,6 @@ pub trait SolveTridiagonalInplace { ) -> Result<&'a mut ArrayBase>; } -/// Represents the LU factorization of a tridiagonal matrix `A` as `A = P*L*U`. -#[derive(Clone)] -pub struct LUFactorizedTridiagonal { - /// A tridiagonal matrix which consists of - /// - l : layout of raw matrix - /// - dl: (n-1) multipliers that define the matrix L. - /// - d : (n) diagonal elements of the upper triangular matrix U. - /// - du: (n-1) elements of the first super-diagonal of U. - pub a: Tridiagonal, - /// (n-2) elements of the second super-diagonal of U. - pub du2: Vec, - /// 1-norm of raw matrix (used in .rcond_tridiagonal()). - pub anom: A::Real, - /// The pivot indices that define the permutation matrix `P`. - pub ipiv: Pivot, -} - impl SolveTridiagonal for LUFactorizedTridiagonal where A: Scalar + Lapack, @@ -661,14 +565,8 @@ impl FactorizeTridiagonalInto for Tridiagonal where A: Scalar + Lapack, { - fn factorize_tridiagonal_into(mut self) -> Result> { - let (du2, anom, ipiv) = unsafe { A::lu_tridiagonal(&mut self)? }; - Ok(LUFactorizedTridiagonal { - a: self, - du2, - anom, - ipiv, - }) + fn factorize_tridiagonal_into(self) -> Result> { + Ok(unsafe { A::lu_tridiagonal(self)? }) } } @@ -677,9 +575,8 @@ where A: Scalar + Lapack, { fn factorize_tridiagonal(&self) -> Result> { - let mut a = self.clone(); - let (du2, anom, ipiv) = unsafe { A::lu_tridiagonal(&mut a)? }; - Ok(LUFactorizedTridiagonal { a, du2, anom, ipiv }) + let a = self.clone(); + Ok(unsafe { A::lu_tridiagonal(a)? }) } } @@ -689,9 +586,8 @@ where S: Data, { fn factorize_tridiagonal(&self) -> Result> { - let mut a = self.extract_tridiagonal()?; - let (du2, anom, ipiv) = unsafe { A::lu_tridiagonal(&mut a)? }; - Ok(LUFactorizedTridiagonal { a, du2, anom, ipiv }) + let a = self.extract_tridiagonal()?; + Ok(unsafe { A::lu_tridiagonal(a)? }) } } @@ -781,7 +677,7 @@ where A: Scalar + Lapack, { fn rcond_tridiagonal(&self) -> Result { - unsafe { A::rcond_tridiagonal(&self) } + unsafe { Ok(A::rcond_tridiagonal(&self)?) } } }