From b74d8f8545715be212e394231b8972a619446844 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sun, 23 Aug 2026 13:05:31 +0200 Subject: [PATCH] Give ParallelFor a home and a test The two shapes - a fixed contiguous split, and work stealing off an atomic - had been copied into whichever file wanted them: an anonymous namespace in RotationScaleMerge.cpp, and another, byte for byte the same, in ShadowFinder.cpp, whose comment claimed it was "the only other file that wants it". The header is taken from the beam-stop GPU commit on the performance branch, which is not otherwise being picked. ShadowFinder now uses it, so the construct is exercised rather than shipped unused, and its own copy is gone. The beam-stop mask is unchanged, which its reference count already pins. Tested for the properties the callers rely on and which are easy to lose in a rewrite: the chunked slices tile the range in order with no empty one, work stealing visits every item exactly once, one thread means the caller's loop in order, an empty or negative count does nothing, an exception in a worker reaches the caller, and - the point of the whole thing - the answer is the serial answer bit for bit at every thread count. RotationScaleMerge.cpp keeps its own copy for now; consolidating it belongs with the scaling work, which is not being touched here. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01VfYvJT5Nb71suJCowRBn5z --- common/ParallelFor.h | 55 ++++++++++ image_analysis/beam_stop/ShadowFinder.cpp | 40 +------ tests/CMakeLists.txt | 1 + tests/ParallelForTest.cpp | 127 ++++++++++++++++++++++ 4 files changed, 184 insertions(+), 39 deletions(-) create mode 100644 common/ParallelFor.h create mode 100644 tests/ParallelForTest.cpp diff --git a/common/ParallelFor.h b/common/ParallelFor.h new file mode 100644 index 00000000..767a070f --- /dev/null +++ b/common/ParallelFor.h @@ -0,0 +1,55 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include +#include + +// Two shapes of "run this over a range on several threads", used by the analysis code. Both take the +// worker count from the caller rather than asking the hardware, so a run that was told how many +// threads to use keeps to it. + +// Chunked: each worker gets one contiguous [lo, hi) range and there is no per-item synchronisation. +// Right for millions of cheap uniform items - the CPU stand-in for a flat CUDA grid-stride kernel. +// The split is fixed and deterministic, so a pass whose per-element work is independent gives the +// same answer as the serial loop, bit for bit. +template +void ParallelChunks(int n, size_t nthreads, Fn fn) { + if (n <= 0) return; + const int nt = static_cast(std::max(1, std::min(nthreads, static_cast(n)))); + if (nt == 1) { fn(0, n); return; } + const int chunk = (n + nt - 1) / nt; + std::vector> futures; + futures.reserve(nt); + for (int t = 0; t < nt; t++) { + const int lo = t * chunk, hi = std::min(n, lo + chunk); + if (lo >= hi) break; + futures.emplace_back(std::async(std::launch::async, [&fn, lo, hi] { fn(lo, hi); })); + } + for (auto &f : futures) f.get(); +} + +// Work-stealing per item, off a shared atomic counter: one atomic per item, so use it only where the +// per-item work is heavy and uneven (per-frame fits, per-ring selections) and the atomic amortises. +// For millions of tiny uniform items a per-item atomic is pure contention - use ParallelChunks. +template +void ParallelFor(int n, size_t nthreads, Fn fn) { + if (n <= 0) return; + if (nthreads <= 1 || n == 1) { + for (int i = 0; i < n; i++) fn(i); + return; + } + const size_t local = std::min(nthreads, static_cast(n)); + std::atomic next = 0; + std::vector> futures; + futures.reserve(local); + for (size_t t = 0; t < local; t++) + futures.emplace_back(std::async(std::launch::async, [&] { + for (int i = next.fetch_add(1); i < n; i = next.fetch_add(1)) + fn(i); + })); + for (auto &f : futures) f.get(); +} diff --git a/image_analysis/beam_stop/ShadowFinder.cpp b/image_analysis/beam_stop/ShadowFinder.cpp index 233a8d74..29eb36bc 100644 --- a/image_analysis/beam_stop/ShadowFinder.cpp +++ b/image_analysis/beam_stop/ShadowFinder.cpp @@ -2,6 +2,7 @@ // SPDX-License-Identifier: GPL-3.0-only #include "ShadowFinder.h" +#include "../../common/ParallelFor.h" #include #include @@ -47,45 +48,6 @@ constexpr int MIN_RING_PIXELS = 32; // at GetMask() time; the BFS forms keep them O(pixels) rather than O(pixels * radius). namespace { -// Run fn(lo, hi) over contiguous slices of [0, n), one per worker. For the per-pixel passes: the -// same self-load-balancing shape the rest of the codebase uses, kept local because this is the only -// other file that wants it. -template -void ParallelChunks(int n, size_t nthreads, Fn fn) { - if (n <= 0) return; - const int nt = static_cast(std::max(1, std::min(nthreads, static_cast(n)))); - if (nt == 1) { fn(0, n); return; } - const int chunk = (n + nt - 1) / nt; - std::vector> futures; - futures.reserve(nt); - for (int t = 0; t < nt; t++) { - const int lo = t * chunk, hi = std::min(n, lo + chunk); - if (lo >= hi) break; - futures.emplace_back(std::async(std::launch::async, [&fn, lo, hi] { fn(lo, hi); })); - } - for (auto &f : futures) f.get(); -} - -// Work-stealing per-item parallel, for items whose cost is very uneven - the rings, whose sizes go -// as the circumference and then fall away at the detector corners. -template -void ParallelFor(int n, size_t nthreads, Fn fn) { - if (n <= 0) return; - if (nthreads <= 1 || n == 1) { - for (int i = 0; i < n; i++) fn(i); - return; - } - const size_t local = std::min(nthreads, static_cast(n)); - std::atomic next = 0; - std::vector> futures; - futures.reserve(local); - for (size_t t = 0; t < local; t++) - futures.emplace_back(std::async(std::launch::async, [&] { - for (int i = next.fetch_add(1); i < n; i = next.fetch_add(1)) fn(i); - })); - for (auto &f : futures) f.get(); -} - // 8-connected dilation by `r` pixels, i.e. every pixel within Chebyshev distance r of a set one. // // This was a multi-source BFS, which is what the distance is defined by - but on a full rectangle diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 790296a7..9892c0b5 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -51,6 +51,7 @@ ADD_EXECUTABLE(jfjoch_test ZMQMetadataSocketTest.cpp JFJochReaderTest.cpp ShadowFinderTest.cpp + ParallelForTest.cpp RugnuxTest.cpp ResultReportTest.cpp RugnuxLargeTest.cpp diff --git a/tests/ParallelForTest.cpp b/tests/ParallelForTest.cpp new file mode 100644 index 00000000..cec2988a --- /dev/null +++ b/tests/ParallelForTest.cpp @@ -0,0 +1,127 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include +#include +#include +#include +#include + +#include "../common/ParallelFor.h" + +namespace { + const std::vector THREAD_COUNTS = {0, 1, 2, 3, 8, 64}; + const std::vector SIZES = {0, 1, 2, 7, 64, 1000}; +} + +// The split is fixed and contiguous, which is what lets a pass whose per-element work is independent +// give the serial answer bit for bit. Every element has to land in exactly one slice, the slices have +// to tile [0, n) in order, and there must never be an empty one. +TEST_CASE("ParallelChunks_SlicesTileTheRange", "[ParallelFor]") { + for (int n: SIZES) { + for (size_t nthreads: THREAD_COUNTS) { + CAPTURE(n, nthreads); + std::mutex m; + std::vector> slices; + std::vector visits(static_cast(std::max(n, 1)), 0); + + ParallelChunks(n, nthreads, [&](int lo, int hi) { + { + std::scoped_lock lock(m); + slices.emplace_back(lo, hi); + } + for (int i = lo; i < hi; i++) + visits[static_cast(i)]++; // slices are disjoint, so no lock needed + }); + + std::sort(slices.begin(), slices.end()); + int expected_lo = 0; + for (const auto &[lo, hi]: slices) { + CHECK(lo == expected_lo); + CHECK(hi > lo); // never an empty slice + expected_lo = hi; + } + CHECK(expected_lo == n); // and they reach the end + CHECK(slices.size() <= static_cast(std::max(n, 0))); + for (int i = 0; i < n; i++) + CHECK(visits[static_cast(i)] == 1); + } + } +} + +// Work stealing, so a worker takes whatever is next rather than a fixed slice - but every item still +// has to be done exactly once, whatever order they come out in. +TEST_CASE("ParallelFor_VisitsEveryItemOnce", "[ParallelFor]") { + for (int n: SIZES) { + for (size_t nthreads: THREAD_COUNTS) { + CAPTURE(n, nthreads); + std::vector> visits(static_cast(std::max(n, 1))); + for (auto &v: visits) + v.store(0); + + ParallelFor(n, nthreads, [&](int i) { visits[static_cast(i)].fetch_add(1); }); + + for (int i = 0; i < n; i++) + CHECK(visits[static_cast(i)].load() == 1); + } + } +} + +// One thread means the caller's loop, in order - the property the serial fallbacks rely on. +TEST_CASE("ParallelFor_IsSerialAndInOrderForOneThread", "[ParallelFor]") { + for (size_t nthreads: {size_t{0}, size_t{1}}) { + CAPTURE(nthreads); + std::vector order; + ParallelFor(16, nthreads, [&](int i) { order.push_back(i); }); + std::vector expected(16); + std::iota(expected.begin(), expected.end(), 0); + CHECK(order == expected); + } +} + +// Splitting the work must not change the answer. Each element is written by exactly one worker, so +// the result has to match the serial loop element for element, at every thread count. +TEST_CASE("ParallelFor_SplitDoesNotChangeTheResult", "[ParallelFor]") { + constexpr int N = 1000; + std::vector serial(N); + for (int i = 0; i < N; i++) + serial[static_cast(i)] = std::sin(i * 0.001) * 1e6 + i; + + for (size_t nthreads: THREAD_COUNTS) { + CAPTURE(nthreads); + std::vector chunked(N, 0.0), stolen(N, 0.0); + ParallelChunks(N, nthreads, [&](int lo, int hi) { + for (int i = lo; i < hi; i++) + chunked[static_cast(i)] = std::sin(i * 0.001) * 1e6 + i; + }); + ParallelFor(N, nthreads, [&](int i) { + stolen[static_cast(i)] = std::sin(i * 0.001) * 1e6 + i; + }); + CHECK(chunked == serial); // bit for bit, not approximately + CHECK(stolen == serial); + } +} + +// A negative or zero count is a no-op rather than an error - callers pass a computed size. +TEST_CASE("ParallelFor_DoesNothingForAnEmptyRange", "[ParallelFor]") { + int calls = 0; + for (int n: {0, -1, -1000}) { + ParallelChunks(n, 8, [&](int, int) { calls++; }); + ParallelFor(n, 8, [&](int) { calls++; }); + } + CHECK(calls == 0); +} + +// An exception thrown in a worker reaches the caller rather than terminating: the futures are waited +// on, so the other workers finish first and only then does it propagate. +TEST_CASE("ParallelFor_PropagatesAnException", "[ParallelFor]") { + CHECK_THROWS_AS(ParallelChunks(64, 4, [](int lo, int) { + if (lo == 0) throw std::runtime_error("from a chunk"); + }), std::runtime_error); + + CHECK_THROWS_AS(ParallelFor(64, 4, [](int i) { + if (i == 0) throw std::runtime_error("from an item"); + }), std::runtime_error); +}