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);