lax/
lib.rs

1//! Safe Rust wrapper for LAPACK without external dependency.
2//!
3//! [Lapack] trait
4//! ----------------
5//!
6//! This crates provides LAPACK wrapper as a traits.
7//! For example, LU decomposition of general matrices is provided like:
8//!
9//! ```ignore
10//! pub trait Lapack {
11//!     fn lu(l: MatrixLayout, a: &mut [Self]) -> Result<Pivot>;
12//! }
13//! ```
14//!
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`:
19//!
20//! ```
21//! use lax::{Lapack, layout::MatrixLayout, Transpose};
22//!
23//! let mut a = vec![
24//!   1.0, 2.0,
25//!   3.0, 4.0
26//! ];
27//! let mut b = vec![1.0, 2.0];
28//! let layout = MatrixLayout::C { row: 2, lda: 2 };
29//! let pivot = f64::lu(layout, &mut a).unwrap();
30//! f64::solve(layout, Transpose::No, &a, &pivot, &mut b).unwrap();
31//! ```
32//!
33//! When you want to write generic algorithm for real and complex matrices,
34//! this trait can be used as a trait bound:
35//!
36//! ```
37//! use lax::{Lapack, layout::MatrixLayout, Transpose};
38//!
39//! fn solve_at_once<T: Lapack>(layout: MatrixLayout, a: &mut [T], b: &mut [T]) -> Result<(), lax::error::Error> {
40//!   let pivot = T::lu(layout, a)?;
41//!   T::solve(layout, Transpose::No, a, &pivot, b)?;
42//!   Ok(())
43//! }
44//! ```
45//!
46//! There are several similar traits as described below to keep development easy.
47//! They are merged into a single trait, [Lapack].
48//!
49//! Linear equation, Inverse matrix, Condition number
50//! --------------------------------------------------
51//!
52//! According to the property input metrix, several types of triangular decomposition are used:
53//!
54//! - [solve] module provides methods for LU-decomposition for general matrix.
55//! - [solveh] module provides methods for Bunch-Kaufman diagonal pivoting method for symmetric/Hermitian indefinite matrix.
56//! - [cholesky] module provides methods for Cholesky decomposition for symmetric/Hermitian positive dinite matrix.
57//!
58//! Eigenvalue Problem
59//! -------------------
60//!
61//! According to the property input metrix,
62//! there are several types of eigenvalue problem API
63//!
64//! - [eig] module for eigenvalue problem for general matrix.
65//! - [eigh] module for eigenvalue problem for symmetric/Hermitian matrix.
66//! - [eigh_generalized] module for generalized eigenvalue problem for symmetric/Hermitian matrix.
67//!
68//! Singular Value Decomposition
69//! -----------------------------
70//!
71//! - [svd] module for singular value decomposition (SVD) for general matrix
72//! - [svddc] module for singular value decomposition (SVD) with divided-and-conquer algorithm for general matrix
73//! - [least_squares] module for solving least square problem using SVD
74//!
75
76#![deny(rustdoc::broken_intra_doc_links, rustdoc::private_intra_doc_links)]
77
78#[cfg(any(feature = "intel-mkl-system", feature = "intel-mkl-static"))]
79extern crate intel_mkl_src as _src;
80
81#[cfg(any(feature = "openblas-system", feature = "openblas-static"))]
82extern crate openblas_src as _src;
83
84#[cfg(any(feature = "netlib-system", feature = "netlib-static"))]
85extern crate netlib_src as _src;
86
87pub mod alloc;
88pub mod cholesky;
89pub mod eig;
90pub mod eigh;
91pub mod eigh_generalized;
92pub mod error;
93pub mod flags;
94pub mod layout;
95pub mod least_squares;
96pub mod opnorm;
97pub mod qr;
98pub mod rcond;
99pub mod solve;
100pub mod solveh;
101pub mod svd;
102pub mod svddc;
103pub mod triangular;
104pub mod tridiagonal;
105
106pub use self::flags::*;
107pub use self::least_squares::LeastSquaresOwned;
108pub use self::svd::{SvdOwned, SvdRef};
109pub use self::tridiagonal::{LUFactorizedTridiagonal, Tridiagonal};
110
111use self::{alloc::*, error::*, layout::*};
112use cauchy::*;
113use std::mem::MaybeUninit;
114
115pub type Pivot = Vec<i32>;
116
117#[cfg_attr(doc, katexit::katexit)]
118/// Trait for primitive types which implements LAPACK subroutines
119pub trait Lapack: Scalar {
120    /// Compute right eigenvalue and eigenvectors for a general matrix
121    fn eig(
122        calc_v: bool,
123        l: MatrixLayout,
124        a: &mut [Self],
125    ) -> Result<(Vec<Self::Complex>, Vec<Self::Complex>)>;
126
127    /// Compute right eigenvalue and eigenvectors for a symmetric or Hermitian matrix
128    fn eigh(
129        calc_eigenvec: bool,
130        layout: MatrixLayout,
131        uplo: UPLO,
132        a: &mut [Self],
133    ) -> Result<Vec<Self::Real>>;
134
135    /// Compute right eigenvalue and eigenvectors for a symmetric or Hermitian matrix
136    fn eigh_generalized(
137        calc_eigenvec: bool,
138        layout: MatrixLayout,
139        uplo: UPLO,
140        a: &mut [Self],
141        b: &mut [Self],
142    ) -> Result<Vec<Self::Real>>;
143
144    /// Execute Householder reflection as the first step of QR-decomposition
145    ///
146    /// For C-continuous array,
147    /// this will call LQ-decomposition of the transposed matrix $ A^T = LQ^T $
148    fn householder(l: MatrixLayout, a: &mut [Self]) -> Result<Vec<Self>>;
149
150    /// Reconstruct Q-matrix from Householder-reflectors
151    fn q(l: MatrixLayout, a: &mut [Self], tau: &[Self]) -> Result<()>;
152
153    /// Execute QR-decomposition at once
154    fn qr(l: MatrixLayout, a: &mut [Self]) -> Result<Vec<Self>>;
155
156    /// Compute singular-value decomposition (SVD)
157    fn svd(l: MatrixLayout, calc_u: bool, calc_vt: bool, a: &mut [Self]) -> Result<SvdOwned<Self>>;
158
159    /// Compute singular value decomposition (SVD) with divide-and-conquer algorithm
160    fn svddc(layout: MatrixLayout, jobz: JobSvd, a: &mut [Self]) -> Result<SvdOwned<Self>>;
161
162    /// Compute a vector $x$ which minimizes Euclidian norm $\| Ax - b\|$
163    /// for a given matrix $A$ and a vector $b$.
164    fn least_squares(
165        a_layout: MatrixLayout,
166        a: &mut [Self],
167        b: &mut [Self],
168    ) -> Result<LeastSquaresOwned<Self>>;
169
170    /// Solve least square problems $\argmin_X \| AX - B\|$
171    fn least_squares_nrhs(
172        a_layout: MatrixLayout,
173        a: &mut [Self],
174        b_layout: MatrixLayout,
175        b: &mut [Self],
176    ) -> Result<LeastSquaresOwned<Self>>;
177
178    /// Computes the LU decomposition of a general $m \times n$ matrix
179    /// with partial pivoting with row interchanges.
180    ///
181    /// For a given matrix $A$, LU decomposition is described as $A = PLU$ where:
182    ///
183    /// - $L$ is lower matrix
184    /// - $U$ is upper matrix
185    /// - $P$ is permutation matrix represented by [Pivot]
186    ///
187    /// This is designed as two step computation according to LAPACK API:
188    ///
189    /// 1. Factorize input matrix $A$ into $L$, $U$, and $P$.
190    /// 2. Solve linear equation $Ax = b$ by [Lapack::solve]
191    ///    or compute inverse matrix $A^{-1}$ by [Lapack::inv] using the output of LU decomposition.
192    ///
193    /// Output
194    /// -------
195    /// - $U$ and $L$ are stored in `a` after LU decomposition has succeeded.
196    /// - $P$ is returned as [Pivot]
197    ///
198    /// Error
199    /// ------
200    /// - if the matrix is singular
201    ///   - On this case, `return_code` in [Error::LapackComputationalFailure] means
202    ///     `return_code`-th diagonal element of $U$ becomes zero.
203    ///
204    fn lu(l: MatrixLayout, a: &mut [Self]) -> Result<Pivot>;
205
206    /// Compute inverse matrix $A^{-1}$ from the output of LU-decomposition
207    fn inv(l: MatrixLayout, a: &mut [Self], p: &Pivot) -> Result<()>;
208
209    /// Solve linear equations $Ax = b$ using the output of LU-decomposition
210    fn solve(l: MatrixLayout, t: Transpose, a: &[Self], p: &Pivot, b: &mut [Self]) -> Result<()>;
211
212    /// Factorize symmetric/Hermitian matrix using Bunch-Kaufman diagonal pivoting method
213    ///
214    /// For a given symmetric matrix $A$,
215    /// this method factorizes $A = U^T D U$ or $A = L D L^T$ where
216    ///
217    /// - $U$ (or $L$) are is a product of permutation and unit upper (lower) triangular matrices
218    /// - $D$ is symmetric and block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
219    ///
220    /// This takes two-step approach based in LAPACK:
221    ///
222    /// 1. Factorize given matrix $A$ into upper ($U$) or lower ($L$) form with diagonal matrix $D$
223    /// 2. Then solve linear equation $Ax = b$, and/or calculate inverse matrix $A^{-1}$
224    ///
225    fn bk(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<Pivot>;
226
227    /// Compute inverse matrix $A^{-1}$ using the result of [Lapack::bk]
228    fn invh(l: MatrixLayout, uplo: UPLO, a: &mut [Self], ipiv: &Pivot) -> Result<()>;
229
230    /// Solve symmetric/Hermitian linear equation $Ax = b$ using the result of [Lapack::bk]
231    fn solveh(l: MatrixLayout, uplo: UPLO, a: &[Self], ipiv: &Pivot, b: &mut [Self]) -> Result<()>;
232
233    /// Solve symmetric/Hermitian positive-definite linear equations using Cholesky decomposition
234    ///
235    /// For a given positive definite matrix $A$,
236    /// Cholesky decomposition is described as $A = U^T U$ or $A = LL^T$ where
237    ///
238    /// - $L$ is lower matrix
239    /// - $U$ is upper matrix
240    ///
241    /// This is designed as two step computation according to LAPACK API
242    ///
243    /// 1. Factorize input matrix $A$ into $L$ or $U$
244    /// 2. Solve linear equation $Ax = b$ by [Lapack::solve_cholesky]
245    ///    or compute inverse matrix $A^{-1}$ by [Lapack::inv_cholesky]
246    ///
247    fn cholesky(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<()>;
248
249    /// Compute inverse matrix $A^{-1}$ using $U$ or $L$ calculated by [Lapack::cholesky]
250    fn inv_cholesky(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<()>;
251
252    /// Solve linear equation $Ax = b$ using $U$ or $L$ calculated by [Lapack::cholesky]
253    fn solve_cholesky(l: MatrixLayout, uplo: UPLO, a: &[Self], b: &mut [Self]) -> Result<()>;
254
255    /// Estimates the the reciprocal of the condition number of the matrix in 1-norm.
256    ///
257    /// `anorm` should be the 1-norm of the matrix `a`.
258    fn rcond(l: MatrixLayout, a: &[Self], anorm: Self::Real) -> Result<Self::Real>;
259
260    /// Compute norm of matrices
261    ///
262    /// For a $n \times m$ matrix
263    /// $$
264    /// A = \begin{pmatrix}
265    ///   a_{11} & \cdots & a_{1m} \\\\
266    ///   \vdots & \ddots & \vdots \\\\
267    ///   a_{n1} & \cdots & a_{nm}
268    /// \end{pmatrix}
269    /// $$
270    /// LAPACK can compute three types of norms:
271    ///
272    /// - Operator norm based on 1-norm for its domain linear space:
273    ///   $$
274    ///   \Vert A \Vert_1 = \sup_{\Vert x \Vert_1 = 1} \Vert Ax \Vert_1
275    ///   = \max_{1 \le j \le m } \sum_{i=1}^n |a_{ij}|
276    ///   $$
277    ///   where
278    ///   $\Vert x\Vert_1 = \sum_{j=1}^m |x_j|$
279    ///   is 1-norm for a vector $x$.
280    ///
281    /// - Operator norm based on $\infty$-norm for its domain linear space:
282    ///   $$
283    ///   \Vert A \Vert_\infty = \sup_{\Vert x \Vert_\infty = 1} \Vert Ax \Vert_\infty
284    ///   = \max_{1 \le i \le n } \sum_{j=1}^m |a_{ij}|
285    ///   $$
286    ///   where
287    ///   $\Vert x\Vert_\infty = \max_{j=1}^m |x_j|$
288    ///   is $\infty$-norm for a vector $x$.
289    ///
290    /// - Frobenious norm
291    ///   $$
292    ///   \Vert A \Vert_F = \sqrt{\mathrm{Tr} \left(AA^\dagger\right)} = \sqrt{\sum_{i=1}^n \sum_{j=1}^m |a_{ij}|^2}
293    ///   $$
294    ///
295    fn opnorm(t: NormType, l: MatrixLayout, a: &[Self]) -> Self::Real;
296
297    fn solve_triangular(
298        al: MatrixLayout,
299        bl: MatrixLayout,
300        uplo: UPLO,
301        d: Diag,
302        a: &[Self],
303        b: &mut [Self],
304    ) -> Result<()>;
305
306    /// Computes the LU factorization of a tridiagonal `m x n` matrix `a` using
307    /// partial pivoting with row interchanges.
308    fn lu_tridiagonal(a: Tridiagonal<Self>) -> Result<LUFactorizedTridiagonal<Self>>;
309
310    fn rcond_tridiagonal(lu: &LUFactorizedTridiagonal<Self>) -> Result<Self::Real>;
311
312    fn solve_tridiagonal(
313        lu: &LUFactorizedTridiagonal<Self>,
314        bl: MatrixLayout,
315        t: Transpose,
316        b: &mut [Self],
317    ) -> Result<()>;
318}
319
320macro_rules! impl_lapack {
321    ($s:ty) => {
322        impl Lapack for $s {
323            fn eig(
324                calc_v: bool,
325                l: MatrixLayout,
326                a: &mut [Self],
327            ) -> Result<(Vec<Self::Complex>, Vec<Self::Complex>)> {
328                use eig::*;
329                let work = EigWork::<$s>::new(calc_v, l)?;
330                let EigOwned { eigs, vr, vl } = work.eval(a)?;
331                Ok((eigs, vr.or(vl).unwrap_or_default()))
332            }
333
334            fn eigh(
335                calc_eigenvec: bool,
336                layout: MatrixLayout,
337                uplo: UPLO,
338                a: &mut [Self],
339            ) -> Result<Vec<Self::Real>> {
340                use eigh::*;
341                let work = EighWork::<$s>::new(calc_eigenvec, layout)?;
342                work.eval(uplo, a)
343            }
344
345            fn eigh_generalized(
346                calc_eigenvec: bool,
347                layout: MatrixLayout,
348                uplo: UPLO,
349                a: &mut [Self],
350                b: &mut [Self],
351            ) -> Result<Vec<Self::Real>> {
352                use eigh_generalized::*;
353                let work = EighGeneralizedWork::<$s>::new(calc_eigenvec, layout)?;
354                work.eval(uplo, a, b)
355            }
356
357            fn householder(l: MatrixLayout, a: &mut [Self]) -> Result<Vec<Self>> {
358                use qr::*;
359                let work = HouseholderWork::<$s>::new(l)?;
360                work.eval(a)
361            }
362
363            fn q(l: MatrixLayout, a: &mut [Self], tau: &[Self]) -> Result<()> {
364                use qr::*;
365                let mut work = QWork::<$s>::new(l)?;
366                work.calc(a, tau)?;
367                Ok(())
368            }
369
370            fn qr(l: MatrixLayout, a: &mut [Self]) -> Result<Vec<Self>> {
371                let tau = Self::householder(l, a)?;
372                let r = Vec::from(&*a);
373                Self::q(l, a, &tau)?;
374                Ok(r)
375            }
376
377            fn svd(
378                l: MatrixLayout,
379                calc_u: bool,
380                calc_vt: bool,
381                a: &mut [Self],
382            ) -> Result<SvdOwned<Self>> {
383                use svd::*;
384                let work = SvdWork::<$s>::new(l, calc_u, calc_vt)?;
385                work.eval(a)
386            }
387
388            fn svddc(layout: MatrixLayout, jobz: JobSvd, a: &mut [Self]) -> Result<SvdOwned<Self>> {
389                use svddc::*;
390                let work = SvdDcWork::<$s>::new(layout, jobz)?;
391                work.eval(a)
392            }
393
394            fn least_squares(
395                l: MatrixLayout,
396                a: &mut [Self],
397                b: &mut [Self],
398            ) -> Result<LeastSquaresOwned<Self>> {
399                let b_layout = l.resized(b.len() as i32, 1);
400                Self::least_squares_nrhs(l, a, b_layout, b)
401            }
402
403            fn least_squares_nrhs(
404                a_layout: MatrixLayout,
405                a: &mut [Self],
406                b_layout: MatrixLayout,
407                b: &mut [Self],
408            ) -> Result<LeastSquaresOwned<Self>> {
409                use least_squares::*;
410                let work = LeastSquaresWork::<$s>::new(a_layout, b_layout)?;
411                work.eval(a, b)
412            }
413
414            fn lu(l: MatrixLayout, a: &mut [Self]) -> Result<Pivot> {
415                use solve::*;
416                LuImpl::lu(l, a)
417            }
418
419            fn inv(l: MatrixLayout, a: &mut [Self], p: &Pivot) -> Result<()> {
420                use solve::*;
421                let mut work = InvWork::<$s>::new(l)?;
422                work.calc(a, p)?;
423                Ok(())
424            }
425
426            fn solve(
427                l: MatrixLayout,
428                t: Transpose,
429                a: &[Self],
430                p: &Pivot,
431                b: &mut [Self],
432            ) -> Result<()> {
433                use solve::*;
434                SolveImpl::solve(l, t, a, p, b)
435            }
436
437            fn bk(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<Pivot> {
438                use solveh::*;
439                let work = BkWork::<$s>::new(l)?;
440                work.eval(uplo, a)
441            }
442
443            fn invh(l: MatrixLayout, uplo: UPLO, a: &mut [Self], ipiv: &Pivot) -> Result<()> {
444                use solveh::*;
445                let mut work = InvhWork::<$s>::new(l)?;
446                work.calc(uplo, a, ipiv)
447            }
448
449            fn solveh(
450                l: MatrixLayout,
451                uplo: UPLO,
452                a: &[Self],
453                ipiv: &Pivot,
454                b: &mut [Self],
455            ) -> Result<()> {
456                use solveh::*;
457                SolvehImpl::solveh(l, uplo, a, ipiv, b)
458            }
459
460            fn cholesky(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<()> {
461                use cholesky::*;
462                CholeskyImpl::cholesky(l, uplo, a)
463            }
464
465            fn inv_cholesky(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<()> {
466                use cholesky::*;
467                InvCholeskyImpl::inv_cholesky(l, uplo, a)
468            }
469
470            fn solve_cholesky(
471                l: MatrixLayout,
472                uplo: UPLO,
473                a: &[Self],
474                b: &mut [Self],
475            ) -> Result<()> {
476                use cholesky::*;
477                SolveCholeskyImpl::solve_cholesky(l, uplo, a, b)
478            }
479
480            fn rcond(l: MatrixLayout, a: &[Self], anorm: Self::Real) -> Result<Self::Real> {
481                use rcond::*;
482                let mut work = RcondWork::<$s>::new(l);
483                work.calc(a, anorm)
484            }
485
486            fn opnorm(t: NormType, l: MatrixLayout, a: &[Self]) -> Self::Real {
487                use opnorm::*;
488                let mut work = OperatorNormWork::<$s>::new(t, l);
489                work.calc(a)
490            }
491
492            fn solve_triangular(
493                al: MatrixLayout,
494                bl: MatrixLayout,
495                uplo: UPLO,
496                d: Diag,
497                a: &[Self],
498                b: &mut [Self],
499            ) -> Result<()> {
500                use triangular::*;
501                SolveTriangularImpl::solve_triangular(al, bl, uplo, d, a, b)
502            }
503
504            fn lu_tridiagonal(a: Tridiagonal<Self>) -> Result<LUFactorizedTridiagonal<Self>> {
505                use tridiagonal::*;
506                let work = LuTridiagonalWork::<$s>::new(a.l);
507                work.eval(a)
508            }
509
510            fn rcond_tridiagonal(lu: &LUFactorizedTridiagonal<Self>) -> Result<Self::Real> {
511                use tridiagonal::*;
512                let mut work = RcondTridiagonalWork::<$s>::new(lu.a.l);
513                work.calc(lu)
514            }
515
516            fn solve_tridiagonal(
517                lu: &LUFactorizedTridiagonal<Self>,
518                bl: MatrixLayout,
519                t: Transpose,
520                b: &mut [Self],
521            ) -> Result<()> {
522                use tridiagonal::*;
523                SolveTridiagonalImpl::solve_tridiagonal(lu, bl, t, b)
524            }
525        }
526    };
527}
528impl_lapack!(c64);
529impl_lapack!(c32);
530impl_lapack!(f64);
531impl_lapack!(f32);