From 88dc157c518d8ca2f5a5d89189b45b5e7ac40fd3 Mon Sep 17 00:00:00 2001 From: tachsin Date: Mon, 5 Oct 2026 11:30:25 +0300 Subject: [PATCH] feat(bench): the bandwidth each iterative solve reaches, its bytes counted by the kernels, against the triad --- ROADMAP.md | 2 +- docs/benchmarks.md | 80 +++++++++++++++++------------------ src/bench.rs | 32 +++++++++++++- src/fdfd/krylov.rs | 33 +++++++++++++++ src/fdfd/three/multigrid.rs | 5 ++- src/lib.rs | 1 + src/traffic.rs | 69 ++++++++++++++++++++++++++++++ studio/src-tauri/src/bench.rs | 28 +++++++++--- 8 files changed, 200 insertions(+), 50 deletions(-) create mode 100644 src/traffic.rs diff --git a/ROADMAP.md b/ROADMAP.md index 4aaa746..2880272 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -206,7 +206,7 @@ Released without the preconditioner, which moved to 0.4.2. Measured on 0.4.0 (Core Ultra 7 265K): most mode solvers run slower on 20 threads than on one (Hadley's corners 13.2 s against 8.0 s), an FDFD job's wavelengths are solved one after another, and the 3D iterative solver reaches about a fifth of the memory bandwidth (an estimate). See [the performance plan](docs/plans/performance.md), Phase A. -- [ ] **A benchmark harness first:** a `bench` command with fixed problems (2D FDFD, the 3D silicon guide, Diel, the strip with ports, the mode-solver examples), recording wall time, iterations, the field's error against a converged reference, peak memory, and bandwidth against the machine's measured roofline (Williams 2009). The same matrices exported and solved by PARDISO (Schenk & Gärtner 2004) and MUMPS (Amestoy 2001) as external programs: the MKL-class baseline, never linked. *(Begun: `photonoxide bench` times the fixed problems of `photonoxide::bench` and the slowest examples, each in a process of its own, with phases, iterations, ns per unknown per iteration, the error against an exact answer or a reference, peak memory and the STREAM triad: [docs/benchmarks.md](docs/benchmarks.md). Still to come: the bandwidth each solve reached, which needs QMR's byte counts per iteration (with its kernels, below), and the exported matrices, after the licence check of the performance plan's decision 4.)* +- [ ] **A benchmark harness first:** a `bench` command with fixed problems (2D FDFD, the 3D silicon guide, Diel, the strip with ports, the mode-solver examples), recording wall time, iterations, the field's error against a converged reference, peak memory, and bandwidth against the machine's measured roofline (Williams 2009). The same matrices exported and solved by PARDISO (Schenk & Gärtner 2004) and MUMPS (Amestoy 2001) as external programs: the MKL-class baseline, never linked. *(Begun: `photonoxide bench` times the fixed problems of `photonoxide::bench` and the slowest examples, each in a process of its own, with phases, iterations, ns per unknown per iteration, the error against an exact answer or a reference, peak memory and the STREAM triad: [docs/benchmarks.md](docs/benchmarks.md). And the bandwidth each iterative solve reached, its bytes counted by the kernels themselves (src/traffic.rs), against the triad: the multigrid at 63 to 68% on Diel and the strip with ports on 20 threads, QMR + ILU(0) at 39% (its sequential triangular solves); on the 40³ guide QMR's vectors fit in the L3 cache and it counts above the roof. Still to come: the exported matrices, after the licence check of the performance plan's decision 4.)* - [x] **Sweeps side by side:** the points of an FDFD wavelength sweep and of the mode solvers' sweeps solved in parallel, each on one thread when there are enough points, collected in order, so events and results are unchanged. *(The jobs' sweeps and `fdfd_s_parameters`, 3 to 6 times faster on 20 threads (#146); the examples that sweep by hand through `parallel::map_in_order`, `hadley_corners` 9.4 s to 2.5 s and `group_index` 13.1 s to 3.0 s, their output unchanged. Inside a point faer shares rayon's pool, so the points don't oversubscribe it.)* - [ ] **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).)* diff --git a/docs/benchmarks.md b/docs/benchmarks.md index 01090dc..a023efd 100644 --- a/docs/benchmarks.md +++ b/docs/benchmarks.md @@ -1,51 +1,51 @@ # Benchmarks -What `photonoxide bench` measured on one machine: the fixed problems of `photonoxide::bench` and the slowest examples, each in a process of its own. Times depend on the machine and on what else ran on it; the errors don't. Peak memory is the process's peak resident set (its peak working set on Windows); ns per unknown per iteration is the QMR phases' time over the unknowns and the iterations. See the performance plan, docs/plans/performance.md. +What `photonoxide bench` measured on one machine: the fixed problems of `photonoxide::bench` and the slowest examples, each in a process of its own. Times depend on the machine and on what else ran on it; the errors don't. Peak memory is the process's peak resident set (its peak working set on Windows); ns per unknown per iteration is the QMR phases' time over the unknowns and the iterations. Bandwidth is the iterating phases' bytes over their time, the bytes counted by the kernels themselves (each nonzero's value and index, each row pointer, each vector read or written once: src/traffic.rs), from memory or from cache; beside it, its share of the triad on the same threads, the roof (Williams et al., Commun. ACM 52(4), 65 (2009)). A problem whose vectors fit in the last-level cache can count above it: the 40³ guide's do. See the performance plan, docs/plans/performance.md. - photonoxide 0.4.1, Intel(R) Core(TM) Ultra 7 265K, 20 logical processors, Windows -- memory bandwidth, the STREAM triad over three arrays of 128 MiB: 33.2 GB/s on 1 thread, 52.5 GB/s on 20 threads +- memory bandwidth, the STREAM triad over three arrays of 128 MiB: 32.9 GB/s on 1 thread, 56.3 GB/s on 20 threads -| Problem | Grid | Unknowns | Threads | Time | Iterations | ns/unknown/iteration | Error | Peak memory | -|---|---|---:|---:|---:|---:|---:|---|---:| -| `fdfd2d/slab-lu` | 440 × 340 cells of 10 nm, PMLs of 20 | 149 600 | 1 | 0.90 s | — | — | 2.6e-13 (S21 = exp(iβL), S11 = 0) | 690 MB | -| `fdfd3d/guide-qmr` | 40 × 40 × 40 cells of 10 nm, PMLs of 10 | 192 000 | 1 | 12.17 s | 2 077 | 27.3 | 1.5e-7 (QMR + ILU(0) to 1e-12, relative) | 320 MB | -| `fdfd3d/guide-ilu` | 40 × 40 × 40 cells of 10 nm, stretched PMLs of 10 | 192 000 | 1 | 5.82 s | 195 | 116.6 | 1.2e-6 (GMRES + multigrid to 1e-12, relative) | 460 MB | -| `fdfd3d/guide-multigrid` | 40 × 40 × 40 cells of 10 nm, stretched PMLs of 10 | 192 000 | 1 | 3.35 s | 16 | 271.4 | 1.6e-6 (GMRES + multigrid to 1e-12, relative) | 1.03 GB | -| `fdfd3d/diel-multigrid` | 40 × 90 × 80 cells of 10 nm, stretched PMLs of 10 | 864 000 | 1 | 22.61 s | 46 | 282.1 | 1.8e-6 (GMRES + multigrid to 1e-12, relative) | 4.64 GB | -| `fdfd3d/strip-ports-multigrid` | 72 × 102 × 82 cells of 20 nm, stretched PMLs of 16 | 1 806 624 | 1 | 52.65 s | 56 | 273.4 | 4.1e-9 (S21 = exp(iβL), S11 = 0) | 9.43 GB | -| `job/strip-modes` | 121 × 111 cells of 20 nm, 11 points | 27 328 | 1 | 4.17 s | — | — | — | 138 MB | -| `job/mmi-fdfd` | 725 × 300 cells of 20 nm, 11 points | 217 500 | 1 | 24.42 s | — | — | — | 1.42 GB | -| `example/strip_waveguide` | examples/strip_waveguide.rs | — | 1 | 6.17 s | — | — | — | 1.45 GB | -| `example/hadley_corners` | examples/hadley_corners.rs | — | 1 | 7.49 s | — | — | — | 235 MB | -| `example/group_index` | examples/group_index.rs | — | 1 | 10.81 s | — | — | — | 258 MB | -| `example/circuit_fit` | examples/circuit_fit.rs | — | 1 | 0.46 s | — | — | — | 11 MB | -| `fdfd2d/slab-lu` | 440 × 340 cells of 10 nm, PMLs of 20 | 149 600 | 20 | 1.88 s | — | — | 2.6e-13 (S21 = exp(iβL), S11 = 0) | 691 MB | -| `fdfd3d/guide-qmr` | 40 × 40 × 40 cells of 10 nm, PMLs of 10 | 192 000 | 20 | 4.12 s | 2 077 | 7.1 | 1.5e-7 (QMR + ILU(0) to 1e-12, relative) | 320 MB | -| `fdfd3d/guide-ilu` | 40 × 40 × 40 cells of 10 nm, stretched PMLs of 10 | 192 000 | 20 | 4.76 s | 195 | 88.3 | 1.2e-6 (GMRES + multigrid to 1e-12, relative) | 460 MB | -| `fdfd3d/guide-multigrid` | 40 × 40 × 40 cells of 10 nm, stretched PMLs of 10 | 192 000 | 20 | 2.49 s | 16 | 166.7 | 1.6e-6 (GMRES + multigrid to 1e-12, relative) | 1.08 GB | -| `fdfd3d/diel-multigrid` | 40 × 90 × 80 cells of 10 nm, stretched PMLs of 10 | 864 000 | 20 | 12.64 s | 46 | 113.8 | 1.8e-6 (GMRES + multigrid to 1e-12, relative) | 4.67 GB | -| `fdfd3d/strip-ports-multigrid` | 72 × 102 × 82 cells of 20 nm, stretched PMLs of 16 | 1 806 624 | 20 | 29.60 s | 56 | 115.4 | 4.1e-9 (S21 = exp(iβL), S11 = 0) | 9.46 GB | -| `job/strip-modes` | 121 × 111 cells of 20 nm, 11 points | 27 328 | 20 | 1.16 s | — | — | — | 1.09 GB | -| `job/mmi-fdfd` | 725 × 300 cells of 20 nm, 11 points | 217 500 | 20 | 8.50 s | — | — | — | 6.90 GB | -| `example/strip_waveguide` | examples/strip_waveguide.rs | — | 20 | 5.60 s | — | — | — | 1.46 GB | -| `example/hadley_corners` | examples/hadley_corners.rs | — | 20 | 2.49 s | — | — | — | 732 MB | -| `example/group_index` | examples/group_index.rs | — | 20 | 2.98 s | — | — | — | 1.01 GB | -| `example/circuit_fit` | examples/circuit_fit.rs | — | 20 | 0.42 s | — | — | — | 11 MB | +| Problem | Grid | Unknowns | Threads | Time | Iterations | ns/unknown/iteration | Bandwidth (of the triad) | Error | Peak memory | +|---|---|---:|---:|---:|---:|---:|---:|---|---:| +| `fdfd2d/slab-lu` | 440 × 340 cells of 10 nm, PMLs of 20 | 149 600 | 1 | 0.92 s | — | — | — | 2.6e-13 (S21 = exp(iβL), S11 = 0) | 690 MB | +| `fdfd3d/guide-qmr` | 40 × 40 × 40 cells of 10 nm, PMLs of 10 | 192 000 | 1 | 13.39 s | 2 077 | 30.3 | 19.4 GB/s (59%) | 1.5e-7 (QMR + ILU(0) to 1e-12, relative) | 320 MB | +| `fdfd3d/guide-ilu` | 40 × 40 × 40 cells of 10 nm, stretched PMLs of 10 | 192 000 | 1 | 6.27 s | 195 | 126.0 | 15.3 GB/s (46%) | 1.2e-6 (GMRES + multigrid to 1e-12, relative) | 460 MB | +| `fdfd3d/guide-multigrid` | 40 × 40 × 40 cells of 10 nm, stretched PMLs of 10 | 192 000 | 1 | 3.52 s | 16 | 288.5 | 12.7 GB/s (39%) | 1.6e-6 (GMRES + multigrid to 1e-12, relative) | 1.03 GB | +| `fdfd3d/diel-multigrid` | 40 × 90 × 80 cells of 10 nm, stretched PMLs of 10 | 864 000 | 1 | 23.72 s | 46 | 303.6 | 13.9 GB/s (42%) | 1.8e-6 (GMRES + multigrid to 1e-12, relative) | 4.65 GB | +| `fdfd3d/strip-ports-multigrid` | 72 × 102 × 82 cells of 20 nm, stretched PMLs of 16 | 1 806 624 | 1 | 55.66 s | 56 | 296.6 | 13.3 GB/s (40%) | 4.1e-9 (S21 = exp(iβL), S11 = 0) | 9.39 GB | +| `job/strip-modes` | 121 × 111 cells of 20 nm, 11 points | 27 328 | 1 | 4.19 s | — | — | — | — | 138 MB | +| `job/mmi-fdfd` | 725 × 300 cells of 20 nm, 11 points | 217 500 | 1 | 24.63 s | — | — | — | — | 1.42 GB | +| `example/strip_waveguide` | examples/strip_waveguide.rs | — | 1 | 6.09 s | — | — | — | — | 1.45 GB | +| `example/hadley_corners` | examples/hadley_corners.rs | — | 1 | 7.51 s | — | — | — | — | 235 MB | +| `example/group_index` | examples/group_index.rs | — | 1 | 10.77 s | — | — | — | — | 258 MB | +| `example/circuit_fit` | examples/circuit_fit.rs | — | 1 | 0.47 s | — | — | — | — | 11 MB | +| `fdfd2d/slab-lu` | 440 × 340 cells of 10 nm, PMLs of 20 | 149 600 | 20 | 1.02 s | — | — | — | 2.6e-13 (S21 = exp(iβL), S11 = 0) | 692 MB | +| `fdfd3d/guide-qmr` | 40 × 40 × 40 cells of 10 nm, PMLs of 10 | 192 000 | 20 | 4.07 s | 2 077 | 7.0 | 84.4 GB/s (150%) | 1.5e-7 (QMR + ILU(0) to 1e-12, relative) | 320 MB | +| `fdfd3d/guide-ilu` | 40 × 40 × 40 cells of 10 nm, stretched PMLs of 10 | 192 000 | 20 | 4.82 s | 195 | 89.2 | 21.6 GB/s (38%) | 1.2e-6 (GMRES + multigrid to 1e-12, relative) | 460 MB | +| `fdfd3d/guide-multigrid` | 40 × 40 × 40 cells of 10 nm, stretched PMLs of 10 | 192 000 | 20 | 2.52 s | 16 | 169.6 | 21.6 GB/s (38%) | 1.6e-6 (GMRES + multigrid to 1e-12, relative) | 1.07 GB | +| `fdfd3d/diel-multigrid` | 40 × 90 × 80 cells of 10 nm, stretched PMLs of 10 | 864 000 | 20 | 12.63 s | 46 | 113.6 | 37.1 GB/s (66%) | 1.8e-6 (GMRES + multigrid to 1e-12, relative) | 4.67 GB | +| `fdfd3d/strip-ports-multigrid` | 72 × 102 × 82 cells of 20 nm, stretched PMLs of 16 | 1 806 624 | 20 | 29.45 s | 56 | 114.6 | 34.3 GB/s (61%) | 4.1e-9 (S21 = exp(iβL), S11 = 0) | 9.51 GB | +| `job/strip-modes` | 121 × 111 cells of 20 nm, 11 points | 27 328 | 20 | 1.27 s | — | — | — | — | 1.09 GB | +| `job/mmi-fdfd` | 725 × 300 cells of 20 nm, 11 points | 217 500 | 20 | 9.28 s | — | — | — | — | 6.02 GB | +| `example/strip_waveguide` | examples/strip_waveguide.rs | — | 20 | 5.70 s | — | — | — | — | 1.46 GB | +| `example/hadley_corners` | examples/hadley_corners.rs | — | 20 | 2.48 s | — | — | — | — | 695 MB | +| `example/group_index` | examples/group_index.rs | — | 20 | 3.01 s | — | — | — | — | 990 MB | +| `example/circuit_fit` | examples/circuit_fit.rs | — | 20 | 0.43 s | — | — | — | — | 11 MB | ## Phases -- `fdfd2d/slab-lu` on 1 thread: assembly and LU 0.75 s; port modes 0.00 s; S-matrix, 2 runs 0.14 s -- `fdfd3d/guide-qmr` on 1 thread: assembly 1.27 s; QMR to 1e-6 10.90 s (2 077 iterations) -- `fdfd3d/guide-ilu` on 1 thread: assembly and ILU(0) 1.45 s; QMR to 1e-8 4.36 s (195 iterations) -- `fdfd3d/guide-multigrid` on 1 thread: assembly and multigrid 2.52 s; GMRES to 1e-8 0.83 s (16 iterations) -- `fdfd3d/diel-multigrid` on 1 thread: assembly and multigrid 11.39 s; GMRES to 1e-8 11.21 s (46 iterations) -- `fdfd3d/strip-ports-multigrid` on 1 thread: assembly and multigrid 24.34 s; port modes 0.65 s; S-matrix, 2 runs of GMRES to 1e-8 27.66 s (56 iterations) -- `fdfd2d/slab-lu` on 20 threads: assembly and LU 1.72 s; port modes 0.00 s; S-matrix, 2 runs 0.16 s -- `fdfd3d/guide-qmr` on 20 threads: assembly 1.30 s; QMR to 1e-6 2.82 s (2 077 iterations) -- `fdfd3d/guide-ilu` on 20 threads: assembly and ILU(0) 1.46 s; QMR to 1e-8 3.31 s (195 iterations) -- `fdfd3d/guide-multigrid` on 20 threads: assembly and multigrid 1.98 s; GMRES to 1e-8 0.51 s (16 iterations) -- `fdfd3d/diel-multigrid` on 20 threads: assembly and multigrid 8.11 s; GMRES to 1e-8 4.52 s (46 iterations) -- `fdfd3d/strip-ports-multigrid` on 20 threads: assembly and multigrid 17.13 s; port modes 0.80 s; S-matrix, 2 runs of GMRES to 1e-8 11.67 s (56 iterations) +- `fdfd2d/slab-lu` on 1 thread: assembly and LU 0.78 s; port modes 0.00 s; S-matrix, 2 runs 0.14 s +- `fdfd3d/guide-qmr` on 1 thread: assembly 1.29 s; QMR to 1e-6 12.10 s (2 077 iterations) +- `fdfd3d/guide-ilu` on 1 thread: assembly and ILU(0) 1.55 s; QMR to 1e-8 4.72 s (195 iterations) +- `fdfd3d/guide-multigrid` on 1 thread: assembly and multigrid 2.63 s; GMRES to 1e-8 0.89 s (16 iterations) +- `fdfd3d/diel-multigrid` on 1 thread: assembly and multigrid 11.65 s; GMRES to 1e-8 12.07 s (46 iterations) +- `fdfd3d/strip-ports-multigrid` on 1 thread: assembly and multigrid 24.99 s; port modes 0.66 s; S-matrix, 2 runs of GMRES to 1e-8 30.01 s (56 iterations) +- `fdfd2d/slab-lu` on 20 threads: assembly and LU 0.86 s; port modes 0.00 s; S-matrix, 2 runs 0.16 s +- `fdfd3d/guide-qmr` on 20 threads: assembly 1.30 s; QMR to 1e-6 2.77 s (2 077 iterations) +- `fdfd3d/guide-ilu` on 20 threads: assembly and ILU(0) 1.48 s; QMR to 1e-8 3.34 s (195 iterations) +- `fdfd3d/guide-multigrid` on 20 threads: assembly and multigrid 2.00 s; GMRES to 1e-8 0.52 s (16 iterations) +- `fdfd3d/diel-multigrid` on 20 threads: assembly and multigrid 8.12 s; GMRES to 1e-8 4.51 s (46 iterations) +- `fdfd3d/strip-ports-multigrid` on 20 threads: assembly and multigrid 16.93 s; port modes 0.93 s; S-matrix, 2 runs of GMRES to 1e-8 11.60 s (56 iterations) ## Problems diff --git a/src/bench.rs b/src/bench.rs index c917573..0e884ab 100644 --- a/src/bench.rs +++ b/src/bench.rs @@ -98,6 +98,23 @@ impl Measurement { .fold((0.0, 0), |(s, n), (t, i)| (s + t, n + i)); (iterations > 0).then(|| seconds * 1e9 / (self.unknowns as f64 * iterations as f64)) } + + /// The memory bandwidth the counted phases reached, in GB/s: their bytes over their time. + /// The kernels' traffic from memory or cache, each byte counted once (src/traffic.rs). + pub fn gigabytes_per_second(&self) -> Option { + let (bytes, seconds) = self + .phases + .iter() + .filter_map(|p| p.bytes.map(|b| (b, p.seconds))) + .fold((0.0, 0.0), |(b, s), (pb, ps)| (b + pb, s + ps)); + (seconds > 0.0).then(|| bytes / seconds / 1e9) + } +} + +/// The bytes the kernels have moved so far in this process ([`Phase::bytes`]): what a phase's +/// bytes are the difference of. +pub fn bytes_moved() -> f64 { + crate::traffic::moved() as f64 } /// One timed phase of a problem. @@ -109,6 +126,10 @@ pub struct Phase { pub seconds: f64, /// QMR's iterations in it, if it iterated. pub iterations: Option, + /// The bytes its kernels moved, by the model of `traffic` (src/traffic.rs: each nonzero's + /// value and index, each row pointer, each vector read or written once), if counted. + #[serde(default)] + pub bytes: Option, } /// The error of what a problem computed. @@ -243,6 +264,7 @@ fn timed(name: &str, phases: &mut Vec, f: impl FnOnce() -> Result) name: name.into(), seconds: t.elapsed().as_secs_f64(), iterations: None, + bytes: None, }); Ok(value) } @@ -439,12 +461,13 @@ fn guide_3d( let mut phases = Vec::new(); let field = { let solver = solve.build(&mut phases, new)?; - let t = Instant::now(); + let (t, moved) = (Instant::now(), bytes_moved()); let (field, convergence) = solver.solve(&source, stopping(tolerance))?; phases.push(Phase { name: format!("{} to {tolerance:e}", solve.label()), seconds: t.elapsed().as_secs_f64(), iterations: Some(convergence.iterations), + bytes: Some(bytes_moved() - moved), }); field }; @@ -541,7 +564,7 @@ fn strip_ports_3d( tolerance, max_iterations: 100_000, }; - let t = Instant::now(); + let (t, moved) = (Instant::now(), bytes_moved()); let (s, runs) = solver.s_matrix_with_convergence(&ports, stopping)?; phases.push(Phase { name: format!( @@ -551,6 +574,7 @@ fn strip_ports_3d( ), seconds: t.elapsed().as_secs_f64(), iterations: Some(runs.iter().map(|c| c.iterations).sum()), + bytes: Some(bytes_moved() - moved), }); timed_done(); let through = (c64::new(0.0, 1.0) * ports[0].mode.beta() * ((right - left) as f64 * h)).exp(); @@ -627,6 +651,7 @@ fn built_in_job(text: &str, timed_done: &mut dyn FnMut()) -> Result name: format!("the job, {points} points"), seconds, iterations: None, + bytes: None, }], accuracy: None, }) @@ -658,16 +683,19 @@ mod tests { name: "assembly".into(), seconds: 1.0, iterations: None, + bytes: None, }, Phase { name: "QMR".into(), seconds: 2.0, iterations: Some(100), + bytes: None, }, Phase { name: "QMR".into(), seconds: 2.0, iterations: Some(300), + bytes: None, }, ], accuracy: None, diff --git a/src/fdfd/krylov.rs b/src/fdfd/krylov.rs index 45b83bb..4b4b4b8 100644 --- a/src/fdfd/krylov.rs +++ b/src/fdfd/krylov.rs @@ -214,6 +214,7 @@ impl Sparse { /// is summed in the same order whatever the threads, so the result is the same bit for bit. fn product(starts: &[usize], columns: &[usize], values: &[c64], v: &[c64], out: &mut [c64]) { use rayon::prelude::*; + crate::traffic::add(crate::traffic::product(out.len(), v.len(), values.len())); let row = |r: usize| -> c64 { (starts[r]..starts[r + 1]) .map(|k| values[k] * v[columns[k]]) @@ -319,6 +320,7 @@ const CHUNK: usize = 16_384; /// chunk summed in order and the chunks' sums added in order. fn dot_conj(u: &[c64], v: &[c64]) -> c64 { use rayon::prelude::*; + crate::traffic::add(crate::traffic::vectors(u.len(), 2, 0)); let parts: Vec = u .par_chunks(CHUNK) .zip(v.par_chunks(CHUNK)) @@ -330,6 +332,7 @@ fn dot_conj(u: &[c64], v: &[c64]) -> c64 { /// w += Σ_j c_j v_j, value by value on rayon's threads, each value's terms in order. fn combine(w: &mut [c64], coefficients: &[c64], vectors: &[Vec]) { use rayon::prelude::*; + crate::traffic::add(crate::traffic::vectors(w.len(), 1 + vectors.len(), 1)); w.par_iter_mut() .with_min_len(CHUNK) .enumerate() @@ -552,6 +555,7 @@ impl Sums { /// Σ a_k b_k, unconjugated (Freund and Nachtigal's bilinear form). fn dot(&mut self, a: &[c64], b: &[c64]) -> c64 { use rayon::prelude::*; + crate::traffic::add(crate::traffic::vectors(a.len(), 2, 0)); a.par_chunks(CHUNK) .zip(b.par_chunks(CHUNK)) .map(|(p, q)| p.iter().zip(q).map(|(x, y)| x * y).sum::()) @@ -571,6 +575,8 @@ impl Sums { left: Option>, ) -> (f64, f64) { use rayon::prelude::*; + let (read, written) = if left.is_some() { (6, 2) } else { (3, 1) }; + crate::traffic::add(crate::traffic::vectors(v_next.len(), read, written)); let right = |k: usize| av[k] - alpha * v[k] - beta_v * v_old[k]; match left { Some(Left { @@ -636,6 +642,9 @@ impl Sums { w_next, r, } = out; + // v, p_{n−1}, p_{n−2} read; x, v_{n+1}, r read and written; p written; and w_{n+1} + let left = usize::from(w_next.is_some()); + crate::traffic::add(crate::traffic::vectors(p.len(), 6 + left, 4 + left)); // the value k's update, but for w; its share of ‖r‖², and v_{n+1}'s value let update = |k: usize, pk: &mut c64, xk: &mut c64, vk: &mut c64, rk: &mut c64| { *pk = (v[k] - epsilon * p_old[k] - theta * p_older[k]) / delta; @@ -1127,7 +1136,13 @@ impl Triangle { }; 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()]; + let diagonal = self.inverse_diagonal.is_some(); for _ in 1..sweeps { + // a product with the factor's off-diagonal part, and b read + crate::traffic::add( + crate::traffic::triangular(b.len(), self.values.len(), diagonal) + + crate::traffic::vectors(b.len(), 1, 0), + ); next.par_iter_mut() .with_min_len(4096) .enumerate() @@ -1146,6 +1161,8 @@ impl Triangle { /// Solves in place: `x` holds the right-hand side on entry. fn solve(&self, x: &mut [c64]) { let n = x.len(); + let diagonal = self.inverse_diagonal.is_some(); + crate::traffic::add(crate::traffic::triangular(n, self.values.len(), diagonal)); let rows: Box> = if self.forward { Box::new(0..n) } else { @@ -1551,6 +1568,22 @@ mod tests { assert!(Sparse::new(n, odd).symmetrized().is_none()); } + #[test] + fn a_product_and_a_triangular_solve_count_their_bytes() { + // other tests' kernels may add to the counter meanwhile, never take from it + let a = tridiagonal(100); + let before = crate::traffic::moved(); + let _ = a.apply(&vec![c64::new(1.0, 0.0); 100]); + assert!(crate::traffic::moved() - before >= crate::traffic::product(100, 100, 298) as u64); + let ilu = Ilu0::new(&a).unwrap(); + let before = crate::traffic::moved(); + let _ = ilu.solve(&vec![c64::new(1.0, 0.0); 100]); + // L (unit) and U (with its inverse diagonal), 99 off-diagonal entries each + let both = + crate::traffic::triangular(100, 99, false) + crate::traffic::triangular(100, 99, true); + assert!(crate::traffic::moved() - before >= both as u64); + } + #[test] fn qmr_is_the_same_bit_for_bit_on_any_number_of_threads() { // plain and preconditioned, on more values than one chunk of its sums, with its history diff --git a/src/fdfd/three/multigrid.rs b/src/fdfd/three/multigrid.rs index f54f00c..4d2aa0b 100644 --- a/src/fdfd/three/multigrid.rs +++ b/src/fdfd/three/multigrid.rs @@ -128,7 +128,9 @@ impl Transfer { fn apply(&self, v: &[c64]) -> Vec { use rayon::prelude::*; - let mut out = vec![c64::new(0.0, 0.0); self.starts.len() - 1]; + let rows = self.starts.len() - 1; + crate::traffic::add(crate::traffic::product(rows, v.len(), self.values.len())); + let mut out = vec![c64::new(0.0, 0.0); rows]; out.par_iter_mut() .with_min_len(4096) .enumerate() @@ -140,6 +142,7 @@ impl Transfer { /// x += d, on rayon's threads. fn add_to(x: &mut [c64], d: &[c64]) { use rayon::prelude::*; + crate::traffic::add(crate::traffic::vectors(x.len(), 2, 1)); x.par_iter_mut() .with_min_len(16_384) .zip(d) diff --git a/src/lib.rs b/src/lib.rs index 00e244f..0285105 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -48,6 +48,7 @@ pub mod raster; pub mod run; mod sparse; pub mod stack; +mod traffic; pub mod units; pub mod validation; diff --git a/src/traffic.rs b/src/traffic.rs new file mode 100644 index 0000000..aeac537 --- /dev/null +++ b/src/traffic.rs @@ -0,0 +1,69 @@ +//! The bytes the iterative solvers' kernels move: what [`crate::bench`] divides by a phase's time +//! for the memory bandwidth a solve reached, against the machine's measured roof (S. Williams, +//! A. Waterman, D. Patterson, Commun. ACM 52(4), 65 (2009), doi:10.1145/1498765.1498785). +//! +//! The model counts what a kernel must read and write at least once: each stored nonzero's value +//! (16 bytes, complex) and column index (8), each row pointer (8), each vector it reads and each +//! it writes (16 bytes a value; one read and written, twice). A gathered vector, x in A x, is +//! counted once, as if it stayed in cache, and nothing is counted for the factorizations' setup. +//! The count is the kernels' traffic from memory or from cache: a problem whose vectors fit in +//! the last-level cache (the 40³ guide's, 3 MB each in a 30 MB L3) moves them faster than memory +//! does, and its bandwidth can exceed the triad's. Each kernel adds its bytes once a call, from +//! all threads into one counter: it counts every solve in the process, which a benchmark runs +//! one at a time. + +use std::sync::atomic::{AtomicU64, Ordering}; + +static MOVED: AtomicU64 = AtomicU64::new(0); + +/// Adds a kernel's bytes. +pub(crate) fn add(bytes: usize) { + MOVED.fetch_add(bytes as u64, Ordering::Relaxed); +} + +/// The bytes moved so far in this process. +pub(crate) fn moved() -> u64 { + MOVED.load(Ordering::Relaxed) +} + +/// A complex value. +const VALUE: usize = 16; +/// A stored index. +const INDEX: usize = 8; + +/// A sparse product by rows, y = A x: A's values and columns, its row pointers, x read once and +/// y written. +pub(crate) fn product(rows: usize, columns: usize, nonzeros: usize) -> usize { + nonzeros * (VALUE + INDEX) + (rows + 1) * INDEX + columns * VALUE + rows * VALUE +} + +/// A triangular solve in place by rows, with or without a stored inverse diagonal: the factor's +/// off-diagonal values and columns, its row pointers, the diagonal, and x read and written. +pub(crate) fn triangular(rows: usize, nonzeros: usize, diagonal: bool) -> usize { + nonzeros * (VALUE + INDEX) + + (rows + 1) * INDEX + + if diagonal { rows * VALUE } else { 0 } + + 2 * rows * VALUE +} + +/// A pass over vectors of `n` values, reading `read` of them and writing `written`. +pub(crate) fn vectors(n: usize, read: usize, written: usize) -> usize { + n * (read + written) * VALUE +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn a_small_products_bytes_by_hand() { + // a 3 x 3 matrix with 7 nonzeros: 7 values (112 B) and 7 columns (56 B), 4 row pointers + // (32 B), x read (48 B) and y written (48 B) + assert_eq!(product(3, 3, 7), 112 + 56 + 32 + 48 + 48); + // a unit lower triangle with 3 off-diagonal entries on 4 rows: 48 + 24 + 40 + 128 + assert_eq!(triangular(4, 3, false), 48 + 24 + 40 + 128); + // with an inverse diagonal: 64 more + assert_eq!(triangular(4, 3, true), 48 + 24 + 40 + 64 + 128); + assert_eq!(vectors(10, 3, 1), 640); + } +} diff --git a/studio/src-tauri/src/bench.rs b/studio/src-tauri/src/bench.rs index 33b8d32..0b0dfe4 100644 --- a/studio/src-tauri/src/bench.rs +++ b/studio/src-tauri/src/bench.rs @@ -277,6 +277,7 @@ fn single(name: &str, seconds: f64, grid: String) -> Measurement { name: name.into(), seconds, iterations: None, + bytes: None, }], accuracy: None, } @@ -291,8 +292,13 @@ fn markdown(chosen: &[Entry], bandwidth: &[(usize, f64)], runs: &[Timed]) -> Str `photonoxide::bench` and the slowest examples, each in a process of its own. Times depend \ on the machine and on what else ran on it; the errors don't. Peak memory is the process's \ peak resident set (its peak working set on Windows); ns per unknown per iteration is the \ - QMR phases' time over the unknowns and the iterations. See the performance plan, \ - docs/plans/performance.md.\n" + QMR phases' time over the unknowns and the iterations. Bandwidth is the iterating phases' \ + bytes over their time, the bytes counted by the kernels themselves (each nonzero's value \ + and index, each row pointer, each vector read or written once: src/traffic.rs), from \ + memory or from cache; beside it, its share of the triad on the same threads, the roof \ + (Williams et al., Commun. ACM 52(4), 65 (2009)). A problem whose vectors fit in the \ + last-level cache can count above it: the 40³ guide's do. See the performance \ + plan, docs/plans/performance.md.\n" ); let _ = writeln!(s, "- photonoxide {}, {}", photonoxide::VERSION, machine()); if !bandwidth.is_empty() { @@ -308,17 +314,27 @@ fn markdown(chosen: &[Entry], bandwidth: &[(usize, f64)], runs: &[Timed]) -> Str } let _ = writeln!( s, - "\n| Problem | Grid | Unknowns | Threads | Time | Iterations | ns/unknown/iteration | Error | Peak memory |" + "\n| Problem | Grid | Unknowns | Threads | Time | Iterations | ns/unknown/iteration | Bandwidth (of the triad) | Error | Peak memory |" ); - let _ = writeln!(s, "|---|---|---:|---:|---:|---:|---:|---|---:|"); + let _ = writeln!(s, "|---|---|---:|---:|---:|---:|---:|---:|---|---:|"); for r in runs { match &r.outcome { Ok(report) => { let m = &report.measurement; let dash = || "—".to_string(); + let roof = bandwidth + .iter() + .find(|(t, _)| *t == r.threads) + .map(|&(_, gb_s)| gb_s); + let reached = m + .gigabytes_per_second() + .map_or_else(dash, |gb_s| match roof { + Some(roof) => format!("{gb_s:.1} GB/s ({:.0}%)", 100.0 * gb_s / roof), + None => format!("{gb_s:.1} GB/s"), + }); let _ = writeln!( s, - "| `{}` | {} | {} | {} | {:.2} s | {} | {} | {} | {} |", + "| `{}` | {} | {} | {} | {:.2} s | {} | {} | {reached} | {} | {} |", r.id, m.grid, if m.unknowns > 0 { @@ -340,7 +356,7 @@ fn markdown(chosen: &[Entry], bandwidth: &[(usize, f64)], runs: &[Timed]) -> Str Err(e) => { let _ = writeln!( s, - "| `{}` | | | {} | failed: {} | | | | |", + "| `{}` | | | {} | failed: {} | | | | | |", r.id, r.threads, e.replace('|', "/")