faer_precond/iluk/
mod.rs

1//! Level-of-fill incomplete LU preconditioner, ILU(k).
2//!
3//! ILU(k) sits between [`crate::Ilu0`] and [`crate::Ilutp`]. Where ILU(0) is
4//! locked to `A`'s own sparsity pattern, ILU(k) admits *structural* fill up to a
5//! chosen level `k`: level-0 entries are those of `A`, and each elimination step
6//! that touches two level-`a`/`b` entries can create a level-`a+b+1` fill entry,
7//! kept only while that level stays `<= k`. Raising `k` gives a more accurate
8//! factor (fewer Krylov iterations) at the cost of more fill and work.
9//!
10//! Unlike [`crate::Ilutp`], the pattern depends only on the *structure* of `A`,
11//! not its values — so the symbolic factor can be built once and reused, and
12//! refactorisation against new values is allocation-free.
13//!
14//! # When to use it
15//!
16//! Reach for ILU(k) on nonsymmetric sparse systems where [`crate::Ilu0`] is too
17//! weak but you would rather control fill *structurally* (a predictable, value-
18//! independent pattern) than through the drop tolerance of [`crate::Ilutp`].
19//! `k = 1` or `k = 2` is the usual range; `k = 0` reproduces [`crate::Ilu0`]
20//! exactly.
21//!
22//! # Repeated factorisation
23//!
24//! Build [`SymbolicIluk`] once, allocate an [`Iluk`] with
25//! [`Iluk::new_with_symbolic`], and call [`Iluk::refactorize`] in the hot loop —
26//! no allocation occurs.
27//!
28//! # Storage
29//!
30//! Identical to [`crate::Ilu0`]: `L` (unit lower) and `U` (upper) are stored CSC
31//! with `L`'s diagonal first and `U`'s diagonal last, and apply uses faer's
32//! sparse triangular solves directly with no heap allocation.
33
34use core::fmt::Debug;
35
36use dyn_stack::{MemStack, StackReq};
37use faer::matrix_free::{BiLinOp, BiPrecond, LinOp, Precond};
38use faer::{Conj, MatMut, MatRef, Par};
39use faer_traits::{ComplexField, Index};
40
41pub mod apply;
42pub mod numeric;
43pub mod symbolic;
44
45pub use numeric::Iluk;
46pub use symbolic::SymbolicIluk;
47
48/// Tuning parameters for [`Iluk`].
49#[derive(Debug, Clone, Copy, PartialEq, Eq)]
50pub struct IlukParams {
51    /// Fill level `k`. `0` reproduces ILU(0); higher admits more fill.
52    pub level: usize,
53}
54
55impl Default for IlukParams {
56    fn default() -> Self {
57        Self { level: 1 }
58    }
59}
60
61/// Error returned by ILU(k) construction or refactorisation.
62#[derive(Debug, Clone, PartialEq, Eq)]
63pub enum IlukError {
64    /// The source matrix was not square.
65    NonSquareMatrix { nrows: usize, ncols: usize },
66    /// Column `col` of the source matrix does not contain its diagonal entry.
67    MissingDiagonal { col: usize },
68    /// Row indices in column `col` are not sorted ascending.
69    UnsortedRowIndices { col: usize },
70    /// A refactorisation was attempted with a matrix whose pattern does not
71    /// match the symbolic factor.
72    PatternMismatch,
73    /// A zero pivot was encountered while eliminating column `col`.
74    ZeroPivot { col: usize },
75}
76
77impl core::fmt::Display for IlukError {
78    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
79        match self {
80            Self::NonSquareMatrix { nrows, ncols } => {
81                write!(f, "matrix must be square but is {nrows}x{ncols}")
82            }
83            Self::MissingDiagonal { col } => write!(f, "column {col} is missing its diagonal entry"),
84            Self::UnsortedRowIndices { col } => write!(f, "column {col} has unsorted row indices"),
85            Self::PatternMismatch => f.write_str("refactorisation pattern does not match symbolic"),
86            Self::ZeroPivot { col } => write!(f, "encountered a zero pivot at column {col}"),
87        }
88    }
89}
90
91impl core::error::Error for IlukError {}
92
93impl<I, T> LinOp<T> for Iluk<I, T>
94where
95    I: Index,
96    T: ComplexField + Debug + Sync,
97{
98    fn apply_scratch(&self, _rhs_ncols: usize, _par: Par) -> StackReq {
99        StackReq::EMPTY
100    }
101
102    fn nrows(&self) -> usize {
103        self.dim()
104    }
105
106    fn ncols(&self) -> usize {
107        self.dim()
108    }
109
110    fn apply(&self, mut out: MatMut<'_, T>, rhs: MatRef<'_, T>, par: Par, _stack: &mut MemStack) {
111        out.copy_from(rhs);
112        apply::solve_in_place(self, Conj::No, out, par);
113    }
114
115    fn conj_apply(
116        &self,
117        mut out: MatMut<'_, T>,
118        rhs: MatRef<'_, T>,
119        par: Par,
120        _stack: &mut MemStack,
121    ) {
122        out.copy_from(rhs);
123        apply::solve_in_place(self, Conj::Yes, out, par);
124    }
125}
126
127impl<I, T> Precond<T> for Iluk<I, T>
128where
129    I: Index,
130    T: ComplexField + Debug + Sync,
131{
132    fn apply_in_place_scratch(&self, _rhs_ncols: usize, _par: Par) -> StackReq {
133        StackReq::EMPTY
134    }
135
136    fn apply_in_place(&self, rhs: MatMut<'_, T>, par: Par, _stack: &mut MemStack) {
137        apply::solve_in_place(self, Conj::No, rhs, par);
138    }
139
140    fn conj_apply_in_place(&self, rhs: MatMut<'_, T>, par: Par, _stack: &mut MemStack) {
141        apply::solve_in_place(self, Conj::Yes, rhs, par);
142    }
143}
144
145impl<I, T> BiLinOp<T> for Iluk<I, T>
146where
147    I: Index,
148    T: ComplexField + Debug + Sync,
149{
150    fn transpose_apply_scratch(&self, _rhs_ncols: usize, _par: Par) -> StackReq {
151        StackReq::EMPTY
152    }
153
154    fn transpose_apply(
155        &self,
156        mut out: MatMut<'_, T>,
157        rhs: MatRef<'_, T>,
158        par: Par,
159        _stack: &mut MemStack,
160    ) {
161        out.copy_from(rhs);
162        apply::solve_transpose_in_place(self, Conj::No, out, par);
163    }
164
165    fn adjoint_apply(
166        &self,
167        mut out: MatMut<'_, T>,
168        rhs: MatRef<'_, T>,
169        par: Par,
170        _stack: &mut MemStack,
171    ) {
172        out.copy_from(rhs);
173        apply::solve_transpose_in_place(self, Conj::Yes, out, par);
174    }
175}
176
177impl<I, T> BiPrecond<T> for Iluk<I, T>
178where
179    I: Index,
180    T: ComplexField + Debug + Sync,
181{
182    fn transpose_apply_in_place_scratch(&self, _rhs_ncols: usize, _par: Par) -> StackReq {
183        StackReq::EMPTY
184    }
185
186    fn transpose_apply_in_place(&self, rhs: MatMut<'_, T>, par: Par, _stack: &mut MemStack) {
187        apply::solve_transpose_in_place(self, Conj::No, rhs, par);
188    }
189
190    fn adjoint_apply_in_place(&self, rhs: MatMut<'_, T>, par: Par, _stack: &mut MemStack) {
191        apply::solve_transpose_in_place(self, Conj::Yes, rhs, par);
192    }
193}
194
195#[cfg(test)]
196mod tests {
197    use super::*;
198    use crate::ilu0::SymbolicIlu0;
199    use faer::sparse::{SparseColMat, Triplet};
200    use faer::{Mat, MatRef, mat};
201
202    fn assert_close(lhs: MatRef<'_, f64>, rhs: MatRef<'_, f64>, tol: f64) {
203        assert_eq!(lhs.nrows(), rhs.nrows());
204        assert_eq!(lhs.ncols(), rhs.ncols());
205        for j in 0..lhs.ncols() {
206            for i in 0..lhs.nrows() {
207                let diff = (*lhs.get(i, j) - *rhs.get(i, j)).abs();
208                assert!(
209                    diff <= tol,
210                    "mismatch at ({i}, {j}): lhs={}, rhs={}, diff={diff}",
211                    *lhs.get(i, j),
212                    *rhs.get(i, j),
213                );
214            }
215        }
216    }
217
218    fn to_dense(a: &SparseColMat<usize, f64>) -> Mat<f64> {
219        let n = a.nrows();
220        let mut out = Mat::<f64>::zeros(n, a.ncols());
221        let a_ref = a.as_ref();
222        for j in 0..a.ncols() {
223            let rows = a_ref.symbolic().row_idx_of_col_raw(j);
224            let vals = a_ref.val_of_col(j);
225            for (r, v) in rows.iter().zip(vals.iter()) {
226                *out.as_mut().get_mut(*r, j) = *v;
227            }
228        }
229        out
230    }
231
232    fn sparse_view_to_dense(a: faer::sparse::SparseColMatRef<'_, usize, f64>) -> Mat<f64> {
233        let mut dense = Mat::<f64>::zeros(a.nrows(), a.ncols());
234        for j in 0..a.ncols() {
235            let rows = a.symbolic().row_idx_of_col_raw(j);
236            let vals = a.val_of_col(j);
237            for (r, v) in rows.iter().zip(vals.iter()) {
238                *dense.as_mut().get_mut(*r, j) = *v;
239            }
240        }
241        dense
242    }
243
244    fn laplacian_2d(grid: usize) -> SparseColMat<usize, f64> {
245        let n = grid * grid;
246        let mut triplets = Vec::new();
247        for gy in 0..grid {
248            for gx in 0..grid {
249                let idx = gy * grid + gx;
250                triplets.push(Triplet::new(idx, idx, 4.0));
251                if gx > 0 {
252                    triplets.push(Triplet::new(idx, idx - 1, -1.0));
253                }
254                if gx + 1 < grid {
255                    triplets.push(Triplet::new(idx, idx + 1, -1.0));
256                }
257                if gy > 0 {
258                    triplets.push(Triplet::new(idx, idx - grid, -1.0));
259                }
260                if gy + 1 < grid {
261                    triplets.push(Triplet::new(idx, idx + grid, -1.0));
262                }
263            }
264        }
265        SparseColMat::try_new_from_triplets(n, n, &triplets).unwrap()
266    }
267
268    fn tridiagonal(n: usize, diag: f64, sub: f64, sup: f64) -> SparseColMat<usize, f64> {
269        let mut triplets = Vec::new();
270        for i in 0..n {
271            triplets.push(Triplet::new(i, i, diag));
272            if i > 0 {
273                triplets.push(Triplet::new(i, i - 1, sub));
274                triplets.push(Triplet::new(i - 1, i, sup));
275            }
276        }
277        SparseColMat::try_new_from_triplets(n, n, &triplets).unwrap()
278    }
279
280    #[test]
281    fn level_zero_matches_ilu0_pattern() {
282        let a = laplacian_2d(5);
283        let sk = SymbolicIluk::<usize>::try_new(a.as_ref().symbolic(), 0).unwrap();
284        let s0 = SymbolicIlu0::<usize>::try_new(a.as_ref().symbolic()).unwrap();
285        assert_eq!(sk.l_col_ptr, s0.l_col_ptr);
286        assert_eq!(sk.l_row_idx, s0.l_row_idx);
287        assert_eq!(sk.u_col_ptr, s0.u_col_ptr);
288        assert_eq!(sk.u_row_idx, s0.u_row_idx);
289    }
290
291    #[test]
292    fn level_one_grows_the_pattern() {
293        let a = laplacian_2d(5);
294        let s0 = SymbolicIluk::<usize>::try_new(a.as_ref().symbolic(), 0).unwrap();
295        let s1 = SymbolicIluk::<usize>::try_new(a.as_ref().symbolic(), 1).unwrap();
296        assert!(
297            s1.l_nnz() + s1.u_nnz() > s0.l_nnz() + s0.u_nnz(),
298            "ILU(1) should introduce fill over ILU(0)"
299        );
300    }
301
302    #[test]
303    fn factor_matches_a_on_its_own_pattern() {
304        // L*U must agree with A at every structural entry of A.
305        let a = laplacian_2d(5);
306        let pc = Iluk::try_new(a.as_ref(), 1).unwrap();
307        let l = sparse_view_to_dense(pc.l_view());
308        let u = sparse_view_to_dense(pc.u_view());
309        let lu = &l * &u;
310        let a_dense = to_dense(&a);
311        let a_ref = a.as_ref();
312        for j in 0..a.ncols() {
313            for r in a_ref.symbolic().row_idx_of_col_raw(j) {
314                let i = *r;
315                let diff = (*lu.as_ref().get(i, j) - *a_dense.as_ref().get(i, j)).abs();
316                assert!(diff <= 1e-10, "L*U disagrees with A at ({i},{j}): {diff}");
317            }
318        }
319    }
320
321    #[test]
322    fn tridiagonal_is_exact() {
323        // No fill is possible in a tridiagonal, so ILU(k) is the exact LU.
324        let a = tridiagonal(6, 4.0, -1.0, -1.0);
325        let pc = Iluk::try_new(a.as_ref(), 2).unwrap();
326        let a_dense = to_dense(&a);
327        let x_true = mat![[1.0], [-2.0], [3.0], [-1.0], [0.5], [2.0_f64]];
328        let mut rhs = (&a_dense * &x_true).to_owned();
329        pc.apply_in_place(rhs.as_mut(), Par::Seq, MemStack::new(&mut []));
330        assert_close(rhs.as_ref(), x_true.as_ref(), 1e-12);
331    }
332
333    #[test]
334    fn refactorize_matches_fresh() {
335        let a1 = laplacian_2d(4);
336        let a2 = {
337            // Same pattern, scaled values.
338            let mut t = Vec::new();
339            let a1_ref = a1.as_ref();
340            for j in 0..a1.ncols() {
341                for (r, v) in a1_ref
342                    .symbolic()
343                    .row_idx_of_col_raw(j)
344                    .iter()
345                    .zip(a1_ref.val_of_col(j))
346                {
347                    t.push(Triplet::new(*r, j, v * 1.5 + if *r == j { 0.5 } else { 0.0 }));
348                }
349            }
350            SparseColMat::try_new_from_triplets(a1.nrows(), a1.ncols(), &t).unwrap()
351        };
352        let fresh = Iluk::try_new(a2.as_ref(), 1).unwrap();
353        let mut reused = Iluk::try_new(a1.as_ref(), 1).unwrap();
354        reused.refactorize(a2.as_ref()).unwrap();
355        for (a, b) in fresh.l_values.iter().zip(reused.l_values.iter()) {
356            assert!((a - b).abs() < 1e-12);
357        }
358        for (a, b) in fresh.u_values.iter().zip(reused.u_values.iter()) {
359            assert!((a - b).abs() < 1e-12);
360        }
361    }
362
363    #[test]
364    fn reduces_residual_on_laplacian() {
365        let a = laplacian_2d(8);
366        let n = a.nrows();
367        let pc = Iluk::try_new(a.as_ref(), 1).unwrap();
368        let a_dense = to_dense(&a);
369        let b = Mat::<f64>::from_fn(n, 1, |i, _| (i % 7) as f64 - 3.0);
370        let mut x = b.clone();
371        pc.apply_in_place(x.as_mut(), Par::Seq, MemStack::new(&mut []));
372        let residual = &a_dense * &x - &b;
373        let b_norm: f64 = b.as_ref().col(0).iter().map(|v| v * v).sum::<f64>().sqrt();
374        let r_norm: f64 = residual
375            .as_ref()
376            .col(0)
377            .iter()
378            .map(|v| v * v)
379            .sum::<f64>()
380            .sqrt();
381        assert!(r_norm / b_norm < 0.5, "ILU(1) residual ratio too large");
382    }
383
384    #[test]
385    fn rejects_non_square() {
386        let mut triplets = Vec::new();
387        for i in 0..3 {
388            triplets.push(Triplet::new(i, i, 1.0));
389        }
390        let a = SparseColMat::<usize, f64>::try_new_from_triplets(3, 4, &triplets).unwrap();
391        assert_eq!(
392            Iluk::try_new(a.as_ref(), 1).unwrap_err(),
393            IlukError::NonSquareMatrix { nrows: 3, ncols: 4 }
394        );
395    }
396}
faer_precond/iluk/mod.rs

faer_precond/iluk/
mod.rs
