r.proj: memory-bounded parallel banding#7627
Conversation
…re Apple M-series
Replace the whole-map RAM buffer (Path A) with a two-level band loop modeled on r.neighbors, adapted for r.proj's CRS-dependent input access: each output band's input footprint is found by back-projecting the band's edges (dense edge walk), so the loaded input strip is sized per band rather than by a fixed neighborhood stencil. Band height adapts to the memory option; if a single output row's footprint exceeds the cap (oblique or large-halo transforms) the module bails, since that case needs the tile cache path, which is not implemented here. Per-thread PROJ contexts (one PJ clone per thread) are retained. Input strips are loaded serially per band because a single fd read path is not thread-safe; each band's output is written in order. Bit-exact against the serial output at 1/2/4/8 threads, for both column-varying and row-varying inputs, with multiple bands exercised. On the test case (105.3M output cells, EPSG:4326 to EPSG:3857, nearest, memory=50) peak RSS was 130 MB versus 763 MB for the whole-map buffer. Developed with assistance from Claude (Anthropic).
krcoder123
changed the title
github-actions
Bot
added
raster
| G_malloc((size_t)band_orows * outcellhd.cols * cell_size); | ||
|
|
||
| double t1 = omp_get_wtime(); | ||
| #pragma omp parallel |
There was a problem hiding this comment.
Why nested parallelism? How is it supposed to work?
https://docs.oracle.com/cd/E19205-01/819-5270/aewbc/index.html
There was a problem hiding this comment.
I didn't use nested parallelism here. The band loop is serial, and each band runs one parallel region. The omp for inside it does not create a second team, it just splits the loop rows among the threads that omp parallel already made. So it is one level of threads, same as writing parallel for, just split into two directives. I split them because each thread has to clone its own PROJ context before the loop and destroy it after, and that setup has to live inside the parallel region but outside the for. I also grepped the functions called inside the region (GPJ_transform, the strip reader) and none of them open their own parallel regions.
Also because the region sits inside the band loop, the team and the PROJ clones get recreated every band instead of once. The runtime reuses its thread pool so this is cheap, and in the benchmark it gets absorbed into the compute phase, which still scales 5.2x. But if you would rather have the region put above the band loop so everything is created all at once I can do that too.
| double t0 = omp_get_wtime(); | ||
| G_switch_env(); /* -> input */ | ||
| for (int r = imin; r <= imax; r++) | ||
| Rast_get_row(fdi, |
There was a problem hiding this comment.
Please see if it's possible to parallelize reading.
There was a problem hiding this comment.
I don't think reading on a single file descriptor is thread safe. The read path keeps descriptor's state so reads happening at the same time can cause a race. But every thread having its own descriptors like r.neighbors should work. I will try that and benchmark it against the serial load.
There was a problem hiding this comment.
I tried the per thread descriptor approach like r.neighbors and it works. Each thread opens its own descriptor and reads a disjoint block of the band's strip. The read phase went from about 0.50s to 0.26s at 8 threads, so about 1.9x. It doesn't scale fully linearly because the reads share one SSD, and the curve flattens like the SSD is saturating, not like a fixed cost. But a better device may scale a bit more linearly. Regardless the total is now about 3.0x vs serial. So the overall speedup has some improvement because the input read is parallelized. I updated the table in the PR description with the new numbers, all measured with the same LZ4 compression including the serial reference. I also added a "Input read speedup" column.
|
@krcoder123 In your result table, please add a new column for |
| /* Inside the map but outside the loaded strip: the span check under-sized | ||
| * the strip. This is a correctness failure, not a NULL. */ |
There was a problem hiding this comment.
I'm not entirely sure what you mean here.
There was a problem hiding this comment.
The check above it already handles coordinates that fall outside the input map, those become NULL like normal. This second check catches an input row that is inside the map but not a part of the rows we preloaded for this band. That should be impossible if the band's footprint estimate was right, so if it ever triggers it means the estimate was wrong. That is a bug in the band sizing, so the module fails with an error instead of writing a NULL and silently producing wrong output. I re wrote the comment to make it clearer now.
Does it mean the current implementation only works for +-1 rows? +-15% of the entire map, you said. Please consider a fallback parallelization method. Row-band then column-wise parallelization. If this performs comparably or even faster, I would only use this method. If this is slower than your current row-wise parallelization, the column-wise method as a fallback. |
Added a compute speedup column in between "Compute (s)" and "Output write (s)". |
I measured the input row span for every output band across a few CRS pairs. The current banding works fine for pairs near 0% like 4326 to 3857 and about 3 to 7% for UTM or LAEA. The problem is polar transforms, where one output row can touch 50 to 80% of the input rows and no band fitting the memory cap covers that. This is why a fallback approach is needed. The column-wise idea sounds interesting. I will benchmark it against the row-wise method and let you know what I find. |
The banding timers call omp_get_wtime(), which is undefined when GRASS is built without OpenMP and breaks the link in the minimum-config build. Wrap omp.h and omp_get_wtime() behind _OPENMP via a small rproj_wtime() helper that returns 0.0 without OpenMP. Also demote the PHASE_TIMERS line from G_message to G_debug, since it is benchmark scaffolding, not user output.
r.proj now calls proj_* functions directly for per-thread PROJ contexts, so its CMake target must link PROJ::proj. Without it the CMake build fails to link with an undefined reference to proj_context_destroy. The Autotools build already links libproj via the module Makefile.
| * read-only shared; the static METERS_in/out race is benign | ||
| * (constant CRS per run). */ | ||
| struct pj_info tproj_local = tproj; | ||
| PJ_CONTEXT *thread_ctx = proj_context_create(); |
There was a problem hiding this comment.
We should probably create wrappers in gproj library instead of calling proj library directly. (That would also make the dependency changes in Makefile and CMakeLists.txt unnecessary).
There was a problem hiding this comment.
Thank you for the suggestion. I implemented this in 98d43be. I moved the context handling for each thread into gproj as GPJ_clone_transform and GPJ_free_transform_clone. The direct PROJ dependency lines in the Makefile and CMakeLists are removed since r.proj now gets PROJ transitively through libgrass_gproj.
PROJ transformation objects are not safe for concurrent use, so a parallel module needs a private clone per thread. Add GPJ_clone_transform() and GPJ_free_transform_clone(), which bundle a cloned transform with its private PROJ context in struct gpj_transform_clone so ownership is a single unit. r.proj's parallel banding created the per-thread context with proj_context_create(), proj_clone(), and proj_destroy() directly; switch it to these helpers so the PROJ calls live in lib/proj and the module makes none.
…proj r.proj no longer calls proj_* directly (the per-thread context handling moved to GPJ_clone_transform/GPJ_free_transform_clone in lib/proj), so it links PROJ transitively through libgrass_gproj. Remove the explicit PROJLIB (Makefile) and PROJ::proj (CMake) dependencies added earlier to satisfy the direct calls.
Each band's input strip is read in parallel: read_nprocs fresh per-thread fds (Rast_open_old), a static block split of the strip rows across threads, each thread reading its disjoint rows through its own fd into its own strip slice. Rast_disable_omp_on_mask gates the parallelism (serial when a mask is present or without OpenMP); fdi remains the serial-fallback path. Env-switch choreography, compute region, write loop, band sizing, and PJ context cloning are unchanged. Experimental, not for merge.
The memory-bounded banding path halves the band height until a full-width input strip fits the memory cap. On oblique projections a single output row can back-project to an input footprint larger than the cap at any height, and that path bailed out. Add column tiling as a second search dimension. Phase 1 is identical to the current banding: halve the band height while the input strip spans the full input width, and use that result whenever a full-width band fits. Phase 2 runs only when a single full-width row still exceeds the cap. It keeps the band as tall as its output buffer allows and halves the tile width instead, so each column tile back-projects to a smaller input row span and the per-band parallel region stays populated rather than collapsing to a single row. The width search runs in two tiers to stay cheap. The upper tier estimates the worst tile strip by probing a bounded, evenly spaced subset of tiles, which is a lower bound on the true worst, and narrows to a candidate width. The lower tier validates that width with the exact per-tile edge walk and narrows further if the estimate was optimistic, so the accepted width is always exact-sized against the cap. Input strips stay full width because the raster API reads whole rows, so a tile strip is its input row span times the full input width, and tiling shrinks the row span rather than the width. Each tile loads one strip whose row span comes from the exact edge walk, one tile at a time, bounding peak memory to the worst tile rather than the whole band. Every output cell is computed once and written in row order, so the result is bit-exact with the serial output. Retain the existing fatal error only for the degenerate tile whose footprint cannot fit the cap at any width, at minimum band height; that footprint needs the tile-cache path, which is not implemented.
I implemented the row band then column-wise parallelization you suggested. The way it decides whether to do column wise or row wise is before loading anything, it checks if the next band of output rows can be processed with its input staying under the memory cap. If yes at full width, then we can go ahead with row wise, which is how it was even before. But now, if even a single full width row needs more input than the cap allows, the band gets split into column tiles until each tile's input fits.This is the column-wise parallelization. The only case that still doesn’t work is a tile whose input footprint contains the pole point, since every longitude converges there, so no tile width brings its footprint under the cap. This would be the North and South Poles. Column-wise works for every other CRS pair though. Output is the exact same against serial r.proj on an easy pair (no tiling), LAEA (250 tiles) and Albers (492 tiles), at 1 to 8 threads with both test input patterns. On easy pairs the speed is comparable to the current row-wise code (back to back comparisons, within noise, same code path when tiling does not engage). Based on what you said initially, would it make sense to make column-wise the primary method? For the tilted CRS pairs I tested with a wide extent LAEA job. Previously the row wise approach didn't work on this at memory=50. With column tiling the same job finishes nicely. For a serial reference I ran the old serial r.proj on the same dataset. It takes 10.24 s but needs the whole input held in RAM, so it ran with memory=2000 while the tiled runs used memory=50.
At 8 threads it is 1.56x faster than the whole map serial using a 1/40th of the memory, with compute scaling 5.38x. The remaining serial cost is the fit search (2.61 s) plus the write. The search can come down more, since adjacent bands almost always need the same tile width. So working on that bring down that sizing time is my next task. Write cannot be parallelized just like before. |
3 participants
This is an early stage draft for parallelizing r.proj which I’m still debating on the method of parallelization for.
This branch parallelizes r.proj with OpenMP using a two level banded structure similar to r.neighbors, but it’s adapted for r.proj. Unlike r.neighbors, where the window size is fixed, the input rows needed by an output band depend on the CRS pair. So the band size is computed for every band by projecting the band's edges back into the input and measuring which input rows it touches. Each thread clones its own PROJ context so that PJ_CONTEXT blocker is now solved. Memory stays within the memory option because only the current band's input strip and output buffer are in RAM at any time.
Correctness:
The output is the exact same against serial r.proj at 1, 2, 4, and 8 threads, with both column-varying and row-varying test inputs. It’s verified with multiple bands (21 to 29) and at 105M output cells. EPSG:4326 to EPSG:3857.
I measured it on an 8-core Apple M-series (results below): peak RSS drops from 763 MB (whole map in RAM) to 130 MB with memory=50 on the same test case. Total speedup at 8 threads is 2.99x. Compute alone scales 5.36x, so the problem is the output write.
All runs including the serial baseline use LZ4 compression so every number shares one compressor. Earlier tables in this PR compared LZ4 runs against a ZSTD compressed serial reference, which inflated the speedups slightly since ZSTD is slower to write. I chose LZ4 for benchmarking because it has the lowest compression cost and adds the least noise to the phase being measured.
Just a note, the totals in this table count only the module's own work, which is the sum of the four phase columns. They leave out about 0.48 s of GRASS process startup, which is the same at every thread count and has nothing to do with r.proj itself. Leaving it in would just mess up every speedup number equally without saying anything about the parallelization.
As the number of threads increases, the compute time in seconds goes down but the output write still stays the same and takes up a big block of time. The output write is serial by design in the raster library (Rast_put_row is sequential and compressed row offsets depend on prior rows), so it becomes the floor as threads increase. r.neighbors has the opposite case. Its work per cell is heavy in time so the serial write is a small fraction of its runtime and barely affects its speedup.
Current scope: This method only works for when the transformations are relatively close by. Transforms between very far projections or polar transforms can scatter one output row's reads across most of the input. To address that I’m thinking of using a tile-cache path, which I didn’t implement yet. For now the code detects this case and just exits with a message.