libgis: Add reentrant drand48 for independent parallel streams - #7884
Draft
petrasovaa wants to merge 16 commits into
Draft
petrasovaa wants to merge 16 commits into
petrasovaa wants to merge 16 commits into
Conversation
As random number generator tracks its internal state as a static variable, it is necessary to guard its updates with mutex to prevent misbehavour when multiple threads change values simultaneously
multiple threads As random number generation is guarded by a mutex, it makes generation a serial operation – always leading to the same sequence of numbers
Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
Replace the pthread mutex approach with a lock-free compare-and-swap on the whole 48-bit generator state kept in a single C11 atomic integer. One successful compare-and-swap is exactly one step of the standard drand48 linear congruential generator, so the sequence of generated values for a given seed is identical to a single-threaded run and to the previous limb-based implementation. Compared to the mutex: the feature is gated by the compiler (__STDC_NO_ATOMICS__), not by --with-pthread, so it is active in every supported build, including the OpenMP-without-pthread configurations (macOS CI, Docker images) where the race reported in OSGeo#5742 actually occurs; no build system changes or extra link dependencies are needed (the Makefile and CMake changes are dropped, which also avoids defining HAVE_PTHREAD from the Win32-threads CMake package on MSVC); and there is no lock contention (in a microbenchmark the mutex made threaded runs slower than single-threaded at every thread count, while the atomic version is about twice as fast as the mutex when several threads generate numbers at the same time). A single shared generator still cannot scale with the number of threads; per-thread state remains future work. Where C11 atomics are unavailable, the code falls back to the previous plain implementation, which is not thread-safe (unchanged behavior). Output compatibility is guarded by the PRNG regression tests merged in OSGeo#7633 and OSGeo#7634 (lib/gis/tests/lib_gis_lrand48_test.py and raster/r.mapcalc/tests/r_mapcalc_rand_seed_test.py) and was additionally verified against the previous implementation for multiple seeds. The 64-bit reformulation of the generator and its verification were AI-assisted (Claude).
…nd_thread-atomics # Conflicts: # lib/gis/CMakeLists.txt # lib/gis/Makefile
G_drand48() advances a single shared state, so parallel callers either race on it or serialize on it. Add G_srand48_r() and G_drand48_r(), which advance a generator owned by the caller. Streams are spaced apart by LCG jump-ahead so they are disjoint by construction, and stream 0 starts where G_srand48() does. The LCG step is now shared with the global generator, replacing an equivalent 16-bit limb implementation. r.sim.water and r.sim.sediment give each walker chunk its own stream. Repeated runs with the same random_seed and nprocs are now identical where they previously were not, and single-threaded results are unchanged. Four threads go from 1.25x to 3.1x on an 8-thread machine. Based on OSGeo#6480. Written with substantial help from Claude Code.
github-actions
Bot
added
raster
nilason
added
the
conflicts/needs rebase
# Conflicts: # lib/gis/lrand48.c # lib/gis/testsuite/test_lrand48.py
Contributor
|
I took the liberty to update to main and resolve the conflicts. I hope I had it right in the end. |
nilason
removed
the
conflicts/needs rebase
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
5 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
TL;DR
G_drand48()keeps one shared state, so parallel code either races on it orserializes on it. This adds a caller-owned generator with independent streams
and uses it in r.sim.water / r.sim.sediment.
random_seedandnprocsare now identical;before they differed by up to ±0.024 m against a ~0.6 m max depth
with no fixture changes
Builds on #6480 and needs it to land first.
Why #6480 is not enough on its own
#6480 makes the shared generator thread-safe with an atomic compare-and-swap,
which is the right fix for callers that want one stream. But it is still one
stream: under contention only one CAS succeeds per cache-line round trip, so
r.sim.water (which draws 2-3 times per walker per timestep) gets slower with
every thread added — 31.4 s at 1 thread, 74.1 s at 8. It also cannot make
results reproducible, as the PR itself documents.
How the streams work
One generator step is the affine map
x -> a*x + c, and composing two such mapsgives another, so the map for n steps can be built by repeated squaring. Stream
s starts
s * 2^36steps along the cycle: 4096 streams of ~6.9e10 draws each,disjoint by construction rather than merely unlikely to collide. Stream 0 starts
exactly where
G_srand48()does, which is what keeps single-threaded resultsidentical.
In r.sim, each walker chunk owns a stream seeded by chunk index rather than
thread number, and copies it into a thread-local so the 8-byte states do not
share a cache line. Walker setup draws from chunk 0's stream, which chunk 0 then
continues — with one chunk that is a single stream used in the order the tool
has always used it.
Measurements
nc_spm
elevationat 10 m (1350x1500), 800k walkers,duration=10,random_seed=42, 145 iterations, on a 4-core / 8-thread i5-10210U.Single-thread also improves because the 64-bit LCG is cheaper than the limb
arithmetic it replaces, which benefits every
G_drand48()caller.Tests: 39 pass in
lib/gis/tests(8 new), and both r.sim gunittest suites passunchanged. The existing sequence-pinning tests are what verify that removing the
limb implementation changed no output.
Not addressed here
nprocsstill changes results above one thread, because the number of walkerchunks is still the number of threads. Making the chunk count a parameter of its
own is a follow-up. Two other r.sim issues are also left alone:
grids->difis araster-sized array used as a per-walker scalar (a race and wasted memory), and
the hand-rolled thread chunking silently drops up to
nthreads/2walkers whennwalk/nthreadsrounds down.🤖 Generated with Claude Code