docs(fdfd): parallel ILU(0) measured on our matrices, both halves, and not used

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)

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

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<sweeps} Mʲ D⁻¹, M = I −
+    /// D⁻¹T strictly triangular (H. Anzt, E. Chow, J. Dongarra, "Iterative sparse triangular
+    /// solves for preconditioning", Euro-Par 2015, LNCS 9233, 650, doi:10.1007/978-3-662-48096-0_50,
+    /// their Eqs. 1–2). On the transposed factor it is exactly the transpose of the same series,
+    /// D⁻¹(I − TᵀD⁻¹)ʲ = (I − D⁻¹Tᵀ)ʲD⁻¹, so QMR's two products stay each other's transposes.
+    fn sweeps(&self, b: &[c64], sweeps: usize) -> Vec<c64> {
+        use rayon::prelude::*;
+        let scale = |i: usize, s: c64| match &self.inverse_diagonal {
+            Some(d) => s * d[i],
+            None => s,
+        };
+        let mut x: Vec<c64> = (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<usize>,
 }
 
 impl Ilu0 {
@@ -1221,19 +1254,40 @@ 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<c64> {
+        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);
         x
     }
 
     fn solve_transpose(&self, v: &[c64]) -> Vec<c64> {
+        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<c64> = (0..n).map(|r| c64::new((r % 7) as f64, 1.0)).collect();
+        let v: Vec<c64> = (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<c64> = (0..64).map(|r| c64::new(1.0, r as f64)).collect();
+        let (x, y) = (exact.solve(&v), swept.solve(&v));
+        let d: Vec<c64> = 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

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<IterativeSolver3d> {
+        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

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::<f64>().sqrt();
+    let error = |x: &[c64]| {
+        let d: Vec<c64> = 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}"),
+        }
+    }
+}