faer_precond/ilu0/
numeric.rs

1//! Numeric ILU(0) factorisation.
2
3use faer::sparse::SparseColMatRef;
4use faer_traits::math_utils::{abs2, copy, mul, one, recip, sub, zero};
5use faer_traits::{ComplexField, Index};
6
7use super::Ilu0Error;
8use super::symbolic::SymbolicIlu0;
9
10/// Numeric ILU(0) factor of a CSC matrix `A`.
11///
12/// Holds the `L` and `U` value buffers together with the symbolic structure
13/// and the dense workspace used during refactorisation. After construction the
14/// factor can be refactored against any matrix with the same sparsity pattern
15/// — [`Ilu0::refactorize`] performs zero heap allocations.
16#[derive(Debug, Clone)]
17pub struct Ilu0<I, T> {
18    pub(crate) symbolic: SymbolicIlu0<I>,
19    pub(crate) l_values: Vec<T>,
20    pub(crate) u_values: Vec<T>,
21    pub(crate) workspace_w: Vec<T>,
22    pub(crate) workspace_marker: Vec<usize>,
23}
24
25impl<I: Index, T: ComplexField> Ilu0<I, T> {
26    /// Allocate the value buffers and dense workspace for a given symbolic factor.
27    ///
28    /// The returned factor is *not* yet populated — call [`Ilu0::refactorize`]
29    /// to fill it with values from a matrix matching `symbolic`'s pattern.
30    pub fn new_with_symbolic(symbolic: SymbolicIlu0<I>) -> Self {
31        let n = symbolic.dim;
32        let l_nnz = symbolic.l_nnz();
33        let u_nnz = symbolic.u_nnz();
34        let l_values = (0..l_nnz).map(|_| zero::<T>()).collect();
35        let u_values = (0..u_nnz).map(|_| zero::<T>()).collect();
36        let workspace_w = (0..n).map(|_| zero::<T>()).collect();
37        let workspace_marker = vec![usize::MAX; n];
38        Self {
39            symbolic,
40            l_values,
41            u_values,
42            workspace_w,
43            workspace_marker,
44        }
45    }
46
47    /// Build a fully populated ILU(0) factor from a CSC matrix `A` in one shot.
48    ///
49    /// Allocates internally for the symbolic structure, the `L`/`U` value
50    /// buffers, and the dense workspace.
51    pub fn try_new(a: SparseColMatRef<'_, I, T>) -> Result<Self, Ilu0Error> {
52        let symbolic = SymbolicIlu0::try_new(a.symbolic())?;
53        let mut me = Self::new_with_symbolic(symbolic);
54        me.refactorize(a)?;
55        Ok(me)
56    }
57
58    /// Dimension `n` of the factored matrix.
59    #[inline]
60    pub fn dim(&self) -> usize {
61        self.symbolic.dim
62    }
63
64    /// Borrow the symbolic factor.
65    #[inline]
66    pub fn symbolic(&self) -> &SymbolicIlu0<I> {
67        &self.symbolic
68    }
69
70    /// Refactorise against a new matrix with the same sparsity pattern.
71    ///
72    /// Performs zero heap allocations.
73    ///
74    /// # Errors
75    ///
76    /// - [`Ilu0Error::PatternMismatch`] if `a`'s shape or column lengths
77    ///   disagree with the symbolic structure, or if the diagonal entries land
78    ///   at unexpected positions.
79    /// - [`Ilu0Error::ZeroPivot`] if a pivot of magnitude zero is encountered
80    ///   during elimination.
81    pub fn refactorize(&mut self, a: SparseColMatRef<'_, I, T>) -> Result<(), Ilu0Error> {
82        let n = self.symbolic.dim;
83        if a.nrows() != n || a.ncols() != n {
84            return Err(Ilu0Error::PatternMismatch);
85        }
86        let a_sym = a.symbolic();
87
88        // Stage 1: copy A's values into L and U slots.
89        for j in 0..n {
90            let a_row_idx = a_sym.row_idx_of_col_raw(j);
91            let a_col_len = a_row_idx.len();
92            let d = self.symbolic.diag_pos[j];
93            if d >= a_col_len || a_row_idx[d].zx() != j {
94                return Err(Ilu0Error::PatternMismatch);
95            }
96            let u_start = self.symbolic.u_col_ptr[j].zx();
97            let u_end = self.symbolic.u_col_ptr[j + 1].zx();
98            let l_start = self.symbolic.l_col_ptr[j].zx();
99            let l_end = self.symbolic.l_col_ptr[j + 1].zx();
100            if u_end - u_start != d + 1 || l_end - l_start != a_col_len - d {
101                return Err(Ilu0Error::PatternMismatch);
102            }
103
104            let a_val = a.val_of_col(j);
105            for (k, val) in a_val[..=d].iter().enumerate() {
106                self.u_values[u_start + k] = copy(val);
107            }
108            self.l_values[l_start] = one::<T>();
109            for (k, val) in a_val[d + 1..].iter().enumerate() {
110                self.l_values[l_start + 1 + k] = copy(val);
111            }
112        }
113
114        // Stage 2: ILU(0) elimination — column-by-column, left-looking.
115        // Borrow the workspaces and value buffers separately to avoid aliasing.
116        let symbolic = &self.symbolic;
117        let w = self.workspace_w.as_mut_slice();
118        let marker = self.workspace_marker.as_mut_slice();
119        let l_values = self.l_values.as_mut_slice();
120        let u_values = self.u_values.as_mut_slice();
121
122        for j in 0..n {
123            let l_start = symbolic.l_col_ptr[j].zx();
124            let l_end = symbolic.l_col_ptr[j + 1].zx();
125            let u_start = symbolic.u_col_ptr[j].zx();
126            let u_end = symbolic.u_col_ptr[j + 1].zx();
127
128            // Scatter U-column and L-column (skip L's unit diagonal slot)
129            // into the dense working vector, marking the touched rows.
130            for (raw_row, val) in symbolic.u_row_idx[u_start..u_end]
131                .iter()
132                .zip(u_values[u_start..u_end].iter())
133            {
134                let i = raw_row.zx();
135                marker[i] = j;
136                w[i] = copy(val);
137            }
138            for (raw_row, val) in symbolic.l_row_idx[l_start + 1..l_end]
139                .iter()
140                .zip(l_values[l_start + 1..l_end].iter())
141            {
142                let i = raw_row.zx();
143                marker[i] = j;
144                w[i] = copy(val);
145            }
146
147            // For each row p < j present in U's column j, subtract L[:,p] * U[p,j]
148            // from w[:] — but only at positions present in pattern(A col j).
149            // U's diagonal sits at position u_end - 1; everything before it is p < j.
150            for raw_p in &symbolic.u_row_idx[u_start..u_end - 1] {
151                let p = raw_p.zx();
152                let u_pj = copy(&w[p]);
153                let lp_start = symbolic.l_col_ptr[p].zx();
154                let lp_end = symbolic.l_col_ptr[p + 1].zx();
155                // Skip the unit diagonal slot at lp_start (row p).
156                for (raw_row, lv) in symbolic.l_row_idx[lp_start + 1..lp_end]
157                    .iter()
158                    .zip(l_values[lp_start + 1..lp_end].iter())
159                {
160                    let i = raw_row.zx();
161                    if marker[i] == j {
162                        let l_ip = copy(lv);
163                        let upd = mul(&l_ip, &u_pj);
164                        w[i] = sub(&w[i], &upd);
165                    }
166                }
167            }
168
169            // Gather updated U-column.
170            for (raw_row, dst) in symbolic.u_row_idx[u_start..u_end]
171                .iter()
172                .zip(u_values[u_start..u_end].iter_mut())
173            {
174                let i = raw_row.zx();
175                *dst = copy(&w[i]);
176            }
177
178            // Pivot is U[j,j] which sits at the last U slot of column j.
179            let pivot = copy(&u_values[u_end - 1]);
180            if abs2(&pivot) == zero::<T::Real>() {
181                return Err(Ilu0Error::ZeroPivot { col: j });
182            }
183            let pivot_inv = recip(&pivot);
184
185            // Gather updated L-column (divided by pivot), keeping the unit
186            // diagonal at l_start.
187            l_values[l_start] = one::<T>();
188            for (raw_row, dst) in symbolic.l_row_idx[l_start + 1..l_end]
189                .iter()
190                .zip(l_values[l_start + 1..l_end].iter_mut())
191            {
192                let i = raw_row.zx();
193                *dst = mul(&w[i], &pivot_inv);
194            }
195        }
196
197        Ok(())
198    }
199
200    /// Construct a [`SparseColMatRef`] view over the `L` factor (unit lower
201    /// triangular, diagonal stored first in each column with value `1`).
202    #[inline]
203    pub fn l_view(&self) -> SparseColMatRef<'_, I, T> {
204        let symbolic = unsafe {
205            faer::sparse::SymbolicSparseColMatRef::<'_, I>::new_unchecked(
206                self.symbolic.dim,
207                self.symbolic.dim,
208                &self.symbolic.l_col_ptr,
209                None,
210                &self.symbolic.l_row_idx,
211            )
212        };
213        SparseColMatRef::new(symbolic, &self.l_values)
214    }
215
216    /// Construct a [`SparseColMatRef`] view over the `U` factor (upper
217    /// triangular, diagonal stored last in each column).
218    #[inline]
219    pub fn u_view(&self) -> SparseColMatRef<'_, I, T> {
220        let symbolic = unsafe {
221            faer::sparse::SymbolicSparseColMatRef::<'_, I>::new_unchecked(
222                self.symbolic.dim,
223                self.symbolic.dim,
224                &self.symbolic.u_col_ptr,
225                None,
226                &self.symbolic.u_row_idx,
227            )
228        };
229        SparseColMatRef::new(symbolic, &self.u_values)
230    }
231}
faer_precond/ilu0/numeric.rs

faer_precond/ilu0/
numeric.rs
