diff --git a/lax/src/lib.rs b/lax/src/lib.rs index 6ed8498c..49ee1702 100644 --- a/lax/src/lib.rs +++ b/lax/src/lib.rs @@ -1,21 +1,24 @@ -//! ndarray-free safe Rust wrapper for LAPACK FFI +//! Safe Rust wrapper for LAPACK without external dependency. //! -//! `Lapack` trait and sub-traits -//! ------------------------------- +//! [Lapack] trait +//! ---------------- //! -//! This crates provides LAPACK wrapper as `impl` of traits to base scalar types. -//! For example, LU decomposition to double-precision matrix is provided like: +//! This crates provides LAPACK wrapper as a traits. +//! For example, LU decomposition of general matrices is provided like: //! -//! ```ignore -//! impl Solve_ for f64 { -//! fn lu(l: MatrixLayout, a: &mut [Self]) -> Result { ... } +//! ``` +//! pub trait Lapack{ +//! fn lu(l: MatrixLayout, a: &mut [Self]) -> Result; //! } //! ``` //! -//! see [Solve_] for detail. You can use it like `f64::lu`: +//! see [Lapack] for detail. +//! This trait is implemented for [f32], [f64], [c32] which is an alias to `num::Complex`, +//! and [c64] which is an alias to `num::Complex`. +//! You can use it like `f64::lu`: //! //! ``` -//! use lax::{Solve_, layout::MatrixLayout, Transpose}; +//! use lax::{Lapack, layout::MatrixLayout, Transpose}; //! //! let mut a = vec![ //! 1.0, 2.0, @@ -31,9 +34,9 @@ //! this trait can be used as a trait bound: //! //! ``` -//! use lax::{Solve_, layout::MatrixLayout, Transpose}; +//! use lax::{Lapack, layout::MatrixLayout, Transpose}; //! -//! fn solve_at_once(layout: MatrixLayout, a: &mut [T], b: &mut [T]) -> Result<(), lax::error::Error> { +//! 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(()) @@ -48,7 +51,7 @@ //! //! According to the property input metrix, several types of triangular decomposition are used: //! -//! - [Solve_] trait provides methods for LU-decomposition for general matrix. +//! - [solve] module 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. //! @@ -184,6 +187,18 @@ pub trait Lapack: /// Computes the LU decomposition of a general $m \times n$ matrix /// with partial pivoting with row interchanges. /// + /// 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$ by [Lapack::solve] + /// or compute inverse matrix $A^{-1}$ by [Lapack::inv] using the output of LU decomposition. + /// /// Output /// ------- /// - $U$ and $L$ are stored in `a` after LU decomposition has succeeded. @@ -195,35 +210,12 @@ pub trait Lapack: /// - 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/solve.rs b/lax/src/solve.rs index 6dfa5433..1b3239f5 100644 --- a/lax/src/solve.rs +++ b/lax/src/solve.rs @@ -1,21 +1,17 @@ +//! Solve linear equations 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] +/// Helper trait to abstract `*getrf` LAPACK routines for implementing [Lapack::lu] /// -/// This is designed as two step computation according to LAPACK API: +/// LAPACK correspondance +/// ---------------------- /// -/// 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. +/// | f32 | f64 | c32 | c64 | +/// |:-------|:-------|:-------|:-------| +/// | sgetrf | dgetrf | cgetrf | zgetrf | /// pub trait LuImpl: Scalar { fn lu(l: MatrixLayout, a: &mut [Self]) -> Result; @@ -57,6 +53,36 @@ impl_lu!(c32, lapack_sys::cgetrf_); impl_lu!(f64, lapack_sys::dgetrf_); impl_lu!(f32, lapack_sys::sgetrf_); +/// Helper trait to abstract `*getrs` LAPACK routines for implementing [Lapack::solve] +/// +/// If the array has C layout, then it needs to be handled +/// specially, since LAPACK expects a Fortran-layout array. +/// Reinterpreting a C layout array as Fortran layout is +/// equivalent to transposing it. So, we can handle the "no +/// transpose" and "transpose" cases by swapping to "transpose" +/// or "no transpose", respectively. For the "Hermite" case, we +/// can take advantage of the following: +/// +/// ```text +/// A^H x = b +/// ⟺ conj(A^T) x = b +/// ⟺ conj(conj(A^T) x) = conj(b) +/// ⟺ conj(conj(A^T)) conj(x) = conj(b) +/// ⟺ A^T conj(x) = conj(b) +/// ``` +/// +/// So, we can handle this case by switching to "no transpose" +/// (which is equivalent to transposing the array since it will +/// be reinterpreted as Fortran layout) and applying the +/// elementwise conjugate to `x` and `b`. +/// +/// LAPACK correspondance +/// ---------------------- +/// +/// | f32 | f64 | c32 | c64 | +/// |:-------|:-------|:-------|:-------| +/// | sgetrs | dgetrs | cgetrs | zgetrs | +/// pub trait SolveImpl: Scalar { fn solve(l: MatrixLayout, t: Transpose, a: &[Self], p: &Pivot, b: &mut [Self]) -> Result<()>; } @@ -71,26 +97,6 @@ macro_rules! impl_solve { ipiv: &Pivot, b: &mut [Self], ) -> Result<()> { - // If the array has C layout, then it needs to be handled - // specially, since LAPACK expects a Fortran-layout array. - // Reinterpreting a C layout array as Fortran layout is - // equivalent to transposing it. So, we can handle the "no - // transpose" and "transpose" cases by swapping to "transpose" - // or "no transpose", respectively. For the "Hermite" case, we - // can take advantage of the following: - // - // ```text - // A^H x = b - // ⟺ conj(A^T) x = b - // ⟺ conj(conj(A^T) x) = conj(b) - // ⟺ conj(conj(A^T)) conj(x) = conj(b) - // ⟺ A^T conj(x) = conj(b) - // ``` - // - // So, we can handle this case by switching to "no transpose" - // (which is equivalent to transposing the array since it will - // be reinterpreted as Fortran layout) and applying the - // elementwise conjugate to `x` and `b`. let (t, conj) = match l { MatrixLayout::C { .. } => match t { Transpose::No => (Transpose::Transpose, false), @@ -138,11 +144,21 @@ impl_solve!(f32, lapack_sys::sgetrs_); impl_solve!(c64, lapack_sys::zgetrs_); impl_solve!(c32, lapack_sys::cgetrs_); +/// Working memory for computing inverse matrix pub struct InvWork { pub layout: MatrixLayout, pub work: Vec>, } +/// Helper trait to abstract `*getri` LAPACK rotuines for implementing [Lapack::inv] +/// +/// LAPACK correspondance +/// ---------------------- +/// +/// | f32 | f64 | c32 | c64 | +/// |:-------|:-------|:-------|:-------| +/// | sgetri | dgetri | cgetri | zgetri | +/// pub trait InvWorkImpl: Sized { type Elem: Scalar; fn new(layout: MatrixLayout) -> Result;