r.proj: memory-bounded parallel banding - #7627
Conversation
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. |
| * 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.
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. The table in this comment is a harder tilted pair compared to the pair in the PR description, EPSG:4326 to EPSG:3035 (lat/lon to Lambert Azimuthal Equal Area over a wide European extent)
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. |
|
I also wanted to clarify how I verified the outputs more clearly. I subtracted the tiled output map from the serial output map in GRASS and check that the difference is 0 everywhere, meaning both min and max of the difference map are 0. I also checked that the null cells are in the exact same places and that both maps have the same number of cells. I did this at 1, 2, 4 and 8 threads with two different test inputs. The md5 part is just for the binaries, I hash them so I always know exactly which build produced which benchmark number since I have to try multiple different maps. Sorry for the confusion. |
|
@krcoder123 Please avoid using CI time for a draft PR, it doesn’t need to be kept updated with every single merged PR, especially when the CI queue is a bit overloaded (that would mean exponential CI usage on each merged PR). We usually don’t mind when it’s combined with actual real work or a real chance of semantic conflicts (logic would break when combined to other changes even if textually it merges fine) |
Understood, sorry about that. The update branch I clicked were on PR 7440, I'll stop doing that until review starts there. The push for this draft PR was fixing a bug and adding some new stuff too. |
|
Just an update on the pole handling and benchmarks. The most recent commit can now handle pole maps with a north or south pole inside the output frame, which previous versions couldn’t do. The reason it couldn’t before is because of a bug in how each tile estimates which input rows it needs. The estimate only checked the tile's edges, but for a tile that contains the pole, the northern most input row it needs sits in the middle of the tile, not on the edge. So the code loaded too few rows and stopped with an error. The fix now finds where the pole lands in the output once per map, and when a tile contains it, extends that tile's row range to cover it. The column-wise parallelization can now handle easy, hard and very tough maps like pole maps. They all match serial r.proj output exactly, tested at memory from 2 MB to 50 MB at 1 and 8 threads. I also brought the sizing portion from 2.6 s to 1.6 s on a hard map by reusing the previous band's size between bands. Here are some of the benchmarks I did: Easy map: EPSG:4326 to EPSG:3857 (lat/lon to Web Mercator), 105M output cells, memory=50, serial baseline 7.503 s
Hard map: EPSG:4326 to EPSG:3035 (lat/lon to Lambert Azimuthal Equal Area, wide European extent), 76M output cells, memory=50, serial baseline 11.537 s
Pole map: EPSG:4326 to EPSG:3413 (lat/lon to Polar Stereographic North, north pole inside the frame), 100M output cells, memory=50, serial baseline 29.252 s
The one thing I want to make clear is that for hard maps and pole maps while the parallelization is good, N=1 becomes slower than serial r.proj. This is because of how the memory option works. This code splits the output into horizontal bands so that only a slice of the input has to be in memory at a time, and neighboring bands need some of the same input rows, so those rows get read from disk more than once. A single thread has to do all of that re-reading alone, while serial r.proj reads the whole input into RAM once and never re-reads. But I did verify the slowdown is purely the memory cap and not the parallel code. With the memory option set high enough, N=1 matches serial (11.28 s vs 11.54 s on the hard map) and every other thread count gets faster too. For example I verified that 8 threads goes from 2.0x to 3.5x. So users can trade memory for speed with the existing memory option if they wanted to. One note on the tables, my original hard map test dataset was lost, so the table in this comment uses a rebuilt 4326 to 3035 dataset of the same dimensions with a fresh serial reference. Absolute times are not comparable to the earlier table but the scaling behavior is the same. I used Claude (Fable) to help with testing and benchmarking throughout, all results verified against serial output. it also helped me find some logic errors which are now fixed. I am away for a few days but I'll keep an eye on this PR. But I wanted to ask if you think a fallback method like tile cache is still needed given that this parallelization method can handle poles now, and if so for which cases? Thanks |
|
I am posting the AI suggestion about the implementation structure. Lot of stuff to sort out, hopefully it will make more sense to you when you are familiar with the code base. Couple comments:
Details1. Keep the input resident across bands (rolling window)Problem it fixes: the N=1 regression on hard maps (0.65x). Each band's strip is malloc'd, read, and freed. Adjacent bands need overlapping input rows — on the 3035 map enough overlap that input read is 5.5 s where old serial r.proj spent ~0.5 s reading each row exactly once. How: bands are processed top-to-bottom, and for any continuous transform the input row span advances mostly monotonically with the output row. So instead of malloc → read [imin, imax] → free, keep the strip as a sliding window: when band N+1 needs [imin', imax'] and the buffer holds [imin, imax], memmove the overlap down and read only rows (imax, imax']. Same memory cap (the window is never larger than one band's strip), no new data structure. Impact: each input row is decompressed roughly once for row-banded cases, which should bring N=1 back to parity with serial 8.4. That matters because a parallel rewrite that's slower at defaults than the code it replaces is a hard sell, and r.proj is typically a run-once import step. Cost: small — a memmove and bookkeeping. Doesn't help the column-tiled case (next item).
Problem it fixes: column tiles have structurally unbounded read amplification, and there's still a hard failure path. Why the amplification: strips are full input width (the raster API reads whole rows — his own comment), so column tiling only narrows the row span per tile. Every one of the 250–500 tiles re-reads complete rows, and neighboring tiles' row spans overlap heavily. The narrower the tiles get (i.e., the harder the transform), the more times the same row is decompressed. The design gets more wasteful exactly when the transform is hardest. The rolling window doesn't help here because within a band, tiles' spans jump around rather than advance. How: this is precisely what the existing readcell structure already does — a 2D LRU cache of decompressed tiles under the memory cap, where re-touching a row is a cache hit. Instead of deleting it, put it behind the strip loads (or use it directly for the tiled path): banding decides the parallel compute schedule, the cache decides what's resident. Those are separate concerns; the current design fuses them, and the fusion is where the waste comes from. One caveat to design for: the old cache faults tiles in lazily, which isn't thread-safe — either pre-fault a tile's needed input serially before the parallel region (the band sizing already knows the footprint), or make it read-only during compute. This also answers his direct question ("is a tile-cache fallback still needed now that poles work?"): yes — no longer for correctness, but as the read-amortization layer. And it removes the remaining G_fatal_error for the over-cap single row: old r.proj handled memory=2 on any CRS pair; the new code must not do worse.
Problem it fixes: the easy-map ceiling. Write is ~0.93 s at every thread count; at 8 threads compute is ~1.03 s, so total is ~2 s of work serialized into 2.76 s, capping the flagship benchmark at 2.72x. How: Rast_put_row must stay sequential (compressed row offsets depend on prior rows) — but nothing says it must wait for compute. Double-buffer the band output: while the worker threads compute band N+1, one thread writes band N's rows in order. OpenMP sections or tasks do this without new dependencies. Total goes from compute + write per band to max(compute, write). Impact: easy map at 8T: ~2.76 s → ~1.8 s, i.e. 2.72x → ~4x. It's the cheapest big win in the PR Impact: sizing drops to a single pass with trivial queries; simpler and faster than the seeded search he's currently tuning. It also naturally produces exact band boundaries instead of the power-of-two lattice, so bands are as tall as the cap truly allows.
Problem it addresses: the thing being parallelized may be mostly avoidable overhead. 5.5 s of co3857 is ~53 ns/cell, and the actual math for that pair is a handful of flops. The rest is per-cell GPJ_transform overhead: a strncmp on the proj string, proj_angular_input(), PJ_COORD marshalling, and a single-point proj_trans call — per cell. How: transform a row (or tile row) at a time with proj_trans_array-style batching — GPJ_transform_array already exists in the library. The angular/units decisions get made once per batch instead of once per cell. Why now and not later: if batching cuts serial compute 2–4x on cheap pipelines, the compute/write/read balance the whole architecture is tuned around shifts — the write overlap (item 3) becomes more important, and the case
Problem it fixes: the PR's central thread-safety claim is currently false. GPJ_transform assignsS_in/METERS_out on every call (lib/proj/do_proj.c:42), so the parallel loop races on them despitethe per-thread PJ clones. It happens to work here because every thread writes identical values, but it's UB, TSan will flag it, and the new GPJ_clone_transform API documents a safety guarantee the library doesn't provide. How: make them locals of GPJ_transform (and check the other entry points that touch them). Trivial diff, but it belongs in this PR, because the entire design — and the new public API — rests on "clone the transform and the Here are some bugs in the current code: Details- All non-nearest methods silently produce nearest-neighbor output. The old dispatch interpolate = menu[method].method was deleted, but the method option still accepts bilinear, bicubic, lanczos, and the _f variants (raster/r.proj/main.c:83), and the compute loop unconditionally calls interpolate_strip() (raster/r.proj/main.c:1289), which is nearest-only. A user asking for bicubic gets nearest with no warning. Either implement the strip-based interpolators (the 2-cell margin is already reserved for their stencils) or G_fatal_error on non-nearest methods until they exist. This alone makes the PR unmergeable as-is; a test with method=bilinear would have caught it. - Functional regression for oblique/large-halo transforms. The old readcell tile cache handled any CRS pair within the memory cap by faulting tiles in from a temp file. That fallback is deleted, and the code now aborts (raster/r.proj/main.c:1167) for cases that previously worked, just slowly. The PR description acknowledges this ("tile-cache path not implemented"), but the fix should be to fall back to the existing readcell path (still in the tree, now dead code) rather than to error out. Merging a speedup that removes working functionality is not acceptable for this tool. - GPJ_transform is not actually thread-safe, even with a cloned PJ. It writes the file-scope globals METERS_in/METERS_out on every call (lib/proj/do_proj.c:42, :891). In this PR all threads write identical values so the output is correct in practice, but it is a formal data race (TSan will flag it) and — more importantly — the comments and the new API's docs claim per-thread safety that the library does not deliver. Making those two variables locals of GPJ_transform looks straightforward and would make the claim true. This needs to be part of the same PR. |
e666323 to
c980ea6
Compare
|
Thanks for the review. I added three major things in my most recent commits. All the non-nearest methods now run in the parallel path. This fixes a bug where bilinear, bicubic, and lanczos silently fell back to nearest. I also found and fixed a tiny rounding mismatch in the output row coordinates, off by one bit. And I added the serial tile cache fallback, so the parallel path never fails on a case the old r.proj completed. But I found and verified that these cases that use the fallback are extremely rare. It only happens when the memory cap cannot hold one output row's input strip. So for those cases the module now warns with the exact memory value that would keep the parallel path, then finishes the run through the old readcell cache. I swept the caps to see how often this fires. It never triggered on any real projection pair down to 1 MB. Only a purposely wide input that’s wider than cap/20 columns hits it. I also went through your review items. For number 5, which was the batching the transform calls, I measured it at about 1.39x per point. This is smaller than what the AI review had mentioned it could be. Also in general that’s a pretty small end-to-end gain for something that needs a large code refactor. So I'm not planning it right now, but I'm happy to do it if you still want it. I will work on number 1 and 3 for this week's work since that's the only thing that can help bring the total speedup up some more. On number 4, I agree the footprint grid is cleaner, but it replaces the sizing that all the results below were tested against, and with the July 24 target Huidae and I set, it seems a bit tight. But I can try if you would like to see how it looks. I’ll finish with the globals fix in a separate small PR like you suggested (with NULL checks added to the clone API), and pytest cases for bilinear, pole maps, and the fallback. Benchmarks for all tested scenarios"easy" EPSG:4326 to EPSG:3857, memory=50, method=nearestSerial: 7.472 s
"easy" EPSG:4326 to EPSG:3857, memory=300, method=nearestSerial: 9.557 s
"laea" EPSG:4326 to EPSG:3035, memory=50, method=nearestSerial: 10.696 s
"laea" EPSG:4326 to EPSG:3035, memory=300, method=nearestSerial: 10.729 s
"pole" EPSG:4326 to EPSG:3413, memory=50, method=nearestSerial: 29.219 s
"pole" EPSG:4326 to EPSG:3413, memory=300, method=nearestSerial: 29.106 s
"easy" EPSG:4326 to EPSG:3857, memory=50, method=bilinearSerial: 7.775 s
"easy" EPSG:4326 to EPSG:3857, memory=300, method=bilinearSerial: 10.058 s
"laea" EPSG:4326 to EPSG:3035, memory=50, method=bilinearSerial: 11.112 s
"laea" EPSG:4326 to EPSG:3035, memory=300, method=bilinearSerial: 11.057 s
"pole" EPSG:4326 to EPSG:3413, memory=50, method=bilinearSerial: 29.498 s
"pole" EPSG:4326 to EPSG:3413, memory=300, method=bilinearSerial: 29.635 s
"easy" EPSG:4326 to EPSG:3857, memory=50, method=lanczosSerial: 14.257 s
"laea" EPSG:4326 to EPSG:3035, memory=50, method=lanczosSerial: 16.393 s
Which memory caps force the serial tile cache fallbackThe module falls back to the serial tile cache only when the memory cap cannot hold one output row's input strip. This table shows which inputs hit that at caps of 1, 2, 5, and 10 MB.
Only the purposely wide input falls back, and only at the two smallest caps. Every normal map stays parallel even at 1 MB. The benchmarks show that speedup improves with the memory cap because a bigger cap means taller bands and fewer repeated input reads. At memory=300 the input read phase almost disappears on the pole map. The pole map scales best overall since its runtime is nearly all compute, which parallelizes well. The laea (hard) map at memory=50 is the weakest case because its sizing phase takes about a second of serial work before any threads start. Compute alone scales 4.8x to 5.7x at 8 threads on every map. The totals land lower than that because the sizing phase and the output write stay serial, which is why I’m trying numbers 1 and 3 from your review to hopefully fix this issue. I also tested bit identity against serial with 33 comparisons, and all of them came out exact. That covers all seven methods on a 100M-cell pole pair at two thread counts, all seven methods again through the forced fallback, and nearest, bilinear, and lanczos on the easy and laea pairs at two thread counts. I wanted to mention a behavioral difference from parallel r.projj and serial p.proj as well. For integer inputs below 2^24 the parallel output is exactly the same as serial, which the 33 comparisons above verify. Above 2^24 with the nearest method the two can differ by 1, and the parallel value is the correct one. This is because the old serial stores inputs as float32, which cannot hold integers that large exactly. I kept the full precision for now since it made sense to have the accurate numbers. If you would rather have byte for byte serial reproduction instead, please let me know and I will add it. I also ran my own AI review over all the code, the benchmarks + evidence with Claude. I addressed most of them and the remaining ones I will address early this week. One last thing, the branch was rebased onto current main to pick up #7730, so history was force updated. Thank you. |
github-actions
Bot
added
Python
|
I created the pytest and removed the testsuite folder gunittest. The pytest is 7 method reference tests plus 4 parallel-correctness tests on a 50x50 generated data. While testing the pole edge case in the pytest I came across a band sizing bug, fixed it, and verfied the output is correct. The NULL checks are also added to the clone API, and the globals fix is up separately as #7764. Please let know if there's anything I should fix in the pytest, thanks. |
|
Just a quick question, did you find a bug in the current implementation or are you referring to a bug in your prior parallel implementation? If it's a bug in non-parallel implementation, I would recommend creating a separate PR for it, so that the fix is not hidden here. |
|
Also the pytest would be better as separate PR and this PR would be adding the parallel tests. This will keep this PR smaller and would catch any regressions. |
|
It was a bug in my parallel implementation's band sizing, but I was able to fix it in this PR. I'll move the method reference tests and replacing the testsuite into a separate PR, and keep the tests that check the threaded output matches the serial output here. |
|
I was able to implement numbers 1 and 3 from the review with Claude. Number 1 keeps the input rows resident across bands instead of re-reading them, and number 3 writes the previous band's output while the next band computes. Both are bitwise identical to serial across all 33 comparison cases. I also added the nprocs option to match the other parallel modules, and a thread scaling benchmark script under raster/r.proj/benchmark. The new run of benchmarks are below. Serial was re-measured in the same session, so the speedup columns are the comparison rather than the raw seconds. At more than one thread, the output write is happening while compute is happening, so the write column only shows the final level. The benchmark script runs for the graphs are still going, I'll post them when they finish. But I think the module is ready for review. If anyone gets a chance to run the benchmark script on their machines I would like to see how it scales with more cores. I'll start looking at r.geomorphon while review is happening here and address any feedback or issues as they come. The method tests are in #7766 and the globals fix is in #7764. Thank you. Benchmarks for all tested scenarios"easy" EPSG:4326 to EPSG:3857, memory=50, method=nearestSerial: 9.991s (measured same session)
"easy" EPSG:4326 to EPSG:3857, memory=300, method=nearestSerial: 9.832s (measured same session)
"easy" EPSG:4326 to EPSG:3857, memory=50, method=bilinearSerial: 10.633s (measured same session)
"easy" EPSG:4326 to EPSG:3857, memory=300, method=bilinearSerial: 9.962s (measured same session)
"laea" EPSG:4326 to EPSG:3035, memory=50, method=nearestSerial: 13.466s (measured same session)
"laea" EPSG:4326 to EPSG:3035, memory=300, method=nearestSerial: 10.792s (measured same session)
"laea" EPSG:4326 to EPSG:3035, memory=50, method=bilinearSerial: 13.648s (measured same session)
"laea" EPSG:4326 to EPSG:3035, memory=300, method=bilinearSerial: 11.146s (measured same session)
"pole" EPSG:4326 to EPSG:3413, memory=50, method=nearestSerial: 29.684s (measured same session)
"pole" EPSG:4326 to EPSG:3413, memory=300, method=nearestSerial: 29.276s (measured same session)
"pole" EPSG:4326 to EPSG:3413, memory=50, method=bilinearSerial: 29.659s (measured same session)
"pole" EPSG:4326 to EPSG:3413, memory=300, method=bilinearSerial: 29.169s (measured same session)
"easy" EPSG:4326 to EPSG:3857, memory=50, method=lanczosSerial: 16.118s (measured same session)
"laea" EPSG:4326 to EPSG:3035, memory=50, method=lanczosSerial: 18.258s (measured same session)
Fallback sweep: memory-cap bail boundary (out_mult=2)
|
|
These are my graphs. Scaling seems to be strongest near the default memory setting. I also noticed with a higher memory value the bands get taller and fewer, so there is less overlap between writing one band and computing the next. This makes the final write at the end a larger serial share of the run, which is exactly how I expected it to behave given the write has to stay sequential. 
|
…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).
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.
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.
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.
The tile sizing search finds, for each output band, the tallest band height and its column tiling whose input strip fits the memory cap. Neighboring bands almost always end up with the same size, since the projection changes gradually from one band to the next. The search now tries the previous band's accepted height and width first instead of restarting the descending scan from the top every time. The previous band's size is checked with the same exact edge walk acceptance test the full search uses, and the next taller height is checked to make sure it does not fit. Together these two checks confirm the reused size is the tallest fitting answer, the same result the full scan would have returned. If either check fails, the code falls back to the full descending search, so the worst case costs the same as before. Because acceptance is decided by the same test in both paths, the resulting bands and tiles are identical to before and the output is bit for bit unchanged. In the common case a band is sized in two edge walks instead of a full descending scan.
The memory-bounded band sizing walks each output tile's perimeter to bound the range of input rows the tile needs, then loads that strip. A tile whose interior contains a geographic pole has its northmost or southmost latitude at the pole, in the tile interior, where the perimeter walk never samples it. The strip was therefore sized too small, and projecting a pole-containing map aborted with a "Band strip under-sized" error at every memory setting, though the projection itself was well defined. This computes, once per map, each geographic pole that lies within the input's latitude coverage: its coordinate in the output projection and its input row. When an output tile's rectangle contains a pole, that pole's input row is folded into the tile's row span, so the height and width search sees the true footprint and shrinks pole tiles until they fit the memory cap. The loaded strip then covers every row the fill reads. Only lat/lon input is handled, where a pole is at latitude 90 or -90. If the pole's coordinate transform fails or returns a non-finite value the pole is skipped and the existing under-size guard stays as the backstop. A map with no pole in the output frame is unaffected: the row spans, the band and tile partition, and the output are byte for byte unchanged.
The banded compute path dispatched a nearest-only strip reader for every resampling method, so bilinear, bicubic, lanczos and their fallback variants silently produced nearest-neighbor output. Add interp_strip.c with strip counterparts of the cache kernels (strip_bilinear, strip_cubic, strip_lanczos, and the three _f fallbacks). They read the in-RAM full-width FCELL band strip the banded path already loads, using the same base index, bounds, weights, and null fallback as the readcell-cache kernels in bilinear.c, cubic.c, and lanczos.c. A strip_kernels[] table, ordered like menu[], resolves each method to its strip counterpart once after option parsing; nearest keeps the existing interpolate_strip reader in slot 0. Output is bitwise identical to serial r.proj for all seven methods across the test datasets.
Serial r.proj computes each output row's northing by subtracting ns_res row by row. The banded path computed it directly as north - ns_res/2 - row * ns_res, which can differ by one ulp when ns_res is not exactly representable. Nearest is unaffected, but for the other methods the shifted interpolation weights changed a few cells (29 of 76M on the EPSG:3035 test). Precompute the row northings once with the serial recurrence and use that array in both the sizing walk and the compute loop.
The banded path aborted when the memory cap could not hold even one output row's input strip. In practice that needs an input wider than about cap/20 columns; poles and oblique projections do not trigger it. Warn with the minimum memory that keeps the parallel path, then finish the run through the old serial readcell cache. Output is bitwise identical to serial r.proj. Failed transforms set the cell NULL like the banded path does. R_PROJ_FORCE_TILECACHE forces the fallback for testing.
GPJ_clone_transform did not check the results of proj_context_create() and proj_clone(). Either can return NULL on failure, and the NULL would otherwise surface later as a crash deep inside PROJ when the cloned transform is first used. Both are now checked and fail with G_fatal_error naming the call. r.proj calls this once per worker thread, so a failure terminates the process from inside the parallel region; that is intended, since a clone failure leaves the thread with no usable transform.
…input A polar-stereographic output frame centered on a pole, reprojecting an input truncated below that pole (e.g. input reaching 89 degrees, not 90), aborted with "Band strip under-sized" (or, before that guard existed, silently read garbage). The banded strip sizing sampled only the output tile's perimeter, so the frame-center-proximal interior cell that reaches the input's northernmost edge row was never seen, and the strip loaded too few input rows. The pole footprint fold already handled a pole lying inside the input map; it now also folds in the input's edge row (0 or rows-1) when a pole outside the input's latitude coverage still projects into the frame. Gated by the existing point-in-rect test, so bands that do not image a pole compute byte-identical spans; verified unchanged on non-pole frames.
Adds pytest tests that verify r.proj's banded parallel output matches its serial output on generated CI-sized data. r.proj has no nprocs option, so each run sets OMP_NUM_THREADS on a per-call environment copy to select the serial (1) or parallel (N) path without mutating shared state. Four tests: bilinear identity (with a nearest-vs-bilinear dispatch-liveness guard so a silent fallback to nearest cannot pass the check vacuously), nearest identity under a constrained memory cap that forces band sizing, nearest identity into a pole-centered frame, and a forced tile-cache fallback (R_PROJ_FORCE_TILECACHE) compared against the banded path to cover both code paths. Inputs are integer CELL below 2^24 so the FCELL readcell cache round-trips losslessly.
Ports the 7 interpolation-method tests from testsuite/test_rproj.py to a parametrized pytest (r_proj_methods_test.py) over generated CI-sized data, with references captured from the serial binary. Removes the gunittest file. Drops its 4 output-format tests (test_list_output_plain/json, test_print_output_plain/json), which asserted NC-SPM-specific -l/-p output; r.proj's -l and -p flags are left without test coverage as a result.
Overlapping input rows are kept between single-tile bands and only new rows are read, instead of re-reading each band's full span. Cuts the input read phase about 70 percent on the wide LAEA benchmark; output is bitwise identical to serial.
The band output buffer is now double-buffered when running with more than one thread: one thread writes the previous band's rows in order while the rest compute the current band. Single-thread runs keep the sequential write. Hides most of the output write time at higher thread counts; output is bitwise identical to serial.
The method reference tests and the gunittest-to-pytest testsuite replacement move to a dedicated PR (branch fix-rproj-tests) so this PR stays focused on the OpenMP parallelization and that PR catches serial regressions on its own. This drops raster/r.proj/tests/r_proj_methods_test.py and restores raster/r.proj/testsuite/test_rproj.py to its main state; the split PR carries the deletion, so this restore is temporary and goes away once that PR merges. The parallel-correctness tests (r_proj_parallel_test.py) and their conftest.py stay here; conftest.py gains a note that it is duplicated on the split branch.
Add the standard G_OPT_M_NPROCS option so the compute thread count can be set with nprocs= instead of only OMP_NUM_THREADS. A value above zero overrides OMP_NUM_THREADS and zero keeps the OpenMP default. The option is read once through compute_nprocs() before the band fit search, so it drives the compute region, the per-thread read fds, and the output double-buffer together. The parallel-correctness tests now pass nprocs= instead of setting OMP_NUM_THREADS. Also shorten the main.c comments to flowing prose, dropping restated design narration and internal shorthand while keeping the load-bearing rationale.
Add a benchmark that sweeps the nprocs= thread count from 1 to 8 at two memory caps and plots the time, speedup, and efficiency metrics, following the r.param.scale benchmark template with grass.benchmark. It builds a source project and reprojects a generated raster from EPSG:4326 into EPSG:3857 in a temporary database, so it is self-contained.
a41dd0f to
003e7ab
Compare
5 participants
This PR parallelizes r.proj with OpenMP while keeping memory use within the amount the user allows.
How it works
The output is split into horizontal bands. For each band, the code projects the band's edges back into the input to see which input rows it touches, then picks a band height where that input strip fits under the memory cap. For easy maps like lat/lon to Web Mercator (EPSG:4326 to EPSG:3857), this is all that's needed, since each output band only touches a small thin part of the input.
Hard maps are pairs where the output rows don't line up with the input rows. For example, lat/lon to Lambert Azimuthal Equal Area over a wide European extent (EPSG:4326 to EPSG:3035). There a single full width output row can need more input rows than the cap allows. When that happens the band is split into column tiles until each tile's input fits. Pole maps like lat/lon to Polar Stereographic North (EPSG:4326 to EPSG:3413) with the pole inside the frame are the extreme case, since every longitude converges at the pole. A tile that has the pole gets its input span extended to cover the pole row.
If no tile width can bring the input under the cap, the module still won't fail. It warns the user with the exact memory value that they would need to allow to keep the parallel path, and still finishes the run through the old serial tile cache. This is so every job the old r.proj could handle can still complete here. In my benchmark sweeps this fallback never activated on any real projection pair down to a 1 MB cap, only on an input that is purposely built wide enough to trigger it. So this case is very rare and most maps would fall under the easy, hard or pole category mentioned above.
Once the bands have their rows, the work is run in parallel. Input rows are read with per thread file descriptors, and each thread computes with its own PROJ context through the new GPJ_clone_transform API in lib/proj. While a band computes, the previous band's output is written at the same time. The write itself has to stay serial because compressed row offsets depend on all prior rows, but it doesn't block the compute threads while it runs. Consecutive bands also need many of the same input rows, so the input strip is kept in memory and slides forward from band to band. It only reads the new rows instead of re-reading the whole strip from disk.
Correctness
Output is bitwise identical to serial r.proj across 33 comparisons. It's all seven interpolation methods on a 100M cell polar pair (4326 to 3413) at two thread counts, all seven methods again through the forced fallback path, and nearest, bilinear, and lanczos on the easy (4326 to 3857) and hard (4326 to 3035) pairs at two thread counts. Null cells match exactly. Parallel correctness pytest tests are included in this PR. The method reference tests and testsuite migration are in #7766, and a fix for a data race on globals in lib/proj is in #7764.
One thing I want to mention regarding behavior. For integer inputs above 2^24 with "nearest" method the parallel output can differ from serial by 1 and the parallel value is the correct one. For example an input cell of 16777217 comes out of old r.proj as 16777216 but stays 16777217 here, because the old code caches inputs as float32 which cannot hold integers that large exactly. Below 2^24 the two are identical.
Performance (The machine I used and tested all this on is an 8 core Apple M3)
At 8 threads I got 4.4x to 4.9x on the easy pair, 4.2x to 4.9x on the polar pair, and 2.1x to 3.5x on the hard pair depending on the memory cap. Single thread runs equally with serial. Peak RSS on the easy pair drops from 763 MB with the whole map in RAM to 130 MB at memory=50. Speedup ratios also grow as the memory cap increases because taller bands mean fewer repeated input reads and less serial sizing work. So, users can trade memory for speed with the existing option. A benchmark script to show thread scaling is a part of this PR under raster/r.proj/benchmark.
Current tables and graphs are in the most recent comments all the way down below.