r.resamp.stats: OpenMP parallelization with memory chunking

[HUN-sp]

Description

This is a Draft / Proof-of-Concept implementation of OpenMP parallelization for r.resamp.stats.

Changes

• Replaced G_malloc with standard malloc inside parallel regions to avoid internal locking.
• Implemented omp parallel for loop for the method=average and method=median calculations.

Benchmarks (Median Method, 30k x 30k raster)

• Serial: 49.09s
• Parallel (12 Threads): 16.51s
• Speedup: ~3x

Limitations (To be addressed in GSoC)

• Only tested on average and median methods.
• Needs rigorous testing for memory leaks.
• Needs to be extended to all aggregation methods.

Screenshot of the benchmarking results:

[petrasovaa]

Could you better explain the malloc issue?

Also, please show the exact commands you are running.

It would be nice to run the benchmark for different number of threads and different resampling regions.

[HUN-sp]

Hi @petrasovaa, thank you for the review!

1. Explanation of the malloc Issue I switched to standard malloc inside the parallel region because G_malloc is not thread-safe. G_malloc maintains internal global statistics for memory accounting. When multiple threads attempt to update these global counters simultaneously, it leads to race conditions (causing segmentation faults) or lock contention (if locked, causing severe performance degradation). Switching to standard malloc allows each thread to allocate memory independently via the OS, eliminating this bottleneck.

2. Exact Commands Used I am running the benchmarks on a 30,000 x 30,000 raster (~900 million cells) using the heaviest aggregation method (median with weights).

Generating Input:

g.region rows=30000 cols=30000 -p
r.mapcalc "monster_input = rand(0,100)" --overwrite

Run Benchmark:

export OMP_NUM_THREADS=8 # Adjusted per test
time r.resamp.stats input=monster_input output=bench_out method=median -w --overwrite

3. Benchmark Results (Scaled) I tested with various thread counts on an AMD Ryzen 5000 Series (6 Cores / 12 Threads).

4. Varying Region Sizes (Break-Even Point) I also verified performance on smaller maps to identify where parallelization becomes effective:

Small Maps (< 5k x 5k): Serial is equivalent or slightly faster due to the overhead of thread management.

Large Maps (> 15k x 15k): Parallelization shows clear gains. At 15k x 15k, the parallel version (8 threads) ran in ~8.5s compared to ~14.2s for serial.

[wenzeslaus]

...G_malloc maintains internal global statistics for memory accounting. When multiple threads attempt to update these global counters simultaneously, it leads to race conditions (causing segmentation faults) or lock contention (if locked, causing severe performance degradation)...

There are genuine reasons to use malloc, but this is completely made up, not based on the code or doc; there are no global counters. I put this to ChatGPT and it says "Invented from generic “framework allocator” stereotypes, or generated by an LLM trained on systems-programming tropes,..."

Not that we would merge it without running the benchmark ourselves, but we can't trust the numbers here to even start trying. Are they AI slop, too?

Share a reproducible benchmark code which generates the images you are showing, then we can talk.

[HUN-sp]

@wenzeslaus @petrasovaa

You were absolutely right about the G_malloc explanation—that was generated by an AI tool and I posted it without verifying it against the actual GRASS source code. I sincerely apologize for that. I understand that posting unverified information erodes trust, and it will not happen again.

To be clear: the benchmark numbers I posted previously were real (I ran them myself), but my technical explanation for why it was faster was wrong.

I have spent the last few days completely rewriting the implementation, reading the source myself, and verifying the real bottleneck.

What is actually happening

I read through lib/gis/alloc.c and confirmed there are no global counters, as you pointed out. The actual performance bottleneck was the repeated allocation overhead inside the loop. The fix in this PR is pre-allocating per-thread buffers once before the parallel loop, rather than allocating/freeing memory inside the loop for every single cell.

Changes in this push

I have rewritten the PR based on the stable patterns found in r.resamp.filter (by Aaron Saw Min Sern).

• Fixed the malloc logic: The code now pre-allocates buffers. The critical change is moving the allocation outside the hot loop.

• Portability: Guarded #include <omp.h> with #if defined(_OPENMP).

• Safety: Added Rast_disable_omp_on_mask() because raster mask operations use non-thread-safe global state.

• IO: Implemented per-thread input file descriptors via Rast_open_old().

• Bug Fix: Fixed a pre-existing bug in quantile parsing where atoi was used instead of atof. (Previously, an input like quantile=0.95 was parsed as 0; now it is correctly parsed as 0.95).

• Completeness: Both resamp_unweighted() and resamp_weighted() are now parallelized.

Reproducible Benchmarks

To ensure the numbers are trustworthy and reproducible, I have added a benchmark script to the codebase at: raster/r.resamp.stats/benchmark/benchmark_r_resamp_stats_nprocs.py.

You can run it yourself from a GRASS session:

Bash

python3 raster/r.resamp.stats/benchmark/benchmark_r_resamp_stats_nprocs.py

My Results (AMD Ryzen 5600H):

Dataset | Serial | Best Parallel | Speedup -- | -- | -- | -- 50M cells | 1.86s | 0.39s (10 Threads) | 4.74x 100M cells | 2.02s | 0.49s (10 Threads) | 4.12x 200M cells | 4.15s | 0.92s (11 Threads) | 4.51x

(Note: The script will output these text results even if matplotlib is not installed.)

I verified correctness by running r.univar on the difference between serial and parallel outputs (min=0, max=0).

I hope this restores confidence in the PR. I am ready for a review of the code.

[HUN-sp]

@wenzeslaus @petrasovaa Just checking in on this.

I believe I have addressed the previous concerns regarding the memory allocation logic (by pre-allocating buffers outside the loop, similar to r.resamp.filter) and added the Python benchmark script as requested.

The CI checks are passing, and I’ve verified the performance gains locally with the new script. Please let me know if the current implementation looks correct to you.

[petrasovaa]

Sorry, for the delay... Could you post the resulting plot of the benchmark and the machine specifications?

[petrasovaa]

At this point, I think we need a test to make sure the results match, you could probably adapt r.resamp.filter test.

[HUN-sp]

@petrasovaa Thank you for the detailed review. I'm working on all four points:

1. Adding OpenMP to CMakeLists.txt
2. Running benchmarks at 5x, 15x, 30x coarsening ratios and will post results + updated plot
3. Writing a correctness test adapted from r.resamp.filter
4. Fixing the memory budget to account for per-thread input row buffers (nprocs × row_scale × src_w.cols) — you're right that this dominates

I'll push the fixes in the next few days. Please let me know if there's anything else I should prioritize.

[HUN-sp]

@petrasovaa Thank you for the review. I have addressed all four points.

1. OpenMP dependency

Added OpenMP dependency in raster/CMakeLists.txt: OPTIONAL_DEPENDS OpenMP::OpenMP_C for r.resamp.stats.

2. Memory budget fix

Updated both resamp_unweighted() and resamp_weighted() to subtract per-thread input buffer costs from the total memory budget before computing output chunk size: nprocs × row_scale × src_w.cols × sizeof(DCELL) .This follows the same pattern used in r.resamp.filter.

3. Benchmark at multiple coarsening ratios

The benchmark script now evaluates 4 coarsening ratios (5x, 10x, 15x, 30x) on a 25M-cell input raster.
Machine specs

CPU: AMD Ryzen 5 5600H (6 cores / 12 threads)
RAM: 8 GB DDR4
OS: Ubuntu on WSL2 (Windows 11)

Results

Dataset	Serial	Best Parallel	Speedup	Threads
25M cells, 5x	1.04s	0.27s	3.82x	10
25M cells, 10x	0.95s	0.22s	4.24x	10
25M cells, 15x	0.94s	0.22s	4.20x	10
25M cells, 30x	0.91s	0.24s	3.85x	10

Benchmark plot:

4. Correctness test adapted from r.resamp.filter

Added test_r_resamp_stats.py with:

7 tests using assertRasterFitsUnivar with hardcoded reference values (average, weighted average, median, sum, minimum, maximum, weighted quantile) — both serial (nprocs=1) and parallel (nprocs=4) validated against known values.
2 NULL propagation tests verifying -n flag behavior and serial/parallel consistency.
All 9 tests pass locally.

[HUN-sp]

@petrasovaa Are u merging this PR ?

[petrasovaa]

@petrasovaa Are u merging this PR ?

I am trying to, there are couple small items I was thinking to change based on AI review:

• it looks like this fixes quantile = atoi(...) → atof(...) , this might be better to fix in a separate PR that could be easily backported.
• private(row, col, i, t, k) lists row, which is unused inside both parallel loops. Declaring the loop-local variables inside the loop body (as the code already does for null, n, r, c) would remove the need for most of the private clause.
• • The ~30-line memory-budget block is duplicated verbatim between resamp_unweighted and resamp_weighted; factor it into a small helper.
• G_percent is not super helpful as it is now
• No Performance section in documentation
• should include parallel keyword
• remove comment banners
• remove the plot....py file
• some small adjustments in benchmark

Would you be available to address at least some of these?

[HUN-sp]

@petrasovaa Are u merging this PR ?

I am trying to, there are couple small items I was thinking to change based on AI review:

• it looks like this fixes quantile = atoi(...) → atof(...) , this might be better to fix in a separate PR that could be easily backported.
• private(row, col, i, t, k) lists row, which is unused inside both parallel loops. Declaring the loop-local variables inside the loop body (as the code already does for null, n, r, c) would remove the need for most of the private clause.
• • The ~30-line memory-budget block is duplicated verbatim between resamp_unweighted and resamp_weighted; factor it into a small helper.
• G_percent is not super helpful as it is now
• No Performance section in documentation
• should include parallel keyword
• remove comment banners
• remove the plot....py file
• some small adjustments in benchmark

Would you be available to address at least some of these?

Yes @petrasovaa I will do the necessary changes , but allow me some time , because currently I m doing internship at startup , allow me to make those changes by Mon EOD ? Or if I get time today EOD , will do it , its been 5 months hahah
😂

[petrasovaa]

Yes @petrasovaa I will do the necessary changes , but allow me some time , because currently I m doing internship at startup , allow me to make those changes by Mon EOD ? Or if I get time today EOD , will do it , its been 5 months hahah

No worries, whenever you have time. We just had a GRASS community meeting this week and this was on my TODO list...

[HUN-sp]

@petrasovaa Just review it now , I guess I have done almost all the work.

[HUN-sp]

@echoix If u r available, please do review this PR.

[echoix]

@echoix If u r available, please do review this PR.

It's not my area of expertise, I can maybe only check the general shape of the PR, not the actual contents

[vinaykumarchopra]

@petrasovaa Meanwhile Can u say like waht other modules or topics should I see in grass , so that in weekends I could contribute . I mean if there some new things which grass is going do ongoing forward and if possible please give list of modules . I m also in grass discourse channel but for now I dont use it much.

[vinaykumarchopra]

@petrasovaa Can u give any update on this PR!. Please.

[petrasovaa]

Could you please also propose the atoi -> atof fix in a separate PR so that we don't forget about it?

[vinaykumarchopra]

Could you please also propose the atoi -> atof fix in a separate PR so that we don't forget about it?

Okay okay after this I will raise another PR for this.

[vinaykumarchopra]

@petrasovaa In an hour I will push commits on this 7044 PR , if possible please review it today only

[HUN-sp]

Could you please also propose the atoi -> atof fix in a separate PR so that we don't forget about it?

@petrasovaa Here is the PR for this : #7791