Skip to content

Commit 23baa44

Browse files
termoshtttermoshtt
committed
Update documents
1 parent db39505 commit 23baa44

File tree

2 files changed

+76
-68
lines changed

2 files changed

+76
-68
lines changed

‎lax/src/lib.rs‎

Lines changed: 28 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -1,21 +1,24 @@
1-
//! ndarray-free safe Rust wrapper for LAPACK FFI
1+
//! Safe Rust wrapper for LAPACK without external dependency.
22
//!
3-
//! `Lapack` trait and sub-traits
4-
//! -------------------------------
3+
//! [Lapack] trait
4+
//! ----------------
55
//!
6-
//! This crates provides LAPACK wrapper as `impl` of traits to base scalar types.
7-
//! For example, LU decomposition to double-precision matrix is provided like:
6+
//! This crates provides LAPACK wrapper as a traits.
7+
//! For example, LU decomposition of general matrices is provided like:
88
//!
9-
//! ```ignore
10-
//! impl Solve_ for f64 {
11-
//! fn lu(l: MatrixLayout, a: &mut [Self]) -> Result<Pivot> { ... }
9+
//! ```
10+
//! pub trait Lapack{
11+
//! fn lu(l: MatrixLayout, a: &mut [Self]) -> Result<Pivot>;
1212
//! }
1313
//! ```
1414
//!
15-
//! see [Solve_] for detail. You can use it like `f64::lu`:
15+
//! see [Lapack] for detail.
16+
//! This trait is implemented for [f32], [f64], [c32] which is an alias to `num::Complex<f32>`,
17+
//! and [c64] which is an alias to `num::Complex<f64>`.
18+
//! You can use it like `f64::lu`:
1619
//!
1720
//! ```
18-
//! use lax::{Solve_, layout::MatrixLayout, Transpose};
21+
//! use lax::{Lapack, layout::MatrixLayout, Transpose};
1922
//!
2023
//! let mut a = vec![
2124
//! 1.0, 2.0,
@@ -31,9 +34,9 @@
3134
//! this trait can be used as a trait bound:
3235
//!
3336
//! ```
34-
//! use lax::{Solve_, layout::MatrixLayout, Transpose};
37+
//! use lax::{Lapack, layout::MatrixLayout, Transpose};
3538
//!
36-
//! fn solve_at_once<T: Solve_>(layout: MatrixLayout, a: &mut [T], b: &mut [T]) -> Result<(), lax::error::Error> {
39+
//! fn solve_at_once<T: Lapack>(layout: MatrixLayout, a: &mut [T], b: &mut [T]) -> Result<(), lax::error::Error> {
3740
//! let pivot = T::lu(layout, a)?;
3841
//! T::solve(layout, Transpose::No, a, &pivot, b)?;
3942
//! Ok(())
@@ -48,7 +51,7 @@
4851
//!
4952
//! According to the property input metrix, several types of triangular decomposition are used:
5053
//!
51-
//! - [Solve_] trait provides methods for LU-decomposition for general matrix.
54+
//! - [solve] module provides methods for LU-decomposition for general matrix.
5255
//! - [Solveh_] triat provides methods for Bunch-Kaufman diagonal pivoting method for symmetric/hermite indefinite matrix.
5356
//! - [Cholesky_] triat provides methods for Cholesky decomposition for symmetric/hermite positive dinite matrix.
5457
//!
@@ -184,6 +187,18 @@ pub trait Lapack:
184187
/// Computes the LU decomposition of a general $m \times n$ matrix
185188
/// with partial pivoting with row interchanges.
186189
///
190+
/// For a given matrix $A$, LU decomposition is described as $A = PLU$ where:
191+
///
192+
/// - $L$ is lower matrix
193+
/// - $U$ is upper matrix
194+
/// - $P$ is permutation matrix represented by [Pivot]
195+
///
196+
/// This is designed as two step computation according to LAPACK API:
197+
///
198+
/// 1. Factorize input matrix $A$ into $L$, $U$, and $P$.
199+
/// 2. Solve linear equation $Ax = b$ by [Lapack::solve]
200+
/// or compute inverse matrix $A^{-1}$ by [Lapack::inv] using the output of LU decomposition.
201+
///
187202
/// Output
188203
/// -------
189204
/// - $U$ and $L$ are stored in `a` after LU decomposition has succeeded.
@@ -195,35 +210,12 @@ pub trait Lapack:
195210
/// - On this case, `return_code` in [Error::LapackComputationalFailure] means
196211
/// `return_code`-th diagonal element of $U$ becomes zero.
197212
///
198-
/// LAPACK correspondance
199-
/// ----------------------
200-
///
201-
/// | f32 | f64 | c32 | c64 |
202-
/// |:-------|:-------|:-------|:-------|
203-
/// | sgetrf | dgetrf | cgetrf | zgetrf |
204-
///
205213
fn lu(l: MatrixLayout, a: &mut [Self]) -> Result<Pivot>;
206214

207215
/// Compute inverse matrix $A^{-1}$ from the output of LU-decomposition
208-
///
209-
/// LAPACK correspondance
210-
/// ----------------------
211-
///
212-
/// | f32 | f64 | c32 | c64 |
213-
/// |:-------|:-------|:-------|:-------|
214-
/// | sgetri | dgetri | cgetri | zgetri |
215-
///
216216
fn inv(l: MatrixLayout, a: &mut [Self], p: &Pivot) -> Result<()>;
217217

218218
/// Solve linear equations $Ax = b$ using the output of LU-decomposition
219-
///
220-
/// LAPACK correspondance
221-
/// ----------------------
222-
///
223-
/// | f32 | f64 | c32 | c64 |
224-
/// |:-------|:-------|:-------|:-------|
225-
/// | sgetrs | dgetrs | cgetrs | zgetrs |
226-
///
227219
fn solve(l: MatrixLayout, t: Transpose, a: &[Self], p: &Pivot, b: &mut [Self]) -> Result<()>;
228220
}
229221

‎lax/src/solve.rs‎

Lines changed: 48 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -1,21 +1,17 @@
1+
//! Solve linear equations using LU-decomposition
2+
13
use crate::{error::*, layout::MatrixLayout, *};
24
use cauchy::*;
35
use num_traits::{ToPrimitive, Zero};
46

5-
#[cfg_attr(doc, katexit::katexit)]
6-
/// Solve linear equations using LU-decomposition
7-
///
8-
/// For a given matrix $A$, LU decomposition is described as $A = PLU$ where:
9-
///
10-
/// - $L$ is lower matrix
11-
/// - $U$ is upper matrix
12-
/// - $P$ is permutation matrix represented by [Pivot]
7+
/// Helper trait to abstract `*getrf` LAPACK routines for implementing [Lapack::lu]
138
///
14-
/// This is designed as two step computation according to LAPACK API:
9+
/// LAPACK correspondance
10+
/// ----------------------
1511
///
16-
/// 1. Factorize input matrix $A$ into $L$, $U$, and $P$.
17-
/// 2. Solve linear equation $Ax = b$ or compute inverse matrix $A^{-1}$
18-
/// using the output of LU decomposition.
12+
/// | f32 | f64 | c32 | c64 |
13+
/// |:-------|:-------|:-------|:-------|
14+
/// | sgetrf | dgetrf | cgetrf | zgetrf |
1915
///
2016
pub trait LuImpl: Scalar {
2117
fn lu(l: MatrixLayout, a: &mut [Self]) -> Result<Pivot>;
@@ -57,6 +53,36 @@ impl_lu!(c32, lapack_sys::cgetrf_);
5753
impl_lu!(f64, lapack_sys::dgetrf_);
5854
impl_lu!(f32, lapack_sys::sgetrf_);
5955

56+
/// Helper trait to abstract `*getrs` LAPACK routines for implementing [Lapack::solve]
57+
///
58+
/// If the array has C layout, then it needs to be handled
59+
/// specially, since LAPACK expects a Fortran-layout array.
60+
/// Reinterpreting a C layout array as Fortran layout is
61+
/// equivalent to transposing it. So, we can handle the "no
62+
/// transpose" and "transpose" cases by swapping to "transpose"
63+
/// or "no transpose", respectively. For the "Hermite" case, we
64+
/// can take advantage of the following:
65+
///
66+
/// ```text
67+
/// A^H x = b
68+
/// ⟺ conj(A^T) x = b
69+
/// ⟺ conj(conj(A^T) x) = conj(b)
70+
/// ⟺ conj(conj(A^T)) conj(x) = conj(b)
71+
/// ⟺ A^T conj(x) = conj(b)
72+
/// ```
73+
///
74+
/// So, we can handle this case by switching to "no transpose"
75+
/// (which is equivalent to transposing the array since it will
76+
/// be reinterpreted as Fortran layout) and applying the
77+
/// elementwise conjugate to `x` and `b`.
78+
///
79+
/// LAPACK correspondance
80+
/// ----------------------
81+
///
82+
/// | f32 | f64 | c32 | c64 |
83+
/// |:-------|:-------|:-------|:-------|
84+
/// | sgetrs | dgetrs | cgetrs | zgetrs |
85+
///
6086
pub trait SolveImpl: Scalar {
6187
fn solve(l: MatrixLayout, t: Transpose, a: &[Self], p: &Pivot, b: &mut [Self]) -> Result<()>;
6288
}
@@ -71,26 +97,6 @@ macro_rules! impl_solve {
7197
ipiv: &Pivot,
7298
b: &mut [Self],
7399
) -> Result<()> {
74-
// If the array has C layout, then it needs to be handled
75-
// specially, since LAPACK expects a Fortran-layout array.
76-
// Reinterpreting a C layout array as Fortran layout is
77-
// equivalent to transposing it. So, we can handle the "no
78-
// transpose" and "transpose" cases by swapping to "transpose"
79-
// or "no transpose", respectively. For the "Hermite" case, we
80-
// can take advantage of the following:
81-
//
82-
// ```text
83-
// A^H x = b
84-
// ⟺ conj(A^T) x = b
85-
// ⟺ conj(conj(A^T) x) = conj(b)
86-
// ⟺ conj(conj(A^T)) conj(x) = conj(b)
87-
// ⟺ A^T conj(x) = conj(b)
88-
// ```
89-
//
90-
// So, we can handle this case by switching to "no transpose"
91-
// (which is equivalent to transposing the array since it will
92-
// be reinterpreted as Fortran layout) and applying the
93-
// elementwise conjugate to `x` and `b`.
94100
let (t, conj) = match l {
95101
MatrixLayout::C { .. } => match t {
96102
Transpose::No => (Transpose::Transpose, false),
@@ -138,11 +144,21 @@ impl_solve!(f32, lapack_sys::sgetrs_);
138144
impl_solve!(c64, lapack_sys::zgetrs_);
139145
impl_solve!(c32, lapack_sys::cgetrs_);
140146

147+
/// Working memory for computing inverse matrix
141148
pub struct InvWork<T: Scalar> {
142149
pub layout: MatrixLayout,
143150
pub work: Vec<MaybeUninit<T>>,
144151
}
145152

153+
/// Helper trait to abstract `*getri` LAPACK rotuines for implementing [Lapack::inv]
154+
///
155+
/// LAPACK correspondance
156+
/// ----------------------
157+
///
158+
/// | f32 | f64 | c32 | c64 |
159+
/// |:-------|:-------|:-------|:-------|
160+
/// | sgetri | dgetri | cgetri | zgetri |
161+
///
146162
pub trait InvWorkImpl: Sized {
147163
type Elem: Scalar;
148164
fn new(layout: MatrixLayout) -> Result<Self>;

0 commit comments

Comments
 (0)