From 249866f9e4c1246a9cb8b3814bc1ea974cee6f1c Mon Sep 17 00:00:00 2001 From: tachsin Date: Mon, 5 Oct 2026 10:38:46 +0300 Subject: [PATCH] docs(fdfd): parallel ILU(0) measured on our matrices, both halves, and not used --- ROADMAP.md | 2 +- docs/methods/fdfd-3d.md | 25 +++++++++ src/fdfd/krylov.rs | 86 +++++++++++++++++++++++++++++- src/fdfd/three/iterative.rs | 15 ++++++ src/fdfd/three/tests.rs | 101 ++++++++++++++++++++++++++++++++++++ 5 files changed, 227 insertions(+), 2 deletions(-) diff --git a/ROADMAP.md b/ROADMAP.md index 372c6f7..39babf9 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -211,7 +211,7 @@ Measured on 0.4.0 (Core Ultra 7 265K): most mode solvers run slower on 20 thread - [ ] **QMR's kernels:** a persistent thread pool, no allocation inside the iteration, fused vector updates, dot products over fixed-size chunks summed in order, then a matrix-free Yee operator (the stencil Malas 2016 optimize). *(Done but the last: the vector work in two fused passes on rayon's threads, into vectors made once, its sums over fixed chunks: 3.2 times faster an iteration on 20 threads, 13 ns per unknown on the 40³ guide, the same bits on any number of threads. The matrix-free operator, if it pays, is still to come.)* - [x] **Symmetric QMR** for the complex-symmetric curl-curl system (Freund 1992; COCG, van der Vorst & Melissen 1990, to compare): one product per iteration and no Aᵀ. *(On the curl-curl operator's diagonal similarity S A S⁻¹, which is complex symmetric: twice as fast on the 40³ guide for the same field, COCG no better. No preconditioner moves to it: Shin and Fan's operator, on which ILU(0) and the multigrid work, isn't symmetrizable by any diagonal (measured), and ILU(0) fails on the curl-curl one, so an incomplete LDLᵀ has nothing to precondition. See [FDFD in 3D](docs/methods/fdfd-3d.md#the-iterative-solver).)* - [x] **Nested dissection** (George 1973) as the ordering of faer's sparse LU on our structured grids. *(The 3D direct solver: twice as fast at 40³, 33 s against 62 s, with a quarter less memory, its separators two steps wide since faer pivots rows. In 2D COLAMD wins, and the 2D solver, the mode solvers and the circuits keep it. See [FDFD in 3D](docs/methods/fdfd-3d.md#cost).)* -- [ ] **Deterministic parallel ILU:** the factorization by fixed synchronous sweeps (Chow & Patel 2015) and the triangular solves by Jacobi sweeps (Anzt 2015), the same bits on any thread count. +- [x] **Deterministic parallel ILU:** the factorization by fixed synchronous sweeps (Chow & Patel 2015) and the triangular solves by Jacobi sweeps (Anzt 2015), the same bits on any thread count. *(Measured, and not used: the factorization isn't where the time goes, and swept triangular solves made QMR + ILU(0) 3.5 times slower at best on the guide, the sweeps' non-normal growth on an indefinite matrix's factors outrunning their parallelism. See [FDFD in 3D](docs/methods/fdfd-3d.md#preconditioning-qmr), "What didn't pay".)* - [ ] **Validation:** every existing case unchanged; each parallel kernel bit-identical on 1, 4 and 20 threads; symmetric QMR against the direct solver to 1e-10; nested dissection's solutions equal to round-off. ### 0.5: Finite-difference time-domain (FDTD) diff --git a/docs/methods/fdfd-3d.md b/docs/methods/fdfd-3d.md index 0370636..6fe4f72 100644 --- a/docs/methods/fdfd-3d.md +++ b/docs/methods/fdfd-3d.md @@ -28,6 +28,10 @@ papers: doi: 10.1109/JLT.2002.800371 - cite: "A. George, SIAM J. Numer. Anal. 10, 345 (1973) (nested dissection)" doi: 10.1137/0710032 + - cite: "E. Chow, A. Patel, SIAM J. Sci. Comput. 37, C169 (2015) (parallel ILU, measured and not used)" + doi: 10.1137/140968896 + - cite: "H. Anzt, E. Chow, J. Dongarra, Euro-Par 2015, LNCS 9233, 650 (iterative triangular solves, measured and not used)" + doi: 10.1007/978-3-662-48096-0_50 - cite: "Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd ed., SIAM (2003) (ILU(0) and GMRES)" doi: 10.1137/1.9780898718003 - cite: "B. Reps, W. Vanroose, H. bin Zubair, J. Comput. Phys. 229, 8384 (2010) (complex-stretched layers, and the multigrid cycle)" @@ -666,6 +670,27 @@ further for the same field. Plain QMR on the guide ran 35.2 s in this run, again **What didn't pay.** Jacobi (the diagonal): 10 798 iterations instead of 1 976 on the curl-curl operator, 1 335 instead of 1 629 on Shin and Fan's (40³ guide, 1e-6). +ILU(0) in parallel, both halves. **Its factorization** by Chow and Patel's fixed-point sweeps +(SIAM J. Sci. Comput. 37, C169 (2015), doi:10.1137/140968896) wasn't built: measured, it isn't +where the time goes. The multigrid's smoothers already factorize slab by slab in parallel, about +1.3 s of Diel's hierarchy, and QMR + ILU(0) on the guide spends 4.0 of its 5.6 s in its +iterations, not its factorization; and Chow and Patel find that matrices far from diagonally +dominant need 3 to 5 synchronous sweeps (their Section 4.2, Table 3). **Its triangular solves** +by k Jacobi sweeps each (Anzt, Chow and Dongarra, Euro-Par 2015, LNCS 9233, 650, +doi:10.1007/978-3-662-48096-0_50: the truncated Neumann series of the factor, all rows at once, +the same bits on any number of threads, and on the transposed factors exactly its transpose, so +QMR stays consistent), on the guide to 1e-8, 20 threads (`ilu_sweeps` in src/fdfd/three/tests.rs): + +| Triangular solves | Exact | 1 sweep | 2 | 3 | 4 | 6 | 8 | +|---|---|---|---|---|---|---|---| +| QMR iterations | 195 | 1 163 | 747 | 730 | 929 | 1 782 | 2 322 | +| Time | 3.3 s | 12.9 s | 10.4 s | 11.5 s | 15.5 s | 32.7 s | 48.1 s | + +Three and a half times slower at best, and worse with more sweeps: the factors of an indefinite +matrix with PMLs are far from diagonally dominant, and the sweeps' non-normal iteration grows +before it converges, as Anzt et al. warn (their Section 1). Neither half is used; the exact +triangular solves stay, sequential. + A first multigrid, before this one: damped Jacobi smoothing and coarse grids rediscretized over the same box, with plain PMLs. As a solver it cut the residual 3 to 6 times per cycle in vacuum or on a silicon guide between walls; with PMLs it stalled at 0.96 per cycle, or diverged even when diff --git a/src/fdfd/krylov.rs b/src/fdfd/krylov.rs index a5c2fde..45b83bb 100644 --- a/src/fdfd/krylov.rs +++ b/src/fdfd/krylov.rs @@ -1112,6 +1112,37 @@ impl Triangle { } } + /// `sweeps` synchronous Jacobi sweeps from x = 0 for T x = `b`, x ← D⁻¹(b − (T − D) x), each + /// row from the previous sweep's values only, all rows at once on rayon's threads: the same + /// bits on any number of them. The truncated Neumann series Σ_{j Vec { + use rayon::prelude::*; + let scale = |i: usize, s: c64| match &self.inverse_diagonal { + Some(d) => s * d[i], + None => s, + }; + let mut x: Vec = (0..b.len()).map(|i| scale(i, b[i])).collect(); + let mut next = vec![c64::new(0.0, 0.0); b.len()]; + for _ in 1..sweeps { + next.par_iter_mut() + .with_min_len(4096) + .enumerate() + .for_each(|(i, out)| { + let mut s = b[i]; + for k in self.starts[i]..self.starts[i + 1] { + s -= self.values[k] * x[self.columns[k]]; + } + *out = scale(i, s); + }); + std::mem::swap(&mut x, &mut next); + } + x + } + /// Solves in place: `x` holds the right-hand side on entry. fn solve(&self, x: &mut [c64]) { let n = x.len(); @@ -1143,6 +1174,8 @@ pub(crate) struct Ilu0 { u: Triangle, ut: Triangle, lt: Triangle, + /// Each triangular solve as this many Jacobi sweeps ([`Triangle::sweeps`]), or exact. + sweeps: Option, } impl Ilu0 { @@ -1221,12 +1254,30 @@ impl Ilu0 { ); let ut = Triangle::new(n, |i| upper_t[i].clone(), Some(inverse), true); let lt = Triangle::new(n, |i| lower_t[i].clone(), None, false); - Ok(Ilu0 { l, u, ut, lt }) + Ok(Ilu0 { + l, + u, + ut, + lt, + sweeps: None, + }) + } + + /// The same factors, each triangular solve taken as `sweeps` Jacobi sweeps instead of + /// exactly: in parallel where the exact solves are sequential, and a fixed linear operator, + /// so QMR can use it as its preconditioner. + #[cfg(test)] + pub(crate) fn with_sweeps(mut self, sweeps: usize) -> Ilu0 { + self.sweeps = Some(sweeps.max(1)); + self } } impl Preconditioner for Ilu0 { fn solve(&self, v: &[c64]) -> Vec { + if let Some(k) = self.sweeps { + return self.u.sweeps(&self.l.sweeps(v, k), k); + } let mut x = v.to_vec(); self.l.solve(&mut x); self.u.solve(&mut x); @@ -1234,6 +1285,9 @@ impl Preconditioner for Ilu0 { } fn solve_transpose(&self, v: &[c64]) -> Vec { + if let Some(k) = self.sweeps { + return self.lt.sweeps(&self.ut.sweeps(v, k), k); + } let mut x = v.to_vec(); self.ut.solve(&mut x); self.lt.solve(&mut x); @@ -1614,6 +1668,36 @@ mod tests { assert!(a.apply(&v) == sequential); } + #[test] + fn jacobi_swept_ilu_is_its_own_transpose_and_converges_to_the_exact_solves() { + // M⁻¹ applied by k sweeps and M⁻ᵀ by k sweeps on the transposed factors are each + // other's transposes, uᵀ(M⁻¹v) = (M⁻ᵀu)ᵀv, for any k; and with enough sweeps (a factor + // of n rows needs at most n) they are the exact solves + let n = 2 * CHUNK + 9; + let (a, _) = gmres_case(n); + let u: Vec = (0..n).map(|r| c64::new((r % 7) as f64, 1.0)).collect(); + let v: Vec = (0..n) + .map(|r| c64::new(1.0, (r % 5) as f64 - 2.0)) + .collect(); + for sweeps in [1, 2, 5] { + let ilu = Ilu0::new(&a).unwrap().with_sweeps(sweeps); + let (left, right) = (dot(&u, &ilu.solve(&v)), dot(&ilu.solve_transpose(&u), &v)); + assert!( + (left - right).norm() < 1e-12 * left.norm(), + "{sweeps}: {left} {right}" + ); + } + // the tridiagonal case's factors are bidiagonal: n sweeps are exact, and far fewer are + // within round-off of it for this diagonally dominant matrix + let t = tridiagonal(64); + let exact = Ilu0::new(&t).unwrap(); + let swept = Ilu0::new(&t).unwrap().with_sweeps(64); + let v: Vec = (0..64).map(|r| c64::new(1.0, r as f64)).collect(); + let (x, y) = (exact.solve(&v), swept.solve(&v)); + let d: Vec = x.iter().zip(&y).map(|(p, q)| p - q).collect(); + assert!(norm(&d) < 1e-12 * norm(&x), "{}", norm(&d) / norm(&x)); + } + #[test] fn ilu0_of_a_tridiagonal_matrix_is_its_lu() { // no fill: ILU(0) is the exact factorization, and its inverse the matrix's diff --git a/src/fdfd/three/iterative.rs b/src/fdfd/three/iterative.rs index 33dec4f..b8623ff 100644 --- a/src/fdfd/three/iterative.rs +++ b/src/fdfd/three/iterative.rs @@ -129,6 +129,21 @@ impl IterativeSolver3d { Ok(self) } + /// As [`IterativeSolver3d::with_ilu`], each of the factors' triangular solves taken as + /// `sweeps` Jacobi sweeps, in parallel, where the exact ones are sequential. + /// + /// # Errors + /// + /// As [`IterativeSolver3d::with_ilu`]. + #[cfg(test)] + pub(crate) fn with_ilu_sweeps(self, sweeps: usize) -> Result { + let mut solver = self.with_ilu()?; + if let Preconditioning::Ilu(ilu) = solver.preconditioner { + solver.preconditioner = Preconditioning::Ilu(Box::new(ilu.with_sweeps(sweeps))); + } + Ok(solver) + } + /// The same problem, solved by GMRES preconditioned from the right by a multigrid cycle on /// Shin and Fan's operator ([`Multigrid`]; GMRES, restarted, as Y. Saad, *Iterative Methods /// for Sparse Linear Systems*, 2nd ed., SIAM (2003), doi:10.1137/1.9780898718003, Algorithms diff --git a/src/fdfd/three/tests.rs b/src/fdfd/three/tests.rs index 2344ce1..cc08927 100644 --- a/src/fdfd/three/tests.rs +++ b/src/fdfd/three/tests.rs @@ -1613,3 +1613,104 @@ fn symmetric_solvers() { error(&x) ); } + +#[test] +#[ignore = "ILU(0) with triangular solves by Jacobi sweeps, for the docs: SWEEP_CASE=guide|diel cargo test --release fdfd::three::tests::ilu_sweeps -- --ignored --nocapture"] +fn ilu_sweeps() { + // the benchmark's guide (40^3 cells of 10 nm, stretched PMLs of 10) or Diel (40 x 90 x 80), + // by QMR + ILU(0) on Shin and Fan's operator to 1e-8, its triangular solves exact or by k + // Jacobi sweeps (Anzt, Chow and Dongarra 2015); the field's error against the exact solves + // to 1e-10 + use crate::fdfd::{Formulation, IterativeSolver3d, Stopping}; + use std::time::Instant; + let case = std::env::var("SWEEP_CASE").unwrap_or_else(|_| "guide".into()); + let (h, pml) = (0.01, 10); + let (n, (wy, wz)): ([usize; 3], (f64, f64)) = if case == "diel" { + ([40, 70 + 2 * pml, 60 + 2 * pml], (0.2, 0.15)) + } else { + ([20 + 2 * pml; 3], (0.05, 0.05)) + }; + let grid = Grid3d { + nx: n[0], + ny: n[1], + nz: n[2], + dx: h, + dy: h, + dz: h, + x0: -(n[0] as f64) * h / 2.0, + y0: -(n[1] as f64) * h / 2.0, + z0: -(n[2] as f64) * h / 2.0, + }; + let eps = move |_: f64, y: f64, z: f64| { + c64::new( + if y.abs() < wy && z.abs() < wz { + 12.09 + } else { + 1.0 + }, + 0.0, + ) + }; + let lam = Wavelength::um(1.55).unwrap(); + let mut source = vec![c64::new(0.0, 0.0); grid.unknowns()]; + if case == "diel" { + for k in 0..grid.nz { + for j in 0..grid.ny { + let [_, y, z] = grid.e_position(Axis::Y, (15, j, k)); + if y.abs() < wy && z.abs() < wz { + source[grid.index(Axis::Y, (15, j, k))] = c64::new(1.0, 0.0); + } + } + } + } else { + source[grid.index(Axis::X, (n[0] / 2 + 3, n[1] / 2 + 3, n[2] / 2 + 3))] = + c64::new(1.0, 0.0); + } + let new = || { + IterativeSolver3d::new( + grid, + lam, + eps, + Boundaries3d::stretched_pml(pml), + Formulation::ShinFan, + ) + .unwrap() + }; + let stop = |tolerance| Stopping { + tolerance, + max_iterations: 20_000, + }; + let exact = new().with_ilu().unwrap(); + let (reference, _) = exact.solve(&source, stop(1e-10)).unwrap(); + let norm = |v: &[c64]| v.iter().map(|z| z.norm_sqr()).sum::().sqrt(); + let error = |x: &[c64]| { + let d: Vec = x + .iter() + .zip(reference.values()) + .map(|(p, q)| p - q) + .collect(); + norm(&d) / norm(reference.values()) + }; + let t = Instant::now(); + let (x, how) = exact.solve(&source, stop(1e-8)).unwrap(); + println!( + "SWEEP {case} exact: {} iterations, {:?}, error {:.1e}", + how.iterations, + t.elapsed(), + error(x.values()) + ); + drop(exact); + for sweeps in [1, 2, 3, 4, 6, 8] { + let solver = new().with_ilu_sweeps(sweeps).unwrap(); + let t = Instant::now(); + match solver.solve(&source, stop(1e-8)) { + Ok((x, how)) => println!( + "SWEEP {case} {sweeps} sweeps: {} iterations, {:?}, error {:.1e}", + how.iterations, + t.elapsed(), + error(x.values()) + ), + Err(e) => println!("SWEEP {case} {sweeps} sweeps: {e}"), + } + } +}