From da367372c6270e2f06d394b8af64531488490746 Mon Sep 17 00:00:00 2001 From: Kaushik Raja Date: Sun, 15 Mar 2026 23:21:16 -0500 Subject: [PATCH 01/11] POC: parallelize r.proj via RAM-resident buffer, 2.5x speedup on 8-core Apple M-series --- raster/r.proj/main.c | 119 +++++++++++++++++++++++++++++-------------- 1 file changed, 80 insertions(+), 39 deletions(-) diff --git a/raster/r.proj/main.c b/raster/r.proj/main.c index 59702cdfc92..2d57c1b52d0 100644 --- a/raster/r.proj/main.c +++ b/raster/r.proj/main.c @@ -66,6 +66,10 @@ #include #include "r.proj.h" +#define PACKAGE "grassmods" + +#include + /* modify this table to add new methods */ struct menu menu[] = { {p_nearest, "nearest", "nearest neighbor"}, @@ -80,6 +84,28 @@ struct menu menu[] = { static char *make_ipol_list(void); static char *make_ipol_desc(void); + +/* Custom Interpolation for RAM Bypass - Lock-Free */ +void interpolate_ram(void *full_map, void *obufptr, int cell_type, + double col_idx, double row_idx, struct Cell_head *incellhd) +{ + int c = (int)floor(col_idx); + int r = (int)floor(row_idx); + int cell_size = Rast_cell_size(cell_type); + + /* Boundary check */ + if (r < 0 || r >= incellhd->rows || c < 0 || c >= incellhd->cols) { + Rast_set_null_value(obufptr, 1, cell_type); + return; + } + + /* Direct memory access - Thread Safe for Reads */ + unsigned char *src = (unsigned char *)full_map + + (((size_t)r * incellhd->cols + c) * cell_size); + memcpy(obufptr, src, cell_size); +} + + int main(int argc, char **argv) { char *mapname, /* ptr to name of output layer */ @@ -671,7 +697,6 @@ int main(int argc, char **argv) cell_type = Rast_get_map_type(fdi); ibuffer = readcell(fdi, memory->answer); Rast_close(fdi); - /* And switch back to original location */ G_switch_env(); Rast_set_output_window(&outcellhd); @@ -702,58 +727,73 @@ int main(int argc, char **argv) xcoord2 = outcellhd.west + (outcellhd.ew_res / 2); ycoord2 = outcellhd.north - (outcellhd.ns_res / 2); - G_important_message(_("Projecting...")); - for (row = 0; row < outcellhd.rows; row++) { - /* obufptr = obuffer */; - G_percent(row, outcellhd.rows - 1, 2); - -#if 0 - /* parallelization does not always work, - * segfaults in the interpolation functions - * can happen */ -#pragma omp parallel for schedule(static) -#endif + /* --- RAM BYPASS START --- */ + /* 1. Allocate a flat buffer for the entire input map in RAM */ + size_t total_cells = (size_t)incellhd.rows * incellhd.cols; - for (col = 0; col < outcellhd.cols; col++) { - void *obufptr = - (void *)((const unsigned char *)obuffer + col * cell_size); + /* Simple memory safety check */ + double memory_mb = (double)(total_cells * cell_size) / (1024.0 * 1024.0); + G_debug(1, "Bypass buffer requires %.2f MB of RAM", memory_mb); - double xcoord1 = xcoord2 + (col)*outcellhd.ew_res; - double ycoord1 = ycoord2; + if (memory_mb > 4000) { + G_warning(_("Input map requires %.2f MB of RAM for parallel processing. " + "If this causes a crash, reduce the region size with g.region."), + memory_mb); + } + /* Proceed with allocation */ + void *full_map_array = G_malloc(total_cells * cell_size); + + if (!full_map_array) + G_fatal_error("Insufficient RAM for bypass. Try a smaller region."); + + G_important_message(_("Loading map into RAM buffer for parallel bypass...")); + + /* 2. Sequential load (Single-threaded, library-safe) */ + /* ibuffer is already loaded by readcell, but we move it to a flat array + so we can access it lock-free without the readcell tile cache logic */ + for (int r = 0; r < incellhd.rows; r++) { + void *row_ptr = (void *)((unsigned char *)full_map_array + ((size_t)r * incellhd.cols * cell_size)); + /* We use the already-loaded cache here, but single-threaded */ + interpolate(ibuffer, row_ptr, cell_type, 0.0, (double)r, &incellhd); + G_percent(r, incellhd.rows - 1, 5); + } - /* project coordinates in output matrix to */ - /* coordinates in input matrix */ - if (GPJ_transform(&oproj, &iproj, &tproj, PJ_FWD, &xcoord1, - &ycoord1, NULL) < 0) { - G_fatal_error(_("Error in %s"), "GPJ_transform()"); - Rast_set_null_value(obufptr, 1, cell_type); - } - else { - /* convert to row/column indices of input matrix */ + G_important_message(_("Projecting (Lock-Free Parallel)...")); + + #pragma omp parallel for private(row, col) schedule(dynamic) + for (row = 0; row < outcellhd.rows; row++) { + void *local_obuffer = Rast_allocate_output_buf(cell_type); + double local_y = outcellhd.north - (outcellhd.ns_res / 2) - (row * outcellhd.ns_res); + double local_x_start = outcellhd.west + (outcellhd.ew_res / 2); - /* column index in input matrix */ - double col_idx = (xcoord1 - incellhd.west) / incellhd.ew_res; + for (col = 0; col < outcellhd.cols; col++) { + void *obufptr = (void *)((unsigned char *)local_obuffer + (size_t)col * cell_size); + double x1 = local_x_start + (col * outcellhd.ew_res); + double y1 = local_y; - /* row index in input matrix */ - double row_idx = (incellhd.north - ycoord1) / incellhd.ns_res; + if (GPJ_transform(&oproj, &iproj, &tproj, 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; - /* and resample data point */ - interpolate(ibuffer, obufptr, cell_type, col_idx, row_idx, - &incellhd); + /* CALL OUR LOCK-FREE RAM INTERPOLATOR */ + interpolate_ram(full_map_array, obufptr, cell_type, c_idx, r_idx, &incellhd); } - - /* obufptr = G_incr_void_ptr(obufptr, cell_size); */ } - Rast_put_row(fdo, obuffer, cell_type); - - xcoord2 = outcellhd.west + (outcellhd.ew_res / 2); - ycoord2 -= outcellhd.ns_res; + #pragma omp critical + { + Rast_put_row(fdo, local_obuffer, cell_type); + } + G_free(local_obuffer); } + /*RAM BYPASS END*/ Rast_close(fdo); release_cache(ibuffer); + G_free(full_map_array); /* Clean up our bypass buffer */ if (have_colors > 0) { Rast_write_colors(mapname, G_mapset(), &colr); @@ -811,3 +851,4 @@ char *make_ipol_desc(void) return buf; } + From df9952297736c6dbbed099c2324d79fdf944f0bf Mon Sep 17 00:00:00 2001 From: Kaushik Raja Date: Mon, 16 Mar 2026 17:27:58 -0500 Subject: [PATCH 02/11] r.proj: use memory option to limit RAM buffer per community feedback --- raster/r.proj/main.c | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/raster/r.proj/main.c b/raster/r.proj/main.c index 2d57c1b52d0..40a3f8be2c2 100644 --- a/raster/r.proj/main.c +++ b/raster/r.proj/main.c @@ -736,10 +736,13 @@ int main(int argc, char **argv) double memory_mb = (double)(total_cells * cell_size) / (1024.0 * 1024.0); G_debug(1, "Bypass buffer requires %.2f MB of RAM", memory_mb); - if (memory_mb > 4000) { - G_warning(_("Input map requires %.2f MB of RAM for parallel processing. " - "If this causes a crash, reduce the region size with g.region."), - memory_mb); + double user_limit_mb = atof(memory->answer); + if (memory_mb > user_limit_mb) { + G_warning(_("Input map requires %.2f MB of RAM for parallel processing, " + "which exceeds the current memory limit (%.0f MB)."), + memory_mb, user_limit_mb); + G_important_message(_("The process may crash if system RAM is insufficient. " + "Increase the 'memory' option or use g.region to reduce the area.")); } /* Proceed with allocation */ void *full_map_array = G_malloc(total_cells * cell_size); From ffed3e2c47e2483d04cda6cfeb370859b5e4c5d7 Mon Sep 17 00:00:00 2001 From: Kaushik Raja Date: Wed, 1 Jul 2026 21:26:58 -0700 Subject: [PATCH 03/11] r.proj: replace full-map RAM buffer with memory-bounded banding 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). --- raster/r.proj/main.c | 331 ++++++++++++++++++++++++++++++------------- 1 file changed, 233 insertions(+), 98 deletions(-) diff --git a/raster/r.proj/main.c b/raster/r.proj/main.c index 40a3f8be2c2..7afc15e60e8 100644 --- a/raster/r.proj/main.c +++ b/raster/r.proj/main.c @@ -84,27 +84,105 @@ struct menu menu[] = { static char *make_ipol_list(void); static char *make_ipol_desc(void); - -/* Custom Interpolation for RAM Bypass - Lock-Free */ -void interpolate_ram(void *full_map, void *obufptr, int cell_type, - double col_idx, double row_idx, struct Cell_head *incellhd) +/* 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); - /* Boundary check */ + /* 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; } - /* Direct memory access - Thread Safe for Reads */ - unsigned char *src = (unsigned char *)full_map + - (((size_t)r * incellhd->cols + c) * cell_size); + /* Inside the map but outside the loaded strip: the span check under-sized + * the strip. This is a correctness failure, not a NULL. */ + 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) { @@ -124,14 +202,7 @@ int main(int argc, char **argv) overwrite, /* Overwrite */ curr_proj; /* output projection (see gis.h) */ - void *obuffer; /* buffer that holds one output row */ - - struct cache *ibuffer; /* buffer that holds the input map */ - func interpolate; /* interpolation routine */ - - double xcoord2, /* temporary x coordinates */ - ycoord2, /* temporary y coordinates */ - onorth, osouth, /* save original border coords */ + 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]; @@ -302,7 +373,6 @@ int main(int argc, char **argv) if (!ipolname) G_fatal_error(_("<%s=%s> unknown %s"), interpol->key, interpol->answer, interpol->key); - interpolate = menu[method].method; mapname = outmap->answer ? outmap->answer : inmap->answer; if (mapname && !list->answer && !overwrite && !print_bounds->answer && @@ -690,18 +760,23 @@ int main(int argc, char **argv) G_message(_("NS-res: %f"), outcellhd.ns_res); G_message(" "); - /* open and read the relevant parts of the input map and close it */ + /* 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); - ibuffer = readcell(fdi, memory->answer); - Rast_close(fdi); - /* And switch back to original location */ + 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); - - /* reproject from output to input */ G_unset_window(); G_set_window(&outcellhd); tproj.def = NULL; @@ -712,91 +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) { + if (strcmp(interpol->answer, "nearest") == 0) fdo = Rast_open_new(mapname, cell_type); - obuffer = (CELL *)Rast_allocate_output_buf(cell_type); - } - else { + else fdo = Rast_open_fp_new(mapname); - cell_type = FCELL_TYPE; - obuffer = (FCELL *)Rast_allocate_output_buf(cell_type); - } - - cell_size = Rast_cell_size(cell_type); - xcoord2 = outcellhd.west + (outcellhd.ew_res / 2); - ycoord2 = outcellhd.north - (outcellhd.ns_res / 2); - - - /* --- RAM BYPASS START --- */ - /* 1. Allocate a flat buffer for the entire input map in RAM */ - size_t total_cells = (size_t)incellhd.rows * incellhd.cols; - - /* Simple memory safety check */ - double memory_mb = (double)(total_cells * cell_size) / (1024.0 * 1024.0); - G_debug(1, "Bypass buffer requires %.2f MB of RAM", memory_mb); + /* 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, + (unsigned char *)strip + + (size_t)(r - imin) * incellhd.cols * cell_size, + r, cell_type); + G_switch_env(); /* -> output */ + t_fill += omp_get_wtime() - t0; + } - double user_limit_mb = atof(memory->answer); - if (memory_mb > user_limit_mb) { - G_warning(_("Input map requires %.2f MB of RAM for parallel processing, " - "which exceeds the current memory limit (%.0f MB)."), - memory_mb, user_limit_mb); - G_important_message(_("The process may crash if system RAM is insufficient. " - "Increase the 'memory' option or use g.region to reduce the area.")); - } - /* Proceed with allocation */ - void *full_map_array = G_malloc(total_cells * cell_size); - - if (!full_map_array) - G_fatal_error("Insufficient RAM for bypass. Try a smaller region."); - - G_important_message(_("Loading map into RAM buffer for parallel bypass...")); - - /* 2. Sequential load (Single-threaded, library-safe) */ - /* ibuffer is already loaded by readcell, but we move it to a flat array - so we can access it lock-free without the readcell tile cache logic */ - for (int r = 0; r < incellhd.rows; r++) { - void *row_ptr = (void *)((unsigned char *)full_map_array + ((size_t)r * incellhd.cols * cell_size)); - /* We use the already-loaded cache here, but single-threaded */ - interpolate(ibuffer, row_ptr, cell_type, 0.0, (double)r, &incellhd); - G_percent(r, incellhd.rows - 1, 5); - } + /* 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); - G_important_message(_("Projecting (Lock-Free Parallel)...")); - - #pragma omp parallel for private(row, col) schedule(dynamic) - for (row = 0; row < outcellhd.rows; row++) { - void *local_obuffer = Rast_allocate_output_buf(cell_type); - 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 = (void *)((unsigned char *)local_obuffer + (size_t)col * cell_size); - double x1 = local_x_start + (col * outcellhd.ew_res); - double y1 = local_y; - - if (GPJ_transform(&oproj, &iproj, &tproj, 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; - - /* CALL OUR LOCK-FREE RAM INTERPOLATOR */ - interpolate_ram(full_map_array, obufptr, cell_type, c_idx, r_idx, &incellhd); + double t1 = omp_get_wtime(); +#pragma omp parallel + { + /* 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(); + 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); + } + } } - } - #pragma omp critical - { - Rast_put_row(fdo, local_obuffer, cell_type); + proj_destroy(tproj_local.pj); + proj_context_destroy(thread_ctx); } - G_free(local_obuffer); + 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; } - /*RAM BYPASS END*/ + 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); - release_cache(ibuffer); - G_free(full_map_array); /* Clean up our bypass buffer */ if (have_colors > 0) { Rast_write_colors(mapname, G_mapset(), &colr); @@ -854,4 +990,3 @@ char *make_ipol_desc(void) return buf; } - From 0fb958e25729bde1700b03eac5ad7fde80295d6b Mon Sep 17 00:00:00 2001 From: Kaushik Raja Date: Thu, 2 Jul 2026 00:09:18 -0700 Subject: [PATCH 04/11] r.proj: link libproj for direct proj_* calls --- raster/r.proj/Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/raster/r.proj/Makefile b/raster/r.proj/Makefile index 083cd95c72f..147b47fe2d8 100644 --- a/raster/r.proj/Makefile +++ b/raster/r.proj/Makefile @@ -2,7 +2,7 @@ MODULE_TOPDIR = ../.. PGM = r.proj -LIBES = $(GPROJLIB) $(RASTERLIB) $(GISLIB) $(MATHLIB) $(PARSONLIB) +LIBES = $(GPROJLIB) $(RASTERLIB) $(GISLIB) $(MATHLIB) $(PARSONLIB) $(PROJLIB) DEPENDENCIES = $(GPROJDEP) $(RASTERDEP) $(GISDEP) EXTRA_LIBS = $(OPENMP_LIBPATH) $(OPENMP_LIB) From 0250f8c175d5f02a1185cbee1acda4ad8c37a944 Mon Sep 17 00:00:00 2001 From: Kaushik Raja Date: Thu, 2 Jul 2026 18:49:36 -0700 Subject: [PATCH 05/11] r.proj: address review comments (remove stray PACKAGE define, clarify comments) --- raster/r.proj/main.c | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/raster/r.proj/main.c b/raster/r.proj/main.c index 7afc15e60e8..d4d181a6936 100644 --- a/raster/r.proj/main.c +++ b/raster/r.proj/main.c @@ -66,8 +66,6 @@ #include #include "r.proj.h" -#define PACKAGE "grassmods" - #include /* modify this table to add new methods */ @@ -104,8 +102,12 @@ static void interpolate_strip(void *strip, void *obufptr, int 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. */ + /* This input row is inside the input map (the check above already handled + * coordinates that fall outside it), but it is not among the rows we + * preloaded into this band's strip. That cannot happen if the band's + * footprint estimate was right, so it means the estimate was wrong: a bug + * in band sizing, not a normal case. Fail loudly rather than write a NULL + * and silently produce wrong output. */ if (r < imin || r > imax) G_fatal_error(_("Band strip under-sized: input row %d outside loaded " "range [%d, %d] at column %d"), @@ -863,6 +865,14 @@ int main(int argc, char **argv) G_malloc((size_t)band_orows * outcellhd.cols * cell_size); double t1 = omp_get_wtime(); + /* Each band runs one parallel region. This is not nested parallelism: + * the "omp for" below does not create a second thread team, it only + * divides the band's output rows among the threads that this "omp + * parallel" created. The two directives are kept separate instead of a + * combined "omp parallel for" because every thread must clone its own + * PROJ context before the row loop starts and destroy it after the loop + * ends, and that per-thread setup has to sit inside the parallel region + * but outside the for. */ #pragma omp parallel { /* Per-thread PROJ context + private transform clone (KEEP: this is From 4aeb88afd4814089e925f51cfb7daaf2445d04b6 Mon Sep 17 00:00:00 2001 From: Kaushik Raja Date: Tue, 7 Jul 2026 20:24:03 -0700 Subject: [PATCH 06/11] r.proj: guard OpenMP timer calls for non-OpenMP builds 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. --- raster/r.proj/main.c | 33 ++++++++++++++++++++++----------- 1 file changed, 22 insertions(+), 11 deletions(-) diff --git a/raster/r.proj/main.c b/raster/r.proj/main.c index d4d181a6936..5c126fd8ae0 100644 --- a/raster/r.proj/main.c +++ b/raster/r.proj/main.c @@ -66,7 +66,18 @@ #include #include "r.proj.h" +#ifdef _OPENMP #include +static inline double rproj_wtime(void) +{ + return omp_get_wtime(); +} +#else +static inline double rproj_wtime(void) +{ + return 0.0; +} +#endif /* modify this table to add new methods */ struct menu menu[] = { @@ -811,7 +822,7 @@ int main(int argc, char **argv) * 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(); + double ts = rproj_wtime(); int band_orows = outcellhd.rows - obr0; int imin = 0, imax = -1; for (;;) { @@ -834,7 +845,7 @@ int main(int argc, char **argv) strip_rows, cap_mb); band_orows = (band_orows + 1) / 2; /* halve (round up), re-sample */ } - t_size += omp_get_wtime() - ts; + t_size += rproj_wtime() - ts; int obr1 = obr0 + band_orows; int strip_rows = imax - imin + 1; @@ -848,7 +859,7 @@ int main(int argc, char **argv) void *strip = NULL; if (strip_rows > 0) { strip = G_malloc((size_t)strip_rows * incellhd.cols * cell_size); - double t0 = omp_get_wtime(); + double t0 = rproj_wtime(); G_switch_env(); /* -> input */ for (int r = imin; r <= imax; r++) Rast_get_row(fdi, @@ -856,7 +867,7 @@ int main(int argc, char **argv) (size_t)(r - imin) * incellhd.cols * cell_size, r, cell_type); G_switch_env(); /* -> output */ - t_fill += omp_get_wtime() - t0; + t_fill += rproj_wtime() - t0; } /* Per-band output buffer, lock-free disjoint row slots (band-relative @@ -864,7 +875,7 @@ int main(int argc, char **argv) void *band_out = G_malloc((size_t)band_orows * outcellhd.cols * cell_size); - double t1 = omp_get_wtime(); + double t1 = rproj_wtime(); /* Each band runs one parallel region. This is not nested parallelism: * the "omp for" below does not create a second thread team, it only * divides the band's output rows among the threads that this "omp @@ -914,17 +925,17 @@ int main(int argc, char **argv) proj_destroy(tproj_local.pj); proj_context_destroy(thread_ctx); } - t_compute += omp_get_wtime() - t1; + t_compute += rproj_wtime() - t1; /* Serial in-order write of the band's rows (Rast_put_row sequential). */ - double t2 = omp_get_wtime(); + double t2 = rproj_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; + t_write += rproj_wtime() - t2; G_percent(obr1, outcellhd.rows, 5); @@ -934,9 +945,9 @@ int main(int argc, char **argv) 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); + G_debug(1, + "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 */ From e71d691441676b17d59d77b5e01666134292f7ec Mon Sep 17 00:00:00 2001 From: Kaushik Raja Date: Tue, 7 Jul 2026 20:25:14 -0700 Subject: [PATCH 07/11] cmake: Add PROJ dependency for r.proj 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. --- raster/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/raster/CMakeLists.txt b/raster/CMakeLists.txt index ff92138b765..91bd13800b0 100644 --- a/raster/CMakeLists.txt +++ b/raster/CMakeLists.txt @@ -429,6 +429,7 @@ build_program_in_subdir( grass_raster grass_gproj grass_parson + PROJ::proj ${LIBM} OPTIONAL_DEPENDS OpenMP::OpenMP_C) From 98d43bec54f4ee3d0d304f90417f2a65a52b0cdd Mon Sep 17 00:00:00 2001 From: Kaushik Raja Date: Wed, 8 Jul 2026 19:26:09 -0700 Subject: [PATCH 08/11] lib/proj: add per-thread transform clone helpers 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. --- include/grass/defs/gprojects.h | 2 ++ include/grass/gprojects.h | 9 +++++++++ lib/proj/do_proj.c | 35 ++++++++++++++++++++++++++++++++++ raster/r.proj/main.c | 12 +++++------- 4 files changed, 51 insertions(+), 7 deletions(-) diff --git a/include/grass/defs/gprojects.h b/include/grass/defs/gprojects.h index 1270ed66c5a..2f21ba91230 100644 --- a/include/grass/defs/gprojects.h +++ b/include/grass/defs/gprojects.h @@ -9,6 +9,8 @@ int GPJ_transform(const struct pj_info *, const struct pj_info *, int GPJ_transform_array(const struct pj_info *, const struct pj_info *, const struct pj_info *, int, double *, double *, double *, int); +void GPJ_clone_transform(const struct pj_info *, struct gpj_transform_clone *); +void GPJ_free_transform_clone(struct gpj_transform_clone *); /* old API, to be removed */ int pj_do_proj(double *, double *, const struct pj_info *, diff --git a/include/grass/gprojects.h b/include/grass/gprojects.h index 803eecf8477..8fc71b8f2a9 100644 --- a/include/grass/gprojects.h +++ b/include/grass/gprojects.h @@ -48,6 +48,15 @@ struct pj_info { char *wkt; }; +/* Per-thread clone of a transform, filled by GPJ_clone_transform() and + * released by GPJ_free_transform_clone(). Bundles the cloned transform with + * the private PROJ context it was cloned into, so ownership is a single unit. + */ +struct gpj_transform_clone { + struct pj_info info; + PJ_CONTEXT *ctx; +}; + struct gpj_datum { char *name, *longname, *ellps; double dx, dy, dz; diff --git a/lib/proj/do_proj.c b/lib/proj/do_proj.c index d297022384d..a086d47cdd7 100644 --- a/lib/proj/do_proj.c +++ b/lib/proj/do_proj.c @@ -1412,3 +1412,38 @@ int pj_do_transform(int count, double *x, double *y, double *h, } return ok; } + +/*! + * \brief Clone a transform into a fresh per-thread PROJ context + * + * PROJ transformation objects are not safe for concurrent use, so each thread + * needs its own. This fills \p clone with a copy of \p src whose PJ is cloned + * into a new private context. Release it with GPJ_free_transform_clone(). + * + * Safe to call concurrently from multiple threads with the same \p src, + * provided \p src is not modified during the calls: each call clones into its + * own new context and touches no shared mutable state. + * + * \param src source transform (as set up by GPJ_init_transform()) + * \param[out] clone receives the per-thread clone (info plus private context) + */ +void GPJ_clone_transform(const struct pj_info *src, + struct gpj_transform_clone *clone) +{ + clone->ctx = proj_context_create(); + clone->info = *src; + clone->info.pj = proj_clone(clone->ctx, src->pj); +} + +/*! + * \brief Free a per-thread transform clone and its context + * + * \param clone clone filled by GPJ_clone_transform(); its cloned PJ is set to + * NULL after release + */ +void GPJ_free_transform_clone(struct gpj_transform_clone *clone) +{ + proj_destroy(clone->info.pj); + proj_context_destroy(clone->ctx); + clone->info.pj = NULL; +} diff --git a/raster/r.proj/main.c b/raster/r.proj/main.c index 5c126fd8ae0..fec3a319ca9 100644 --- a/raster/r.proj/main.c +++ b/raster/r.proj/main.c @@ -890,9 +890,8 @@ int main(int argc, char **argv) * 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(); - tproj_local.pj = proj_clone(thread_ctx, tproj.pj); + struct gpj_transform_clone tproj_local; + GPJ_clone_transform(&tproj, &tproj_local); #pragma omp for private(row, col) schedule(dynamic) for (row = obr0; row < obr1; row++) { @@ -909,8 +908,8 @@ int main(int argc, char **argv) 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) { + if (GPJ_transform(&oproj, &iproj, &tproj_local.info, PJ_FWD, + &x1, &y1, NULL) < 0) { Rast_set_null_value(obufptr, 1, cell_type); } else { @@ -922,8 +921,7 @@ int main(int argc, char **argv) } } - proj_destroy(tproj_local.pj); - proj_context_destroy(thread_ctx); + GPJ_free_transform_clone(&tproj_local); } t_compute += rproj_wtime() - t1; From 9dd408fe84a4ce32ae547b8eaf4a7772877e27d6 Mon Sep 17 00:00:00 2001 From: Kaushik Raja Date: Wed, 8 Jul 2026 19:32:52 -0700 Subject: [PATCH 09/11] r.proj: drop direct PROJ dependency now that proj calls moved to lib/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. --- raster/CMakeLists.txt | 1 - raster/r.proj/Makefile | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/raster/CMakeLists.txt b/raster/CMakeLists.txt index 91bd13800b0..ff92138b765 100644 --- a/raster/CMakeLists.txt +++ b/raster/CMakeLists.txt @@ -429,7 +429,6 @@ build_program_in_subdir( grass_raster grass_gproj grass_parson - PROJ::proj ${LIBM} OPTIONAL_DEPENDS OpenMP::OpenMP_C) diff --git a/raster/r.proj/Makefile b/raster/r.proj/Makefile index 147b47fe2d8..083cd95c72f 100644 --- a/raster/r.proj/Makefile +++ b/raster/r.proj/Makefile @@ -2,7 +2,7 @@ MODULE_TOPDIR = ../.. PGM = r.proj -LIBES = $(GPROJLIB) $(RASTERLIB) $(GISLIB) $(MATHLIB) $(PARSONLIB) $(PROJLIB) +LIBES = $(GPROJLIB) $(RASTERLIB) $(GISLIB) $(MATHLIB) $(PARSONLIB) DEPENDENCIES = $(GPROJDEP) $(RASTERDEP) $(GISDEP) EXTRA_LIBS = $(OPENMP_LIBPATH) $(OPENMP_LIB) From 3a48d213fec77c91ef67c2c1255709ba0b28481a Mon Sep 17 00:00:00 2001 From: Kaushik Raja Date: Wed, 8 Jul 2026 15:38:09 -0700 Subject: [PATCH 10/11] r.proj: experimental per-thread-fd parallel strip reads 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. --- raster/r.proj/main.c | 65 ++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 60 insertions(+), 5 deletions(-) diff --git a/raster/r.proj/main.c b/raster/r.proj/main.c index fec3a319ca9..13d90126c7a 100644 --- a/raster/r.proj/main.c +++ b/raster/r.proj/main.c @@ -784,6 +784,32 @@ int main(int argc, char **argv) cell_type = FCELL_TYPE; cell_size = Rast_cell_size(cell_type); + /* Parallel input reads: decide the read-thread count here, in the INPUT + * env, so the mask guard checks the source mapset's mask (the mask that + * would apply to Rast_get_row on fdi). Rast_disable_omp_on_mask returns 1 + * (serial) if a mask is present or without OpenMP, and does NOT touch the + * thread count when no mask exists (lib/raster/mask_info.c:226-231), so the + * compute region's threads are unperturbed in the common case. When + * read_nprocs > 1 we open that many FRESH read fds (one per thread); fdi is + * used only by the serial fallback. + * INFERRED-safe (not yet runtime-verified; the gate converts it): + * concurrent Rast_open_old fds on the same map across locations is sound + * from the r.neighbors same-location precedent (in_fd[t]) plus the Stage 1 + * fcb analysis (each fd carries its own cur_row/data/data_fd; reads depend + * only on the fcb and R__.rd_window). */ +#ifdef _OPENMP + int want_nprocs = omp_get_max_threads(); +#else + int want_nprocs = 1; +#endif + int read_nprocs = Rast_disable_omp_on_mask(want_nprocs); + int *fd_read = NULL; + if (read_nprocs > 1) { + fd_read = G_malloc((size_t)read_nprocs * sizeof(int)); + for (int t = 0; t < read_nprocs; t++) + fd_read[t] = Rast_open_old(inmap->answer, setname); + } + /* 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. @@ -861,11 +887,35 @@ int main(int argc, char **argv) strip = G_malloc((size_t)strip_rows * incellhd.cols * cell_size); double t0 = rproj_wtime(); G_switch_env(); /* -> input */ - for (int r = imin; r <= imax; r++) - Rast_get_row(fdi, - (unsigned char *)strip + - (size_t)(r - imin) * incellhd.cols * cell_size, - r, cell_type); + if (read_nprocs > 1) { +#ifdef _OPENMP + /* Parallel read: each thread reads a contiguous, disjoint + * block of strip rows (schedule(static)) through its OWN fd + * into its own disjoint strip slice (slice = row r - imin). + * No two threads share an fd or a strip row. */ +#pragma omp parallel num_threads(read_nprocs) + { + int t = omp_get_thread_num(); +#pragma omp for schedule(static) + for (int r = imin; r <= imax; r++) + Rast_get_row(fd_read[t], + (unsigned char *)strip + + (size_t)(r - imin) * incellhd.cols * + cell_size, + r, cell_type); + } +#endif + } + else { + /* Serial fallback (nprocs==1, mask present, or no OpenMP): + * original loop, unchanged, through fdi. */ + for (int r = imin; r <= imax; r++) + Rast_get_row(fdi, + (unsigned char *)strip + (size_t)(r - imin) * + incellhd.cols * + cell_size, + r, cell_type); + } G_switch_env(); /* -> output */ t_fill += rproj_wtime() - t0; } @@ -950,6 +1000,11 @@ int main(int argc, char **argv) /* Close input map in its own env, then the output map. */ G_switch_env(); /* -> input */ Rast_close(fdi); + if (fd_read) { + for (int t = 0; t < read_nprocs; t++) + Rast_close(fd_read[t]); + G_free(fd_read); + } G_switch_env(); /* -> output */ Rast_close(fdo); From b46ee6ca18122ece67e91b83e840c266e96f1c81 Mon Sep 17 00:00:00 2001 From: Kaushik Raja Date: Thu, 9 Jul 2026 23:24:42 -0700 Subject: [PATCH 11/11] r.proj: add adaptive column tiling for oblique reprojections 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. --- raster/r.proj/main.c | 388 +++++++++++++++++++++++++++++++------------ 1 file changed, 278 insertions(+), 110 deletions(-) diff --git a/raster/r.proj/main.c b/raster/r.proj/main.c index 13d90126c7a..e968b79c740 100644 --- a/raster/r.proj/main.c +++ b/raster/r.proj/main.c @@ -130,29 +130,30 @@ static void interpolate_strip(void *strip, void *obufptr, int cell_type, 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. */ +/* Dense edge-walk of an output tile's rectangle [obr0, obr1) x [obc0, obc1) + * projected into input space; returns the min/max INPUT ROW touched, plus a + * 2-cell margin, clamped to the input map. Samples the tile's top and bottom + * rows across its columns [obc0, obc1) and its left and right columns across + * its rows (bordwalk-style), so a curved transform's interior-edge extremum is + * caught -- corner-only sampling can under-size the strip. A full-width band is + * the case obc0=0, obc1=cols. Called serially, before the parallel region, so + * the shared tproj is safe here. Returns imax < imin for a tile 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) + int obc0, int obc1, 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 */ + /* top edge (row obr0) and bottom edge (row obr1-1), tile 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++) { + for (c = obc0; c < obc1; 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) @@ -164,9 +165,9 @@ static void band_input_row_span(const struct Cell_head *ohd, rmax = ri; } } - /* left edge (col 0) and right edge (col cols-1), all band rows */ + /* left edge (col obc0) and right edge (col obc1-1), all band rows */ for (e = 0; e < 2; e++) { - int ocol = (e == 0) ? 0 : (ohd->cols - 1); + int ocol = (e == 0) ? obc0 : (obc1 - 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; @@ -197,6 +198,75 @@ static void band_input_row_span(const struct Cell_head *ohd, *imax = hi; } +/* Largest input-row strip (in rows) among the column tiles of width tilew that + * partition output columns [0, ohd->cols) for the band [obr0, obr1). Tiles are + * loaded one at a time, so peak strip memory is set by the worst tile, not the + * union of the band's tiles; the fit search sizes this against the cap. Every + * call is a full tile edge-walk, so this is O(tiles * perimeter) -- paid in the + * serial size phase, and only when column splitting is actually entered. + * Returns 0 if every tile projects entirely outside the input. */ +static int worst_tile_strip_rows(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 tilew) +{ + int worst = 0, obc0; + + for (obc0 = 0; obc0 < ohd->cols; obc0 += tilew) { + int obc1 = obc0 + tilew; + int imin, imax, rows; + + if (obc1 > ohd->cols) + obc1 = ohd->cols; + band_input_row_span(ohd, ihd, oproj, iproj, tproj, obr0, obr1, obc0, + obc1, &imin, &imax); + rows = imax - imin + 1; /* imax < imin (empty) -> <= 0, ignored */ + if (rows > worst) + worst = rows; + } + return worst; +} + +#define TILE_PROBE 16 /* tiles sampled by the Phase-2 width-search estimate */ + +/* Cheap estimate of worst_tile_strip_rows: the largest input-row strip among + * at most `probe` column tiles, evenly spaced across the band width and always + * including the first and last. A subset max is a LOWER bound on the true + * worst, so it only PRUNES the Phase-2 search; the chosen width is exact- + * validated by worst_tile_strip_rows before use. */ +static int est_worst_tile_strip_rows(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 tilew, int probe) +{ + int ntiles = (ohd->cols + tilew - 1) / tilew; + int worst = 0, k; + + if (probe < 1) + probe = 1; + if (probe > ntiles) + probe = ntiles; + for (k = 0; k < probe; k++) { + int ti = (probe == 1) ? 0 : (int)((long)k * (ntiles - 1) / (probe - 1)); + int obc0 = ti * tilew; + int obc1 = obc0 + tilew; + int imin, imax, rows; + + if (obc1 > ohd->cols) + obc1 = ohd->cols; + band_input_row_span(ohd, ihd, oproj, iproj, tproj, obr0, obr1, obc0, + obc1, &imin, &imax); + rows = imax - imin + 1; + if (rows > worst) + worst = rows; + } + return worst; +} + int main(int argc, char **argv) { char *mapname, /* ptr to name of output layer */ @@ -839,21 +909,28 @@ int main(int argc, char **argv) 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; + int max_tiles = 1; /* most column tiles used by any single band */ 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. */ + /* Fit search. Phase 1 (fast path, unchanged): halve the band height + * until the FULL-WIDTH strip plus the band output buffer fit the cap; + * the span is re-run per candidate height. Phase 2 (oblique fallback): + * only if a single full-width output row still busts the cap, split the + * row into column tiles and halve tile WIDTH until the worst tile's + * strip fits. Strips are full input width (the raster API reads whole + * rows), so width splitting shrinks a tile's input ROW span, not its + * width. Easy pairs never leave Phase 1. */ double ts = rproj_wtime(); int band_orows = outcellhd.rows - obr0; + int tilew = outcellhd.cols; int imin = 0, imax = -1; for (;;) { band_input_row_span(&outcellhd, &incellhd, &oproj, &iproj, &tproj, - obr0, obr0 + band_orows, &imin, &imax); + obr0, obr0 + band_orows, 0, outcellhd.cols, + &imin, &imax); int strip_rows = imax - imin + 1; size_t strip_bytes = strip_rows > 0 ? (size_t)strip_rows * incellhd.cols * cell_size @@ -862,121 +939,213 @@ int main(int argc, char **argv) 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); + break; /* height exhausted: fall through to column splitting */ band_orows = (band_orows + 1) / 2; /* halve (round up), re-sample */ } + if (band_orows == 1) { + /* Phase 2 (oblique only): Phase 1 could not fit even a single + * full-width row, so prefer a TALL tiled band instead of a 1-row + * one. Rescan from the full remaining height downward; at the + * tallest height whose output buffer fits the cap, halve tile WIDTH + * until the worst column tile's strip fits, and only reduce height + * when no width fits. Keeping the band tall gives the per-band + * parallel region many output rows. Runs only on this path, so the + * easy-pair size phase (Phase 1) is unaffected. */ + band_orows = outcellhd.rows - obr0; + for (;;) { + size_t out_bytes = + (size_t)band_orows * outcellhd.cols * cell_size; + if (out_bytes <= cap_bytes) { + /* Upper tier: cheap probe estimate narrows to a candidate + * width, checking down to tilew==1. The estimate is a lower + * bound, so est-no-fit at tilew==1 implies exact-no-fit -> + * the exact validation below is skipped entirely at heights + * where no width can fit (this is what keeps the search + * cheap; scanning every tile there was the cost). */ + tilew = outcellhd.cols; + int est_fit = 0; + for (;;) { + int est = est_worst_tile_strip_rows( + &outcellhd, &incellhd, &oproj, &iproj, &tproj, obr0, + obr0 + band_orows, tilew, TILE_PROBE); + size_t est_bytes = + est > 0 ? (size_t)est * incellhd.cols * cell_size + : 0; + if (est_bytes + out_bytes <= cap_bytes) { + est_fit = 1; + break; + } + if (tilew == 1) + break; + tilew = (tilew + 1) / 2; + } + /* Lower tier: EXACT validation, only when the estimate + * found a candidate. Narrow and re-validate if the estimate + * was optimistic; this exact-sizes the accepted width so + * the cap is honored. */ + int fit = 0; + if (est_fit) { + for (;;) { + int worst = worst_tile_strip_rows( + &outcellhd, &incellhd, &oproj, &iproj, &tproj, + obr0, obr0 + band_orows, tilew); + size_t strip_bytes = + worst > 0 + ? (size_t)worst * incellhd.cols * cell_size + : 0; + if (strip_bytes + out_bytes <= cap_bytes) { + fit = 1; + break; + } + if (tilew == 1) + break; /* exact: no width fits at this height */ + tilew = (tilew + 1) / 2; + } + } + if (fit) + break; + } + if (band_orows == 1) { + /* Single output row at minimum width still over cap = + * singular/large-halo; needs the tile-cache path. */ + size_t out1 = (size_t)outcellhd.cols * cell_size; + int worst = worst_tile_strip_rows(&outcellhd, &incellhd, + &oproj, &iproj, &tproj, + obr0, obr0 + 1, 1); + size_t strip_bytes = + worst > 0 ? (size_t)worst * incellhd.cols * cell_size + : 0; + 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 + out1) / (1024.0 * 1024.0), worst, + cap_mb); + } + band_orows = (band_orows + 1) / 2; /* shrink height, retry */ + } + } t_size += rproj_wtime() - ts; int obr1 = obr0 + band_orows; - int strip_rows = imax - imin + 1; n_bands++; + int n_tiles = (outcellhd.cols + tilew - 1) / tilew; + if (n_tiles > max_tiles) + max_tiles = n_tiles; - /* 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 = rproj_wtime(); - G_switch_env(); /* -> input */ - if (read_nprocs > 1) { + /* Per-band output buffer, lock-free disjoint row slots, filled column + * tile by column tile and written once after all tiles. Full width + * regardless of tiling. */ + void *band_out = + G_malloc((size_t)band_orows * outcellhd.cols * cell_size); + + /* Column tiles processed one at a time: only the current tile's strip + * is resident, so peak strip memory is the worst tile, not the band's + * union. tilew == cols is the single-tile fast path (obc0=0, + * obc1=cols), identical to un-tiled banding. */ + for (int obc0 = 0; obc0 < outcellhd.cols; obc0 += tilew) { + int obc1 = obc0 + tilew; + if (obc1 > outcellhd.cols) + obc1 = outcellhd.cols; + + /* Per-tile input row span (full-width strip: the raster API reads + * whole rows, so columns are not cropped). */ + band_input_row_span(&outcellhd, &incellhd, &oproj, &iproj, &tproj, + obr0, obr1, obc0, obc1, &imin, &imax); + int strip_rows = imax - imin + 1; + + /* Serial strip load (single fd -> get_row not thread-safe). EMPTY + * TILE: strip_rows <= 0 -> projects outside input, no read; 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 = rproj_wtime(); + G_switch_env(); /* -> input */ + if (read_nprocs > 1) { #ifdef _OPENMP - /* Parallel read: each thread reads a contiguous, disjoint - * block of strip rows (schedule(static)) through its OWN fd - * into its own disjoint strip slice (slice = row r - imin). - * No two threads share an fd or a strip row. */ + /* Parallel read: each thread reads a contiguous, disjoint + * block of strip rows through its OWN fd into its own + * disjoint strip slice. No two threads share an fd/row. */ #pragma omp parallel num_threads(read_nprocs) - { - int t = omp_get_thread_num(); + { + int t = omp_get_thread_num(); #pragma omp for schedule(static) + for (int r = imin; r <= imax; r++) + Rast_get_row(fd_read[t], + (unsigned char *)strip + + (size_t)(r - imin) * + incellhd.cols * cell_size, + r, cell_type); + } +#endif + } + else { + /* Serial fallback (nprocs==1, mask, or no OpenMP). */ for (int r = imin; r <= imax; r++) - Rast_get_row(fd_read[t], + Rast_get_row(fdi, (unsigned char *)strip + (size_t)(r - imin) * incellhd.cols * cell_size, r, cell_type); } -#endif + G_switch_env(); /* -> output */ + t_fill += rproj_wtime() - t0; } - else { - /* Serial fallback (nprocs==1, mask present, or no OpenMP): - * original loop, unchanged, through fdi. */ - for (int r = imin; r <= imax; r++) - Rast_get_row(fdi, - (unsigned char *)strip + (size_t)(r - imin) * - incellhd.cols * - cell_size, - r, cell_type); - } - G_switch_env(); /* -> output */ - t_fill += rproj_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 = rproj_wtime(); - /* Each band runs one parallel region. This is not nested parallelism: - * the "omp for" below does not create a second thread team, it only - * divides the band's output rows among the threads that this "omp - * parallel" created. The two directives are kept separate instead of a - * combined "omp parallel for" because every thread must clone its own - * PROJ context before the row loop starts and destroy it after the loop - * ends, and that per-thread setup has to sit inside the parallel region - * but outside the for. */ + double t1 = rproj_wtime(); + /* One parallel region per tile. Not nested: the "omp for" divides + * the band's output rows among this region's threads. Separate + * directives so each thread clones its PROJ context before the row + * loop and destroys it after. */ #pragma omp parallel - { - /* 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 gpj_transform_clone tproj_local; - GPJ_clone_transform(&tproj, &tproj_local); + { + struct gpj_transform_clone tproj_local; + GPJ_clone_transform(&tproj, &tproj_local); #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.info, 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); + 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 = obc0; col < obc1; 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.info, + 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); + } } } + + GPJ_free_transform_clone(&tproj_local); } + t_compute += rproj_wtime() - t1; - GPJ_free_transform_clone(&tproj_local); + if (strip) + G_free(strip); } - t_compute += rproj_wtime() - t1; - /* Serial in-order write of the band's rows (Rast_put_row sequential). - */ + /* Serial in-order write of the band's rows once all tiles filled + * band_out (Rast_put_row sequential). */ double t2 = rproj_wtime(); for (row = obr0; row < obr1; row++) Rast_put_row(fdo, @@ -987,15 +1156,14 @@ int main(int argc, char **argv) G_percent(obr1, outcellhd.rows, 5); - if (strip) - G_free(strip); G_free(band_out); obr0 = obr1; } G_debug(1, - "PHASE_TIMERS size=%.4f fill=%.4f compute=%.4f write=%.4f bands=%d", - t_size, t_fill, t_compute, t_write, n_bands); + "PHASE_TIMERS size=%.4f fill=%.4f compute=%.4f write=%.4f bands=%d " + "tiles=%d", + t_size, t_fill, t_compute, t_write, n_bands, max_tiles); /* Close input map in its own env, then the output map. */ G_switch_env(); /* -> input */