math_audio_solvers/preconditioners/
ilu.rs

1//! ILU(0) preconditioner
2//!
3//! Incomplete LU factorization with no fill-in.
4//! Effective for many sparse systems, especially those arising from discretization.
5
6use crate::sparse::CsrMatrix;
7use crate::traits::{ComplexField, Preconditioner};
8use ndarray::Array1;
9use num_traits::FromPrimitive;
10
11/// ILU(0) preconditioner
12///
13/// Computes an incomplete LU factorization where the sparsity pattern
14/// of L and U is the same as the original matrix.
15#[derive(Debug, Clone)]
16pub struct IluPreconditioner<T: ComplexField> {
17    /// Lower triangular factor (stored in CSR format conceptually)
18    l_values: Vec<T>,
19    l_col_indices: Vec<usize>,
20    l_row_ptrs: Vec<usize>,
21    /// Upper triangular factor
22    u_values: Vec<T>,
23    u_col_indices: Vec<usize>,
24    u_row_ptrs: Vec<usize>,
25    /// Diagonal of U (for fast access)
26    u_diag: Vec<T>,
27    /// Matrix dimension
28    n: usize,
29    #[allow(dead_code)]
30    /// Precomputed diagonal indices for fast lookup during factorization
31    diag_indices: Vec<usize>,
32}
33
34impl<T: ComplexField> IluPreconditioner<T> {
35    /// Create ILU(0) preconditioner from a CSR matrix
36    pub fn from_csr(matrix: &CsrMatrix<T>) -> Self {
37        let n = matrix.num_rows;
38
39        let col_indices = &matrix.col_indices;
40        let row_ptrs = &matrix.row_ptrs;
41
42        let mut diag_indices = vec![usize::MAX; n];
43        for i in 0..n {
44            #[allow(clippy::needless_range_loop)]
45            for idx in row_ptrs[i]..row_ptrs[i + 1] {
46                if col_indices[idx] == i {
47                    diag_indices[i] = idx;
48                    break;
49                }
50            }
51        }
52
53        let mut values = matrix.values.clone();
54
55        for i in 0..n {
56            for idx in row_ptrs[i]..row_ptrs[i + 1] {
57                let k = col_indices[idx];
58                if k >= i {
59                    break;
60                }
61
62                let u_kk_idx = diag_indices[k];
63                if u_kk_idx == usize::MAX {
64                    continue;
65                }
66
67                let u_kk = values[u_kk_idx];
68
69                if u_kk.norm() < T::Real::from_f64(1e-30).unwrap() {
70                    continue;
71                }
72
73                let l_ik = values[idx] * u_kk.inv();
74                values[idx] = l_ik;
75
76                for j_idx in row_ptrs[i]..row_ptrs[i + 1] {
77                    let j = col_indices[j_idx];
78                    if j <= k {
79                        continue;
80                    }
81
82                    let u_kj_idx = diag_indices[k] + 1;
83                    if u_kj_idx < row_ptrs[k + 1] && col_indices[u_kj_idx] == j {
84                        values[j_idx] = values[j_idx] - l_ik * values[u_kj_idx];
85                    } else {
86                        for search_idx in (row_ptrs[k] + 1)..row_ptrs[k + 1] {
87                            if col_indices[search_idx] == j {
88                                values[j_idx] = values[j_idx] - l_ik * values[search_idx];
89                                break;
90                            }
91                        }
92                    }
93                }
94            }
95        }
96
97        let mut l_values = Vec::new();
98        let mut l_col_indices = Vec::new();
99        let mut l_row_ptrs = vec![0];
100
101        let mut u_values = Vec::new();
102        let mut u_col_indices = Vec::new();
103        let mut u_row_ptrs = vec![0];
104        let mut u_diag = vec![T::one(); n];
105
106        for i in 0..n {
107            for idx in row_ptrs[i]..row_ptrs[i + 1] {
108                let j = col_indices[idx];
109                let val = values[idx];
110
111                if j < i {
112                    l_values.push(val);
113                    l_col_indices.push(j);
114                } else {
115                    u_values.push(val);
116                    u_col_indices.push(j);
117                    if j == i {
118                        u_diag[i] = val;
119                    }
120                }
121            }
122            l_row_ptrs.push(l_values.len());
123            u_row_ptrs.push(u_values.len());
124        }
125
126        Self {
127            l_values,
128            l_col_indices,
129            l_row_ptrs,
130            u_values,
131            u_col_indices,
132            u_row_ptrs,
133            u_diag,
134            n,
135            diag_indices,
136        }
137    }
138}
139
140impl<T: ComplexField> Preconditioner<T> for IluPreconditioner<T> {
141    fn apply(&self, r: &Array1<T>) -> Array1<T> {
142        let mut y = r.clone();
143
144        // Forward substitution: Ly = r (L has unit diagonal)
145        for i in 0..self.n {
146            for idx in self.l_row_ptrs[i]..self.l_row_ptrs[i + 1] {
147                let j = self.l_col_indices[idx];
148                let l_ij = self.l_values[idx];
149                y[i] = y[i] - l_ij * y[j];
150            }
151        }
152
153        // Backward substitution: Ux = y
154        let mut x = y;
155        for i in (0..self.n).rev() {
156            for idx in self.u_row_ptrs[i]..self.u_row_ptrs[i + 1] {
157                let j = self.u_col_indices[idx];
158                if j > i {
159                    let u_ij = self.u_values[idx];
160                    x[i] = x[i] - u_ij * x[j];
161                }
162            }
163
164            let u_ii = self.u_diag[i];
165            if u_ii.norm() > T::Real::from_f64(1e-30).unwrap() {
166                x[i] *= u_ii.inv();
167            }
168        }
169
170        x
171    }
172}
173
174#[cfg(test)]
175mod tests {
176    use super::*;
177    use crate::iterative::{GmresConfig, gmres, gmres_preconditioned};
178    use ndarray::array;
179    use num_complex::Complex64;
180
181    #[test]
182    fn test_ilu_preconditioner() {
183        // Simple tridiagonal matrix
184        let dense = array![
185            [
186                Complex64::new(4.0, 0.0),
187                Complex64::new(-1.0, 0.0),
188                Complex64::new(0.0, 0.0)
189            ],
190            [
191                Complex64::new(-1.0, 0.0),
192                Complex64::new(4.0, 0.0),
193                Complex64::new(-1.0, 0.0)
194            ],
195            [
196                Complex64::new(0.0, 0.0),
197                Complex64::new(-1.0, 0.0),
198                Complex64::new(4.0, 0.0)
199            ],
200        ];
201
202        let matrix = CsrMatrix::from_dense(&dense, 1e-15);
203        let precond = IluPreconditioner::from_csr(&matrix);
204
205        let r = array![
206            Complex64::new(1.0, 0.0),
207            Complex64::new(2.0, 0.0),
208            Complex64::new(3.0, 0.0)
209        ];
210
211        let result = precond.apply(&r);
212
213        // Verify that M*result ≈ r (where M = L*U)
214        // For ILU(0) on this matrix, M should be close to A
215        let check = matrix.matvec(&result);
216        for i in 0..3 {
217            // Allow some error since ILU is approximate
218            assert!(
219                (check[i] - r[i]).norm() < 0.5,
220                "ILU should approximately invert: got {:?} expected {:?}",
221                check[i],
222                r[i]
223            );
224        }
225    }
226
227    #[test]
228    fn test_ilu_with_gmres() {
229        // Larger system where preconditioning helps
230        let n = 10;
231        let mut dense = ndarray::Array2::from_elem((n, n), Complex64::new(0.0, 0.0));
232
233        // Tridiagonal matrix
234        for i in 0..n {
235            dense[[i, i]] = Complex64::new(4.0, 0.0);
236            if i > 0 {
237                dense[[i, i - 1]] = Complex64::new(-1.0, 0.0);
238            }
239            if i < n - 1 {
240                dense[[i, i + 1]] = Complex64::new(-1.0, 0.0);
241            }
242        }
243
244        let matrix = CsrMatrix::from_dense(&dense, 1e-15);
245        let precond = IluPreconditioner::from_csr(&matrix);
246
247        let b = Array1::from_iter((0..n).map(|i| Complex64::new((i as f64).sin(), 0.0)));
248
249        let config = GmresConfig {
250            max_iterations: 50,
251            restart: 10,
252            tolerance: 1e-10,
253            print_interval: 0,
254        };
255
256        // Solve without preconditioning
257        let sol_no_precond = gmres(&matrix, &b, &config);
258
259        // Solve with preconditioning
260        let sol_precond = gmres_preconditioned(&matrix, &precond, &b, &config);
261
262        // Both should converge
263        assert!(sol_no_precond.converged);
264        assert!(sol_precond.converged);
265
266        // Preconditioned should need fewer iterations (usually)
267        // For this well-conditioned problem, difference may be small
268        assert!(
269            sol_precond.iterations <= sol_no_precond.iterations + 5,
270            "Preconditioning should not hurt: {} vs {}",
271            sol_precond.iterations,
272            sol_no_precond.iterations
273        );
274    }
275}
math_audio_solvers/preconditioners/ilu.rs

math_audio_solvers/preconditioners/
ilu.rs
