faer_precond/iluk/
numeric.rs

1//! Numeric ILU(k) factorisation.
2//!
3//! The elimination is identical to ILU(0)'s left-looking sweep — incomplete LU
4//! restricted to a prescribed pattern — only here the pattern is the larger
5//! level-of-fill pattern from [`SymbolicIluk`]. `A`'s values are scattered into
6//! the (wider) `L`/`U` slots, fill positions starting at zero.
7
8use faer::sparse::{SparseColMatRef, SymbolicSparseColMatRef};
9use faer_traits::math_utils::{abs2, copy, mul, one, recip, sub, zero};
10use faer_traits::{ComplexField, Index};
11
12use super::IlukError;
13use super::symbolic::SymbolicIluk;
14
15/// Numeric ILU(k) factor of a CSC matrix `A`.
16///
17/// Holds the `L`/`U` value buffers, the symbolic structure and the dense
18/// workspace used during refactorisation. [`Iluk::refactorize`] performs zero
19/// heap allocations.
20#[derive(Debug, Clone)]
21pub struct Iluk<I, T> {
22    pub(crate) symbolic: SymbolicIluk<I>,
23    pub(crate) l_values: Vec<T>,
24    pub(crate) u_values: Vec<T>,
25    pub(crate) workspace_w: Vec<T>,
26    pub(crate) workspace_marker: Vec<usize>,
27}
28
29impl<I: Index, T: ComplexField> Iluk<I, T> {
30    /// Allocate the value buffers and dense workspace for a symbolic factor.
31    pub fn new_with_symbolic(symbolic: SymbolicIluk<I>) -> Self {
32        let n = symbolic.dim;
33        let l_values = (0..symbolic.l_nnz()).map(|_| zero::<T>()).collect();
34        let u_values = (0..symbolic.u_nnz()).map(|_| zero::<T>()).collect();
35        let workspace_w = (0..n).map(|_| zero::<T>()).collect();
36        let workspace_marker = vec![usize::MAX; n];
37        Self {
38            symbolic,
39            l_values,
40            u_values,
41            workspace_w,
42            workspace_marker,
43        }
44    }
45
46    /// Build a fully populated ILU(k) factor from a CSC matrix `A`.
47    pub fn try_new(a: SparseColMatRef<'_, I, T>, level: usize) -> Result<Self, IlukError> {
48        let symbolic = SymbolicIluk::try_new(a.symbolic(), level)?;
49        let mut me = Self::new_with_symbolic(symbolic);
50        me.refactorize(a)?;
51        Ok(me)
52    }
53
54    /// Dimension `n` of the factored matrix.
55    #[inline]
56    pub fn dim(&self) -> usize {
57        self.symbolic.dim
58    }
59
60    /// Borrow the symbolic factor.
61    #[inline]
62    pub fn symbolic(&self) -> &SymbolicIluk<I> {
63        &self.symbolic
64    }
65
66    /// Refactorise against a new matrix with the same sparsity pattern.
67    ///
68    /// Performs zero heap allocations.
69    ///
70    /// # Errors
71    ///
72    /// - [`IlukError::PatternMismatch`] if `a`'s shape disagrees with the
73    ///   symbolic structure.
74    /// - [`IlukError::ZeroPivot`] if a pivot of magnitude zero is encountered.
75    pub fn refactorize(&mut self, a: SparseColMatRef<'_, I, T>) -> Result<(), IlukError> {
76        let n = self.symbolic.dim;
77        if a.nrows() != n || a.ncols() != n {
78            return Err(IlukError::PatternMismatch);
79        }
80
81        let symbolic = &self.symbolic;
82        let w = self.workspace_w.as_mut_slice();
83        let marker = self.workspace_marker.as_mut_slice();
84        let l_values = self.l_values.as_mut_slice();
85        let u_values = self.u_values.as_mut_slice();
86
87        for j in 0..n {
88            let u_start = symbolic.u_col_ptr[j].zx();
89            let u_end = symbolic.u_col_ptr[j + 1].zx();
90            let l_start = symbolic.l_col_ptr[j].zx();
91            let l_end = symbolic.l_col_ptr[j + 1].zx();
92
93            // Mark every pattern position in column j and zero it.
94            for raw in &symbolic.u_row_idx[u_start..u_end] {
95                let i = raw.zx();
96                marker[i] = j;
97                w[i] = zero::<T>();
98            }
99            for raw in &symbolic.l_row_idx[l_start..l_end] {
100                let i = raw.zx();
101                marker[i] = j;
102                w[i] = zero::<T>();
103            }
104            // Overlay A's column j (its pattern is a subset of the fill pattern).
105            let a_rows = a.symbolic().row_idx_of_col_raw(j);
106            let a_vals = a.val_of_col(j);
107            for (raw, val) in a_rows.iter().zip(a_vals.iter()) {
108                let i = raw.zx();
109                if marker[i] == j {
110                    w[i] = copy(val);
111                }
112            }
113
114            // Left-looking elimination: for each p < j in U's column j, subtract
115            // L[:,p] * U[p,j] from w at pattern positions.
116            for raw_p in &symbolic.u_row_idx[u_start..u_end - 1] {
117                let p = raw_p.zx();
118                let u_pj = copy(&w[p]);
119                let lp_start = symbolic.l_col_ptr[p].zx();
120                let lp_end = symbolic.l_col_ptr[p + 1].zx();
121                for (raw_row, lv) in symbolic.l_row_idx[lp_start + 1..lp_end]
122                    .iter()
123                    .zip(l_values[lp_start + 1..lp_end].iter())
124                {
125                    let i = raw_row.zx();
126                    if marker[i] == j {
127                        let upd = mul(lv, &u_pj);
128                        w[i] = sub(&w[i], &upd);
129                    }
130                }
131            }
132
133            // Gather U column j.
134            for (raw_row, dst) in symbolic.u_row_idx[u_start..u_end]
135                .iter()
136                .zip(u_values[u_start..u_end].iter_mut())
137            {
138                let i = raw_row.zx();
139                *dst = copy(&w[i]);
140            }
141
142            let pivot = copy(&u_values[u_end - 1]);
143            if abs2(&pivot) == zero::<T::Real>() {
144                return Err(IlukError::ZeroPivot { col: j });
145            }
146            let pivot_inv = recip(&pivot);
147
148            // Gather L column j (off-diagonal entries divided by the pivot).
149            l_values[l_start] = one::<T>();
150            for (raw_row, dst) in symbolic.l_row_idx[l_start + 1..l_end]
151                .iter()
152                .zip(l_values[l_start + 1..l_end].iter_mut())
153            {
154                let i = raw_row.zx();
155                *dst = mul(&w[i], &pivot_inv);
156            }
157        }
158
159        Ok(())
160    }
161
162    /// View over the `L` factor (unit lower triangular, diagonal first).
163    #[inline]
164    pub fn l_view(&self) -> SparseColMatRef<'_, I, T> {
165        let symbolic = unsafe {
166            SymbolicSparseColMatRef::<'_, I>::new_unchecked(
167                self.symbolic.dim,
168                self.symbolic.dim,
169                &self.symbolic.l_col_ptr,
170                None,
171                &self.symbolic.l_row_idx,
172            )
173        };
174        SparseColMatRef::new(symbolic, &self.l_values)
175    }
176
177    /// View over the `U` factor (upper triangular, diagonal last).
178    #[inline]
179    pub fn u_view(&self) -> SparseColMatRef<'_, I, T> {
180        let symbolic = unsafe {
181            SymbolicSparseColMatRef::<'_, I>::new_unchecked(
182                self.symbolic.dim,
183                self.symbolic.dim,
184                &self.symbolic.u_col_ptr,
185                None,
186                &self.symbolic.u_row_idx,
187            )
188        };
189        SparseColMatRef::new(symbolic, &self.u_values)
190    }
191}
faer_precond/iluk/numeric.rs

faer_precond/iluk/
numeric.rs
