pub fn lu_decompose(a: &Matrix) -> Option<(Matrix, Matrix, Vec<usize>)>
LU decomposition with partial pivoting. Returns (L, U, pivot) or None if singular.