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]

So that we don't forget about it, we need to add this to raster/CMakeLists.txt

[petrasovaa]

I would also like to see results for different values (like 5, 15, 30 or something like that), I assume that would influence the parallel efficiency.

[petrasovaa]

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

[petrasovaa]

If I am understanding this correctly, the memory budget calculation only accounts for the output chunk buffer but not the per-thread input row buffers (nprocs × row_scale × src_w.cols), which can dominate memory usage at high resolution ratios.

[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.