Repository navigation
r.mapcalc: Enable parallel rand() with reproducible results - #8114
Merged
Merged
Conversation
rand() drew from the generator shared by the whole program, which hands out values to whichever thread asks first, so r.mapcalc evaluated expressions with rand() on one thread only, whatever nprocs was. Every rand() call draws one value per cell of a row, and every row evaluates every call once, so a row draws columns times calls values. r.mapcalc now builds an exact random number layout from the seed with a row as its unit (in r3.mapcalc, the rows of all depths one after another), and each thread places a generator state of its own at the row's stream before evaluating the row. The streams of an exact layout are the parts of the sequence which the rows drew when they were evaluated one after another from the shared generator, so the result for a seed stays the same as in earlier versions and is the same for any nprocs, and rand() no longer turns off parallel evaluation. A warning is given when the whole computation draws more than 2^46 values, after which the values repeat plus a constant. The seed is parsed as a long long, since a long does not hold the range of seeds on Windows, and checked right after parsing. A seed outside -2147483648 to 4294967295 is now an error rather than being silently reduced to its low 32 bits, and a seed with trailing characters is an error rather than being cut off. The automatic seed comes from G_random_generate_seed(), so it can be set with GRASS_RANDOM_SEED or SOURCE_DATE_EPOCH, and is still written to the history. rand() stays in the calc library. The library gains calc_set_random_state(), through which a tool gives rand() a function returning the calling thread's generator state, and r.mapcalc sets it to the state it places for the row. When no function is set, rand() draws from the generator shared by the whole program as before, so other users of the library get the same values as before. Written with the help of Claude Code.
Add the SPDX lines to the new test file and describe a row's random numbers as its stream, the term the library documentation uses. Written with the help of Claude Code.
Call G_random_state_from_seed() and G_random_layout_batches(), the new names of G_random_seed() and G_random_layout_runs() in libgis. The values of rand() do not change. Written with the help of Claude Code.
The libgis functions take and return the seed as int64_t, so the tool keeps it in that type, printed with PRId64, and forms the counts it passes to the layout in it. Written with the help of Claude Code.
The type comes with grass/gis.h, but code which uses it includes the header itself rather than relying on that. Written with the help of Claude Code.
The seed range, the environment variables behind a generated seed and the comparison with earlier versions go: a bad seed is refused with a message, the variables are documented on the environment variables page, and the comparison belongs to the release notes. What stays is that the result does not depend on nprocs and depends on all rand() calls of the expression. Written with the help of Claude Code.
Written with the help of Claude Code.
The function parses and checks the seed with the rules of the library, in place of the tool's own strtoll() and the throwaway state it built to have the range checked before the output maps are created.
The warning joined two complete clauses with a semicolon and had no final period. The style guide asks for periods after complete sentences, and most messages with two sentences use them. Written with the help of Claude Code.
The per-thread states were adjacent elements of one array, 8 bytes each, so up to eight threads wrote to the same cache line with every draw. Each row now places a state local to its evaluation, and the array only holds a pointer to it, written once per row. With 48 rand() calls per cell on a region of 2000 by 8000 cells, the best of three runs took: nprocs before after 1 5.16 s 5.17 s 2 3.03 s 2.80 s 4 1.96 s 1.67 s 8 2.47 s 1.40 s 16 2.04 s 1.56 s The results are the same as before for every nprocs.
calc_set_random_state() and rand() are public, so they get Doxygen comments in the source file, and the header keeps only the prototype. The static parts get short comments, and the comment no longer describes what rand() did in earlier versions.
The comment on the layout now says what the layout is instead of what rand() drew in earlier versions. The comments on the order of evaluation are dropped. The arguments of a function and the expressions of a row are evaluated one after another on the thread of the row, since G_init_workers() is never called, so the worker threads of G_begin_execute() do not exist. The tests of the order of the rand() calls fail if that changes. The reason the size check is only a warning is in the PR description instead.
The comment now says why following a variable would count the rand() calls of its value again for every use, while they draw once per cell. The function counting the calls of all expressions gets a short comment too.
Saying that the values depend on the number of rand() calls is enough, and it makes the explanation after it unnecessary.
A test that a feature no longer turns off parallel evaluation would not be written for new code, and the message it checked is not about rand(). The test file docstring no longer describes earlier versions; that history is in the PR description.
One test of three expressions of the three types replaces the tests of each type, of the order of expressions and of expressions from a file. The tests that the result does not depend on nprocs or that a seed gives the same result are dropped, since every test compares exact values for nprocs 1 and 4. The tests of different seeds and of the value range are dropped, since this change does not affect them, and the two tests of invalid seeds are merged.
r3.mapcalc always evaluates on one thread, so the run with four threads repeated the run with one.
The sentence on threads started with the negative. The values now come first, with the number of threads in parentheses.
AGENTS.md asks for no decorative dividers.
petrasovaa
reviewed
Oct 27, 2026
petrasovaa
left a comment
petrasovaa
left a comment
Contributor
There was a problem hiding this comment.
raster/r.mapcalc/tests/r_mapcalc_nprocs_test.py has a skipped test
@pytest.mark.skip(reason="Waiting for a fix to work on all platforms for any nprocs")
The file tests error states, and rand() without a seed is not one since automatic seeding. The test was skipped because parallel rand() failed differently per platform; rand() without a seed is covered by the automatic seed tests, which run with the default nprocs.
Member
Author
|
So, that skipped test can now be removed. It is not an error state anymore because |
wenzeslaus
changed the title
Member
Author
|
This is actually enabling the parallelization itself, not just making it reproducible, so I updated the title. My original title was underselling the impact of this. So the parallelization is all new for rand(), although technically, we could have enabled parallelization with irreproducible results and worse performance after #6480 (that's the gray dashed line in the plots). So, now, rand() no longer slows down r.mapcalc expressions, so expressions with rand() run 2-6 times faster with 8-16 threads on 52-420M cells. An expression with rand() alone runs 1.2 times faster on 52 cells. |
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
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.
rand() drew from the generator shared by the whole program, so r.mapcalc evaluated any expression with rand() on one thread, whatever nprocs was (to avoid the threads hitting the same library generator which was at that time not thread-safe).
Each row now gets a generator state of its own from an exact random number layout of the seed, so the result for a seed is the same as before and the same for any nprocs, and rand() no longer turns off parallel evaluation.
rand() in the calc library needs a modification since that's the place where the RNG library is called. The old global RNG is kept as a fallback, but a new function allows to set a single callback which can return a state to be used with the new RNG API. This is used to return the relevant state for each thread.
The picture shows wall time against the number of threads for two expressions on 52.5 and 420 million cells. This PR is blue and current main, where rand() runs on one thread, is orange. The dashed gray line is current main without the one-thread rule, where the threads share the atomic generator, which costs time when rand() is a large part of the work and many threads draw at the same time (an implementation with no race, but also no reproducibility). The dotted gray line is main from before the atomic generator (#6480) without the rule, where the threads share the generator with no synchronization, which is fast but gives values that change with the number of threads. Only the dotted line still has the shared floating point flag removed in #8031, so its stop at about 110 s on the costlier expression comes from the flag, not from the data race.
Builds on the random number layouts (#7979), the seed option parsing function (#8033), and the integer functions (#8079), all merged.
Details (AI generated)
Design
struct G_random_statelocal to its evaluation at the row's stream before evaluating the row. The streams of an exact layout are the positions of a serial run, which is why the values for a seed do not change. rand() reaches the state through a per-thread pointer. The states are not kept in one array of states, because the draws of up to eight threads would then write to the same cache line (see Performance of the state per row).lib/calc. The library gainscalc_set_random_state(), through which a tool gives rand() a function returning the calling thread's state; r.mapcalc sets it. When nothing is set, rand() draws from the shared generator as before, so other users of the library are unaffected.Alternative layouts
The layout here has a row as its unit, and all rand() calls of a row share its stream. So the values of a seed depend on the number of rand() calls, and adding a call changes the values of the others. The alternative is a layout in which each rand() call has a stream of its own, for example with the call as a batch of the layout and the row as the unit. Then a single rand() would give the same values as the first rand() of an expression with any number of calls. It was not chosen for three reasons.
Seed handling
The seed is read with
G_random_seed_from_option(), so a seed outside -2147483648 to 4294967295 is an error instead of being silently reduced to its low 32 bits, and a seed with trailing characters is an error instead of being cut off. When rand() is used without a seed, the automatic seed comes fromG_random_generate_seed(), soGRASS_RANDOM_SEEDandSOURCE_DATE_EPOCHare honored, and it is still written to the history.nprocsis no longer forced to 1 when the expression contains rand().Documentation
The .md files are the source and the .html files follow them in the random number generator section, which is the only part this PR changes in the .html files. Elsewhere the .html files were already out of sync with the .md files on main, and this PR leaves that as is. Most notably, r.mapcalc.html has no Performance section, so the removed note that rand() disables parallelization has no counterpart there.
Verification
raster/r.mapcalc/tests/r_mapcalc_rand_seed_test.pywere produced by earlier versions, which evaluated expressions with rand() on one thread, so they also check that the values stay the same. They now pass for nprocs 1 and 4. The test that rand() turned off parallel evaluation is removed.raster/r.mapcalc/tests/r_mapcalc_nprocs_test.pyis removed. The file tests error states, and rand() without a seed is not one since automatic seeding. rand() without a seed is covered by the automatic seed tests, which run with the default nprocs, so on all threads.raster/r.mapcalc/tests/r_mapcalc_rand_layout_test.py, 15 tests in about 3 s, comparing exact values with the values of the library's layout for nprocs 1 and 4 (1 for r3.mapcalc, which runs on one thread): expressions of the three types in order, several rand() calls in one expression, if() with rand() in both branches, a variable of eval(), and r3.mapcalc over several depths. It also checks that an automatic seed gives the result of the seed in the history and that invalid seeds are errors.Benchmarks
Two benchmarks of rand() on 1 to 16 threads, run on 2026-10-21 and 2026-10-22 with this PR at 40a3b8a merged with main f557dea. The later commits of this PR change only comments, documentation and tests, so the timed code is the same as at 7d3339c. Both were first posted as comments with the figures below.
Protocol
Shared Intel Xeon W-2295 @ 3.00GHz (18 cores, 36 logical CPUs), Linux 6.8. All builds are Autotools with the same configure line, gcc, CFLAGS
-g -O2, OpenMP on. Onegrass --tmp-mapsetsession on the SECREF North Carolina dataset drives the runs; each run uses the binary and libraries of its build, checked through/proc/<pid>/maps. Builds are interleaved within each number of threads; the figures show the median and the range of the runs.seed=1for all runs.Expressions, each with one
rand(0.0, 1.0):x = rand(0.0, 1.0)on 52.5 million cells (the SECREF extent at 0.1 m)Runs:
x = rand(0.0, 1.0), whose runs take about 1.5 s, and 3 for the others, 630 more runs, load 0.8 to 17.9 (median 8.1).This PR, main, earlier versions and alternatives
This PR is the only version which gives the same values for any number of threads and also gets faster with more threads. The versions which force one thread for rand() stay flat. Current main without the one-thread rule is as fast as this PR on the heaviest expression, but on the cheaper ones it gets slower above 4 to 8 threads, as the threads compete for the shared generator (4 math functions + rand(): 1.55 s at 4 threads, 3.13 s at 16; this PR 1.70 s). The cost depends on how much of each row's work is rand(), not on the number of draws. With one rand() alone, the threads mostly wait for the rows to be written one after another, so few draw at the same time, and at 16 threads it uses as much CPU time as this PR (21 s), although single runs vary from 1.56 s to 2.65 s. With 4 math functions + rand(), many threads draw at once, and 79% of the CPU time is spent in G_drand48(), which advances the shared generator with a compare-and-swap (47 s of CPU time against 24 s for this PR). On the heaviest expression, rand() is too small a part of the work for the draws of the threads to collide often. The version with the data race keeps up with this PR on the cheap expressions, but it is built on main before #8031, so on the costlier ones it stops at 5.1 s and 108 s, like this PR on main before #8031 in the next figure. The values of both versions without the rule change with the number of threads.
31 math functions, 2 graph(), 2 if() + rand(0.0, 1.0), 420 million cells; wall time in seconds, median (min-max) of 3 runs:
4 math functions + rand(0.0, 1.0), 52.5 million cells; wall time in seconds, median (min-max) of 3 runs (6 at 1 thread before OpenMP):
x = rand(0.0, 1.0), 52.5 million cells; wall time in seconds, median (min-max) of 5 runs:
With one thread, this PR is faster than current main (one rand(): 1.48 s against 1.76 s), most likely because it does not advance a shared generator with atomics.
This PR before and after #8031
#8031 removed a shared floating point flag written by all threads in r.mapcalc. With main before #8031, this PR gets the costlier expressions down to 2 to 3 times faster than main, and then stops. With #8031, the same PR goes on: 6.2 times faster for 14 math functions + rand() and 4.4 times for the heaviest expression. The cheap expressions are the same on both. #8031 alone does not change main, since rand() still forces one thread there. The values are the same for any number of threads in all builds.
Main is f557dea; the builds before #8031 have its commit bf7d8b6 reverted, so each pair differs only by #8031.
14 math functions + rand(0.0, 1.0), 52.5 million cells; wall time in seconds, median (min-max) of 3 runs:
31 math functions, 2 graph(), 2 if() + rand(0.0, 1.0), 420 million cells; wall time in seconds, median (min-max) of 3 runs (5 at 5 and 12 threads for this PR on main before #8031):
For one rand() and for 4 math functions + rand(), this PR gives the same times on both mains (1.34 s at 2 threads for one rand(), about 1.6 to 1.7 s from 5 threads on for 4 math functions), so #8031 matters only when there is more work per cell.
The cheaper expressions get slower after their fastest thread count (e.g. 4 math functions + rand(): 1.56 s at 3 threads, 1.70 s at 16). The same expressions without rand() do the same on main with #8031 and on this PR (4 math functions: 1.54 s at 3 threads, 1.70 s at 16;
x = 1.0: 0.62 s to 0.69 s), so this comes from r.mapcalc, not from the random numbers.Values
Values of
c = rand(0.0, 1.0)on 52.5 million cells, compared cell by cell: every build at nprocs 1 gives the same map as this PR at nprocs 1, including r.mapcalc before OpenMP and the builds before #8031. Between nprocs 1 and 16, the map differs in 0 cells for this PR and the builds which force one thread, in 128,062 cells for main without the rule and in all 52,500,000 cells for the version with the data race.Performance of the state per row
An expression with 48 rand() calls per cell on a region of 2000 by 8000 cells, best of three runs on a 36-core machine, with the same results in all variants. Adjacent states is an earlier version of this PR with one array of per-thread states. States 64 bytes apart is the same array with the states a cache line apart, tried only to check the cause. State per row is the code of this PR.
The states a cache line apart are as fast as the state per row, so the slowdown came from the threads writing to the same cache line. Above 8 threads the time no longer drops, most likely because the rows are written to the output one after another.
On one thread, the state per row costs nothing measurable compared with the array of states: +0.1% with
-O2and +0.2% with-O3for 14 math functions + rand() on 52.5 million cells, medians of 12 runs per build with the build order reversed halfway, between two builds which differ only in this change. Separately linked builds can differ by about 5% on one thread without any change in work. CallingG_random_uint32()from the calc library adds one imported function, which moves all code of the library 32 bytes later; with the same 118 billion instructions, the same expression then took 4% to 5% more cycles, all of it in the calc library. This is code placement, and other changes to the library can move it either way.Suggested squash commit message
This description and the code were written with Claude Code (Claude Opus 5.5 and earlier Claude models) and reviewed by the author of the PR.