mdarray_linalg/
lu.rs

1//! LU, Cholesky, matrix inversion, and determinant computation utilities
2//!
3//! ```rust,ignore
4//! use mdarray_linalg::prelude::*;
5//! use mdarray_linalg_backend::Backend;
6//!
7//! let bd = Backend::default();
8//!
9//! let a = Array::from_fn([3, 3], |i| {
10//!     (i[0] + 1) as f64 + 2. * (i[1] + 1) as f64 + if i[0] == i[1] { 3. } else { 0. }
11//! }); // invertible matrix
12//!
13//! // ----- LU decomposition -----
14//! // P * A = L * U  where P is a permutation matrix.
15//! let (l, u, p) = bd.lu(&mut a.clone());
16//!
17//! // ----- Determinant and inverse -----
18//! let d = bd.det(&mut a.clone());
19//! let a_inv = bd.inv(&mut a.clone()).expect("Can't compute inverse");
20//!
21//! // ----- Cholesky decomposition -----
22//! // For a symmetric positive-definite matrix: A = L * L^T
23//! let s = a.clone() + a.permute([1, 0]); // symmetric matrix
24//! let l = bd.cholesky(&mut s).unwrap();
25//! // Reconstruct: A ≈ L * L^T
26//! let a_reconstructed = l.dot(&l.transpose());
27//! ```
28
29use mdarray::{Array, Dim, Layout, Slice};
30use thiserror::Error;
31
32/// Error types related to matrix inversion
33#[derive(Debug, Error)]
34pub enum InvError {
35    /// The input or output matrix is not square
36    #[error("Matrix must be square: got {rows}x{cols}")]
37    NotSquare { rows: i32, cols: i32 },
38
39    /// Backend returned a non-zero error code
40    #[error("Backend error code: {0}")]
41    BackendError(i32),
42
43    /// Matrix is singular: U(i,i) is exactly zero
44    #[error("Matrix is singular: zero pivot at position {pivot}")]
45    Singular { pivot: i32 },
46
47    /// The leading principal minor is not positive (Cholesky decomp)
48    #[error("The leading principal minor is not positive")]
49    NotPositiveDefinite { lpm: i32 },
50}
51
52///  LU decomposition and matrix inversion
53pub trait LU<T, D0: Dim, D1: Dim> {
54    /// Computes LU decomposition overwriting existing matrices
55    fn lu_write<L: Layout, Ll: Layout, Lu: Layout, Lp: Layout>(
56        &self,
57        a: &mut Slice<T, (D0, D1), L>,
58        l: &mut Slice<T, (D0, D0), Ll>,
59        u: &mut Slice<T, (D0, D1), Lu>,
60        p: &mut Slice<T, (D0, D0), Lp>,
61    );
62
63    /// Computes LU decomposition with new allocated matrices: L, U, P (permutation matrix)
64    fn lu<L: Layout>(
65        &self,
66        a: &mut Slice<T, (D0, D1), L>,
67    ) -> (Array<T, (D0, D0)>, Array<T, (D0, D1)>, Array<T, (D0, D0)>);
68
69    /// Computes inverse overwriting the input matrix
70    fn inv_write<L: Layout>(&self, a: &mut Slice<T, (D0, D1), L>) -> Result<(), InvError>;
71
72    /// Computes inverse with new allocated matrix
73    fn inv<L: Layout>(
74        &self,
75        a: &mut Slice<T, (D0, D1), L>,
76    ) -> Result<Array<T, (D0, D1)>, InvError>;
77
78    /// Computes the determinant of a square matrix. Panics if the
79    /// matrix is non-square.
80    fn det<L: Layout>(&self, a: &mut Slice<T, (D0, D1), L>) -> T;
81
82    /// Computes the Cholesky decomposition, returning a lower-triangular matrix
83    fn cholesky<L: Layout>(
84        &self,
85        a: &mut Slice<T, (D0, D1), L>,
86    ) -> Result<Array<T, (D0, D1)>, InvError>;
87
88    /// Computes the Cholesky decomposition in-place, overwriting the input matrix
89    fn cholesky_write<L: Layout>(&self, a: &mut Slice<T, (D0, D1), L>) -> Result<(), InvError>;
90}
mdarray_linalg/lu.rs

mdarray_linalg/
lu.rs
