diff --git a/.github/workflows/gh-pages.yml b/.github/workflows/gh-pages.yml index 90ad7b90..f8509459 100644 --- a/.github/workflows/gh-pages.yml +++ b/.github/workflows/gh-pages.yml @@ -40,8 +40,6 @@ jobs: with: command: doc args: --no-deps - env: - RUSTDOCFLAGS: "--html-in-header ndarray-linalg/katex-header.html" - name: Setup Pages uses: actions/configure-pages@v2 diff --git a/lax/Cargo.toml b/lax/Cargo.toml index a36c34ab..f46bdbb8 100644 --- a/lax/Cargo.toml +++ b/lax/Cargo.toml @@ -33,6 +33,7 @@ thiserror = "1.0.24" cauchy = "0.4.0" num-traits = "0.2.14" lapack-sys = "0.14.0" +katexit = "0.1.2" [dependencies.intel-mkl-src] version = "0.7.0" @@ -50,6 +51,3 @@ version = "0.10.4" optional = true default-features = false features = ["cblas"] - -[package.metadata.docs.rs] -rustdoc-args = ["--html-in-header", "katex-header.html"] diff --git a/lax/katex-header.html b/lax/katex-header.html deleted file mode 100644 index 6e10c052..00000000 --- a/lax/katex-header.html +++ /dev/null @@ -1,16 +0,0 @@ - - - - - diff --git a/lax/src/cholesky.rs b/lax/src/cholesky.rs index 94e683ff..9b213246 100644 --- a/lax/src/cholesky.rs +++ b/lax/src/cholesky.rs @@ -1,21 +1,53 @@ -//! Cholesky decomposition - use super::*; use crate::{error::*, layout::*}; use cauchy::*; +#[cfg_attr(doc, katexit::katexit)] +/// Solve symmetric/hermite positive-definite linear equations using Cholesky decomposition +/// +/// For a given positive definite matrix $A$, +/// Cholesky decomposition is described as $A = U^T U$ or $A = LL^T$ where +/// +/// - $L$ is lower matrix +/// - $U$ is upper matrix +/// +/// This is designed as two step computation according to LAPACK API +/// +/// 1. Factorize input matrix $A$ into $L$ or $U$ +/// 2. Solve linear equation $Ax = b$ or compute inverse matrix $A^{-1}$ +/// using $U$ or $L$. pub trait Cholesky_: Sized { - /// Cholesky: wrapper of `*potrf` + /// Compute Cholesky decomposition $A = U^T U$ or $A = L L^T$ according to [UPLO] + /// + /// LAPACK correspondance + /// ---------------------- + /// + /// | f32 | f64 | c32 | c64 | + /// |:-------|:-------|:-------|:-------| + /// | spotrf | dpotrf | cpotrf | zpotrf | /// - /// **Warning: Only the portion of `a` corresponding to `UPLO` is written.** fn cholesky(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<()>; - /// Wrapper of `*potri` + /// Compute inverse matrix $A^{-1}$ using $U$ or $L$ + /// + /// LAPACK correspondance + /// ---------------------- + /// + /// | f32 | f64 | c32 | c64 | + /// |:-------|:-------|:-------|:-------| + /// | spotri | dpotri | cpotri | zpotri | /// - /// **Warning: Only the portion of `a` corresponding to `UPLO` is written.** fn inv_cholesky(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<()>; - /// Wrapper of `*potrs` + /// Solve linear equation $Ax = b$ using $U$ or $L$ + /// + /// LAPACK correspondance + /// ---------------------- + /// + /// | f32 | f64 | c32 | c64 | + /// |:-------|:-------|:-------|:-------| + /// | spotrs | dpotrs | cpotrs | zpotrs | + /// fn solve_cholesky(l: MatrixLayout, uplo: UPLO, a: &[Self], b: &mut [Self]) -> Result<()>; } diff --git a/lax/src/eig.rs b/lax/src/eig.rs index d3c07222..deafba2f 100644 --- a/lax/src/eig.rs +++ b/lax/src/eig.rs @@ -1,12 +1,19 @@ -//! Eigenvalue decomposition for general matrices - use crate::{error::*, layout::MatrixLayout, *}; use cauchy::*; use num_traits::{ToPrimitive, Zero}; -/// Wraps `*geev` for general matrices +#[cfg_attr(doc, katexit::katexit)] +/// Eigenvalue problem for general matrix pub trait Eig_: Scalar { - /// Calculate Right eigenvalue + /// Compute right eigenvalue and eigenvectors $Ax = \lambda x$ + /// + /// LAPACK correspondance + /// ---------------------- + /// + /// | f32 | f64 | c32 | c64 | + /// |:------|:------|:------|:------| + /// | sgeev | dgeev | cgeev | zgeev | + /// fn eig( calc_v: bool, l: MatrixLayout, diff --git a/lax/src/eigh.rs b/lax/src/eigh.rs index a9406ee6..055b2f47 100644 --- a/lax/src/eigh.rs +++ b/lax/src/eigh.rs @@ -1,12 +1,20 @@ -//! Eigenvalue decomposition for Symmetric/Hermite matrices - use super::*; use crate::{error::*, layout::MatrixLayout}; use cauchy::*; use num_traits::{ToPrimitive, Zero}; +#[cfg_attr(doc, katexit::katexit)] +/// Eigenvalue problem for symmetric/hermite matrix pub trait Eigh_: Scalar { - /// Wraps `*syev` for real and `*heev` for complex + /// Compute right eigenvalue and eigenvectors $Ax = \lambda x$ + /// + /// LAPACK correspondance + /// ---------------------- + /// + /// | f32 | f64 | c32 | c64 | + /// |:------|:------|:------|:------| + /// | ssyev | dsyev | cheev | zheev | + /// fn eigh( calc_eigenvec: bool, layout: MatrixLayout, @@ -14,7 +22,15 @@ pub trait Eigh_: Scalar { a: &mut [Self], ) -> Result>; - /// Wraps `*sygv` for real and `*hegv` for complex + /// Compute generalized right eigenvalue and eigenvectors $Ax = \lambda B x$ + /// + /// LAPACK correspondance + /// ---------------------- + /// + /// | f32 | f64 | c32 | c64 | + /// |:------|:------|:------|:------| + /// | ssygv | dsygv | chegv | zhegv | + /// fn eigh_generalized( calc_eigenvec: bool, layout: MatrixLayout, diff --git a/lax/src/flags.rs b/lax/src/flags.rs index 37a11b3c..dd9dde3d 100644 --- a/lax/src/flags.rs +++ b/lax/src/flags.rs @@ -92,17 +92,19 @@ impl JobEv { } } -/// Specifies how many of the columns of *U* and rows of *V*ᵀ are computed and returned. +/// Specifies how many singular vectors are computed /// -/// For an input array of shape *m*×*n*, the following are computed: +/// For an input matrix $A$ of shape $m \times n$, +/// the following are computed on the singular value decomposition $A = U\Sigma V^T$: +#[cfg_attr(doc, katexit::katexit)] #[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)] #[repr(u8)] pub enum JobSvd { - /// All *m* columns of *U* and all *n* rows of *V*ᵀ. + /// All $m$ columns of $U$, and/or all $n$ rows of $V^T$. All = b'A', - /// The first min(*m*,*n*) columns of *U* and the first min(*m*,*n*) rows of *V*ᵀ. + /// The first $\min(m, n)$ columns of $U$ and/or the first $\min(m, n)$ rows of $V^T$. Some = b'S', - /// No columns of *U* or rows of *V*ᵀ. + /// No columns of $U$ and/or rows of $V^T$. None = b'N', } diff --git a/lax/src/least_squares.rs b/lax/src/least_squares.rs index 6be44f33..5699257d 100644 --- a/lax/src/least_squares.rs +++ b/lax/src/least_squares.rs @@ -12,14 +12,18 @@ pub struct LeastSquaresOutput { pub rank: i32, } -/// Wraps `*gelsd` +#[cfg_attr(doc, katexit::katexit)] +/// Solve least square problem pub trait LeastSquaresSvdDivideConquer_: Scalar { + /// Compute a vector $x$ which minimizes Euclidian norm $\| Ax - b\|$ + /// for a given matrix $A$ and a vector $b$. fn least_squares( a_layout: MatrixLayout, a: &mut [Self], b: &mut [Self], ) -> Result>; + /// Solve least square problems $\argmin_X \| AX - B\|$ fn least_squares_nrhs( a_layout: MatrixLayout, a: &mut [Self], diff --git a/lax/src/lib.rs b/lax/src/lib.rs index 83cf3658..47f3eedc 100644 --- a/lax/src/lib.rs +++ b/lax/src/lib.rs @@ -1,63 +1,75 @@ -//! Linear Algebra eXtension (LAX) -//! =============================== -//! //! ndarray-free safe Rust wrapper for LAPACK FFI //! -//! Linear equation, Inverse matrix, Condition number -//! -------------------------------------------------- +//! `Lapack` trait and sub-traits +//! ------------------------------- //! -//! As the property of $A$, several types of triangular factorization are used: +//! This crates provides LAPACK wrapper as `impl` of traits to base scalar types. +//! For example, LU decomposition to double-precision matrix is provided like: //! -//! - 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. +//! ```ignore +//! impl Solve_ for f64 { +//! fn lu(l: MatrixLayout, a: &mut [Self]) -> Result { ... } +//! } +//! ``` //! -//! | 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] | - | +//! see [Solve_] for detail. You can use it like `f64::lu`: //! -//! [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 +//! ``` +//! use lax::{Solve_, layout::MatrixLayout, Transpose}; //! -//! [bk]: solveh/trait.Solveh_.html#tymethod.bk -//! [solveh]: solveh/trait.Solveh_.html#tymethod.solveh -//! [invh]: solveh/trait.Solveh_.html#tymethod.invh +//! let mut a = vec![ +//! 1.0, 2.0, +//! 3.0, 4.0 +//! ]; +//! let mut b = vec![1.0, 2.0]; +//! let layout = MatrixLayout::C { row: 2, lda: 2 }; +//! let pivot = f64::lu(layout, &mut a).unwrap(); +//! f64::solve(layout, Transpose::No, &a, &pivot, &mut b).unwrap(); +//! ``` //! -//! Eigenvalue Problem -//! ------------------- +//! When you want to write generic algorithm for real and complex matrices, +//! this trait can be used as a trait bound: +//! +//! ``` +//! use lax::{Solve_, layout::MatrixLayout, Transpose}; //! -//! Solve eigenvalue problem for a matrix $A$ +//! fn solve_at_once(layout: MatrixLayout, a: &mut [T], b: &mut [T]) -> Result<(), lax::error::Error> { +//! let pivot = T::lu(layout, a)?; +//! T::solve(layout, Transpose::No, a, &pivot, b)?; +//! Ok(()) +//! } +//! ``` //! -//! $$ Av_i = \lambda_i v_i $$ +//! There are several similar traits as described below to keep development easy. +//! They are merged into a single trait, [Lapack]. //! -//! or generalized eigenvalue problem +//! Linear equation, Inverse matrix, Condition number +//! -------------------------------------------------- +//! +//! According to the property input metrix, several types of triangular decomposition are used: //! -//! $$ Av_i = \lambda_i B v_i $$ +//! - [Solve_] trait provides methods for LU-decomposition for general matrix. +//! - [Solveh_] triat provides methods for Bunch-Kaufman diagonal pivoting method for symmetric/hermite indefinite matrix. +//! - [Cholesky_] triat provides methods for Cholesky decomposition for symmetric/hermite positive dinite matrix. +//! +//! Eigenvalue Problem +//! ------------------- //! -//! | matrix type | Eigenvalue (EV) | Generalized Eigenvalue Problem (EG) | -//! |:--------------------------------|:----------------|:------------------------------------| -//! | General (GE) |[eig] | - | -//! | Symmetric (SY) / Hermitian (HE) |[eigh] |[eigh_generalized] | +//! According to the property input metrix, +//! there are several types of eigenvalue problem API //! -//! [eig]: eig/trait.Eig_.html#tymethod.eig -//! [eigh]: eigh/trait.Eigh_.html#tymethod.eigh -//! [eigh_generalized]: eigh/trait.Eigh_.html#tymethod.eigh_generalized +//! - [Eig_] trait provides methods for eigenvalue problem for general matrix. +//! - [Eigh_] trait provides methods for eigenvalue problem for symmetric/hermite matrix. //! -//! Singular Value Decomposition (SVD), Least square problem -//! ---------------------------------------------------------- +//! Singular Value Decomposition +//! ----------------------------- //! -//! | matrix type | Singular Value Decomposition (SVD) | SVD with divided-and-conquer (SDD) | Least square problem (LSD) | -//! |:-------------|:-----------------------------------|:-----------------------------------|:---------------------------| -//! | General (GE) | [svd] | [svddc] | [least_squares] | +//! - [SVD_] trait provides methods for singular value decomposition for general matrix +//! - [SVDDC_] trait provides methods for singular value decomposition for general matrix +//! with divided-and-conquer algorithm +//! - [LeastSquaresSvdDivideConquer_] trait provides methods +//! for solving least square problem by SVD //! -//! [svd]: svd/trait.SVD_.html#tymethod.svd -//! [svddc]: svddck/trait.SVDDC_.html#tymethod.svddc -//! [least_squares]: least_squares/trait.LeastSquaresSvdDivideConquer_.html#tymethod.least_squares #[cfg(any(feature = "intel-mkl-system", feature = "intel-mkl-static"))] extern crate intel_mkl_src as _src; diff --git a/lax/src/solve.rs b/lax/src/solve.rs index ae76f190..9f25c9bf 100644 --- a/lax/src/solve.rs +++ b/lax/src/solve.rs @@ -1,24 +1,66 @@ -//! Solve linear problem using LU decomposition - use crate::{error::*, layout::MatrixLayout, *}; use cauchy::*; use num_traits::{ToPrimitive, Zero}; +#[cfg_attr(doc, katexit::katexit)] +/// Solve linear equations using LU-decomposition +/// +/// For a given matrix $A$, LU decomposition is described as $A = PLU$ where: +/// +/// - $L$ is lower matrix +/// - $U$ is upper matrix +/// - $P$ is permutation matrix represented by [Pivot] +/// +/// This is designed as two step computation according to LAPACK API: +/// +/// 1. Factorize input matrix $A$ into $L$, $U$, and $P$. +/// 2. Solve linear equation $Ax = b$ or compute inverse matrix $A^{-1}$ +/// using the output of LU decomposition. +/// pub trait Solve_: Scalar + Sized { - /// Computes the LU factorization of a general `m x n` matrix `a` using - /// partial pivoting with row interchanges. + /// Computes the LU decomposition of a general $m \times n$ matrix + /// with partial pivoting with row interchanges. /// - /// $ PA = LU $ + /// Output + /// ------- + /// - $U$ and $L$ are stored in `a` after LU decomposition has succeeded. + /// - $P$ is returned as [Pivot] /// /// Error /// ------ - /// - `LapackComputationalFailure { return_code }` when the matrix is singular - /// - Division by zero will occur if it is used to solve a system of equations - /// because `U[(return_code-1, return_code-1)]` is exactly zero. + /// - if the matrix is singular + /// - On this case, `return_code` in [Error::LapackComputationalFailure] means + /// `return_code`-th diagonal element of $U$ becomes zero. + /// + /// LAPACK correspondance + /// ---------------------- + /// + /// | f32 | f64 | c32 | c64 | + /// |:-------|:-------|:-------|:-------| + /// | sgetrf | dgetrf | cgetrf | zgetrf | + /// fn lu(l: MatrixLayout, a: &mut [Self]) -> Result; + /// Compute inverse matrix $A^{-1}$ from the output of LU-decomposition + /// + /// LAPACK correspondance + /// ---------------------- + /// + /// | f32 | f64 | c32 | c64 | + /// |:-------|:-------|:-------|:-------| + /// | sgetri | dgetri | cgetri | zgetri | + /// fn inv(l: MatrixLayout, a: &mut [Self], p: &Pivot) -> Result<()>; + /// Solve linear equations $Ax = b$ using the output of LU-decomposition + /// + /// LAPACK correspondance + /// ---------------------- + /// + /// | f32 | f64 | c32 | c64 | + /// |:-------|:-------|:-------|:-------| + /// | sgetrs | dgetrs | cgetrs | zgetrs | + /// fn solve(l: MatrixLayout, t: Transpose, a: &[Self], p: &Pivot, b: &mut [Self]) -> Result<()>; } diff --git a/lax/src/solveh.rs b/lax/src/solveh.rs index 9f65978d..5d1a3d02 100644 --- a/lax/src/solveh.rs +++ b/lax/src/solveh.rs @@ -1,17 +1,55 @@ -//! Solve symmetric linear problem using the Bunch-Kaufman diagonal pivoting method. -//! -//! See also [the manual of dsytrf](http://www.netlib.org/lapack/lapack-3.1.1/html/dsytrf.f.html) - use crate::{error::*, layout::MatrixLayout, *}; use cauchy::*; use num_traits::{ToPrimitive, Zero}; +#[cfg_attr(doc, katexit::katexit)] +/// Solve symmetric/hermite indefinite linear problem using the [Bunch-Kaufman diagonal pivoting method][BK]. +/// +/// For a given symmetric matrix $A$, +/// this method factorizes $A = U^T D U$ or $A = L D L^T$ where +/// +/// - $U$ (or $L$) are is a product of permutation and unit upper (lower) triangular matrices +/// - $D$ is symmetric and block diagonal with 1-by-1 and 2-by-2 diagonal blocks. +/// +/// This takes two-step approach based in LAPACK: +/// +/// 1. Factorize given matrix $A$ into upper ($U$) or lower ($L$) form with diagonal matrix $D$ +/// 2. Then solve linear equation $Ax = b$, and/or calculate inverse matrix $A^{-1}$ +/// +/// [BK]: https://doi.org/10.2307/2005787 +/// pub trait Solveh_: Sized { - /// Bunch-Kaufman: wrapper of `*sytrf` and `*hetrf` + /// Factorize input matrix using Bunch-Kaufman diagonal pivoting method + /// + /// LAPACK correspondance + /// ---------------------- + /// + /// | f32 | f64 | c32 | c64 | + /// |:-------|:-------|:-------|:-------| + /// | ssytrf | dsytrf | chetrf | zhetrf | + /// fn bk(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result; - /// Wrapper of `*sytri` and `*hetri` + + /// Compute inverse matrix $A^{-1}$ from factroized result + /// + /// LAPACK correspondance + /// ---------------------- + /// + /// | f32 | f64 | c32 | c64 | + /// |:-------|:-------|:-------|:-------| + /// | ssytri | dsytri | chetri | zhetri | + /// fn invh(l: MatrixLayout, uplo: UPLO, a: &mut [Self], ipiv: &Pivot) -> Result<()>; - /// Wrapper of `*sytrs` and `*hetrs` + + /// Solve linear equation $Ax = b$ using factroized result + /// + /// LAPACK correspondance + /// ---------------------- + /// + /// | f32 | f64 | c32 | c64 | + /// |:-------|:-------|:-------|:-------| + /// | ssytrs | dsytrs | chetrs | zhetrs | + /// fn solveh(l: MatrixLayout, uplo: UPLO, a: &[Self], ipiv: &Pivot, b: &mut [Self]) -> Result<()>; } diff --git a/lax/src/svd.rs b/lax/src/svd.rs index 0a509a0e..0f0f81ab 100644 --- a/lax/src/svd.rs +++ b/lax/src/svd.rs @@ -14,9 +14,18 @@ pub struct SVDOutput { pub vt: Option>, } -/// Wraps `*gesvd` +#[cfg_attr(doc, katexit::katexit)] +/// Singular value decomposition pub trait SVD_: Scalar { - /// Calculate singular value decomposition $ A = U \Sigma V^T $ + /// Compute singular value decomposition $A = U \Sigma V^T$ + /// + /// LAPACK correspondance + /// ---------------------- + /// + /// | f32 | f64 | c32 | c64 | + /// |:-------|:-------|:-------|:-------| + /// | sgesvd | dgesvd | cgesvd | zgesvd | + /// fn svd(l: MatrixLayout, calc_u: bool, calc_vt: bool, a: &mut [Self]) -> Result>; } diff --git a/lax/src/svddc.rs b/lax/src/svddc.rs index bb59348f..55ddfd46 100644 --- a/lax/src/svddc.rs +++ b/lax/src/svddc.rs @@ -2,7 +2,18 @@ use crate::{error::*, layout::MatrixLayout, *}; use cauchy::*; use num_traits::{ToPrimitive, Zero}; +#[cfg_attr(doc, katexit::katexit)] +/// Singular value decomposition with divide-and-conquer method pub trait SVDDC_: Scalar { + /// Compute singular value decomposition $A = U \Sigma V^T$ + /// + /// LAPACK correspondance + /// ---------------------- + /// + /// | f32 | f64 | c32 | c64 | + /// |:-------|:-------|:-------|:-------| + /// | sgesdd | dgesdd | cgesdd | zgesdd | + /// fn svddc(l: MatrixLayout, jobz: JobSvd, a: &mut [Self]) -> Result>; } diff --git a/ndarray-linalg/Cargo.toml b/ndarray-linalg/Cargo.toml index 6147a2eb..bfab605f 100644 --- a/ndarray-linalg/Cargo.toml +++ b/ndarray-linalg/Cargo.toml @@ -30,6 +30,7 @@ intel-mkl-system = ["lax/intel-mkl-system"] [dependencies] cauchy = "0.4.0" +katexit = "0.1.2" num-complex = "0.4.0" num-traits = "0.2.14" rand = "0.8.3" @@ -78,6 +79,3 @@ harness = false [[bench]] name = "solveh" harness = false - -[package.metadata.docs.rs] -rustdoc-args = ["--html-in-header", "katex-header.html"] diff --git a/ndarray-linalg/katex-header.html b/ndarray-linalg/katex-header.html deleted file mode 100644 index 6e10c052..00000000 --- a/ndarray-linalg/katex-header.html +++ /dev/null @@ -1,16 +0,0 @@ - - - - - diff --git a/ndarray-linalg/src/eig.rs b/ndarray-linalg/src/eig.rs index 833f6d49..bde5bd7d 100644 --- a/ndarray-linalg/src/eig.rs +++ b/ndarray-linalg/src/eig.rs @@ -5,6 +5,7 @@ use crate::layout::*; use crate::types::*; use ndarray::*; +#[cfg_attr(doc, katexit::katexit)] /// Eigenvalue decomposition of general matrix reference pub trait Eig { /// EigVec is the right eivenvector