mdarray_linalg/
lu.rs

1//! LU, Cholesky, matrix inversion, and determinant computation utilities
2use mdarray::{DSlice, DTensor, Layout};
3use thiserror::Error;
4
5/// Error types related to matrix inversion
6#[derive(Debug, Error)]
7pub enum InvError {
8    /// The input or output matrix is not square
9    #[error("Matrix must be square: got {rows}x{cols}")]
10    NotSquare { rows: i32, cols: i32 },
11
12    /// Backend returned a non-zero error code
13    #[error("Backend error code: {0}")]
14    BackendError(i32),
15
16    /// Matrix is singular: U(i,i) is exactly zero
17    #[error("Matrix is singular: zero pivot at position {pivot}")]
18    Singular { pivot: i32 },
19
20    /// The leading principal minor is not positive (Cholesky decomp)
21    #[error("The leading principal minor is not positive")]
22    NotPositiveDefinite { lpm: i32 },
23}
24
25/// Result type for matrix inversion
26pub type InvResult<T> = Result<DTensor<T, 2>, InvError>;
27
28///  LU decomposition and matrix inversion
29pub trait LU<T> {
30    /// Computes LU decomposition overwriting existing matrices
31    fn lu_overwrite<L: Layout, Ll: Layout, Lu: Layout, Lp: Layout>(
32        &self,
33        a: &mut DSlice<T, 2, L>,
34        l: &mut DSlice<T, 2, Ll>,
35        u: &mut DSlice<T, 2, Lu>,
36        p: &mut DSlice<T, 2, Lp>,
37    );
38
39    /// Computes LU decomposition with new allocated matrices: L, U, P (permutation matrix)
40    fn lu<L: Layout>(
41        &self,
42        a: &mut DSlice<T, 2, L>,
43    ) -> (DTensor<T, 2>, DTensor<T, 2>, DTensor<T, 2>);
44
45    /// Computes inverse overwriting the input matrix
46    fn inv_overwrite<L: Layout>(&self, a: &mut DSlice<T, 2, L>) -> Result<(), InvError>;
47
48    /// Computes inverse with new allocated matrix
49    fn inv<L: Layout>(&self, a: &mut DSlice<T, 2, L>) -> InvResult<T>;
50
51    /// Computes the determinant of a square matrix. Panics if the
52    /// matrix is non-square.
53    fn det<L: Layout>(&self, a: &mut DSlice<T, 2, L>) -> T;
54
55    /// Computes the Cholesky decomposition, returning a lower-triangular matrix
56    fn choleski<L: Layout>(&self, a: &mut DSlice<T, 2, L>) -> InvResult<T>;
57
58    /// Computes the Cholesky decomposition in-place, overwriting the input matrix
59    fn choleski_overwrite<L: Layout>(&self, a: &mut DSlice<T, 2, L>) -> Result<(), InvError>;
60}