r.proj: memory-bounded parallel banding #7627
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -66,6 +66,10 @@ | |
| #include <grass/gjson.h> | ||
| #include "r.proj.h" | ||
|
|
||
| #define PACKAGE "grassmods" | ||
|
|
||
| #include <omp.h> | ||
|
|
||
| /* modify this table to add new methods */ | ||
| struct menu menu[] = { | ||
| {p_nearest, "nearest", "nearest neighbor"}, | ||
|
|
@@ -80,6 +84,106 @@ struct menu menu[] = { | |
| static char *make_ipol_list(void); | ||
| static char *make_ipol_desc(void); | ||
|
|
||
| /* Nearest read from an in-RAM input STRIP holding input rows [imin, imax]. | ||
| * col_idx/row_idx are full-map input indices; the strip is addressed relative | ||
| * to imin. A sample inside the full input map but outside the loaded strip | ||
| * means the band footprint was under-sized: this is the stop-on-divergence | ||
| * trip (must never fire if band_input_row_span is correct). Lock-free: reads | ||
| * only, disjoint output slots per thread. */ | ||
| static void interpolate_strip(void *strip, void *obufptr, int cell_type, | ||
| double col_idx, double row_idx, | ||
| struct Cell_head *incellhd, int imin, int imax) | ||
| { | ||
| int c = (int)floor(col_idx); | ||
| int r = (int)floor(row_idx); | ||
| int cell_size = Rast_cell_size(cell_type); | ||
|
|
||
| /* Outside the full input map: legitimate NULL (same as p_nearest). */ | ||
| if (r < 0 || r >= incellhd->rows || c < 0 || c >= incellhd->cols) { | ||
| Rast_set_null_value(obufptr, 1, cell_type); | ||
| return; | ||
| } | ||
|
|
||
| /* Inside the map but outside the loaded strip: the span check under-sized | ||
| * the strip. This is a correctness failure, not a NULL. */ | ||
|
Member
There was a problem hiding this comment. I'm not entirely sure what you mean here.
Contributor
Author
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. |
||
| if (r < imin || r > imax) | ||
| G_fatal_error(_("Band strip under-sized: input row %d outside loaded " | ||
| "range [%d, %d] at column %d"), | ||
| r, imin, imax, c); | ||
|
|
||
| unsigned char *src = | ||
| (unsigned char *)strip + | ||
| (((size_t)(r - imin) * incellhd->cols + c) * cell_size); | ||
| memcpy(obufptr, src, cell_size); | ||
| } | ||
|
|
||
| /* Dense edge-walk of an output band's rectangle [obr0, obr1) projected into | ||
| * input space; returns the min/max INPUT ROW touched, plus a 2-cell margin, | ||
| * clamped to the input map. Samples the band's top and bottom rows across all | ||
| * columns and its left and right columns across all band rows (bordwalk-style), | ||
| * so a curved transform's interior-edge extremum is caught -- corner-only | ||
| * sampling can under-size the strip. Called serially, before the parallel | ||
| * region, so the shared tproj is safe here. Returns imax < imin for a band | ||
| * that projects entirely outside the input. */ | ||
| static void band_input_row_span(const struct Cell_head *ohd, | ||
| const struct Cell_head *ihd, | ||
| const struct pj_info *oproj, | ||
| const struct pj_info *iproj, | ||
| const struct pj_info *tproj, int obr0, int obr1, | ||
| int *imin, int *imax) | ||
| { | ||
| double rmin = 1e300, rmax = -1e300; | ||
| int e, r, c; | ||
|
|
||
| /* top edge (row obr0) and bottom edge (row obr1-1), all columns */ | ||
| for (e = 0; e < 2; e++) { | ||
| int orow = (e == 0) ? obr0 : (obr1 - 1); | ||
| double y = ohd->north - (orow + 0.5) * ohd->ns_res; | ||
| for (c = 0; c < ohd->cols; c++) { | ||
| double x = ohd->west + (c + 0.5) * ohd->ew_res; | ||
| double xx = x, yy = y; | ||
| if (GPJ_transform(oproj, iproj, tproj, PJ_FWD, &xx, &yy, NULL) < 0) | ||
| continue; | ||
| double ri = (ihd->north - yy) / ihd->ns_res; | ||
| if (ri < rmin) | ||
| rmin = ri; | ||
| if (ri > rmax) | ||
| rmax = ri; | ||
| } | ||
| } | ||
| /* left edge (col 0) and right edge (col cols-1), all band rows */ | ||
| for (e = 0; e < 2; e++) { | ||
| int ocol = (e == 0) ? 0 : (ohd->cols - 1); | ||
| double x = ohd->west + (ocol + 0.5) * ohd->ew_res; | ||
| for (r = obr0; r < obr1; r++) { | ||
| double y = ohd->north - (r + 0.5) * ohd->ns_res; | ||
| double xx = x, yy = y; | ||
| if (GPJ_transform(oproj, iproj, tproj, PJ_FWD, &xx, &yy, NULL) < 0) | ||
| continue; | ||
| double ri = (ihd->north - yy) / ihd->ns_res; | ||
| if (ri < rmin) | ||
| rmin = ri; | ||
| if (ri > rmax) | ||
| rmax = ri; | ||
| } | ||
| } | ||
|
|
||
| if (rmax < rmin) { /* band projects entirely outside the input */ | ||
| *imin = 0; | ||
| *imax = -1; | ||
| return; | ||
| } | ||
|
|
||
| int lo = (int)floor(rmin) - 2; /* 2-cell margin for interp stencils */ | ||
| int hi = (int)floor(rmax) + 2; | ||
| if (lo < 0) | ||
| lo = 0; | ||
| if (hi > ihd->rows - 1) | ||
| hi = ihd->rows - 1; | ||
| *imin = lo; | ||
| *imax = hi; | ||
| } | ||
|
|
||
| int main(int argc, char **argv) | ||
| { | ||
| char *mapname, /* ptr to name of output layer */ | ||
|
|
@@ -98,14 +202,7 @@ int main(int argc, char **argv) | |
| overwrite, /* Overwrite */ | ||
| curr_proj; /* output projection (see gis.h) */ | ||
|
|
||
| double onorth, osouth, /* save original border coords */ | ||
| oeast, owest, inorth, isouth, ieast, iwest; | ||
| char north_str[30], south_str[30], east_str[30], west_str[30]; | ||
|
|
||
|
|
@@ -276,7 +373,6 @@ int main(int argc, char **argv) | |
| if (!ipolname) | ||
| G_fatal_error(_("<%s=%s> unknown %s"), interpol->key, interpol->answer, | ||
| interpol->key); | ||
|
|
||
| mapname = outmap->answer ? outmap->answer : inmap->answer; | ||
| if (mapname && !list->answer && !overwrite && !print_bounds->answer && | ||
|
|
@@ -664,19 +760,23 @@ int main(int argc, char **argv) | |
| G_message(_("NS-res: %f"), outcellhd.ns_res); | ||
| G_message(" "); | ||
|
|
||
| /* Open the input map (input location env). Banding loads only per-band | ||
| * input strips, not the whole map, so fdi stays open across the band loop. | ||
| */ | ||
| G_switch_env(); | ||
| Rast_set_input_window(&incellhd); | ||
| fdi = Rast_open_old(inmap->answer, setname); | ||
| cell_type = Rast_get_map_type(fdi); | ||
| if (strcmp(interpol->answer, "nearest") != 0) | ||
| cell_type = FCELL_TYPE; | ||
| cell_size = Rast_cell_size(cell_type); | ||
|
|
||
| /* Back to the output location: set output window, init transform, open | ||
| * output map. Both fds now stay open; rd_window/wr_window are set and | ||
| * survive env switches, so reads/writes use the right windows throughout. | ||
| */ | ||
| G_switch_env(); | ||
| Rast_set_output_window(&outcellhd); | ||
| G_unset_window(); | ||
| G_set_window(&outcellhd); | ||
| tproj.def = NULL; | ||
|
|
@@ -687,73 +787,152 @@ int main(int argc, char **argv) | |
| if (GPJ_init_transform(&oproj, &iproj, &tproj) < 0) | ||
| G_fatal_error(_("Unable to initialize coordinate transformation")); | ||
|
|
||
| if (strcmp(interpol->answer, "nearest") == 0) | ||
| fdo = Rast_open_new(mapname, cell_type); | ||
| else | ||
| fdo = Rast_open_fp_new(mapname); | ||
|
|
||
| /* Banding (r.neighbors two-level structure): outer serial band loop -> | ||
| * serial strip load -> parallel compute into a per-band buffer -> serial | ||
| * in-order per-band write -> next band. Bounds peak memory by the cap | ||
| * instead of the whole input map (Path A). */ | ||
| double cap_mb = atof(memory->answer); | ||
| size_t cap_bytes = (size_t)(cap_mb * 1024.0 * 1024.0); | ||
| double t_size = 0.0, t_fill = 0.0, t_compute = 0.0, t_write = 0.0; | ||
| int n_bands = 0; | ||
|
|
||
| G_important_message(_("Projecting (banded, per-thread PROJ context)...")); | ||
|
|
||
| int obr0 = 0; | ||
| while (obr0 < outcellhd.rows) { | ||
| /* Band height: start from all remaining rows, halve until the input | ||
| * strip plus the band's output buffer fit the cap. band_input_row_span | ||
| * is RE-RUN for every candidate height -- the previous band's span is | ||
| * never reused. */ | ||
| double ts = omp_get_wtime(); | ||
| int band_orows = outcellhd.rows - obr0; | ||
| int imin = 0, imax = -1; | ||
| for (;;) { | ||
| band_input_row_span(&outcellhd, &incellhd, &oproj, &iproj, &tproj, | ||
| obr0, obr0 + band_orows, &imin, &imax); | ||
| int strip_rows = imax - imin + 1; | ||
| size_t strip_bytes = | ||
| strip_rows > 0 ? (size_t)strip_rows * incellhd.cols * cell_size | ||
| : 0; | ||
| size_t out_bytes = (size_t)band_orows * outcellhd.cols * cell_size; | ||
| if (strip_bytes + out_bytes <= cap_bytes) | ||
| break; | ||
| if (band_orows == 1) | ||
| G_fatal_error( | ||
| _("A single output row needs %.1f MB (input footprint %d " | ||
| "rows), exceeding the memory cap (%.1f MB). This " | ||
| "large-halo/oblique case needs the tile-cache path, " | ||
| "which is not implemented."), | ||
| (double)(strip_bytes + out_bytes) / (1024.0 * 1024.0), | ||
| strip_rows, cap_mb); | ||
| band_orows = (band_orows + 1) / 2; /* halve (round up), re-sample */ | ||
| } | ||
| t_size += omp_get_wtime() - ts; | ||
|
|
||
| int obr1 = obr0 + band_orows; | ||
| int strip_rows = imax - imin + 1; | ||
| n_bands++; | ||
|
|
||
| /* Serial strip load (single fd -> get_row not thread-safe). Reads in | ||
| * the INPUT env (matching the serial code's invariant), then back to | ||
| * OUTPUT for compute+write. EMPTY BAND: strip_rows <= 0 means the band | ||
| * projects entirely outside the input -> no malloc, no read; its cells | ||
| * become NULL via interpolate_strip's out-of-map path. */ | ||
| void *strip = NULL; | ||
| if (strip_rows > 0) { | ||
| strip = G_malloc((size_t)strip_rows * incellhd.cols * cell_size); | ||
| double t0 = omp_get_wtime(); | ||
| G_switch_env(); /* -> input */ | ||
| for (int r = imin; r <= imax; r++) | ||
| Rast_get_row(fdi, | ||
|
Member
There was a problem hiding this comment. Please see if it's possible to parallelize reading.
Contributor
Author
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.
Contributor
Author
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. |
||
| (unsigned char *)strip + | ||
| (size_t)(r - imin) * incellhd.cols * cell_size, | ||
| r, cell_type); | ||
| G_switch_env(); /* -> output */ | ||
| t_fill += omp_get_wtime() - t0; | ||
| } | ||
|
|
||
| /* Per-band output buffer, lock-free disjoint row slots (band-relative | ||
| * index), mirroring r.neighbors' outputs[i].buf. */ | ||
| void *band_out = | ||
| G_malloc((size_t)band_orows * outcellhd.cols * cell_size); | ||
|
|
||
| double t1 = omp_get_wtime(); | ||
| #pragma omp parallel | ||
|
Member
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
Contributor
Author
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. |
||
| { | ||
| /* Per-thread PROJ context + private transform clone (KEEP: this is | ||
| * the bit-exact-verified, banding-agnostic part). oproj/iproj are | ||
| * 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(); | ||
|
Contributor
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).
Contributor
Author
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. |
||
| tproj_local.pj = proj_clone(thread_ctx, tproj.pj); | ||
|
|
||
| #pragma omp for private(row, col) schedule(dynamic) | ||
| for (row = obr0; row < obr1; row++) { | ||
| void *out_row = | ||
| (unsigned char *)band_out + | ||
| (size_t)(row - obr0) * outcellhd.cols * cell_size; | ||
| double local_y = outcellhd.north - (outcellhd.ns_res / 2) - | ||
| (row * outcellhd.ns_res); | ||
| double local_x_start = outcellhd.west + (outcellhd.ew_res / 2); | ||
|
|
||
| for (col = 0; col < outcellhd.cols; col++) { | ||
| void *obufptr = | ||
| (unsigned char *)out_row + (size_t)col * cell_size; | ||
| double x1 = local_x_start + (col * outcellhd.ew_res); | ||
| double y1 = local_y; | ||
|
|
||
| if (GPJ_transform(&oproj, &iproj, &tproj_local, PJ_FWD, &x1, | ||
| &y1, NULL) < 0) { | ||
| Rast_set_null_value(obufptr, 1, cell_type); | ||
| } | ||
| else { | ||
| double c_idx = (x1 - incellhd.west) / incellhd.ew_res; | ||
| double r_idx = (incellhd.north - y1) / incellhd.ns_res; | ||
| interpolate_strip(strip, obufptr, cell_type, c_idx, | ||
| r_idx, &incellhd, imin, imax); | ||
| } | ||
| } | ||
| } | ||
|
|
||
| proj_destroy(tproj_local.pj); | ||
| proj_context_destroy(thread_ctx); | ||
| } | ||
| t_compute += omp_get_wtime() - t1; | ||
|
|
||
| /* Serial in-order write of the band's rows (Rast_put_row sequential). | ||
| */ | ||
| double t2 = omp_get_wtime(); | ||
| for (row = obr0; row < obr1; row++) | ||
| Rast_put_row(fdo, | ||
| (unsigned char *)band_out + | ||
| (size_t)(row - obr0) * outcellhd.cols * cell_size, | ||
| cell_type); | ||
| t_write += omp_get_wtime() - t2; | ||
|
|
||
| G_percent(obr1, outcellhd.rows, 5); | ||
|
|
||
| if (strip) | ||
| G_free(strip); | ||
| G_free(band_out); | ||
| obr0 = obr1; | ||
| } | ||
|
|
||
| G_message("PHASE_TIMERS size=%.4f fill=%.4f compute=%.4f write=%.4f " | ||
| "bands=%d", | ||
| t_size, t_fill, t_compute, t_write, n_bands); | ||
|
|
||
| /* Close input map in its own env, then the output map. */ | ||
| G_switch_env(); /* -> input */ | ||
| Rast_close(fdi); | ||
| G_switch_env(); /* -> output */ | ||
| Rast_close(fdo); | ||
|
|
||
| if (have_colors > 0) { | ||
| Rast_write_colors(mapname, G_mapset(), &colr); | ||
|
|
||
Uh oh!
There was an error while loading. Please reload this page.