Skip to content

Commit

Permalink
Update documents
Browse files Browse the repository at this point in the history
  • Loading branch information
@termoshtt
termoshtt committed Sep 30, 2022
1 parent db39505 commit 23baa44
Show file tree
Hide file tree
Showing 2 changed files with 76 additions and 68 deletions.
64 changes: 28 additions & 36 deletions lax/src/lib.rs
View file
@@ -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<Pivot> { ... }
//! ```
//! pub trait Lapack{
//! fn lu(l: MatrixLayout, a: &mut [Self]) -> Result<Pivot>;
//! }
//! ```
//!
//! 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<f32>`,
//! and [c64] which is an alias to `num::Complex<f64>`.
//! 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,
Expand All @@ -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<T: Solve_>(layout: MatrixLayout, a: &mut [T], b: &mut [T]) -> Result<(), lax::error::Error> {
//! fn solve_at_once<T: Lapack>(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(())
Expand All @@ -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.
//!
Expand Down Expand Up @@ -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.
Expand All @@ -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<Pivot>;

/// 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<()>;
}

Expand Down
80 changes: 48 additions & 32 deletions lax/src/solve.rs
View file
@@ -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<Pivot>;
Expand Down Expand Up @@ -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<()>;
}
Expand All @@ -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),
Expand Down Expand Up @@ -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<T: Scalar> {
pub layout: MatrixLayout,
pub work: Vec<MaybeUninit<T>>,
}

/// 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<Self>;
Expand Down

0 comments on commit 23baa44

Please sign in to comment.