diff --git a/crates/foundation/apsis-math/src/lib.rs b/crates/foundation/apsis-math/src/lib.rs index 86782ba..1422afc 100644 --- a/crates/foundation/apsis-math/src/lib.rs +++ b/crates/foundation/apsis-math/src/lib.rs @@ -29,6 +29,7 @@ #![deny(missing_docs)] mod mat3; +pub mod parallel; mod quat; mod vec3; diff --git a/crates/foundation/apsis-math/src/parallel.rs b/crates/foundation/apsis-math/src/parallel.rs new file mode 100644 index 0000000..1b5aef7 --- /dev/null +++ b/crates/foundation/apsis-math/src/parallel.rs @@ -0,0 +1,150 @@ +//! Scoped parallelism built on the standard library alone. +//! +//! The workspace forbids third party crates, so a threading pool would have to +//! be written here. It is not: the operations in this workspace that benefit +//! from several cores are all embarrassingly parallel loops over independent +//! items, and the standard library's scoped threads express those directly. +//! `std::thread::scope` borrows its data, so a caller can share a slice or a +//! reader without an [`Arc`], and joins every thread before returning, so no +//! result escapes the scope still being written. +//! +//! The entry point is [`for_each_index`], which splits a range into one chunk +//! per thread and runs a closure over each element. The closure must be +//! callable from another thread, so anything it captures has to be `Sync`; +//! in practice that means reading shared input and writing to a disjoint +//! slot of the output, which is the shape every parallel loop here has. +//! +//! Because there is no pool, the cost of a parallel call includes starting +//! the threads. That is worth paying only when each thread has real work, so +//! [`thread_count`] reports what the machine offers and the callers in this +//! workspace use it to skip the parallel path for small inputs. +//! +//! ``` +//! use apsis_math::parallel::map_index; +//! +//! // Square every element, one value per index. +//! let input: Vec = (0..1000).map(f64::from).collect(); +//! let output = map_index(input.len(), |index| input[index] * input[index]); +//! assert_eq!(output[10], 100.0); +//! ``` + +use core::hint; + +/// The number of threads to use for a piece of work. +/// +/// This is the machine's parallelism as reported by the standard library, at +/// least one. It is the number the callers use to decide whether a parallel +/// run is worth starting. +pub fn thread_count() -> usize { + std::thread::available_parallelism() + .map(|count| count.get()) + .unwrap_or(1) +} + +/// The number of threads that should run a loop over `len` items. +/// +/// A loop shorter than the thread count would leave most threads with nothing +/// to do while still paying to start them, so this returns one thread until +/// there is enough work for a second. The threshold is a per item cost of +/// roughly a microsecond, below which starting a thread cannot pay for +/// itself. +pub fn threads_for(len: usize) -> usize { + const ITEMS_PER_THREAD: usize = 1_024; + thread_count().min(len.div_ceil(ITEMS_PER_THREAD)).max(1) +} + +/// Runs a closure once for each index of a range, spread over several threads. +/// +/// The range is split into one contiguous chunk per thread and each thread +/// runs the closure for the indices in its chunk. The closure may be called +/// from any thread, so it must be `Send` and `Sync`; anything it captures is +/// shared rather than copied. The call returns only after every thread has +/// finished, so the data the closure wrote is safe to read afterwards. +/// +/// The closure is handed the index rather than the element, which keeps the +/// signature independent of the collection: the caller reads the input and +/// writes the output itself, and the two never overlap. +/// +/// # Panics +/// +/// If the closure panics, the panic is resumed on the calling thread after +/// the other threads have been joined, so a failure surfaces at the call +/// rather than in a detached thread. +pub fn for_each_index(len: usize, operation: F) +where + F: Fn(usize) + Sync, +{ + let threads = threads_for(len); + if threads <= 1 { + // A single thread would spawn and join one thread to do exactly what + // this loop does, so it runs directly. + for index in 0..len { + operation(index); + } + return; + } + let chunk = len.div_ceil(threads); + std::thread::scope(|scope| { + for worker in 0..threads { + let first = worker * chunk; + if first >= len { + break; + } + let last = hint::black_box((first + chunk).min(len)); + let operation = &operation; + scope.spawn(move || { + for index in first..last { + operation(index); + } + }); + } + }); +} + +/// Collects one value for each index of a range, spread over several threads. +/// +/// This is [`for_each_index`] for the common case where the result of the +/// work is a value rather than a write into a shared array. The output is +/// returned in index order regardless of which thread produced which element, +/// so it does not depend on the scheduling. +/// +/// The closure is handed the index and returns the value for that index. As +/// with [`for_each_index`], it may run on any thread, so captured data must be +/// `Sync`. +pub fn map_index(len: usize, operation: F) -> Vec +where + T: Send, + F: Fn(usize) -> T + Sync, +{ + let threads = threads_for(len); + if threads <= 1 { + return (0..len).map(operation).collect(); + } + let chunk = len.div_ceil(threads); + // Each thread writes a contiguous run of the output, so the buffer is + // split into disjoint mutable slices up front and no slot is written + // twice. `Vec>` is filled in parallel and unwrapped after the + // scope, which is where the values are moved into their final order. + let mut slots: Vec> = (0..len).map(|_| None).collect(); + let operation = &operation; + std::thread::scope(|scope| { + let mut rest = &mut slots[..]; + let mut offset = 0; + while offset < len { + let take = chunk.min(len - offset); + let (this, tail) = core::mem::take(&mut rest).split_at_mut(take); + let base = offset; + scope.spawn(move || { + for (index, slot) in this.iter_mut().enumerate() { + *slot = Some(operation(base + index)); + } + }); + rest = tail; + offset += take; + } + }); + slots + .into_iter() + .map(|slot| slot.expect("every index was assigned")) + .collect() +} diff --git a/crates/foundation/apsis-math/tests/parallel.rs b/crates/foundation/apsis-math/tests/parallel.rs new file mode 100644 index 0000000..929dd45 --- /dev/null +++ b/crates/foundation/apsis-math/tests/parallel.rs @@ -0,0 +1,136 @@ +//! Tests for the scoped parallelism helpers. +//! +//! The central property is that a parallel run and a serial run produce the +//! same values, because the whole point of the helpers is to be a drop-in for +//! the sequential loop they replace. Each test therefore compares against the +//! sequential answer rather than against a constant, and forces both the +//! single and multi thread paths by choosing lengths that fall on each side +//! of the threshold. + +use std::sync::atomic::Ordering; + +use apsis_math::parallel::{for_each_index, map_index, thread_count, threads_for}; + +#[test] +fn the_thread_count_is_at_least_one() { + assert!(thread_count() >= 1); +} + +#[test] +fn small_work_uses_a_single_thread() { + // Below the per thread threshold there is nothing to gain from starting + // a thread, so the helper says so. + assert_eq!(threads_for(0), 1); + assert_eq!(threads_for(1), 1); + assert_eq!(threads_for(64), 1); +} + +#[test] +fn large_work_uses_more_than_one_thread() { + let count = thread_count(); + if count > 1 { + assert!(threads_for(1_000_000) > 1); + // Never more threads than there are items. + assert!(threads_for(10) <= 10); + } +} + +#[test] +fn every_index_is_visited_exactly_once() { + use std::sync::atomic::AtomicUsize; + for &len in &[0_usize, 1, 7, 1_000, 5_000, 100_000] { + let visits = (0..len).map(|_| AtomicUsize::new(0)).collect::>(); + for_each_index(len, |index| { + visits[index].fetch_add(1, Ordering::Relaxed); + }); + assert!( + visits + .iter() + .all(|count| count.load(Ordering::Relaxed) == 1), + "length {len} visited an index more than once" + ); + } +} + +#[test] +fn the_parallel_result_matches_the_sequential_one() { + for &len in &[0_usize, 1, 1_000, 50_000, 300_000] { + let sequential: Vec = (0..len) + .map(|index| (index as f64).sin() * 1_000.0) + .collect(); + let parallel = map_index(len, |index| (index as f64).sin() * 1_000.0); + assert_eq!(parallel, sequential, "length {len} differs"); + } +} + +#[test] +fn the_output_is_in_index_order() { + // Threads finish out of order, so this checks that the helper restores + // index order rather than completion order. + let len = 200_000; + let values = map_index(len, |index| index); + for (position, value) in values.iter().enumerate() { + assert_eq!(*value, position); + } +} + +#[test] +fn a_closure_may_capture_shared_input() { + // The shape every caller in the workspace uses: read a shared slice and + // produce one value per element. `map_index` owns the output buffer, so + // no shared mutable state is needed. + let input: Vec = (0..100_000).map(f64::from).collect(); + let output = map_index(input.len(), |index| input[index] * input[index]); + for index in [0, 1, 999, 50_000, 99_999] { + assert_eq!(output[index], input[index] * input[index]); + } +} + +#[test] +fn a_closure_may_accumulate_atomically() { + // Writing to a shared counter needs atomics; this checks the helper + // tolerates that, since summing a large range is a common pattern. + use std::sync::atomic::AtomicU64; + let total = AtomicU64::new(0); + let len = 200_000_u64; + for_each_index(len as usize, |index| { + total.fetch_add(index as u64, Ordering::Relaxed); + }); + let expected: u64 = (0..len).sum(); + assert_eq!(total.load(Ordering::Relaxed), expected); +} + +#[test] +fn an_empty_range_does_nothing() { + let empty: Vec = map_index(0, |_| 1.0); + assert!(empty.is_empty()); + // The index form accepts an empty range without spawning anything. + for_each_index(0, |_| unreachable!("an empty range runs no closure")); +} + +#[test] +fn a_range_that_does_not_divide_evenly_covers_everything() { + // A length that is not a multiple of the thread count is where an + // off by one in the chunking would drop the tail. + use std::sync::atomic::AtomicUsize; + for &len in &[1_001_usize, 4_097, 65_537, 100_003] { + let marks = (0..len).map(|_| AtomicUsize::new(0)).collect::>(); + for_each_index(len, |index| { + marks[index].store(1, Ordering::Relaxed); + }); + assert!( + marks.iter().all(|mark| mark.load(Ordering::Relaxed) == 1), + "length {len} missed an index" + ); + } +} + +#[test] +fn values_of_a_non_copy_type_survive_the_collect() { + // A string is not `Copy`, so this exercises the unwrap-and-move path that + // the parallel branch uses. + let values = map_index(10_000, |index| format!("item-{index}")); + assert_eq!(values[0], "item-0"); + assert_eq!(values[9_999], "item-9999"); + assert_eq!(values.len(), 10_000); +}