faer_precond/ict/
numeric.rs

1//! Numeric threshold incomplete Cholesky (ICT).
2//!
3//! Left-looking column Cholesky with dual dropping (a drop tolerance and a
4//! per-column fill budget), restricted to the lower triangle. Because the fill
5//! pattern depends on the matrix *values*, there is no symbolic phase — the
6//! factor is rebuilt in full each time (reusing buffer capacity).
7//!
8//! To find, when forming column `j`, every earlier column `k` that has a stored
9//! entry in row `j`, we keep the standard linked-list of active columns: bucket
10//! `head[j]` chains the columns whose current cursor points at row `j`.
11
12use faer::sparse::{SparseColMatRef, SymbolicSparseColMatRef};
13use faer_traits::math_utils::{
14    abs, abs2, add, conj, copy, from_f64, from_real, mul, neg, real, recip, sqrt, sub, zero,
15};
16use faer_traits::{ComplexField, Index};
17
18use super::{IctError, IctParams};
19use crate::ilutp::{FillControl, RowNorm};
20
21/// Numeric ICT factor of a Hermitian positive-definite CSC matrix `A`.
22///
23/// Stores the lower-triangular factor `L` (diagonal first per column) such that
24/// `L L^H ~= A`. Apply uses the same two triangular solves as [`crate::Ic0`].
25#[derive(Debug, Clone)]
26pub struct Ict<I, T> {
27    pub(crate) dim: usize,
28    pub(crate) params: IctParams,
29    pub(crate) l_col_ptr: Vec<I>,
30    pub(crate) l_row_idx: Vec<I>,
31    pub(crate) l_values: Vec<T>,
32}
33
34impl<I: Index, T: ComplexField> Ict<I, T> {
35    /// Build an ICT factor with the default parameters.
36    pub fn try_new(a: SparseColMatRef<'_, I, T>) -> Result<Self, IctError> {
37        Self::try_new_with_params(a, IctParams::default())
38    }
39
40    /// Build an ICT factor with explicit parameters.
41    ///
42    /// Only the lower triangle of `A` is read.
43    ///
44    /// # Errors
45    ///
46    /// - [`IctError::NonSquareMatrix`] if `A` is not square.
47    /// - [`IctError::InvalidDropTol`] / [`IctError::InvalidFillControl`] for bad
48    ///   parameters.
49    /// - [`IctError::NotPositiveDefinite`] if a non-positive pivot appears.
50    pub fn try_new_with_params(
51        a: SparseColMatRef<'_, I, T>,
52        params: IctParams,
53    ) -> Result<Self, IctError> {
54        params.validate()?;
55        if a.nrows() != a.ncols() {
56            return Err(IctError::NonSquareMatrix {
57                nrows: a.nrows(),
58                ncols: a.ncols(),
59            });
60        }
61        let n = a.nrows();
62        let mut me = Self {
63            dim: n,
64            params,
65            l_col_ptr: Vec::with_capacity(n + 1),
66            l_row_idx: Vec::new(),
67            l_values: Vec::new(),
68        };
69        me.factor(a)?;
70        Ok(me)
71    }
72
73    /// Refactorise against a new matrix with the same dimension.
74    ///
75    /// Reuses the existing buffer capacity but, because the pattern is value-
76    /// dependent, recomputes it from scratch — this is **not** allocation-free
77    /// in the strict sense (it may grow the buffers).
78    ///
79    /// # Errors
80    ///
81    /// As [`Ict::try_new_with_params`], plus [`IctError::PatternMismatch`] if
82    /// `a`'s dimension differs.
83    pub fn refactorize(&mut self, a: SparseColMatRef<'_, I, T>) -> Result<(), IctError> {
84        if a.nrows() != self.dim || a.ncols() != self.dim {
85            return Err(IctError::PatternMismatch);
86        }
87        self.factor(a)
88    }
89
90    /// Dimension `n` of the factored matrix.
91    #[inline]
92    pub fn dim(&self) -> usize {
93        self.dim
94    }
95
96    fn per_column_budget(&self, a: SparseColMatRef<'_, I, T>) -> usize {
97        match self.params.fill {
98            FillControl::PerRow(p) => p,
99            FillControl::Factor(f) => {
100                let nnz: usize = (0..self.dim)
101                    .map(|j| a.symbolic().row_idx_of_col_raw(j).len())
102                    .sum();
103                let per = (f * nnz as f64 / self.dim.max(1) as f64).round() as usize;
104                per.max(1)
105            }
106        }
107    }
108
109    fn factor(&mut self, a: SparseColMatRef<'_, I, T>) -> Result<(), IctError> {
110        let n = self.dim;
111        let budget = self.per_column_budget(a);
112        let drop_tol = from_f64::<T::Real>(self.params.drop_tol);
113
114        self.l_col_ptr.clear();
115        self.l_row_idx.clear();
116        self.l_values.clear();
117        self.l_col_ptr.push(I::truncate(0));
118
119        // Dense accumulator and pattern marker.
120        let mut w = (0..n).map(|_| zero::<T>()).collect::<Vec<T>>();
121        let mut marker = vec![usize::MAX; n];
122        let mut touched: Vec<usize> = Vec::new();
123
124        // Linked list of active columns: `head[r]` is the first column whose
125        // cursor points at row r; `next_col` chains the rest. `cursor[k]` is the
126        // position in l_row_idx of column k's next-to-use entry.
127        let mut head = vec![NONE; n];
128        let mut next_col = vec![NONE; n];
129        let mut cursor = vec![0usize; n];
130
131        for j in 0..n {
132            touched.clear();
133            // Diagonal slot is always part of the pattern.
134            w[j] = zero::<T>();
135            marker[j] = j;
136            touched.push(j);
137
138            // Gather A's lower column j (rows >= j) and accumulate its norm.
139            let a_rows = a.symbolic().row_idx_of_col_raw(j);
140            let a_vals = a.val_of_col(j);
141            for (raw, val) in a_rows.iter().zip(a_vals.iter()) {
142                let i = raw.zx();
143                if i < j {
144                    continue;
145                }
146                if marker[i] != j {
147                    marker[i] = j;
148                    touched.push(i);
149                    w[i] = copy(val);
150                } else {
151                    w[i] = copy(val);
152                }
153            }
154
155            // Left-looking updates from every column k with an entry in row j.
156            let mut k = head[j];
157            while k != NONE {
158                let kc = k as usize;
159                let nk = next_col[kc];
160                let pos = cursor[kc];
161                let l_jk = copy(&self.l_values[pos]);
162                let conj_l_jk = conj(&l_jk);
163                let k_end = self.l_col_ptr[kc + 1].zx();
164                for p in pos..k_end {
165                    let i = self.l_row_idx[p].zx();
166                    let upd = mul(&self.l_values[p], &conj_l_jk);
167                    if marker[i] == j {
168                        w[i] = sub(&w[i], &upd);
169                    } else {
170                        marker[i] = j;
171                        touched.push(i);
172                        w[i] = neg(&upd);
173                    }
174                }
175                // Advance column k's cursor and re-bucket it.
176                let new_pos = pos + 1;
177                cursor[kc] = new_pos;
178                if new_pos < k_end {
179                    let nr = self.l_row_idx[new_pos].zx();
180                    next_col[kc] = head[nr];
181                    head[nr] = kc as isize;
182                }
183                k = nk;
184            }
185            head[j] = NONE;
186
187            // Pivot.
188            let pivot_real = real(&w[j]);
189            if pivot_real.partial_cmp(&zero::<T::Real>()) != Some(core::cmp::Ordering::Greater) {
190                return Err(IctError::NotPositiveDefinite { col: j });
191            }
192            let pivot = from_real::<T>(&sqrt::<T::Real>(&pivot_real));
193            let pivot_inv = recip(&pivot);
194
195            // Candidate off-diagonal entries (i > j) and their relative-drop norm.
196            let mut col_norm = zero::<T::Real>();
197            let mut candidates: Vec<usize> = Vec::new();
198            for &i in &touched {
199                if i > j {
200                    candidates.push(i);
201                    col_norm = match self.params.norm {
202                        RowNorm::One => add(&col_norm, &abs(&w[i])),
203                        RowNorm::Two => add(&col_norm, &abs2(&w[i])),
204                    };
205                }
206            }
207            if matches!(self.params.norm, RowNorm::Two) {
208                col_norm = sqrt::<T::Real>(&col_norm);
209            }
210            let tau = mul(&drop_tol, &col_norm);
211
212            // Drop sub-threshold entries, then keep the `budget` largest.
213            candidates.retain(|&i| abs(&w[i]) >= tau);
214            if candidates.len() > budget {
215                candidates
216                    .sort_by(|&a_i, &b_i| abs(&w[b_i]).partial_cmp(&abs(&w[a_i])).unwrap_or(core::cmp::Ordering::Equal));
217                candidates.truncate(budget);
218            }
219            candidates.sort_unstable();
220
221            // Append column j: diagonal first, then kept entries scaled by 1/pivot.
222            self.l_row_idx.push(I::truncate(j));
223            self.l_values.push(copy(&pivot));
224            for &i in &candidates {
225                self.l_row_idx.push(I::truncate(i));
226                self.l_values.push(mul(&w[i], &pivot_inv));
227            }
228            self.l_col_ptr.push(I::truncate(self.l_row_idx.len()));
229
230            // Insert column j into the linked list at its first off-diagonal row.
231            let col_start = self.l_col_ptr[j].zx();
232            let col_end = self.l_col_ptr[j + 1].zx();
233            if col_end > col_start + 1 {
234                cursor[j] = col_start + 1;
235                let nr = self.l_row_idx[col_start + 1].zx();
236                next_col[j] = head[nr];
237                head[nr] = j as isize;
238            }
239        }
240
241        Ok(())
242    }
243
244    /// View over the lower-triangular factor `L` (diagonal first per column).
245    #[inline]
246    pub fn l_view(&self) -> SparseColMatRef<'_, I, T> {
247        let symbolic = unsafe {
248            SymbolicSparseColMatRef::<'_, I>::new_unchecked(
249                self.dim,
250                self.dim,
251                &self.l_col_ptr,
252                None,
253                &self.l_row_idx,
254            )
255        };
256        SparseColMatRef::new(symbolic, &self.l_values)
257    }
258}
259
260const NONE: isize = -1;
faer_precond/ict/numeric.rs

faer_precond/ict/
numeric.rs
