faer_precond/ict/
mod.rs

1//! Threshold incomplete Cholesky preconditioner (ICT).
2//!
3//! ICT is to [`crate::Ic0`] what [`crate::Ilutp`] is to [`crate::Ilu0`]: the
4//! symmetric-positive-definite factorisation made *adaptive*. Instead of being
5//! locked to `A`'s sparsity pattern, it keeps the most significant fill each
6//! column produces — governed by a relative drop tolerance and a per-column fill
7//! budget — and discards the rest. That makes it markedly more effective than
8//! IC(0) on ill-conditioned SPD systems where zero-fill is too weak, at the cost
9//! of a value-dependent pattern.
10//!
11//! Tune it through [`IctParams`]:
12//!
13//! - `drop_tol` — relative threshold; an entry is dropped when its magnitude
14//!   falls below `drop_tol * ||column||`.
15//! - `fill` — how many off-diagonal entries to keep per column (see
16//!   [`FillControl`]).
17//!
18//! # When to use it
19//!
20//! Reach for ICT on symmetric positive-definite problems where [`crate::Ic0`]
21//! stalls — strongly anisotropic or high-contrast diffusion, ill-conditioned
22//! elasticity — and you are solving with conjugate gradient. For well-behaved
23//! SPD operators [`crate::Ic0`] is cheaper; for nonsymmetric systems use
24//! [`crate::Ilutp`].
25//!
26//! # Value-dependent pattern
27//!
28//! Like [`crate::Ilutp`], ICT's fill pattern depends on the matrix values, so
29//! there is no separate symbolic phase and [`Ict::refactorize`] rebuilds the
30//! pattern (reusing buffer capacity but **not** allocation-free). Apply itself
31//! allocates nothing.
32//!
33//! # Storage
34//!
35//! The factor `L` is stored CSC with its diagonal *first* in each column,
36//! matching [`faer::sparse::linalg::triangular_solve`]. Apply is the same two
37//! triangular solves as [`crate::Ic0`] (`M^{-1} = L^{-H} L^{-1}`).
38
39use core::fmt::Debug;
40
41use dyn_stack::{MemStack, StackReq};
42use faer::matrix_free::{BiLinOp, BiPrecond, LinOp, Precond};
43use faer::{Conj, MatMut, MatRef, Par};
44use faer_traits::{ComplexField, Index};
45
46pub mod apply;
47pub mod numeric;
48
49pub use numeric::Ict;
50
51// ICT shares the fill-budget and norm knobs with ILUTP.
52pub use crate::ilutp::{FillControl, RowNorm};
53
54/// Tuning parameters for [`Ict`].
55///
56/// Construct with [`IctParams::default`] and override fields as needed.
57#[derive(Debug, Clone, Copy, PartialEq)]
58pub struct IctParams {
59    /// Relative drop tolerance. An entry is dropped when its magnitude is below
60    /// `drop_tol * ||column||`. Default `1e-3`.
61    pub drop_tol: f64,
62    /// Per-column fill budget. Default `FillControl::Factor(5.0)`.
63    pub fill: FillControl,
64    /// Norm used to scale `drop_tol`. Default `RowNorm::Two`.
65    pub norm: RowNorm,
66}
67
68impl Default for IctParams {
69    fn default() -> Self {
70        Self {
71            drop_tol: 1e-3,
72            fill: FillControl::Factor(5.0),
73            norm: RowNorm::Two,
74        }
75    }
76}
77
78impl IctParams {
79    pub(crate) fn validate(&self) -> Result<(), IctError> {
80        if !self.drop_tol.is_finite() || self.drop_tol < 0.0 {
81            return Err(IctError::InvalidDropTol);
82        }
83        if let FillControl::Factor(f) = self.fill
84            && (!f.is_finite() || f <= 0.0)
85        {
86            return Err(IctError::InvalidFillControl);
87        }
88        Ok(())
89    }
90}
91
92/// Error returned by ICT construction or refactorisation.
93#[derive(Debug, Clone, PartialEq, Eq)]
94pub enum IctError {
95    /// The source matrix was not square.
96    NonSquareMatrix { nrows: usize, ncols: usize },
97    /// A non-positive pivot was encountered at column `col` — the matrix is not
98    /// positive definite, or ICT has broken down for these parameters.
99    NotPositiveDefinite { col: usize },
100    /// A refactorisation was attempted with a matrix of a different dimension.
101    PatternMismatch,
102    /// `drop_tol` was negative or NaN.
103    InvalidDropTol,
104    /// `FillControl::Factor` was non-positive or NaN.
105    InvalidFillControl,
106}
107
108impl core::fmt::Display for IctError {
109    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
110        match self {
111            Self::NonSquareMatrix { nrows, ncols } => {
112                write!(f, "matrix must be square but is {nrows}x{ncols}")
113            }
114            Self::NotPositiveDefinite { col } => {
115                write!(f, "encountered a non-positive pivot at column {col}")
116            }
117            Self::PatternMismatch => f.write_str("refactorisation dimension does not match"),
118            Self::InvalidDropTol => f.write_str("drop_tol must be finite and non-negative"),
119            Self::InvalidFillControl => f.write_str("fill factor must be finite and positive"),
120        }
121    }
122}
123
124impl core::error::Error for IctError {}
125
126impl<I, T> LinOp<T> for Ict<I, T>
127where
128    I: Index,
129    T: ComplexField + Debug + Sync,
130{
131    fn apply_scratch(&self, _rhs_ncols: usize, _par: Par) -> StackReq {
132        StackReq::EMPTY
133    }
134
135    fn nrows(&self) -> usize {
136        self.dim()
137    }
138
139    fn ncols(&self) -> usize {
140        self.dim()
141    }
142
143    fn apply(&self, mut out: MatMut<'_, T>, rhs: MatRef<'_, T>, par: Par, _stack: &mut MemStack) {
144        out.copy_from(rhs);
145        apply::solve_in_place(self, Conj::No, out, par);
146    }
147
148    fn conj_apply(
149        &self,
150        mut out: MatMut<'_, T>,
151        rhs: MatRef<'_, T>,
152        par: Par,
153        _stack: &mut MemStack,
154    ) {
155        out.copy_from(rhs);
156        apply::solve_in_place(self, Conj::Yes, out, par);
157    }
158}
159
160impl<I, T> Precond<T> for Ict<I, T>
161where
162    I: Index,
163    T: ComplexField + Debug + Sync,
164{
165    fn apply_in_place_scratch(&self, _rhs_ncols: usize, _par: Par) -> StackReq {
166        StackReq::EMPTY
167    }
168
169    fn apply_in_place(&self, rhs: MatMut<'_, T>, par: Par, _stack: &mut MemStack) {
170        apply::solve_in_place(self, Conj::No, rhs, par);
171    }
172
173    fn conj_apply_in_place(&self, rhs: MatMut<'_, T>, par: Par, _stack: &mut MemStack) {
174        apply::solve_in_place(self, Conj::Yes, rhs, par);
175    }
176}
177
178impl<I, T> BiLinOp<T> for Ict<I, T>
179where
180    I: Index,
181    T: ComplexField + Debug + Sync,
182{
183    fn transpose_apply_scratch(&self, _rhs_ncols: usize, _par: Par) -> StackReq {
184        StackReq::EMPTY
185    }
186
187    fn transpose_apply(
188        &self,
189        mut out: MatMut<'_, T>,
190        rhs: MatRef<'_, T>,
191        par: Par,
192        _stack: &mut MemStack,
193    ) {
194        // M = L L^H is Hermitian, so M^{-T} = conj(M^{-1}).
195        out.copy_from(rhs);
196        apply::solve_in_place(self, Conj::Yes, out, par);
197    }
198
199    fn adjoint_apply(
200        &self,
201        mut out: MatMut<'_, T>,
202        rhs: MatRef<'_, T>,
203        par: Par,
204        _stack: &mut MemStack,
205    ) {
206        // M^{-H} = M^{-1} for Hermitian M.
207        out.copy_from(rhs);
208        apply::solve_in_place(self, Conj::No, out, par);
209    }
210}
211
212impl<I, T> BiPrecond<T> for Ict<I, T>
213where
214    I: Index,
215    T: ComplexField + Debug + Sync,
216{
217    fn transpose_apply_in_place_scratch(&self, _rhs_ncols: usize, _par: Par) -> StackReq {
218        StackReq::EMPTY
219    }
220
221    fn transpose_apply_in_place(&self, rhs: MatMut<'_, T>, par: Par, _stack: &mut MemStack) {
222        apply::solve_in_place(self, Conj::Yes, rhs, par);
223    }
224
225    fn adjoint_apply_in_place(&self, rhs: MatMut<'_, T>, par: Par, _stack: &mut MemStack) {
226        apply::solve_in_place(self, Conj::No, rhs, par);
227    }
228}
229
230#[cfg(test)]
231mod tests {
232    use super::*;
233    use faer::sparse::{SparseColMat, Triplet};
234    use faer::{Mat, MatRef, mat};
235
236    fn assert_close(lhs: MatRef<'_, f64>, rhs: MatRef<'_, f64>, tol: f64) {
237        assert_eq!(lhs.nrows(), rhs.nrows());
238        assert_eq!(lhs.ncols(), rhs.ncols());
239        for j in 0..lhs.ncols() {
240            for i in 0..lhs.nrows() {
241                let diff = (*lhs.get(i, j) - *rhs.get(i, j)).abs();
242                assert!(
243                    diff <= tol,
244                    "mismatch at ({i}, {j}): lhs={}, rhs={}, diff={diff}",
245                    *lhs.get(i, j),
246                    *rhs.get(i, j),
247                );
248            }
249        }
250    }
251
252    fn to_dense(a: &SparseColMat<usize, f64>) -> Mat<f64> {
253        let n = a.nrows();
254        let mut out = Mat::<f64>::zeros(n, a.ncols());
255        let a_ref = a.as_ref();
256        for j in 0..a.ncols() {
257            let rows = a_ref.symbolic().row_idx_of_col_raw(j);
258            let vals = a_ref.val_of_col(j);
259            for (r, v) in rows.iter().zip(vals.iter()) {
260                *out.as_mut().get_mut(*r, j) = *v;
261            }
262        }
263        out
264    }
265
266    fn tridiagonal(n: usize, diag: f64, off: f64) -> SparseColMat<usize, f64> {
267        let mut triplets = Vec::new();
268        for i in 0..n {
269            triplets.push(Triplet::new(i, i, diag));
270            if i > 0 {
271                triplets.push(Triplet::new(i, i - 1, off));
272                triplets.push(Triplet::new(i - 1, i, off));
273            }
274        }
275        SparseColMat::try_new_from_triplets(n, n, &triplets).unwrap()
276    }
277
278    fn laplacian_2d(grid: usize) -> SparseColMat<usize, f64> {
279        let n = grid * grid;
280        let mut triplets = Vec::new();
281        for gy in 0..grid {
282            for gx in 0..grid {
283                let idx = gy * grid + gx;
284                triplets.push(Triplet::new(idx, idx, 4.0));
285                if gx > 0 {
286                    triplets.push(Triplet::new(idx, idx - 1, -1.0));
287                }
288                if gx + 1 < grid {
289                    triplets.push(Triplet::new(idx, idx + 1, -1.0));
290                }
291                if gy > 0 {
292                    triplets.push(Triplet::new(idx, idx - grid, -1.0));
293                }
294                if gy + 1 < grid {
295                    triplets.push(Triplet::new(idx, idx + grid, -1.0));
296                }
297            }
298        }
299        SparseColMat::try_new_from_triplets(n, n, &triplets).unwrap()
300    }
301
302    /// drop_tol = 0 and a large fill budget makes ICT an exact Cholesky, so
303    /// M^{-1} A x = x.
304    fn exact_params(n: usize) -> IctParams {
305        IctParams {
306            drop_tol: 0.0,
307            fill: FillControl::PerRow(n),
308            norm: RowNorm::Two,
309        }
310    }
311
312    #[test]
313    fn exact_keep_inverts_tridiagonal() {
314        let a = tridiagonal(7, 4.0, -1.0);
315        let pc = Ict::try_new_with_params(a.as_ref(), exact_params(7)).unwrap();
316        let a_dense = to_dense(&a);
317        let x_true = mat![[1.0], [-2.0], [3.0], [-1.0], [0.5], [2.0], [-1.5_f64]];
318        let mut rhs = (&a_dense * &x_true).to_owned();
319        pc.apply_in_place(rhs.as_mut(), Par::Seq, MemStack::new(&mut []));
320        assert_close(rhs.as_ref(), x_true.as_ref(), 1e-10);
321    }
322
323    #[test]
324    fn exact_keep_reconstructs_laplacian() {
325        // With no dropping, L L^T == A exactly (full Cholesky of the lower part).
326        let a = laplacian_2d(4);
327        let n = a.nrows();
328        let pc = Ict::try_new_with_params(a.as_ref(), exact_params(n)).unwrap();
329        let a_dense = to_dense(&a);
330        let x_true = Mat::<f64>::from_fn(n, 1, |i, _| (i % 5) as f64 - 2.0);
331        let mut rhs = (&a_dense * &x_true).to_owned();
332        pc.apply_in_place(rhs.as_mut(), Par::Seq, MemStack::new(&mut []));
333        assert_close(rhs.as_ref(), x_true.as_ref(), 1e-8);
334    }
335
336    #[test]
337    fn reduces_residual_on_laplacian() {
338        let a = laplacian_2d(8);
339        let n = a.nrows();
340        let pc = Ict::try_new(a.as_ref()).unwrap();
341        let a_dense = to_dense(&a);
342        let b = Mat::<f64>::from_fn(n, 1, |i, _| (i % 7) as f64 - 3.0);
343        let mut x = b.clone();
344        pc.apply_in_place(x.as_mut(), Par::Seq, MemStack::new(&mut []));
345        let residual = &a_dense * &x - &b;
346        let b_norm: f64 = b.as_ref().col(0).iter().map(|v| v * v).sum::<f64>().sqrt();
347        let r_norm: f64 = residual
348            .as_ref()
349            .col(0)
350            .iter()
351            .map(|v| v * v)
352            .sum::<f64>()
353            .sqrt();
354        assert!(r_norm / b_norm < 0.5, "ICT residual ratio too large");
355    }
356
357    #[test]
358    fn symmetric_transpose_equals_apply() {
359        let a = tridiagonal(6, 4.0, -1.0);
360        let pc = Ict::try_new(a.as_ref()).unwrap();
361        let rhs = mat![[1.0_f64], [-2.0], [3.0], [0.5], [-1.0], [2.0]];
362        let mut fwd = rhs.clone();
363        pc.apply_in_place(fwd.as_mut(), Par::Seq, MemStack::new(&mut []));
364        let mut tr = rhs.clone();
365        pc.transpose_apply_in_place(tr.as_mut(), Par::Seq, MemStack::new(&mut []));
366        assert_close(fwd.as_ref(), tr.as_ref(), 1e-12);
367    }
368
369    #[test]
370    fn refactorize_matches_fresh() {
371        let a1 = tridiagonal(7, 4.0, -1.0);
372        let a2 = tridiagonal(7, 5.0, -1.5);
373        let fresh = Ict::try_new(a2.as_ref()).unwrap();
374        let mut reused = Ict::try_new(a1.as_ref()).unwrap();
375        reused.refactorize(a2.as_ref()).unwrap();
376        assert_eq!(fresh.l_values.len(), reused.l_values.len());
377        for (a, b) in fresh.l_values.iter().zip(reused.l_values.iter()) {
378            assert!((a - b).abs() < 1e-12);
379        }
380    }
381
382    #[test]
383    fn rejects_non_positive_definite() {
384        // Indefinite matrix: negative diagonal.
385        let a = mat_to_sparse(&[&[1.0, 2.0], &[2.0, 1.0]]);
386        assert_eq!(
387            Ict::try_new(a.as_ref()).unwrap_err(),
388            IctError::NotPositiveDefinite { col: 1 }
389        );
390    }
391
392    #[test]
393    fn rejects_invalid_params() {
394        let a = tridiagonal(3, 4.0, -1.0);
395        let bad = IctParams {
396            drop_tol: -1.0,
397            ..Default::default()
398        };
399        assert_eq!(
400            Ict::try_new_with_params(a.as_ref(), bad).unwrap_err(),
401            IctError::InvalidDropTol
402        );
403    }
404
405    fn mat_to_sparse(rows: &[&[f64]]) -> SparseColMat<usize, f64> {
406        let n = rows.len();
407        let mut triplets = Vec::new();
408        for (i, row) in rows.iter().enumerate() {
409            for (j, &v) in row.iter().enumerate() {
410                if v != 0.0 {
411                    triplets.push(Triplet::new(i, j, v));
412                }
413            }
414        }
415        SparseColMat::try_new_from_triplets(n, n, &triplets).unwrap()
416    }
417}
faer_precond/ict/mod.rs

faer_precond/ict/
mod.rs
