qualia_core_db/solvers/linear_algebra/
lu.rs

1//! Dynamic LU decomposition with partial pivoting (`P·A = L·U`) and determinant.
2//!
3//! The engine's canonical **dynamic** LU. The fixed-size [`super::StaticLuDecomposition`]
4//! handles only 4×4; this is the general `n×n` routine that the specialized libraries
5//! call (they keep only a thin error-mapping facade). Row-major, fails closed on a
6//! shape mismatch; a zero pivot is recorded in [`Lu::singular`] (not an error) so the
7//! determinant correctly comes out 0.
8//!
9//! Doolittle elimination with partial pivoting — `O(n³)`, numerically robust.
10
11use crate::solvers::SolversError;
12
13/// An in-place LU decomposition with partial pivoting (Doolittle), `P·A = L·U`.
14#[derive(Debug, Clone)]
15pub struct Lu {
16    /// Combined factors, row-major `n×n`: `U` on/above the diagonal, the strictly-lower
17    /// part of `L` below it (L's unit diagonal is implicit).
18    pub lu: Vec<f64>,
19    /// Row permutation: `pivots[i]` is the original row now in position `i`.
20    pub pivots: Vec<usize>,
21    /// Sign of the permutation (`+1`/`-1`), i.e. `det(P)`.
22    pub sign: f64,
23    /// `true` if a zero pivot was encountered (matrix is singular).
24    pub singular: bool,
25    pub n: usize,
26}
27
28impl Lu {
29    /// Solve `A x = b` using this factorization (apply `P`, forward-substitute `L`,
30    /// back-substitute `U`). `None` if the matrix is singular or `b` has the wrong length.
31    /// `O(n²)` given the existing `O(n³)` factorization — the reusable solve behind any
32    /// dense linear system (BVP Newton steps, implicit ODE stages, least-squares normal
33    /// equations).
34    pub fn solve(&self, b: &[f64]) -> Option<Vec<f64>> {
35        if self.singular || b.len() != self.n {
36            return None;
37        }
38        let n = self.n;
39        // Pb — apply the row permutation.
40        let mut y: Vec<f64> = (0..n).map(|i| b[self.pivots[i]]).collect();
41        // Forward substitution: L has an implicit unit diagonal.
42        for i in 0..n {
43            let mut s = y[i];
44            for j in 0..i {
45                s -= self.lu[i * n + j] * y[j];
46            }
47            y[i] = s;
48        }
49        // Back substitution against U.
50        for i in (0..n).rev() {
51            let mut s = y[i];
52            for j in (i + 1)..n {
53                s -= self.lu[i * n + j] * y[j];
54            }
55            let diag = self.lu[i * n + i];
56            if diag == 0.0 {
57                return None;
58            }
59            y[i] = s / diag;
60        }
61        Some(y)
62    }
63
64    /// `det(A) = sign · Π U[i][i]`.
65    pub fn determinant(&self) -> f64 {
66        if self.singular {
67            return 0.0;
68        }
69        let mut det = self.sign;
70        for i in 0..self.n {
71            det *= self.lu[i * self.n + i];
72        }
73        det
74    }
75}
76
77/// LU-decompose a row-major `n×n` matrix with partial pivoting. The reusable primitive
78/// behind [`determinant`] (and a building block for solves / condition estimates). O(n³).
79/// Returns [`SolversError::InvalidDimension`] for an empty or non-square input.
80pub fn lu_decompose(n: usize, data: &[f64]) -> Result<Lu, SolversError> {
81    if n == 0 || data.len() != n * n {
82        return Err(SolversError::InvalidDimension);
83    }
84    let mut a = data.to_vec();
85    let mut pivots: Vec<usize> = (0..n).collect();
86    let mut sign = 1.0_f64;
87    let mut singular = false;
88
89    for col in 0..n {
90        // Partial pivot: largest magnitude in this column at/below the diagonal.
91        let mut pivot = col;
92        let mut maxv = a[col * n + col].abs();
93        for r in (col + 1)..n {
94            let v = a[r * n + col].abs();
95            if v > maxv {
96                maxv = v;
97                pivot = r;
98            }
99        }
100        if maxv == 0.0 {
101            singular = true;
102            continue; // leave a zero on the diagonal; det → 0
103        }
104        if pivot != col {
105            for k in 0..n {
106                a.swap(col * n + k, pivot * n + k);
107            }
108            pivots.swap(col, pivot);
109            sign = -sign;
110        }
111        let diag = a[col * n + col];
112        for r in (col + 1)..n {
113            let factor = a[r * n + col] / diag;
114            a[r * n + col] = factor; // store L's multiplier in the lower triangle
115            for k in (col + 1)..n {
116                a[r * n + k] -= factor * a[col * n + k];
117            }
118        }
119    }
120
121    Ok(Lu {
122        lu: a,
123        pivots,
124        sign,
125        singular,
126        n,
127    })
128}
129
130/// Determinant of a row-major `n×n` matrix via LU decomposition with partial pivoting.
131/// O(n³), numerically robust; returns 0.0 for a singular matrix.
132pub fn determinant(n: usize, data: &[f64]) -> Result<f64, SolversError> {
133    Ok(lu_decompose(n, data)?.determinant())
134}
135
136/// Solve a general row-major `n×n` system `A x = b` via LU with partial pivoting. `None`
137/// on a shape mismatch or a singular matrix. The canonical dense solve for the engine.
138pub fn lu_solve(n: usize, a: &[f64], b: &[f64]) -> Option<Vec<f64>> {
139    lu_decompose(n, a).ok()?.solve(b)
140}
141
142#[cfg(test)]
143mod tests {
144    use super::*;
145
146    #[test]
147    fn determinant_2x2_and_3x3() {
148        // det[[1,2],[3,4]] = -2
149        assert!((determinant(2, &[1.0, 2.0, 3.0, 4.0]).unwrap() + 2.0).abs() < 1e-12);
150        // det[[6,1,1],[4,-2,5],[2,8,7]] = -306
151        let d = determinant(3, &[6.0, 1.0, 1.0, 4.0, -2.0, 5.0, 2.0, 8.0, 7.0]).unwrap();
152        assert!((d + 306.0).abs() < 1e-9, "det = {d}");
153    }
154
155    #[test]
156    fn singular_has_zero_determinant() {
157        let sing = lu_decompose(2, &[1.0, 2.0, 2.0, 4.0]).unwrap();
158        assert!(sing.singular && sing.determinant() == 0.0);
159    }
160
161    #[test]
162    fn reconstructs_permuted_a() {
163        // Rebuild L and U from the factors and verify L·U == P·A.
164        let n = 3;
165        let a = [4.0, 3.0, 2.0, 2.0, 1.0, 3.0, 3.0, 2.0, 1.0];
166        let f = lu_decompose(n, &a).unwrap();
167        assert!(!f.singular);
168        let mut l = vec![0.0; n * n];
169        let mut u = vec![0.0; n * n];
170        for i in 0..n {
171            l[i * n + i] = 1.0;
172            for j in 0..n {
173                if j < i {
174                    l[i * n + j] = f.lu[i * n + j];
175                } else {
176                    u[i * n + j] = f.lu[i * n + j];
177                }
178            }
179        }
180        // P·A
181        let mut pa = vec![0.0; n * n];
182        for i in 0..n {
183            for j in 0..n {
184                pa[i * n + j] = a[f.pivots[i] * n + j];
185            }
186        }
187        // L·U
188        for i in 0..n {
189            for j in 0..n {
190                let mut s = 0.0;
191                for k in 0..n {
192                    s += l[i * n + k] * u[k * n + j];
193                }
194                assert!((s - pa[i * n + j]).abs() < 1e-9);
195            }
196        }
197    }
198
199    #[test]
200    fn rejects_bad_dims() {
201        assert_eq!(
202            determinant(2, &[1.0, 2.0, 3.0]),
203            Err(SolversError::InvalidDimension)
204        );
205    }
206
207    #[test]
208    fn lu_solve_recovers_known_solution() {
209        // [[2,1,1],[1,3,2],[1,0,0]] x = [4,6,1] → x = [1,1,1].
210        let a = [2.0, 1.0, 1.0, 1.0, 3.0, 2.0, 1.0, 0.0, 0.0];
211        let x = lu_solve(3, &a, &[4.0, 6.0, 1.0]).unwrap();
212        for xi in &x {
213            assert!((xi - 1.0).abs() < 1e-9, "x = {x:?}");
214        }
215        // Residual A·x − b ≈ 0 on a second system.
216        let a2 = [4.0, 3.0, 6.0, 3.0];
217        let b2 = [10.0, 12.0];
218        let x2 = lu_solve(2, &a2, &b2).unwrap();
219        assert!((4.0 * x2[0] + 3.0 * x2[1] - 10.0).abs() < 1e-9);
220        assert!((6.0 * x2[0] + 3.0 * x2[1] - 12.0).abs() < 1e-9);
221        // Singular → None.
222        assert!(lu_solve(2, &[1.0, 2.0, 2.0, 4.0], &[1.0, 2.0]).is_none());
223    }
224}
qualia_core_db/solvers/linear_algebra/lu.rs

qualia_core_db/solvers/linear_algebra/
lu.rs
