From f033870a9c3f71625cf1965fa711d5398deac2b0 Mon Sep 17 00:00:00 2001 From: Kaushik Raja Date: Thu, 28 May 2026 18:13:05 -0700 Subject: [PATCH 01/13] WIP r.param.scale per-thread file descriptor parallelization --- raster/r.param.scale/Makefile | 7 +- raster/r.param.scale/interface.c | 5 +- raster/r.param.scale/main.c | 3 +- raster/r.param.scale/param.h | 3 +- raster/r.param.scale/process.c | 423 +++++++++++++++---------------- 5 files changed, 209 insertions(+), 232 deletions(-) diff --git a/raster/r.param.scale/Makefile b/raster/r.param.scale/Makefile index 6cebd8a5d93..fc5779d5413 100644 --- a/raster/r.param.scale/Makefile +++ b/raster/r.param.scale/Makefile @@ -1,10 +1,9 @@ MODULE_TOPDIR = ../.. - PGM = r.param.scale - LIBES = $(GMATHLIB) $(RASTERLIB) $(GISLIB) $(MATHLIB) +EXTRA_LIBS = $(OPENMP_LIBPATH) $(OPENMP_LIB) DEPENDENCIES = $(GMATHDEP) $(RASTERDEP) $(GISDEP) - +EXTRA_CFLAGS = $(OPENMP_CFLAGS) +EXTRA_INC = $(OPENMP_INCPATH) include $(MODULE_TOPDIR)/include/Make/Module.make - default: cmd diff --git a/raster/r.param.scale/interface.c b/raster/r.param.scale/interface.c index 2dac926caad..4a923c3d81c 100644 --- a/raster/r.param.scale/interface.c +++ b/raster/r.param.scale/interface.c @@ -31,7 +31,8 @@ void interface(int argc, char **argv) *win_size, /* Size of side of local window. */ *parameter, /* Morphometric parameter to calculate. */ *expon, /* Inverse distance exponent for weight. */ - *vert_sc; /* Vertical scaling factor. */ + *vert_sc, /* Vertical scaling factor. */ + *nprocs_opt; struct Flag *constr; /* Forces quadratic through the central */ @@ -67,6 +68,7 @@ void interface(int argc, char **argv) parameter = G_define_option(); expon = G_define_option(); vert_sc = G_define_option(); + nprocs_opt = G_define_standard_option(G_OPT_M_NPROCS); constr = G_define_flag(); @@ -130,6 +132,7 @@ void interface(int argc, char **argv) rast_in_name = rast_in->answer; /* Now keyboard input has been parsed, */ rast_out_name = rast_out->answer; /* can place the contents into strings */ wsize = atoi(win_size->answer); + nprocs = G_set_omp_num_threads(nprocs_opt); constrained = constr->answer; sscanf(expon->answer, "%lf", &exponent); sscanf(vert_sc->answer, "%lf", &zscale); diff --git a/raster/r.param.scale/main.c b/raster/r.param.scale/main.c index 8bb76b3ec67..9b13b6d4e3a 100644 --- a/raster/r.param.scale/main.c +++ b/raster/r.param.scale/main.c @@ -22,7 +22,8 @@ int constrained; /* Flag that forces quadtratic through */ int fd_in, /* File descriptor for input and */ fd_out, /* output raster files. */ wsize, /* Size of local processing window. */ - mparam; /* Morphometric parameter to calculate. */ + mparam, /* Morphometric parameter to calculate. */ + nprocs; /* NUmber of threads in process. */ double resoln, /* Planimetric resolution. */ exponent, /* Distance weighting exponent. */ diff --git a/raster/r.param.scale/param.h b/raster/r.param.scale/param.h index f9f6d5e80f5..a0d806b7031 100644 --- a/raster/r.param.scale/param.h +++ b/raster/r.param.scale/param.h @@ -88,7 +88,8 @@ extern int constrained; /* Flag that forces quadtratic through */ extern int fd_in, /* File descriptor for input and */ fd_out, /* output raster files. */ wsize, /* Size of local processing window. */ - mparam; /* Morphometric parameter to calculate. */ + mparam, /* Morphometric parameter to calculate. */ + nprocs; extern double resoln, /* Planimetric resolution. */ exponent, /* Distance weighting exponent. */ diff --git a/raster/r.param.scale/process.c b/raster/r.param.scale/process.c index 6a0da1dbd87..a5d28a2b740 100644 --- a/raster/r.param.scale/process.c +++ b/raster/r.param.scale/process.c @@ -1,15 +1,10 @@ /*****************************************************************************/ - -/*** ***/ - /*** process() ***/ - -/*** Reads in a raster map row by row for processing. ***/ - -/*** Jo Wood, Project ASSIST, V1.0 7th February 1993 ***/ - -/*** ***/ - +/*** Parallel (per-thread file descriptors): each thread opens its own ***/ +/*** file descriptor and reads its strip + halo directly from disk ***/ +/*** into a private buffer. No full map in RAM. Memory-light path for ***/ +/*** maps that do not fit in RAM. ***/ +/*** Based on the serial version by Jo Wood, Project ASSIST, 1993. ***/ /*****************************************************************************/ #include @@ -20,56 +15,21 @@ #include "param.h" #include "nrutil.h" +#if defined(_OPENMP) +#include +#endif + void process(void) { - - /*-----------------------------------------------------------------------*/ - /* INITIALISE */ - - /*-----------------------------------------------------------------------*/ - - DCELL *row_in, /* Buffer large enough to hold `wsize' */ - *row_out = NULL, /* raster rows. When GRASS reads in a */ - /* raster row, each element is of type */ - /* DCELL */ - *window_ptr, /* Stores local terrain window. */ - centre, /* Elevation of central cell in window. */ - *window_cell; /* Stores last read cell in window. */ - - CELL *featrow_out = NULL; /* store features in CELL */ - struct Cell_head region; /* Structure to hold region information */ + int nrows, ncols, row; - int nrows, /* Will store the current number of */ - ncols, /* rows and columns in the raster. */ - row, col, /* Counts through each row and column */ - /* of the input raster. */ - wind_row, /* Counts through each row and column */ - wind_col, /* of the local neighbourhood window. */ - *index_ptr; /* Row permutation vector for LU decomp. */ - - double **normal_ptr, /* Cross-products matrix. */ - *obs_ptr, /* Observed vector. */ - temp; /* Unused */ - - double *weight_ptr; /* Weighting matrix for observed values. */ - int found_null; /* If null was found among window cells */ - - /*-----------------------------------------------------------------------*/ - /* GET RASTER AND WINDOW DETAILS */ - - /*-----------------------------------------------------------------------*/ - - G_get_window(®ion); /* Fill out the region structure (the */ - /* geographical limits etc.) */ + G_get_window(®ion); + nrows = Rast_window_rows(); + ncols = Rast_window_cols(); - nrows = Rast_window_rows(); /* Find out the number of rows and */ - ncols = Rast_window_cols(); /* columns of the raster. */ - - if ((region.ew_res / region.ns_res >= - 1.01) || /* If EW and NS resolns are */ - (region.ns_res / region.ew_res >= - 1.01)) { /* >1% different, warn user. */ + if ((region.ew_res / region.ns_res >= 1.01) || + (region.ns_res / region.ew_res >= 1.01)) { G_warning( _("E-W and N-S grid resolutions are different. Taking average.")); resoln = (region.ns_res + region.ew_res) / 2; @@ -78,209 +38,222 @@ void process(void) resoln = region.ns_res; /*-----------------------------------------------------------------------*/ - /* RESERVE MEMORY TO HOLD Z VALUES AND MATRICES */ - + /* CALCULATE LEAST SQUARES COEFFICIENTS (ONCE) */ + /* Same as serial. Weights, normal matrix, and LU decomposition depend */ + /* only on window size, not on data. Computed once, read-only by all */ + /* threads. */ /*-----------------------------------------------------------------------*/ - row_in = (DCELL *)G_malloc(ncols * sizeof(DCELL) * wsize); - /* Reserve `wsize' rows of memory. */ - - if (mparam != FEATURE) { - row_out = - Rast_allocate_buf(DCELL_TYPE); /* Initialise output row buffer. */ - Rast_set_d_null_value(row_out, ncols); - } - else { - featrow_out = - Rast_allocate_buf(CELL_TYPE); /* Initialise output row buffer. */ - Rast_set_c_null_value(featrow_out, ncols); - } - - window_ptr = (DCELL *)G_malloc(SQR(wsize) * sizeof(DCELL)); - /* Reserve enough memory for local wind. */ - - weight_ptr = (double *)G_malloc(SQR(wsize) * sizeof(double)); - /* Reserve enough memory weights matrix. */ - - normal_ptr = dmatrix(0, 5, 0, 5); /* Allocate memory for 6*6 matrix */ - index_ptr = ivector(0, 5); /* and for 1D vector holding indices */ - obs_ptr = dvector(0, 5); /* and for 1D vector holding observed z */ - - /* ---------------------------------------------------------------- */ - /* - CALCULATE LEAST SQUARES COEFFICIENTS - */ - /* ---------------------------------------------------------------- */ - - /*--- Calculate weighting matrix. ---*/ + double *weight_ptr = (double *)G_malloc(SQR(wsize) * sizeof(double)); + double **normal_ptr = dmatrix(0, 5, 0, 5); + int *index_ptr = ivector(0, 5); + double temp; find_weight(weight_ptr); + find_normal(normal_ptr, weight_ptr); - /* Initial coefficients need only be found once since they are - constant for any given window size. The only element that - changes is the observed vector (RHS of normal equations). */ - - /*--- Find normal equations in matrix form. ---*/ + if (constrained) + G_ludcmp(normal_ptr, 5, index_ptr, &temp); + else + G_ludcmp(normal_ptr, 6, index_ptr, &temp); - find_normal(normal_ptr, weight_ptr); + /*-----------------------------------------------------------------------*/ + /* CHANGE 1: Per-thread input file descriptors. */ + /* Each thread opens its own descriptor to the input raster. Each thread */ + /* will read its strip+halo directly from disk into its private buffer. */ + /* This is the r.neighbors pattern. No full_map array exists. */ + /*-----------------------------------------------------------------------*/ - /*--- Apply LU decomposition to normal equations. ---*/ + int *fd_thread = (int *)G_malloc((size_t)nprocs * sizeof(int)); + for (int i = 0; i < nprocs; i++) + fd_thread[i] = Rast_open_old(rast_in_name, ""); - if (constrained) { - G_ludcmp(normal_ptr, 5, index_ptr, &temp); - /* To constrain the quadtratic - through the central cell, ignore - the calculations involving the - coefficient f. Since these are - all in the last row and column of - the matrix, simply redimension. */ - /* disp_matrix(normal_ptr,obs_ptr,obs_ptr,5); - */ - } + /*-----------------------------------------------------------------------*/ + /* CHANGE 2: Output slots, one pointer per row. */ + /* Threads fill their rows into slots, a serial */ + /* loop writes them in order at the end (Rast_put_row is not */ + /* thread-safe and must be called in row order). */ + /*-----------------------------------------------------------------------*/ - else { - G_ludcmp(normal_ptr, 6, index_ptr, &temp); - /* disp_matrix(normal_ptr,obs_ptr,obs_ptr,6); - */ + DCELL **row_out = (DCELL **)G_malloc((size_t)nrows * sizeof(DCELL *)); + CELL **featrow_out = (CELL **)G_malloc((size_t)nrows * sizeof(CELL *)); + for (row = 0; row < nrows; row++) { + row_out[row] = NULL; + featrow_out[row] = NULL; } /*-----------------------------------------------------------------------*/ - /* PROCESS INPUT RASTER AND WRITE OUT RASTER LINE BY LINE */ - + /* CHANGE 3: Parallel region. */ + /* Each thread claims its output row range, then reads its strip + halo */ + /* from disk via its own descriptor into a private buffer, and does */ + /* all window math reading from that private buffer. */ /*-----------------------------------------------------------------------*/ - if (mparam != FEATURE) - for (wind_row = 0; wind_row < EDGE; wind_row++) - Rast_put_row(fd_out, row_out, - DCELL_TYPE); /* Write out the edge cells as NULL. */ - else - for (wind_row = 0; wind_row < EDGE; wind_row++) - Rast_put_row(fd_out, featrow_out, - CELL_TYPE); /* Write out the edge cells as NULL. */ - - for (wind_row = 0; wind_row < wsize - 1; wind_row++) - Rast_get_row(fd_in, row_in + (wind_row * ncols), wind_row, DCELL_TYPE); - /* Read in enough of the first rows to */ - /* allow window to be examined. */ - - for (row = EDGE; row < (nrows - EDGE); row++) { - G_percent(row + 1, nrows - EDGE, 2); - - Rast_get_row(fd_in, row_in + ((wsize - 1) * ncols), row + EDGE, - DCELL_TYPE); - - for (col = EDGE; col < (ncols - EDGE); col++) { - /* Find central z value */ - centre = *(row_in + EDGE * ncols + col); - /* Test for no data and propagate */ - if (Rast_is_d_null_value(¢re)) { - if (mparam == FEATURE) - Rast_set_c_null_value(featrow_out + col, 1); +#pragma omp parallel if (nprocs > 1) + { + int nthreads = 1; + int tid = 0; +#if defined(_OPENMP) + nthreads = omp_get_num_threads(); + tid = omp_get_thread_num(); +#endif + + /* Split processable rows [EDGE .. nrows-EDGE) evenly across threads */ + int proc_first = EDGE; + int proc_last = nrows - EDGE; /* exclusive */ + int proc_count = proc_last - proc_first; + + int my_start = proc_first + (int)((long)proc_count * tid / nthreads); + int my_end = + proc_first + (int)((long)proc_count * (tid + 1) / nthreads); + + if (my_end > my_start) { + /* Input rows this thread needs: my_start-EDGE .. my_end-1+EDGE */ + int buf_first = my_start - EDGE; + int buf_last = my_end - 1 + EDGE; + int buf_rows = buf_last - buf_first + 1; + + /* Allocate private strip and read it from disk via THIS */ + /* thread's own descriptor. No shared RAM map; each thread */ + /* pulls only what it needs. */ + DCELL *strip = + (DCELL *)G_malloc((size_t)buf_rows * ncols * sizeof(DCELL)); + for (int r = 0; r < buf_rows; r++) + Rast_get_d_row(fd_thread[tid], strip + (size_t)r * ncols, + buf_first + r); + + /* Per-thread math scratch */ + DCELL *window_ptr = (DCELL *)G_malloc(SQR(wsize) * sizeof(DCELL)); + double *obs_ptr = dvector(0, 5); + + for (int orow = my_start; orow < my_end; orow++) { + DCELL *out_d = NULL; + CELL *out_c = NULL; + if (mparam != FEATURE) { + out_d = Rast_allocate_buf(DCELL_TYPE); + Rast_set_d_null_value(out_d, ncols); + } else { - Rast_set_d_null_value(row_out + col, 1); + out_c = Rast_allocate_buf(CELL_TYPE); + Rast_set_c_null_value(out_c, ncols); } - continue; - } - found_null = FALSE; - for (wind_row = 0; wind_row < wsize; wind_row++) { - if (found_null) - break; - for (wind_col = 0; wind_col < wsize; wind_col++) { - - /* Express all window values relative */ - /* to the central elevation. */ - window_cell = - row_in + (wind_row * ncols) + col + wind_col - EDGE; - /* Test for no data and propagate */ - if (Rast_is_d_null_value(window_cell)) { + + /* Position of this row inside the private strip */ + int srow = orow - buf_first; + + for (int col = EDGE; col < (ncols - EDGE); col++) { + DCELL centre = *(strip + (size_t)srow * ncols + col); + if (Rast_is_d_null_value(¢re)) { if (mparam == FEATURE) - Rast_set_c_null_value(featrow_out + col, 1); - else { - Rast_set_d_null_value(row_out + col, 1); + Rast_set_c_null_value(out_c + col, 1); + else + Rast_set_d_null_value(out_d + col, 1); + continue; + } + + int found_null = FALSE; + for (int wind_row = 0; wind_row < wsize; wind_row++) { + if (found_null) + break; + for (int wind_col = 0; wind_col < wsize; wind_col++) { + DCELL *window_cell = + strip + + (size_t)(srow - EDGE + wind_row) * ncols + col + + wind_col - EDGE; + if (Rast_is_d_null_value(window_cell)) { + if (mparam == FEATURE) + Rast_set_c_null_value(out_c + col, 1); + else + Rast_set_d_null_value(out_d + col, 1); + found_null = TRUE; + break; + } + *(window_ptr + (wind_row * wsize) + wind_col) = + *(window_cell)-centre; } - found_null = TRUE; - break; } - *(window_ptr + (wind_row * wsize) + wind_col) = - *(window_cell)-centre; + if (found_null) + continue; + + /* Math block: unchanged from serial */ + find_obs(window_ptr, obs_ptr, weight_ptr); + if (constrained) + G_lubksb(normal_ptr, 5, index_ptr, obs_ptr); + else + G_lubksb(normal_ptr, 6, index_ptr, obs_ptr); + + if (mparam == FEATURE) + *(out_c + col) = (CELL)feature(obs_ptr); + else + *(out_d + col) = param(mparam, obs_ptr); + + if (mparam == ELEV) + *(out_d + col) += centre; } + + if (mparam != FEATURE) + row_out[orow] = out_d; + else + featrow_out[orow] = out_c; } - if (found_null) - continue; + G_free(strip); + G_free(window_ptr); + free_dvector(obs_ptr, 0, 5); + } + } /* end parallel */ - /*--- Use LU back substitution to solve normal equations. ---*/ - find_obs(window_ptr, obs_ptr, weight_ptr); + /*-----------------------------------------------------------------------*/ + /* CHANGE 4: Close per-thread descriptors. */ + /*-----------------------------------------------------------------------*/ - /* disp_wind(window_ptr); - disp_matrix(normal_ptr,obs_ptr,obs_ptr,6); - */ + for (int i = 0; i < nprocs; i++) + Rast_close(fd_thread[i]); + G_free(fd_thread); - if (constrained) { - G_lubksb(normal_ptr, 5, index_ptr, obs_ptr); - /* - disp_matrix(normal_ptr,obs_ptr,obs_ptr,5); - */ - } + /*-----------------------------------------------------------------------*/ + /* CHANGE 5: Serial in-order write to disk. */ + /* Top and bottom EDGE rows are written as NULL, then the computed rows */ + /* are written from row_out[]/featrow_out[] in order. */ + /*-----------------------------------------------------------------------*/ - else { - G_lubksb(normal_ptr, 6, index_ptr, obs_ptr); - /* - disp_matrix(normal_ptr,obs_ptr,obs_ptr,6); - */ - } + DCELL *null_d = NULL; + CELL *null_c = NULL; + if (mparam != FEATURE) { + null_d = Rast_allocate_buf(DCELL_TYPE); + Rast_set_d_null_value(null_d, ncols); + } + else { + null_c = Rast_allocate_buf(CELL_TYPE); + Rast_set_c_null_value(null_c, ncols); + } - /*--- Calculate terrain parameter based on quad. coefficients. ---*/ - if (mparam == FEATURE) - *(featrow_out + col) = (CELL)feature(obs_ptr); + for (row = 0; row < nrows; row++) { + if (row < EDGE || row >= nrows - EDGE) { + if (mparam != FEATURE) + Rast_put_row(fd_out, null_d, DCELL_TYPE); else - *(row_out + col) = param(mparam, obs_ptr); - - if (mparam == ELEV) - *(row_out + col) += centre; /* Add central elevation back */ + Rast_put_row(fd_out, null_c, CELL_TYPE); + } + else { + if (mparam != FEATURE) { + Rast_put_row(fd_out, row_out[row], DCELL_TYPE); + G_free(row_out[row]); + } + else { + Rast_put_row(fd_out, featrow_out[row], CELL_TYPE); + G_free(featrow_out[row]); + } } - - if (mparam != FEATURE) - Rast_put_row(fd_out, row_out, - DCELL_TYPE); /* Write the row buffer to the output */ - /* raster. */ - else /* write FEATURE to CELL */ - Rast_put_row(fd_out, featrow_out, - CELL_TYPE); /* Write the row buffer to the output */ - /* raster. */ - - /* 'Shuffle' rows down one, and read in */ - /* one new row. */ - for (wind_row = 0; wind_row < wsize - 1; wind_row++) - for (col = 0; col < ncols; col++) - *(row_in + (wind_row * ncols) + col) = - *(row_in + ((wind_row + 1) * ncols) + col); - } - - if (mparam != FEATURE) - Rast_set_d_null_value(row_out, ncols); - else - Rast_set_c_null_value(featrow_out, ncols); - for (wind_row = 0; wind_row < EDGE; wind_row++) { - if (mparam != FEATURE) - Rast_put_row(fd_out, row_out, - DCELL_TYPE); /* Write out the edge cells as NULL. */ - else - Rast_put_row(fd_out, featrow_out, - CELL_TYPE); /* Write out the edge cells as NULL. */ } - /*-----------------------------------------------------------------------*/ - /* FREE MEMORY USED TO STORE RASTER ROWS, LOCAL WINDOW AND MATRICES */ - - /*-----------------------------------------------------------------------*/ - - G_free(row_in); if (mparam != FEATURE) - G_free(row_out); + G_free(null_d); else - G_free(featrow_out); - - G_free(window_ptr); + G_free(null_c); + G_free(row_out); + G_free(featrow_out); + G_free(weight_ptr); free_dmatrix(normal_ptr, 0, 5, 0, 5); - free_dvector(obs_ptr, 0, 5); free_ivector(index_ptr, 0, 5); } From 64c5f01ae8a71ade7238d9dc6c3eb4da0ada42cd Mon Sep 17 00:00:00 2001 From: Kaushik Raja Date: Fri, 5 Jun 2026 01:44:47 -0700 Subject: [PATCH 02/13] r.param.scale: rework OpenMP parallelization to match r.neighbors MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This reworks my previous disk based parallel version so it follows the same structure r.neighbors uses. Instead of giving each thread a strip with the whole map, there's now an outer loop over a band of rows, and the work inside each band is split across the threads. Each thread keeps a small ring buffer of wsize rows and rotates it as it goes, reading a new row from disk per output row. So, a thread's input memory only depends on the window height, not the size of the map. Memory is now limited by the usual memory= option, which sets how many rows are in a band. A single shared output buffer for the band replaces each row having an allocation and row pointer arrays for the whole map. I also reworked find_obs() to precompute the window coordinates once and reuse w*z across the six weighted sums. This does less floating-point work per cell and the output is bit-identical to before. A few fixes that came out of the review: - When the window is taller than the region, the borders used to write more rows than the map has. The row counts are now clamped so exactly nrows rows get written. (This was a latent bug in the serial module too.) - Also brought back the nprocs < 1 error guard after Rast_disable_omp_on_mask, to match r.neighbors. - The each thread’s rings, window, and obs buffers are now allocated once before the band loop and indexed by thread id, instead of being reallocated on every band. Output is exactly identical to the serial module at one thread and with a raster mask (which forces a single thread). At more than one thread there's a pre-existing around 1 ULP difference from the OpenMP floating-point environment, tracked separately. --- raster/r.param.scale/close_down.c | 3 - raster/r.param.scale/find_normal.c | 27 ++- raster/r.param.scale/interface.c | 7 +- raster/r.param.scale/main.c | 6 +- raster/r.param.scale/open_files.c | 4 - raster/r.param.scale/param.h | 10 +- raster/r.param.scale/process.c | 360 +++++++++++++++++------------ 7 files changed, 246 insertions(+), 171 deletions(-) diff --git a/raster/r.param.scale/close_down.c b/raster/r.param.scale/close_down.c index 70e0f5af2f7..4f3a44e105e 100644 --- a/raster/r.param.scale/close_down.c +++ b/raster/r.param.scale/close_down.c @@ -21,9 +21,6 @@ void close_down(void) struct History history; char map_title[80], map_type[32]; - /* Close connection with existing input raster. */ - Rast_unopen(fd_in); - /* Write output raster map and close connection. */ Rast_close(fd_out); diff --git a/raster/r.param.scale/find_normal.c b/raster/r.param.scale/find_normal.c index 29b953cb0d8..812383d7dbb 100644 --- a/raster/r.param.scale/find_normal.c +++ b/raster/r.param.scale/find_normal.c @@ -103,25 +103,34 @@ void find_obs(DCELL *z, /* Local window of elevs. */ edge = EDGE, /* EDGE = (wsize-1)/2. */ offset; /* Array offset for weights & z */ - double x, y; /* Local window coordinates. */ + double xv[MAX_WSIZE], /* Precomputed window coordinates, */ + yv[MAX_WSIZE]; /* constant for every output cell. */ for (row = 0; row < 6; row++) /* Initialise column vector. */ obs[row] = 0.0; + for (col = 0; col < wsize; col++) /* x depends only on column. */ + xv[col] = resoln * (col - edge); + for (row = 0; row < wsize; row++) /* y depends only on row. */ + yv[row] = resoln * (row - edge); + for (row = 0; row < wsize; row++) for (col = 0; col < wsize; col++) { - x = resoln * (col - edge); - y = resoln * (row - edge); + /* w*z reused across all six sums; multiply order is kept */ + /* identical to the original so the result is unchanged. */ + double wz; + offset = row * wsize + col; + wz = *(w + offset) * *(z + offset); - obs[0] += *(w + offset) * *(z + offset) * x * x; - obs[1] += *(w + offset) * *(z + offset) * y * y; - obs[2] += *(w + offset) * *(z + offset) * x * y; - obs[3] += *(w + offset) * *(z + offset) * x; - obs[4] += *(w + offset) * *(z + offset) * y; + obs[0] += wz * xv[col] * xv[col]; + obs[1] += wz * yv[row] * yv[row]; + obs[2] += wz * xv[col] * yv[row]; + obs[3] += wz * xv[col]; + obs[4] += wz * yv[row]; if (!constrained) /* If constrained, should remain 0.0 */ - obs[5] += *(w + offset) * *(z + offset); + obs[5] += wz; } } diff --git a/raster/r.param.scale/interface.c b/raster/r.param.scale/interface.c index 4a923c3d81c..5f33831aacd 100644 --- a/raster/r.param.scale/interface.c +++ b/raster/r.param.scale/interface.c @@ -32,7 +32,7 @@ void interface(int argc, char **argv) *parameter, /* Morphometric parameter to calculate. */ *expon, /* Inverse distance exponent for weight. */ *vert_sc, /* Vertical scaling factor. */ - *nprocs_opt; + *nprocs_opt, *mem_opt; struct Flag *constr; /* Forces quadratic through the central */ @@ -69,6 +69,7 @@ void interface(int argc, char **argv) expon = G_define_option(); vert_sc = G_define_option(); nprocs_opt = G_define_standard_option(G_OPT_M_NPROCS); + mem_opt = G_define_standard_option(G_OPT_MEMORYMB); constr = G_define_flag(); @@ -133,6 +134,10 @@ void interface(int argc, char **argv) rast_out_name = rast_out->answer; /* can place the contents into strings */ wsize = atoi(win_size->answer); nprocs = G_set_omp_num_threads(nprocs_opt); + nprocs = Rast_disable_omp_on_mask(nprocs); + if (nprocs < 1) + G_fatal_error(_("<%d> is not valid number of nprocs."), nprocs); + memory = atoi(mem_opt->answer); constrained = constr->answer; sscanf(expon->answer, "%lf", &exponent); sscanf(vert_sc->answer, "%lf", &zscale); diff --git a/raster/r.param.scale/main.c b/raster/r.param.scale/main.c index 9b13b6d4e3a..29b8b6ca977 100644 --- a/raster/r.param.scale/main.c +++ b/raster/r.param.scale/main.c @@ -19,11 +19,11 @@ const char *rast_out_name; /* Name of the raster output file. */ int constrained; /* Flag that forces quadtratic through */ /* the central cell of the window. */ -int fd_in, /* File descriptor for input and */ - fd_out, /* output raster files. */ +int fd_out, /* File descriptor for output raster. */ wsize, /* Size of local processing window. */ mparam, /* Morphometric parameter to calculate. */ - nprocs; /* NUmber of threads in process. */ + nprocs, /* Number of threads in process. */ + memory; /* Memory cap in MB for row buffering. */ double resoln, /* Planimetric resolution. */ exponent, /* Distance weighting exponent. */ diff --git a/raster/r.param.scale/open_files.c b/raster/r.param.scale/open_files.c index ad4e6463c8c..efab56816b9 100644 --- a/raster/r.param.scale/open_files.c +++ b/raster/r.param.scale/open_files.c @@ -19,10 +19,6 @@ void open_files(void) { - /* Open existing file and set the input file descriptor. */ - - fd_in = Rast_open_old(rast_in_name, ""); - /* Open new file and set the output file descriptor. */ if (mparam != FEATURE) diff --git a/raster/r.param.scale/param.h b/raster/r.param.scale/param.h index a0d806b7031..4816e830c45 100644 --- a/raster/r.param.scale/param.h +++ b/raster/r.param.scale/param.h @@ -85,11 +85,11 @@ extern int constrained; /* Flag that forces quadtratic through */ /* the central cell of the window. */ -extern int fd_in, /* File descriptor for input and */ - fd_out, /* output raster files. */ - wsize, /* Size of local processing window. */ - mparam, /* Morphometric parameter to calculate. */ - nprocs; +extern int fd_out, /* File descriptor for output raster. */ + wsize, /* Size of local processing window. */ + mparam, /* Morphometric parameter to calculate. */ + nprocs, /* Number of threads in process. */ + memory; /* Memory cap in MB for row buffering. */ extern double resoln, /* Planimetric resolution. */ exponent, /* Distance weighting exponent. */ diff --git a/raster/r.param.scale/process.c b/raster/r.param.scale/process.c index a5d28a2b740..5a17e542c24 100644 --- a/raster/r.param.scale/process.c +++ b/raster/r.param.scale/process.c @@ -19,6 +19,17 @@ #include #endif +/* Rotate the per-thread ring: the oldest slot (0, smallest row index) is */ +/* recycled to become the newest slot (n-1). Mirrors r.neighbors rotate_bufs */ +/* (bufs.c:32-45). Pointer shuffle only, no data copy. */ +static void rotate_ring(DCELL **ring, int n) +{ + DCELL *temp = ring[0]; + for (int i = 1; i < n; i++) + ring[i - 1] = ring[i]; + ring[n - 1] = temp; +} + void process(void) { struct Cell_head region; /* Structure to hold region information */ @@ -69,139 +80,203 @@ void process(void) fd_thread[i] = Rast_open_old(rast_in_name, ""); /*-----------------------------------------------------------------------*/ - /* CHANGE 2: Output slots, one pointer per row. */ - /* Threads fill their rows into slots, a serial */ - /* loop writes them in order at the end (Rast_put_row is not */ - /* thread-safe and must be called in row order). */ + /* Band height from the memory cap. Reserve wsize input rows per thread */ + /* (the eventual ring-buffer cost), spend the rest on the band. */ + /*-----------------------------------------------------------------------*/ + + size_t in_buf_size = + (size_t)ncols * sizeof(DCELL) * (size_t)wsize * (size_t)nprocs; + size_t out_buf_size = (size_t)memory * (1 << 20); + if (out_buf_size <= in_buf_size) + out_buf_size = 0; + else + out_buf_size -= in_buf_size; + int brows = (int)(out_buf_size / (sizeof(DCELL) * (size_t)ncols)); + int proc_rows = nrows - 2 * EDGE; /* processable rows */ + if (brows > proc_rows) + brows = proc_rows; + if (brows < nprocs) + brows = nprocs; + if (brows < 1) + brows = 1; + + /*-----------------------------------------------------------------------*/ + /* CHANGE 2: One shared band buffer (brows rows), reused each band, plus */ + /* a single NULL row for the EDGE borders. Threads write disjoint */ + /* band-relative rows; a serial loop writes each band in row order. */ + /*-----------------------------------------------------------------------*/ + + DCELL *band_d = NULL, *null_d = NULL; + CELL *band_c = NULL, *null_c = NULL; + if (mparam != FEATURE) { + band_d = (DCELL *)G_malloc((size_t)brows * ncols * sizeof(DCELL)); + null_d = Rast_allocate_buf(DCELL_TYPE); + Rast_set_d_null_value(null_d, ncols); + } + else { + band_c = (CELL *)G_malloc((size_t)brows * ncols * sizeof(CELL)); + null_c = Rast_allocate_buf(CELL_TYPE); + Rast_set_c_null_value(null_c, ncols); + } + + /* Top EDGE border rows are written as NULL. top_border is clamped to nrows + */ + /* so a region with fewer than 2*EDGE+1 rows cannot overflow (see bottom). + */ + int top_border = (EDGE < nrows) ? EDGE : nrows; + for (row = 0; row < top_border; row++) { + if (mparam != FEATURE) + Rast_put_row(fd_out, null_d, DCELL_TYPE); + else + Rast_put_row(fd_out, null_c, CELL_TYPE); + } + + /*-----------------------------------------------------------------------*/ + /* Per-thread scratch, allocated ONCE before the band loop and indexed by */ + /* tid (like r.neighbors allocate_bufs), instead of churning the allocator + */ + /* on every band: a wsize-row ring, the window matrix, and the obs vector. + */ /*-----------------------------------------------------------------------*/ - DCELL **row_out = (DCELL **)G_malloc((size_t)nrows * sizeof(DCELL *)); - CELL **featrow_out = (CELL **)G_malloc((size_t)nrows * sizeof(CELL *)); - for (row = 0; row < nrows; row++) { - row_out[row] = NULL; - featrow_out[row] = NULL; + DCELL ***rings = (DCELL ***)G_malloc((size_t)nprocs * sizeof(DCELL **)); + DCELL **window_ptrs = (DCELL **)G_malloc((size_t)nprocs * sizeof(DCELL *)); + double **obs_ptrs = (double **)G_malloc((size_t)nprocs * sizeof(double *)); + for (int t = 0; t < nprocs; t++) { + rings[t] = (DCELL **)G_malloc((size_t)wsize * sizeof(DCELL *)); + for (int k = 0; k < wsize; k++) + rings[t][k] = (DCELL *)G_malloc((size_t)ncols * sizeof(DCELL)); + window_ptrs[t] = (DCELL *)G_malloc(SQR(wsize) * sizeof(DCELL)); + obs_ptrs[t] = dvector(0, 5); } /*-----------------------------------------------------------------------*/ - /* CHANGE 3: Parallel region. */ - /* Each thread claims its output row range, then reads its strip + halo */ - /* from disk via its own descriptor into a private buffer, and does */ - /* all window math reading from that private buffer. */ + /* CHANGE 3: Two-level structure. Outer serial loop over bands of brows */ + /* rows; inner parallel region splits the band across threads. Each */ + /* thread reads its sub-band + halo from disk into a private strip. */ /*-----------------------------------------------------------------------*/ + for (int band_start = EDGE; band_start < nrows - EDGE; + band_start += brows) { + int band_end = band_start + brows; + if (band_end > nrows - EDGE) + band_end = nrows - EDGE; /* exclusive */ + int band_count = band_end - band_start; + #pragma omp parallel if (nprocs > 1) - { - int nthreads = 1; - int tid = 0; + { + int nthreads = 1; + int tid = 0; #if defined(_OPENMP) - nthreads = omp_get_num_threads(); - tid = omp_get_thread_num(); + nthreads = omp_get_num_threads(); + tid = omp_get_thread_num(); #endif - /* Split processable rows [EDGE .. nrows-EDGE) evenly across threads */ - int proc_first = EDGE; - int proc_last = nrows - EDGE; /* exclusive */ - int proc_count = proc_last - proc_first; - - int my_start = proc_first + (int)((long)proc_count * tid / nthreads); - int my_end = - proc_first + (int)((long)proc_count * (tid + 1) / nthreads); - - if (my_end > my_start) { - /* Input rows this thread needs: my_start-EDGE .. my_end-1+EDGE */ - int buf_first = my_start - EDGE; - int buf_last = my_end - 1 + EDGE; - int buf_rows = buf_last - buf_first + 1; - - /* Allocate private strip and read it from disk via THIS */ - /* thread's own descriptor. No shared RAM map; each thread */ - /* pulls only what it needs. */ - DCELL *strip = - (DCELL *)G_malloc((size_t)buf_rows * ncols * sizeof(DCELL)); - for (int r = 0; r < buf_rows; r++) - Rast_get_d_row(fd_thread[tid], strip + (size_t)r * ncols, - buf_first + r); - - /* Per-thread math scratch */ - DCELL *window_ptr = (DCELL *)G_malloc(SQR(wsize) * sizeof(DCELL)); - double *obs_ptr = dvector(0, 5); - - for (int orow = my_start; orow < my_end; orow++) { - DCELL *out_d = NULL; - CELL *out_c = NULL; - if (mparam != FEATURE) { - out_d = Rast_allocate_buf(DCELL_TYPE); - Rast_set_d_null_value(out_d, ncols); + /* Split this band's rows [band_start .. band_end) across threads */ + int my_start = + band_start + (int)((long)band_count * tid / nthreads); + int my_end = + band_start + (int)((long)band_count * (tid + 1) / nthreads); + + if (my_end > my_start) { + /* Per-thread scratch, allocated once before the band loop. */ + DCELL **ring = rings[tid]; + DCELL *window_ptr = window_ptrs[tid]; + double *obs_ptr = obs_ptrs[tid]; + + /* Seed the ring with the first wsize-1 window rows; one new */ + /* row is read per output row below (read-once, mirrors */ + /* r.neighbors main.c:522-528). readrow = next disk row. */ + int readrow = my_start - EDGE; + for (int k = 0; k < wsize - 1; k++) { + rotate_ring(ring, wsize); + Rast_get_d_row(fd_thread[tid], ring[wsize - 1], readrow++); } - else { - out_c = Rast_allocate_buf(CELL_TYPE); - Rast_set_c_null_value(out_c, ncols); - } - - /* Position of this row inside the private strip */ - int srow = orow - buf_first; - for (int col = EDGE; col < (ncols - EDGE); col++) { - DCELL centre = *(strip + (size_t)srow * ncols + col); - if (Rast_is_d_null_value(¢re)) { - if (mparam == FEATURE) - Rast_set_c_null_value(out_c + col, 1); - else - Rast_set_d_null_value(out_d + col, 1); - continue; + for (int orow = my_start; orow < my_end; orow++) { + size_t boff = (size_t)(orow - band_start) * ncols; + DCELL *out_d = NULL; + CELL *out_c = NULL; + if (mparam != FEATURE) { + out_d = band_d + boff; + Rast_set_d_null_value(out_d, ncols); + } + else { + out_c = band_c + boff; + Rast_set_c_null_value(out_c, ncols); } - int found_null = FALSE; - for (int wind_row = 0; wind_row < wsize; wind_row++) { - if (found_null) - break; - for (int wind_col = 0; wind_col < wsize; wind_col++) { - DCELL *window_cell = - strip + - (size_t)(srow - EDGE + wind_row) * ncols + col + - wind_col - EDGE; - if (Rast_is_d_null_value(window_cell)) { - if (mparam == FEATURE) - Rast_set_c_null_value(out_c + col, 1); - else - Rast_set_d_null_value(out_d + col, 1); - found_null = TRUE; + /* Advance the ring by one row: rotate out the oldest, */ + /* read row orow+EDGE into the recycled newest slot; + * ring[k]*/ + /* then holds map row orow-EDGE+k (same layout as 3.1). */ + rotate_ring(ring, wsize); + Rast_get_d_row(fd_thread[tid], ring[wsize - 1], readrow++); + + for (int col = EDGE; col < (ncols - EDGE); col++) { + DCELL centre = *(ring[EDGE] + col); + if (Rast_is_d_null_value(¢re)) { + if (mparam == FEATURE) + Rast_set_c_null_value(out_c + col, 1); + else + Rast_set_d_null_value(out_d + col, 1); + continue; + } + + int found_null = FALSE; + for (int wind_row = 0; wind_row < wsize; wind_row++) { + if (found_null) break; + for (int wind_col = 0; wind_col < wsize; + wind_col++) { + DCELL *window_cell = + ring[wind_row] + col + wind_col - EDGE; + if (Rast_is_d_null_value(window_cell)) { + if (mparam == FEATURE) + Rast_set_c_null_value(out_c + col, 1); + else + Rast_set_d_null_value(out_d + col, 1); + found_null = TRUE; + break; + } + *(window_ptr + (wind_row * wsize) + wind_col) = + *(window_cell)-centre; } - *(window_ptr + (wind_row * wsize) + wind_col) = - *(window_cell)-centre; } + if (found_null) + continue; + + /* Math block: unchanged from serial */ + find_obs(window_ptr, obs_ptr, weight_ptr); + if (constrained) + G_lubksb(normal_ptr, 5, index_ptr, obs_ptr); + else + G_lubksb(normal_ptr, 6, index_ptr, obs_ptr); + + if (mparam == FEATURE) + *(out_c + col) = (CELL)feature(obs_ptr); + else + *(out_d + col) = param(mparam, obs_ptr); + + if (mparam == ELEV) + *(out_d + col) += centre; } - if (found_null) - continue; - - /* Math block: unchanged from serial */ - find_obs(window_ptr, obs_ptr, weight_ptr); - if (constrained) - G_lubksb(normal_ptr, 5, index_ptr, obs_ptr); - else - G_lubksb(normal_ptr, 6, index_ptr, obs_ptr); - - if (mparam == FEATURE) - *(out_c + col) = (CELL)feature(obs_ptr); - else - *(out_d + col) = param(mparam, obs_ptr); - - if (mparam == ELEV) - *(out_d + col) += centre; } - - if (mparam != FEATURE) - row_out[orow] = out_d; - else - featrow_out[orow] = out_c; } - - G_free(strip); - G_free(window_ptr); - free_dvector(obs_ptr, 0, 5); - } - } /* end parallel */ + } /* end parallel */ + + /* Serial in-order write of this band's computed rows. */ + if (mparam != FEATURE) + for (int orow = band_start; orow < band_end; orow++) + Rast_put_row(fd_out, + band_d + (size_t)(orow - band_start) * ncols, + DCELL_TYPE); + else + for (int orow = band_start; orow < band_end; orow++) + Rast_put_row(fd_out, + band_c + (size_t)(orow - band_start) * ncols, + CELL_TYPE); + } /* end band loop */ /*-----------------------------------------------------------------------*/ /* CHANGE 4: Close per-thread descriptors. */ @@ -212,47 +287,40 @@ void process(void) G_free(fd_thread); /*-----------------------------------------------------------------------*/ - /* CHANGE 5: Serial in-order write to disk. */ - /* Top and bottom EDGE rows are written as NULL, then the computed rows */ - /* are written from row_out[]/featrow_out[] in order. */ + /* CHANGE 5: Bottom EDGE border rows are written as NULL. The computed */ + /* rows were written per band inside the loop above. */ /*-----------------------------------------------------------------------*/ - DCELL *null_d = NULL; - CELL *null_c = NULL; + /* Bottom border: the remaining rows are NULL. With band_rows = max(0, */ + /* nrows-2*EDGE) computed rows, top_border + band_rows + bottom always sums + */ + /* to exactly nrows, even when the window is taller than the region. */ + int band_rows = (proc_rows > 0) ? proc_rows : 0; + for (row = top_border + band_rows; row < nrows; row++) { + if (mparam != FEATURE) + Rast_put_row(fd_out, null_d, DCELL_TYPE); + else + Rast_put_row(fd_out, null_c, CELL_TYPE); + } + if (mparam != FEATURE) { - null_d = Rast_allocate_buf(DCELL_TYPE); - Rast_set_d_null_value(null_d, ncols); + G_free(band_d); + G_free(null_d); } else { - null_c = Rast_allocate_buf(CELL_TYPE); - Rast_set_c_null_value(null_c, ncols); + G_free(band_c); + G_free(null_c); } - - for (row = 0; row < nrows; row++) { - if (row < EDGE || row >= nrows - EDGE) { - if (mparam != FEATURE) - Rast_put_row(fd_out, null_d, DCELL_TYPE); - else - Rast_put_row(fd_out, null_c, CELL_TYPE); - } - else { - if (mparam != FEATURE) { - Rast_put_row(fd_out, row_out[row], DCELL_TYPE); - G_free(row_out[row]); - } - else { - Rast_put_row(fd_out, featrow_out[row], CELL_TYPE); - G_free(featrow_out[row]); - } - } + for (int t = 0; t < nprocs; t++) { + for (int k = 0; k < wsize; k++) + G_free(rings[t][k]); + G_free(rings[t]); + G_free(window_ptrs[t]); + free_dvector(obs_ptrs[t], 0, 5); } - - if (mparam != FEATURE) - G_free(null_d); - else - G_free(null_c); - G_free(row_out); - G_free(featrow_out); + G_free(rings); + G_free(window_ptrs); + G_free(obs_ptrs); G_free(weight_ptr); free_dmatrix(normal_ptr, 0, 5, 0, 5); free_ivector(index_ptr, 0, 5); From 7b98052bff66739d0b271b2cb27effe978e56ed8 Mon Sep 17 00:00:00 2001 From: Kaushik Raja Date: Fri, 5 Jun 2026 11:04:28 -0700 Subject: [PATCH 03/13] r.param.scale: add benchmark and synthetic DEM generator scripts --- .../benchmark/benchmark_rparamscale.py | 167 ++++++++++++++++++ .../r.param.scale/benchmark/make_synth_dem.py | 59 +++++++ 2 files changed, 226 insertions(+) create mode 100644 raster/r.param.scale/benchmark/benchmark_rparamscale.py create mode 100644 raster/r.param.scale/benchmark/make_synth_dem.py diff --git a/raster/r.param.scale/benchmark/benchmark_rparamscale.py b/raster/r.param.scale/benchmark/benchmark_rparamscale.py new file mode 100644 index 00000000000..8272f64f5df --- /dev/null +++ b/raster/r.param.scale/benchmark/benchmark_rparamscale.py @@ -0,0 +1,167 @@ +#!/usr/bin/env python3 +# +# benchmark.py - parallel scaling benchmark for r.param.scale +# +# WHAT THIS MEASURES +# Wall-clock parallel SPEEDUP (median time at 1 thread / median time at N threads) +# for r.param.scale on a fixed input raster. For each window size it runs every +# thread count, takes the median of the timed runs (one warmup is discarded), and +# reports the speedup of each thread count against that size's own 1-thread run. +# Defaults: +# input = synth_dem_100m (generate it with the companion DEM script) +# method = elev +# memory = 300 MB +# sizes = [31] (window side; add more to sweep window size) +# nprocs = [1, 2, 4, 8] +# repeat = 4 (1 warmup discarded + 3 timed, median reported) +# A config whose timed runs spread more than 5% is flagged. Occasional flags are +# EXPECTED on machines with variable thermals or scheduling; a flag just means that +# config should be re-run, not that anything is wrong. +# +# The script creates one temporary output raster +# (tmp_benchmark_rparamscale_delete_me) and removes it at the end, even on error. +# +# RESULTS ARE HARDWARE-DEPENDENT. Absolute times and the speedup curve depend on +# core count and type (performance vs efficiency cores), memory bandwidth, cooling, +# and power source. Your numbers will differ from anyone else's. Use AC power. +# +# THREAD PLACEMENT ON HYBRID CPUs +# On CPUs with both performance and efficiency cores (e.g. Apple M-series), thread +# counts below the performance-core count can show large run-to-run variance +# because the scheduler places threads on P or E cores inconsistently. Pin the +# threads to cores to stabilize those configs: +# OMP_PROC_BIND=true OMP_PLACES=cores grass <...> --exec python3 benchmark.py +# +# RUN +# grass /path/to/grassdata// --exec python3 benchmark.py +# (the input raster must exist in that mapset) + +import io +import os +import platform +import statistics +import subprocess +from contextlib import redirect_stdout + +import grass.script as gs +from grass.pygrass.modules import Module +from grass.benchmark import benchmark_single + +# ---- configuration ------------------------------------------------------- +INPUT = "synth_dem_100m" +METHOD = "elev" +MEMORY = 300 # MB, passed to memory= +SIZES = [31] # window sides to test (e.g. [5, 15, 31, 51] to sweep) +NPROCS = [1, 2, 4, 8] # thread counts to test (1 is the baseline) +REPEAT = 4 # 1 warmup (discarded) + (REPEAT - 1) timed runs +SPREAD_FLAG = 5.0 # percent; flag a config whose timed runs spread above this +TMP_OUTPUT = "tmp_benchmark_rparamscale_delete_me" # created then removed +# -------------------------------------------------------------------------- + + +def keep_awake(): + """On macOS, hold an idle-sleep assertion for this process's lifetime.""" + if platform.system() == "Darwin": + try: + return subprocess.Popen(["caffeinate", "-i", "-w", str(os.getpid())]) + except FileNotFoundError: + pass + return None + + +def report_power(): + if platform.system() == "Darwin": + try: + ps = subprocess.check_output(["pmset", "-g", "ps"], text=True) + if "Battery Power" in ps: + print("Power: BATTERY -- WARNING: plug into AC, results will be noisy") + else: + print("Power: AC") + except Exception: + print("Power: could not determine") + else: + print("Power: not checked (non-macOS); ensure AC power and no power saving") + + +def timed_median(all_times): + """Discard the warmup (first run); return (median, spread_percent) of the rest.""" + timed = all_times[1:] + med = statistics.median(timed) + spread = (max(timed) - min(timed)) / med * 100.0 if med else 0.0 + return med, spread + + +def bench(size, nprocs): + mod = Module( + "r.param.scale", + input=INPUT, + output=TMP_OUTPUT, + method=METHOD, + size=size, + memory=MEMORY, + nprocs=nprocs, + overwrite=True, + quiet=True, + stderr_=subprocess.PIPE, + run_=False, + ) + # benchmark_single prints raw per-run times; silence it and keep only the median. + with redirect_stdout(io.StringIO()): + res = benchmark_single(mod, label="size%d n%d" % (size, nprocs), repeat=REPEAT) + return timed_median(res.all_times) + + +def main(): + awake = keep_awake() + try: + print("GRASS %s" % gs.version()["version"]) + print("Machine: %s; %d logical CPUs" % (platform.platform(), os.cpu_count())) + report_power() + print( + "Workload: %s, method=%s, memory=%d MB, sizes=%s, nprocs=%s, " + "1 warmup + %d timed runs (median)" + % (INPUT, METHOD, MEMORY, SIZES, NPROCS, REPEAT - 1) + ) + print( + "Pinning: OMP_PROC_BIND=%s OMP_PLACES=%s" + % ( + os.environ.get("OMP_PROC_BIND", "unset"), + os.environ.get("OMP_PLACES", "unset"), + ) + ) + gs.run_command("g.region", raster=INPUT) + + for size in SIZES: + results = {p: bench(size, p) for p in NPROCS} + base = results[NPROCS[0]][0] + print( + "\nr.param.scale size=%d (%s, memory=%d MB, baseline = %d thread)" + % (size, INPUT, MEMORY, NPROCS[0]) + ) + print(" nprocs | median time (s) | speedup") + print(" -------+-----------------+---------") + for p in NPROCS: + med, spread = results[p] + flag = " *" if spread > SPREAD_FLAG else "" + print(" %d | %11.3f | %5.2fx%s" % (p, med, base / med, flag)) + flagged = [str(p) for p in NPROCS if results[p][1] > SPREAD_FLAG] + if flagged: + print( + " * timed-run spread > %.0f%% (expected on thermally/" + "scheduling-variable machines; just re-run that config): " + "nprocs %s" % (SPREAD_FLAG, ", ".join(flagged)) + ) + finally: + # Always remove the temporary output raster, even on error. + try: + gs.run_command( + "g.remove", type="raster", name=TMP_OUTPUT, flags="f", quiet=True + ) + except Exception: + pass + if awake: + awake.terminate() + + +if __name__ == "__main__": + main() diff --git a/raster/r.param.scale/benchmark/make_synth_dem.py b/raster/r.param.scale/benchmark/make_synth_dem.py new file mode 100644 index 00000000000..4a6dcfa4cec --- /dev/null +++ b/raster/r.param.scale/benchmark/make_synth_dem.py @@ -0,0 +1,59 @@ +#!/usr/bin/env python3 +# +# make_synth_dem.py - generate the synthetic DEM used to benchmark r.param.scale +# +# WHAT THIS MAKES +# A synthetic terrain raster (DCELL) of ROWS x COLS cells, default 10000 x 10000 +# (100 million cells, about 800 MB), named "synth_dem_100m" -- the default input +# of the companion benchmark script. The surface is a deterministic sum of +# sinusoids at several frequencies, built with r.mapcalc from col()/row() only, so +# it is bit-reproducible on every machine and GRASS version (no random seed, no FFT). +# +# r.surf.fractal would give more realistic fractal terrain, but it depends on that +# module being built and on FFT-friendly dimensions; this analytic surface needs +# only r.mapcalc and works at any size. r.param.scale's run time depends on the +# raster size, not the cell values, so this is a valid and stable benchmark input. +# +# OVERWRITE +# This OVERWRITES any existing raster named OUTPUT (it passes --overwrite). Do not +# run it in a mapset where that name holds data you want to keep. +# +# RUN +# grass /path/to/grassdata// --exec python3 make_synth_dem.py +# (use a projected location; r.param.scale does not support lat/long) + +import grass.script as gs + +# ---- configuration ------------------------------------------------------- +ROWS = 10000 # lower ROWS/COLS for a quicker, smaller test (e.g. 2000 x 2000) +COLS = 10000 +OUTPUT = "synth_dem_100m" +# -------------------------------------------------------------------------- + +# Deterministic multi-frequency surface (a sum of sinusoids of col()/row(), with a +# +500 baseline so values read like positive elevations). DCELL because the constants +# and sin()/cos() are double precision. +EXPR = ( + "double(" + "500.0 + " + "256.0 * sin(col() * 0.9) * cos(row() * 0.9) + " + "128.0 * sin(col() * 2.1) * cos(row() * 1.7) + " + " 64.0 * sin(col() * 4.3) * cos(row() * 3.7) + " + " 32.0 * sin((col() + row()) * 6.1)" + ")" +) + + +def main(): + # Region: ROWS x COLS at resolution 1 (extent 0..COLS by 0..ROWS). + gs.run_command("g.region", n=ROWS, s=0, e=COLS, w=0, res=1, quiet=True) + gs.mapcalc("%s = %s" % (OUTPUT, EXPR), overwrite=True, quiet=True) + info = gs.raster_info(OUTPUT) + print( + "created %s: %d x %d cells (%s)" + % (OUTPUT, int(info["rows"]), int(info["cols"]), info["datatype"]) + ) + + +if __name__ == "__main__": + main() From add24cbf26f1da7f0b0f82dd8d8c87049420026e Mon Sep 17 00:00:00 2001 From: Kaushik Raja Date: Sun, 7 Jun 2026 14:04:54 -0700 Subject: [PATCH 04/13] r.param.scale: use DSQR to avoid float overflow in allocation size (CodeQL) --- raster/r.param.scale/process.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/raster/r.param.scale/process.c b/raster/r.param.scale/process.c index 5a17e542c24..792ba495917 100644 --- a/raster/r.param.scale/process.c +++ b/raster/r.param.scale/process.c @@ -55,7 +55,7 @@ void process(void) /* threads. */ /*-----------------------------------------------------------------------*/ - double *weight_ptr = (double *)G_malloc(SQR(wsize) * sizeof(double)); + double *weight_ptr = (double *)G_malloc(DSQR(wsize) * sizeof(double)); double **normal_ptr = dmatrix(0, 5, 0, 5); int *index_ptr = ivector(0, 5); double temp; @@ -146,7 +146,7 @@ void process(void) rings[t] = (DCELL **)G_malloc((size_t)wsize * sizeof(DCELL *)); for (int k = 0; k < wsize; k++) rings[t][k] = (DCELL *)G_malloc((size_t)ncols * sizeof(DCELL)); - window_ptrs[t] = (DCELL *)G_malloc(SQR(wsize) * sizeof(DCELL)); + window_ptrs[t] = (DCELL *)G_malloc(DSQR(wsize) * sizeof(DCELL)); obs_ptrs[t] = dvector(0, 5); } From 6b7a30e018b910d50e5a29a8dd33d14274c41406 Mon Sep 17 00:00:00 2001 From: Kaushik Raja Date: Tue, 9 Jun 2026 18:43:44 -0700 Subject: [PATCH 05/13] r.param.scale: run the LU decomposition single-threaded The LU solve happens once, before the parallel part starts. Running it on a single thread keeps the pivot choice the same every time, so the output now matches the serial version exactly at any thread count. Everything else still runs in parallel, so the speed doesn't change. --- raster/r.param.scale/process.c | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/raster/r.param.scale/process.c b/raster/r.param.scale/process.c index 792ba495917..f8a5e7432bf 100644 --- a/raster/r.param.scale/process.c +++ b/raster/r.param.scale/process.c @@ -63,10 +63,22 @@ void process(void) find_weight(weight_ptr); find_normal(normal_ptr, weight_ptr); + /* The one-time LU decomposition is run single-threaded because G_ludcmp's + * internal pivot search is not thread-safe (it has an unsynchronized shared + * reduction in lib/gmath/lu.c). Parallelizing this tiny one-time 6x6 solve + * gains nothing, and at nprocs>1 the race shifts the factorization by about + * 1 ULP, which then perturbs every per-cell solve. */ +#ifdef _OPENMP + int saved_threads = omp_get_max_threads(); + omp_set_num_threads(1); +#endif if (constrained) G_ludcmp(normal_ptr, 5, index_ptr, &temp); else G_ludcmp(normal_ptr, 6, index_ptr, &temp); +#ifdef _OPENMP + omp_set_num_threads(saved_threads); +#endif /*-----------------------------------------------------------------------*/ /* CHANGE 1: Per-thread input file descriptors. */ From 61f515d0d1cfbb5124ec903b13072d1f74efcd52 Mon Sep 17 00:00:00 2001 From: Kaushik Raja Date: Wed, 10 Jun 2026 21:21:05 -0700 Subject: [PATCH 06/13] r.param.scale: add parallel correctness testsuite gunittest test on a generated deterministic DEM. For each method, window size and the -c flag it compares nprocs=1 against nprocs=4 with a null-aware bit-exact check and against single-thread reference values. Covers the multi-band disk-strip path and the single-band path, asserts output is invariant to band count, and includes a masked case. --- .../testsuite/test_r_param_scale.py | 671 ++++++++++++++++++ 1 file changed, 671 insertions(+) create mode 100644 raster/r.param.scale/testsuite/test_r_param_scale.py diff --git a/raster/r.param.scale/testsuite/test_r_param_scale.py b/raster/r.param.scale/testsuite/test_r_param_scale.py new file mode 100644 index 00000000000..ba5bf948bfd --- /dev/null +++ b/raster/r.param.scale/testsuite/test_r_param_scale.py @@ -0,0 +1,671 @@ +"""Test of r.param.scale parallelization (disk-strip path). + +Verifies that r.param.scale produces the same output when run threaded +(nprocs=4) as it does serially (nprocs=1), and locks in reference +univariate statistics for a sweep of methods, window sizes and the +constrained (-c) flag. + +Self-contained: the input DEM is generated in setUpClass with r.mapcalc +from a deterministic polynomial in row()/col() (pure +, -, * with decimal +constants, no trig and no rand), so it is bit-identical on every platform +and GRASS version. A small cubic term makes local curvature vary across +the raster while keeping the fit well conditioned. + +The region (400 x 400) together with memory=1 forces the module past its +memory budget, so it uses the multi-band per-thread file-descriptor +disk-strip path (brows < proc_rows) rather than holding the map in RAM. +The same memory setting is used for the serial and threaded runs. + +The reference values below were generated from a genuine single-threaded +run (OMP_NUM_THREADS=1, nprocs=1): the one-time LU decomposition relies on +a library pivot search that is only deterministic on a single thread, so +the reference must come from a single-threaded run. + +Used dataset: nc_spm_08_grass7 (only as a host location; the DEM is +generated, the NC rasters are not used). + +Author: Kaushik Raja +""" + +from grass.gunittest.case import TestCase +import grass.script as gs + + +REFERENCE = { + "elev_s3": { + "n": 158404, + "null_cells": 1596, + "cells": 160000, + "min": -564.762512, + "max": 1491.000666, + "range": 2055.763178, + "mean": 426.151076, + "mean_of_abs": 480.263031516931, + "stddev": 442.861043442243, + "variance": 196125.903798752, + "coeff_var": 103.921136982473, + "sum": 67504035.042704, + }, + "elev_s3_c": { + "n": 158404, + "null_cells": 1596, + "cells": 160000, + "min": -564.762512, + "max": 1491.000666, + "range": 2055.763178, + "mean": 426.151076, + "mean_of_abs": 480.263031516931, + "stddev": 442.861043442243, + "variance": 196125.903798752, + "coeff_var": 103.921136982473, + "sum": 67504035.042704, + }, + "elev_s9": { + "n": 153664, + "null_cells": 6336, + "cells": 160000, + "min": -546.478838, + "max": 1467.508536, + "range": 2013.987374, + "mean": 425.280891, + "mean_of_abs": 476.390592098136, + "stddev": 434.633318407327, + "variance": 188906.121469765, + "coeff_var": 102.199117713786, + "sum": 65350362.834624, + }, + "elev_s9_c": { + "n": 153664, + "null_cells": 6336, + "cells": 160000, + "min": -546.478838, + "max": 1467.508536, + "range": 2013.987374, + "mean": 425.280891, + "mean_of_abs": 476.390592098136, + "stddev": 434.633318407327, + "variance": 188906.121469765, + "coeff_var": 102.199117713786, + "sum": 65350362.834624, + }, + "elev_s15": { + "n": 148996, + "null_cells": 11004, + "cells": 160000, + "min": -528.372536, + "max": 1444.287922, + "range": 1972.660458, + "mean": 424.423924, + "mean_of_abs": 472.603709161615, + "stddev": 426.473749373164, + "variance": 181879.858904405, + "coeff_var": 100.482966500532, + "sum": 63237466.980304, + }, + "elev_s15_c": { + "n": 148996, + "null_cells": 11004, + "cells": 160000, + "min": -528.372536, + "max": 1444.287922, + "range": 1972.660458, + "mean": 424.423924, + "mean_of_abs": 472.603709161615, + "stddev": 426.473749373164, + "variance": 181879.858904405, + "coeff_var": 100.482966500532, + "sum": 63237466.980304, + }, + "slope_s3": { + "n": 158404, + "null_cells": 1596, + "cells": 160000, + "min": 31.143096450782, + "max": 83.9114645026956, + "range": 52.7683680519136, + "mean": 72.4250151761828, + "mean_of_abs": 72.4250151761828, + "stddev": 10.9080822146547, + "variance": 118.986257601667, + "coeff_var": 15.0612080482406, + "sum": 11472412.1039681, + }, + "slope_s3_c": { + "n": 158404, + "null_cells": 1596, + "cells": 160000, + "min": 31.143096450782, + "max": 83.9114645026956, + "range": 52.7683680519136, + "mean": 72.4250151761828, + "mean_of_abs": 72.4250151761828, + "stddev": 10.9080822146547, + "variance": 118.986257601667, + "coeff_var": 15.0612080482406, + "sum": 11472412.1039681, + }, + "slope_s9": { + "n": 153664, + "null_cells": 6336, + "cells": 160000, + "min": 31.944182807223, + "max": 83.8375409801301, + "range": 51.8933581729071, + "mean": 72.457518347825, + "mean_of_abs": 72.457518347825, + "stddev": 10.7142059615324, + "variance": 114.794209386136, + "coeff_var": 14.7868795479579, + "sum": 11134112.0994002, + }, + "slope_s9_c": { + "n": 153664, + "null_cells": 6336, + "cells": 160000, + "min": 31.944182807223, + "max": 83.8375409801301, + "range": 51.8933581729071, + "mean": 72.457518347825, + "mean_of_abs": 72.457518347825, + "stddev": 10.7142059615324, + "variance": 114.794209386136, + "coeff_var": 14.7868795479579, + "sum": 11134112.0994002, + }, + "slope_s15": { + "n": 148996, + "null_cells": 11004, + "cells": 160000, + "min": 32.7563367059122, + "max": 83.7623712622949, + "range": 51.0060345563827, + "mean": 72.4906681285151, + "mean_of_abs": 72.4906681285151, + "stddev": 10.5190698109164, + "variance": 110.650829686932, + "coeff_var": 14.5109295892647, + "sum": 10800819.5884762, + }, + "slope_s15_c": { + "n": 148996, + "null_cells": 11004, + "cells": 160000, + "min": 32.7563367059122, + "max": 83.7623712622949, + "range": 51.0060345563827, + "mean": 72.4906681285151, + "mean_of_abs": 72.4906681285151, + "stddev": 10.5190698109164, + "variance": 110.650829686932, + "coeff_var": 14.5109295892647, + "sum": 10800819.5884762, + }, + "profc_s3": { + "n": 158404, + "null_cells": 1596, + "cells": 160000, + "min": -0.00393749267654076, + "max": 0.00164711282024063, + "range": 0.00558460549678139, + "mean": -0.000342552958886785, + "mean_of_abs": 0.000491699464986382, + "stddev": 0.000829231153728635, + "variance": 6.87624306314123e-07, + "coeff_var": -242.073855214516, + "sum": -54.2617588995023, + }, + "profc_s3_c": { + "n": 158404, + "null_cells": 1596, + "cells": 160000, + "min": -0.00393749267654189, + "max": 0.00164711282024106, + "range": 0.00558460549678296, + "mean": -0.000342552958886785, + "mean_of_abs": 0.000491699464986382, + "stddev": 0.000829231153728635, + "variance": 6.87624306314123e-07, + "coeff_var": -242.073855214516, + "sum": -54.2617588995023, + }, + "profc_s9": { + "n": 153664, + "null_cells": 6336, + "cells": 160000, + "min": -0.00393703319719625, + "max": 0.00156307905481015, + "range": 0.0055001122520064, + "mean": -0.000343203532196393, + "mean_of_abs": 0.000487504566970183, + "stddev": 0.000819692530679842, + "variance": 6.71895844852323e-07, + "coeff_var": -238.835691880579, + "sum": -52.7380275714266, + }, + "profc_s9_c": { + "n": 153664, + "null_cells": 6336, + "cells": 160000, + "min": -0.00393703319719643, + "max": 0.00156307905480931, + "range": 0.00550011225200574, + "mean": -0.000343203532196394, + "mean_of_abs": 0.000487504566970183, + "stddev": 0.000819692530679842, + "variance": 6.71895844852324e-07, + "coeff_var": -238.835691880579, + "sum": -52.7380275714266, + }, + "profc_s15": { + "n": 148996, + "null_cells": 11004, + "cells": 160000, + "min": -0.00393547758337053, + "max": 0.0014817471036016, + "range": 0.00541722468697213, + "mean": -0.000343511486125908, + "mean_of_abs": 0.000483002036366709, + "stddev": 0.000809295470184593, + "variance": 6.54959158061301e-07, + "coeff_var": -235.594879027701, + "sum": -51.1818373868158, + }, + "profc_s15_c": { + "n": 148996, + "null_cells": 11004, + "cells": 160000, + "min": -0.00393547758337018, + "max": 0.00148174710360166, + "range": 0.00541722468697183, + "mean": -0.000343511486125908, + "mean_of_abs": 0.000483002036366709, + "stddev": 0.000809295470184592, + "variance": 6.54959158061301e-07, + "coeff_var": -235.594879027701, + "sum": -51.1818373868158, + }, + "crosc_s3": { + "n": 158404, + "null_cells": 1596, + "cells": 160000, + "min": -0.00798326486875052, + "max": 0.0140198635034536, + "range": 0.0220031283722041, + "mean": 0.00523063372892867, + "mean_of_abs": 0.00700064364287339, + "stddev": 0.0057384315690963, + "variance": 3.2929596873201e-05, + "coeff_var": 109.708151372924, + "sum": 828.553305197217, + }, + "crosc_s3_c": { + "n": 158404, + "null_cells": 1596, + "cells": 160000, + "min": -0.00798326486872778, + "max": 0.0140198635036052, + "range": 0.0220031283723329, + "mean": 0.00523063372892867, + "mean_of_abs": 0.00700064364287339, + "stddev": 0.0057384315690963, + "variance": 3.2929596873201e-05, + "coeff_var": 109.708151372924, + "sum": 828.553305197217, + }, + "crosc_s9": { + "n": 153664, + "null_cells": 6336, + "cells": 160000, + "min": -0.00790623053867785, + "max": 0.0139600812849915, + "range": 0.0218663118236694, + "mean": 0.00530129926236255, + "mean_of_abs": 0.00703005215756205, + "stddev": 0.00570183973099957, + "variance": 3.25109763180052e-05, + "coeff_var": 107.555515144763, + "sum": 814.618849851679, + }, + "crosc_s9_c": { + "n": 153664, + "null_cells": 6336, + "cells": 160000, + "min": -0.00790623053867993, + "max": 0.0139600812849679, + "range": 0.0218663118236478, + "mean": 0.00530129926236254, + "mean_of_abs": 0.00703005215756204, + "stddev": 0.00570183973099955, + "variance": 3.25109763180051e-05, + "coeff_var": 107.555515144763, + "sum": 814.618849851678, + }, + "crosc_s15": { + "n": 148996, + "null_cells": 11004, + "cells": 160000, + "min": -0.00782884301345861, + "max": 0.0139001591796006, + "range": 0.0217290021930592, + "mean": 0.00537273496983023, + "mean_of_abs": 0.0070588784164276, + "stddev": 0.00566360869239518, + "variance": 3.20764634205742e-05, + "coeff_var": 105.413885557324, + "sum": 800.516019564825, + }, + "crosc_s15_c": { + "n": 148996, + "null_cells": 11004, + "cells": 160000, + "min": -0.00782884301345516, + "max": 0.0139001591796018, + "range": 0.021729002193057, + "mean": 0.00537273496983024, + "mean_of_abs": 0.00705887841642761, + "stddev": 0.00566360869239518, + "variance": 3.20764634205743e-05, + "coeff_var": 105.413885557324, + "sum": 800.516019564826, + }, + "feature_s3": { + "n": 158404, + "null_cells": 1596, + "cells": 160000, + "min": 1.0, + "max": 5.0, + "range": 4.0, + "mean": 4.56869776015757, + "mean_of_abs": 4.56869776015757, + "stddev": 0.841755158740511, + "variance": 0.708551747266264, + "coeff_var": 18.4244001886323, + "sum": 723700.0, + }, + "feature_s3_c": { + "n": 158404, + "null_cells": 1596, + "cells": 160000, + "min": 1.0, + "max": 5.0, + "range": 4.0, + "mean": 4.56869776015757, + "mean_of_abs": 4.56869776015757, + "stddev": 0.841755158740511, + "variance": 0.708551747266264, + "coeff_var": 18.4244001886323, + "sum": 723700.0, + }, + "feature_s9": { + "n": 153664, + "null_cells": 6336, + "cells": 160000, + "min": 1.0, + "max": 5.0, + "range": 4.0, + "mean": 4.57504685547688, + "mean_of_abs": 4.57504685547688, + "stddev": 0.837584186678913, + "variance": 0.701547269774577, + "coeff_var": 18.3076635745539, + "sum": 703020.0, + }, + "feature_s9_c": { + "n": 153664, + "null_cells": 6336, + "cells": 160000, + "min": 1.0, + "max": 5.0, + "range": 4.0, + "mean": 4.57504685547688, + "mean_of_abs": 4.57504685547688, + "stddev": 0.837584186678913, + "variance": 0.701547269774577, + "coeff_var": 18.3076635745539, + "sum": 703020.0, + }, + "feature_s15": { + "n": 148996, + "null_cells": 11004, + "cells": 160000, + "min": 1.0, + "max": 5.0, + "range": 4.0, + "mean": 4.58155923648957, + "mean_of_abs": 4.58155923648957, + "stddev": 0.833230429104926, + "variance": 0.69427294798638, + "coeff_var": 18.1866125939988, + "sum": 682634.0, + }, + "feature_s15_c": { + "n": 148996, + "null_cells": 11004, + "cells": 160000, + "min": 1.0, + "max": 5.0, + "range": 4.0, + "mean": 4.58155923648957, + "mean_of_abs": 4.58155923648957, + "stddev": 0.833230429104926, + "variance": 0.69427294798638, + "coeff_var": 18.1866125939988, + "sum": 682634.0, + }, +} + + +class TestParamScale(TestCase): + """Parallel-vs-serial and reference-value regression for r.param.scale + on a generated DEM, exercising the disk-strip (multi-band) path. + """ + + # Generated input DEM and the fixed, reproducible region. + input_dem = "rps_test_dem" + region = {"n": 400, "s": 0, "e": 400, "w": 0, "res": 1} + dem_expr = ( + "double(100.0 + 0.5*col() + 0.4*row() " + "+ 0.003*col()*col() - 0.002*row()*row() + 0.0015*col()*row() " + "+ 0.00001*col()*col()*col() - 0.000008*row()*row()*row())" + ) + # Low memory budget forces the multi-band disk-strip path; identical for + # the serial and threaded runs so the comparison stays valid. + memory = 1 + + test_options = REFERENCE + to_remove = [] + + @classmethod + def setUpClass(cls): + cls.use_temp_region() + cls.runModule("g.region", **cls.region) + cls.runModule( + "r.mapcalc", + expression="{} = {}".format(cls.input_dem, cls.dem_expr), + overwrite=True, + ) + cls.to_remove.append(cls.input_dem) + + @classmethod + def tearDownClass(cls): + cls.del_temp_region() + if cls.to_remove: + cls.runModule( + "g.remove", flags="f", type="raster", name=",".join(cls.to_remove) + ) + + @staticmethod + def _key(method, size, constrained): + return "{}_s{}{}".format(method, size, "_c" if constrained else "") + + def _assert_identical(self, a, b, label): + """Null-aware bit-exact equality of two rasters. + + The diff is 1 wherever the null pattern differs OR (both non-null) + the values differ, and 0 otherwise. It is never null, so r.univar + does not skip cells. The differing-cell sum must be exactly 0, and + the two rasters must report the same non-null count. + """ + diff = "rps_diff_{}".format(label) + self.to_remove.append(diff) + expr = ( + "{d} = if(isnull({a}) != isnull({b}), 1, " + "if(isnull({a}), 0, if({a} != {b}, 1, 0)))" + ).format(d=diff, a=a, b=b) + self.runModule("r.mapcalc", expression=expr, overwrite=True) + dstats = gs.parse_command("r.univar", map=diff, flags="g") + differing = int(float(dstats.get("sum", 0) or 0)) + self.assertEqual( + differing, + 0, + msg="{} and {} differ in {} cells ({})".format(a, b, differing, label), + ) + a_n = int(gs.parse_command("r.univar", map=a, flags="g")["n"]) + b_n = int(gs.parse_command("r.univar", map=b, flags="g")["n"]) + self.assertEqual( + a_n, + b_n, + msg="non-null count mismatch ({} vs {}) for {}".format(a_n, b_n, label), + ) + + def _check_combo(self, method, size, constrained): + key = self._key(method, size, constrained) + serial = "rps_{}_serial".format(key) + threaded = "rps_{}_threaded".format(key) + self.to_remove.extend([serial, threaded]) + + common = { + "input": self.input_dem, + "method": method, + "size": size, + "memory": self.memory, + "overwrite": True, + } + serial_kw = dict(common, output=serial, nprocs=1) + threaded_kw = dict(common, output=threaded, nprocs=4) + if constrained: + serial_kw["flags"] = "c" + threaded_kw["flags"] = "c" + + self.assertModule("r.param.scale", **serial_kw) + self.assertModule("r.param.scale", **threaded_kw) + + # (a) reference-value regression guard (includes cell/null counts). + reference = self.test_options[key] + self.assertRasterFitsUnivar(raster=serial, reference=reference, precision=1e-5) + self.assertRasterFitsUnivar( + raster=threaded, reference=reference, precision=1e-5 + ) + # (b) null-aware bit-exact serial vs threaded. + self._assert_identical(serial, threaded, key) + + def _check_method(self, method): + for size in (3, 9, 15): + for constrained in (False, True): + self._check_combo(method, size, constrained) + + def test_elev(self): + """Method elev at sizes 3, 9, 15, with and without -c.""" + self._check_method("elev") + + def test_slope(self): + """Method slope at sizes 3, 9, 15, with and without -c.""" + self._check_method("slope") + + def test_profc(self): + """Method profc (profile curvature) at sizes 3, 9, 15, +/- -c.""" + self._check_method("profc") + + def test_crosc(self): + """Method crosc (cross-sectional curvature) at sizes 3, 9, 15, +/- -c.""" + self._check_method("crosc") + + def test_feature(self): + """Method feature (morphometric features) at sizes 3, 9, 15, +/- -c.""" + self._check_method("feature") + + def test_single_band_and_memory_invariance(self): + """Single-band path (memory=300) and invariance to band count. + + At size 15 the whole 400x400 map fits in memory=300 (one band), + whereas memory=1 forces multiple bands. For every method (default + and -c): the single-band serial and threaded outputs must be + identical, and the single-band serial output must be identical to + the multi-band (memory=1) serial output, i.e. the result is + invariant to how many bands the module uses. + """ + for method in ("elev", "slope", "profc", "crosc", "feature"): + for constrained in (False, True): + tag = "{}{}".format(method, "_c" if constrained else "") + s_m1 = "rps_inv_{}_s_m1".format(tag) + s_m300 = "rps_inv_{}_s_m300".format(tag) + t_m300 = "rps_inv_{}_t_m300".format(tag) + self.to_remove.extend([s_m1, s_m300, t_m300]) + + base = { + "input": self.input_dem, + "method": method, + "size": 15, + "overwrite": True, + } + runs = [ + dict(base, output=s_m1, nprocs=1, memory=1), + dict(base, output=s_m300, nprocs=1, memory=300), + dict(base, output=t_m300, nprocs=4, memory=300), + ] + if constrained: + for r in runs: + r["flags"] = "c" + for r in runs: + self.assertModule("r.param.scale", **r) + + # single-band serial == single-band threaded + self._assert_identical(s_m300, t_m300, "inv_{}_par".format(tag)) + # band-count invariance: one band (mem300) == multi band (mem1) + self._assert_identical(s_m300, s_m1, "inv_{}_mem".format(tag)) + # lock single-band output to the stored size-15 reference + reference = self.test_options[self._key(method, 15, constrained)] + self.assertRasterFitsUnivar( + raster=s_m300, reference=reference, precision=1e-5 + ) + self.assertRasterFitsUnivar( + raster=t_m300, reference=reference, precision=1e-5 + ) + + def test_mask(self): + """Masked sub-region: serial and threaded must match with a mask + active (exercises Rast_disable_omp_on_mask); masked cells are null + in both. Method elev, size 9, memory=1. + """ + maskmap = "rps_maskmap" + serial = "rps_mask_serial" + threaded = "rps_mask_threaded" + self.to_remove.extend([maskmap, serial, threaded]) + self.runModule( + "r.mapcalc", + expression=( + "{} = if(row() > 100 && row() < 300 " + "&& col() > 100 && col() < 300, 1, null())".format(maskmap) + ), + overwrite=True, + ) + self.runModule("r.mask", raster=maskmap) + try: + base = { + "input": self.input_dem, + "method": "elev", + "size": 9, + "overwrite": True, + "memory": 1, + } + self.assertModule("r.param.scale", output=serial, nprocs=1, **base) + self.assertModule("r.param.scale", output=threaded, nprocs=4, **base) + self._assert_identical(serial, threaded, "mask") + finally: + self.runModule("r.mask", flags="r") + + +if __name__ == "__main__": + from grass.gunittest.main import test + + test() From 6f6db7fbbdcc21e060e1b230ec4f038d112a8fac Mon Sep 17 00:00:00 2001 From: Kaushik Raja Date: Thu, 11 Jun 2026 14:27:33 -0700 Subject: [PATCH 07/13] r.param.scale: enable OpenMP in the CMake build --- raster/CMakeLists.txt | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/raster/CMakeLists.txt b/raster/CMakeLists.txt index 82efc265aa6..37c66e08575 100644 --- a/raster/CMakeLists.txt +++ b/raster/CMakeLists.txt @@ -413,7 +413,9 @@ build_program_in_subdir( grass_interpfl grass_interpdata grass_qtree - ${LIBM}) + ${LIBM} + OPTIONAL_DEPENDS + OpenMP::OpenMP_C) build_program_in_subdir(r.patch DEPENDS grass_gis grass_raster OPTIONAL_DEPENDS OpenMP::OpenMP_C) From 8f025aa4f041d868c1995636adda57c3851e53f4 Mon Sep 17 00:00:00 2001 From: Kaushik Raja Date: Thu, 11 Jun 2026 14:27:51 -0700 Subject: [PATCH 08/13] r.param.scale: use full executable paths in benchmark script --- raster/r.param.scale/benchmark/benchmark_rparamscale.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/raster/r.param.scale/benchmark/benchmark_rparamscale.py b/raster/r.param.scale/benchmark/benchmark_rparamscale.py index 8272f64f5df..576310453bf 100644 --- a/raster/r.param.scale/benchmark/benchmark_rparamscale.py +++ b/raster/r.param.scale/benchmark/benchmark_rparamscale.py @@ -63,7 +63,9 @@ def keep_awake(): """On macOS, hold an idle-sleep assertion for this process's lifetime.""" if platform.system() == "Darwin": try: - return subprocess.Popen(["caffeinate", "-i", "-w", str(os.getpid())]) + return subprocess.Popen( + ["/usr/bin/caffeinate", "-i", "-w", str(os.getpid())] + ) except FileNotFoundError: pass return None @@ -72,7 +74,7 @@ def keep_awake(): def report_power(): if platform.system() == "Darwin": try: - ps = subprocess.check_output(["pmset", "-g", "ps"], text=True) + ps = subprocess.check_output(["/usr/bin/pmset", "-g", "ps"], text=True) if "Battery Power" in ps: print("Power: BATTERY -- WARNING: plug into AC, results will be noisy") else: From 5661e31b2d690e9b73a00f27b1ec525ba545bb89 Mon Sep 17 00:00:00 2001 From: Kaushik Raja Date: Mon, 15 Jun 2026 17:50:55 -0700 Subject: [PATCH 09/13] r.param.scale: simplify allocations and linearize ring buffers Simplified G_malloc calls to the sizeof(*ptr) form, reorder in_buf_size so the multiplication happens in size_t, document the MiB conversion, and replaced the per-thread per-row ring allocations with a single contiguous block indexed by per-thread pointer tables. --- raster/r.param.scale/process.c | 27 +++++++++++++++------------ 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/raster/r.param.scale/process.c b/raster/r.param.scale/process.c index f8a5e7432bf..b466426a40c 100644 --- a/raster/r.param.scale/process.c +++ b/raster/r.param.scale/process.c @@ -87,7 +87,7 @@ void process(void) /* This is the r.neighbors pattern. No full_map array exists. */ /*-----------------------------------------------------------------------*/ - int *fd_thread = (int *)G_malloc((size_t)nprocs * sizeof(int)); + int *fd_thread = G_malloc(sizeof(*fd_thread) * nprocs); for (int i = 0; i < nprocs; i++) fd_thread[i] = Rast_open_old(rast_in_name, ""); @@ -96,8 +96,8 @@ void process(void) /* (the eventual ring-buffer cost), spend the rest on the band. */ /*-----------------------------------------------------------------------*/ - size_t in_buf_size = - (size_t)ncols * sizeof(DCELL) * (size_t)wsize * (size_t)nprocs; + size_t in_buf_size = sizeof(DCELL) * ncols * wsize * nprocs; + /* memory= is in MiB; convert to bytes */ size_t out_buf_size = (size_t)memory * (1 << 20); if (out_buf_size <= in_buf_size) out_buf_size = 0; @@ -121,12 +121,12 @@ void process(void) DCELL *band_d = NULL, *null_d = NULL; CELL *band_c = NULL, *null_c = NULL; if (mparam != FEATURE) { - band_d = (DCELL *)G_malloc((size_t)brows * ncols * sizeof(DCELL)); + band_d = G_malloc(sizeof(*band_d) * brows * ncols); null_d = Rast_allocate_buf(DCELL_TYPE); Rast_set_d_null_value(null_d, ncols); } else { - band_c = (CELL *)G_malloc((size_t)brows * ncols * sizeof(CELL)); + band_c = G_malloc(sizeof(*band_c) * brows * ncols); null_c = Rast_allocate_buf(CELL_TYPE); Rast_set_c_null_value(null_c, ncols); } @@ -151,13 +151,17 @@ void process(void) */ /*-----------------------------------------------------------------------*/ - DCELL ***rings = (DCELL ***)G_malloc((size_t)nprocs * sizeof(DCELL **)); - DCELL **window_ptrs = (DCELL **)G_malloc((size_t)nprocs * sizeof(DCELL *)); - double **obs_ptrs = (double **)G_malloc((size_t)nprocs * sizeof(double *)); + DCELL ***rings = G_malloc(sizeof(*rings) * nprocs); + DCELL **window_ptrs = G_malloc(sizeof(*window_ptrs) * nprocs); + double **obs_ptrs = G_malloc(sizeof(*obs_ptrs) * nprocs); + /* All ring rows live in one contiguous block; each thread's pointer */ + /* table (rings[t][k]) indexes into it, so rotate_ring only shuffles */ + /* pointers and there are no per-row allocations. */ + DCELL *ring_block = G_malloc(sizeof(*ring_block) * nprocs * wsize * ncols); for (int t = 0; t < nprocs; t++) { - rings[t] = (DCELL **)G_malloc((size_t)wsize * sizeof(DCELL *)); + rings[t] = G_malloc(sizeof(*rings[t]) * wsize); for (int k = 0; k < wsize; k++) - rings[t][k] = (DCELL *)G_malloc((size_t)ncols * sizeof(DCELL)); + rings[t][k] = ring_block + (size_t)(t * wsize + k) * ncols; window_ptrs[t] = (DCELL *)G_malloc(DSQR(wsize) * sizeof(DCELL)); obs_ptrs[t] = dvector(0, 5); } @@ -324,12 +328,11 @@ void process(void) G_free(null_c); } for (int t = 0; t < nprocs; t++) { - for (int k = 0; k < wsize; k++) - G_free(rings[t][k]); G_free(rings[t]); G_free(window_ptrs[t]); free_dvector(obs_ptrs[t], 0, 5); } + G_free(ring_block); G_free(rings); G_free(window_ptrs); G_free(obs_ptrs); From f9a27d83bccfbac8b775fe16f2f72dfd97e6a050 Mon Sep 17 00:00:00 2001 From: Kaushik Raja Date: Wed, 17 Jun 2026 00:20:34 -0700 Subject: [PATCH 10/13] r.param.scale: rework benchmark to nprocs convention and document parallelization --- .../benchmark_r_param_scale_nprocs.py | 57 ++++++ .../benchmark/benchmark_rparamscale.py | 169 ------------------ .../r.param.scale/benchmark/make_synth_dem.py | 59 ------ raster/r.param.scale/r.param.scale.md | 11 ++ 4 files changed, 68 insertions(+), 228 deletions(-) create mode 100644 raster/r.param.scale/benchmark/benchmark_r_param_scale_nprocs.py delete mode 100644 raster/r.param.scale/benchmark/benchmark_rparamscale.py delete mode 100644 raster/r.param.scale/benchmark/make_synth_dem.py diff --git a/raster/r.param.scale/benchmark/benchmark_r_param_scale_nprocs.py b/raster/r.param.scale/benchmark/benchmark_r_param_scale_nprocs.py new file mode 100644 index 00000000000..48d833907c5 --- /dev/null +++ b/raster/r.param.scale/benchmark/benchmark_r_param_scale_nprocs.py @@ -0,0 +1,57 @@ +"""Benchmarking of r.param.scale + +@author Kaushik Raja +""" + +import os +from subprocess import DEVNULL + +import grass.benchmark as bm +from grass.exceptions import CalledModuleError, GrassError +from grass.pygrass.modules import Module + + +def main(): + results = [] + + # Users can add more or modify existing reference maps + benchmark(10000, "r.param.scale_100M", results) + + bm.nprocs_plot(results, filename="r_param_scale_benchmark_nprocs.svg") + + +def benchmark(size, label, results): + reference = "r_param_scale_reference_map" + output = "benchmark_r_param_scale_nprocs" + + generate_map(rows=size, cols=size, fname=reference) + module = Module( + "r.param.scale", + input=reference, + output=output, + method="elev", + size=31, + memory=300, + run_=False, + stdout_=DEVNULL, + overwrite=True, + ) + results.append( + bm.benchmark_nprocs(module, label=label, max_nprocs=os.cpu_count(), repeat=3) + ) + Module("g.remove", quiet=True, flags="f", type="raster", name=(reference, output)) + + +def generate_map(rows, cols, fname): + Module("g.region", flags="p", s=0, n=rows, w=0, e=cols, res=1) + # Generate using r.random.surface if r.surf.fractal fails + try: + print("Generating reference map using r.surf.fractal...") + Module("r.surf.fractal", output=fname) + except (CalledModuleError, GrassError): + print("r.surf.fractal fails, using r.random.surface instead...") + Module("r.random.surface", output=fname) + + +if __name__ == "__main__": + main() diff --git a/raster/r.param.scale/benchmark/benchmark_rparamscale.py b/raster/r.param.scale/benchmark/benchmark_rparamscale.py deleted file mode 100644 index 576310453bf..00000000000 --- a/raster/r.param.scale/benchmark/benchmark_rparamscale.py +++ /dev/null @@ -1,169 +0,0 @@ -#!/usr/bin/env python3 -# -# benchmark.py - parallel scaling benchmark for r.param.scale -# -# WHAT THIS MEASURES -# Wall-clock parallel SPEEDUP (median time at 1 thread / median time at N threads) -# for r.param.scale on a fixed input raster. For each window size it runs every -# thread count, takes the median of the timed runs (one warmup is discarded), and -# reports the speedup of each thread count against that size's own 1-thread run. -# Defaults: -# input = synth_dem_100m (generate it with the companion DEM script) -# method = elev -# memory = 300 MB -# sizes = [31] (window side; add more to sweep window size) -# nprocs = [1, 2, 4, 8] -# repeat = 4 (1 warmup discarded + 3 timed, median reported) -# A config whose timed runs spread more than 5% is flagged. Occasional flags are -# EXPECTED on machines with variable thermals or scheduling; a flag just means that -# config should be re-run, not that anything is wrong. -# -# The script creates one temporary output raster -# (tmp_benchmark_rparamscale_delete_me) and removes it at the end, even on error. -# -# RESULTS ARE HARDWARE-DEPENDENT. Absolute times and the speedup curve depend on -# core count and type (performance vs efficiency cores), memory bandwidth, cooling, -# and power source. Your numbers will differ from anyone else's. Use AC power. -# -# THREAD PLACEMENT ON HYBRID CPUs -# On CPUs with both performance and efficiency cores (e.g. Apple M-series), thread -# counts below the performance-core count can show large run-to-run variance -# because the scheduler places threads on P or E cores inconsistently. Pin the -# threads to cores to stabilize those configs: -# OMP_PROC_BIND=true OMP_PLACES=cores grass <...> --exec python3 benchmark.py -# -# RUN -# grass /path/to/grassdata// --exec python3 benchmark.py -# (the input raster must exist in that mapset) - -import io -import os -import platform -import statistics -import subprocess -from contextlib import redirect_stdout - -import grass.script as gs -from grass.pygrass.modules import Module -from grass.benchmark import benchmark_single - -# ---- configuration ------------------------------------------------------- -INPUT = "synth_dem_100m" -METHOD = "elev" -MEMORY = 300 # MB, passed to memory= -SIZES = [31] # window sides to test (e.g. [5, 15, 31, 51] to sweep) -NPROCS = [1, 2, 4, 8] # thread counts to test (1 is the baseline) -REPEAT = 4 # 1 warmup (discarded) + (REPEAT - 1) timed runs -SPREAD_FLAG = 5.0 # percent; flag a config whose timed runs spread above this -TMP_OUTPUT = "tmp_benchmark_rparamscale_delete_me" # created then removed -# -------------------------------------------------------------------------- - - -def keep_awake(): - """On macOS, hold an idle-sleep assertion for this process's lifetime.""" - if platform.system() == "Darwin": - try: - return subprocess.Popen( - ["/usr/bin/caffeinate", "-i", "-w", str(os.getpid())] - ) - except FileNotFoundError: - pass - return None - - -def report_power(): - if platform.system() == "Darwin": - try: - ps = subprocess.check_output(["/usr/bin/pmset", "-g", "ps"], text=True) - if "Battery Power" in ps: - print("Power: BATTERY -- WARNING: plug into AC, results will be noisy") - else: - print("Power: AC") - except Exception: - print("Power: could not determine") - else: - print("Power: not checked (non-macOS); ensure AC power and no power saving") - - -def timed_median(all_times): - """Discard the warmup (first run); return (median, spread_percent) of the rest.""" - timed = all_times[1:] - med = statistics.median(timed) - spread = (max(timed) - min(timed)) / med * 100.0 if med else 0.0 - return med, spread - - -def bench(size, nprocs): - mod = Module( - "r.param.scale", - input=INPUT, - output=TMP_OUTPUT, - method=METHOD, - size=size, - memory=MEMORY, - nprocs=nprocs, - overwrite=True, - quiet=True, - stderr_=subprocess.PIPE, - run_=False, - ) - # benchmark_single prints raw per-run times; silence it and keep only the median. - with redirect_stdout(io.StringIO()): - res = benchmark_single(mod, label="size%d n%d" % (size, nprocs), repeat=REPEAT) - return timed_median(res.all_times) - - -def main(): - awake = keep_awake() - try: - print("GRASS %s" % gs.version()["version"]) - print("Machine: %s; %d logical CPUs" % (platform.platform(), os.cpu_count())) - report_power() - print( - "Workload: %s, method=%s, memory=%d MB, sizes=%s, nprocs=%s, " - "1 warmup + %d timed runs (median)" - % (INPUT, METHOD, MEMORY, SIZES, NPROCS, REPEAT - 1) - ) - print( - "Pinning: OMP_PROC_BIND=%s OMP_PLACES=%s" - % ( - os.environ.get("OMP_PROC_BIND", "unset"), - os.environ.get("OMP_PLACES", "unset"), - ) - ) - gs.run_command("g.region", raster=INPUT) - - for size in SIZES: - results = {p: bench(size, p) for p in NPROCS} - base = results[NPROCS[0]][0] - print( - "\nr.param.scale size=%d (%s, memory=%d MB, baseline = %d thread)" - % (size, INPUT, MEMORY, NPROCS[0]) - ) - print(" nprocs | median time (s) | speedup") - print(" -------+-----------------+---------") - for p in NPROCS: - med, spread = results[p] - flag = " *" if spread > SPREAD_FLAG else "" - print(" %d | %11.3f | %5.2fx%s" % (p, med, base / med, flag)) - flagged = [str(p) for p in NPROCS if results[p][1] > SPREAD_FLAG] - if flagged: - print( - " * timed-run spread > %.0f%% (expected on thermally/" - "scheduling-variable machines; just re-run that config): " - "nprocs %s" % (SPREAD_FLAG, ", ".join(flagged)) - ) - finally: - # Always remove the temporary output raster, even on error. - try: - gs.run_command( - "g.remove", type="raster", name=TMP_OUTPUT, flags="f", quiet=True - ) - except Exception: - pass - if awake: - awake.terminate() - - -if __name__ == "__main__": - main() diff --git a/raster/r.param.scale/benchmark/make_synth_dem.py b/raster/r.param.scale/benchmark/make_synth_dem.py deleted file mode 100644 index 4a6dcfa4cec..00000000000 --- a/raster/r.param.scale/benchmark/make_synth_dem.py +++ /dev/null @@ -1,59 +0,0 @@ -#!/usr/bin/env python3 -# -# make_synth_dem.py - generate the synthetic DEM used to benchmark r.param.scale -# -# WHAT THIS MAKES -# A synthetic terrain raster (DCELL) of ROWS x COLS cells, default 10000 x 10000 -# (100 million cells, about 800 MB), named "synth_dem_100m" -- the default input -# of the companion benchmark script. The surface is a deterministic sum of -# sinusoids at several frequencies, built with r.mapcalc from col()/row() only, so -# it is bit-reproducible on every machine and GRASS version (no random seed, no FFT). -# -# r.surf.fractal would give more realistic fractal terrain, but it depends on that -# module being built and on FFT-friendly dimensions; this analytic surface needs -# only r.mapcalc and works at any size. r.param.scale's run time depends on the -# raster size, not the cell values, so this is a valid and stable benchmark input. -# -# OVERWRITE -# This OVERWRITES any existing raster named OUTPUT (it passes --overwrite). Do not -# run it in a mapset where that name holds data you want to keep. -# -# RUN -# grass /path/to/grassdata// --exec python3 make_synth_dem.py -# (use a projected location; r.param.scale does not support lat/long) - -import grass.script as gs - -# ---- configuration ------------------------------------------------------- -ROWS = 10000 # lower ROWS/COLS for a quicker, smaller test (e.g. 2000 x 2000) -COLS = 10000 -OUTPUT = "synth_dem_100m" -# -------------------------------------------------------------------------- - -# Deterministic multi-frequency surface (a sum of sinusoids of col()/row(), with a -# +500 baseline so values read like positive elevations). DCELL because the constants -# and sin()/cos() are double precision. -EXPR = ( - "double(" - "500.0 + " - "256.0 * sin(col() * 0.9) * cos(row() * 0.9) + " - "128.0 * sin(col() * 2.1) * cos(row() * 1.7) + " - " 64.0 * sin(col() * 4.3) * cos(row() * 3.7) + " - " 32.0 * sin((col() + row()) * 6.1)" - ")" -) - - -def main(): - # Region: ROWS x COLS at resolution 1 (extent 0..COLS by 0..ROWS). - gs.run_command("g.region", n=ROWS, s=0, e=COLS, w=0, res=1, quiet=True) - gs.mapcalc("%s = %s" % (OUTPUT, EXPR), overwrite=True, quiet=True) - info = gs.raster_info(OUTPUT) - print( - "created %s: %d x %d cells (%s)" - % (OUTPUT, int(info["rows"]), int(info["cols"]), info["datatype"]) - ) - - -if __name__ == "__main__": - main() diff --git a/raster/r.param.scale/r.param.scale.md b/raster/r.param.scale/r.param.scale.md index 737b685c4a0..6cd0d48acad 100644 --- a/raster/r.param.scale/r.param.scale.md +++ b/raster/r.param.scale/r.param.scale.md @@ -42,6 +42,17 @@ downslope) is stored as (West is 0 degree, East is +/- 180 degree): Note that the aspect map is calculated differently from *[r.slope.aspect](r.slope.aspect.md)*. +### Performance + +To enable parallel processing, the user can specify the number of +threads to be used with the **nprocs** parameter (default 0, which uses +all available cores). The **memory** parameter (default 300) can also +be provided to determine the size of the buffer for computation. + +To reduce the memory requirements to minimum, set option **memory** to +zero. To take advantage of the parallelization, GRASS needs to be +compiled with OpenMP enabled. + ## EXAMPLE The next commands will create a geomorphological map of the Spearfish From 95a1404c567c42f05a0be4cf75d44846cd4dad15 Mon Sep 17 00:00:00 2001 From: Kaushik Raja Date: Sat, 27 Jun 2026 00:13:57 -0700 Subject: [PATCH 11/13] r.param.scale: replace testsuite with unskipped pytest and update benchmark --- .../benchmark_r_param_scale_nprocs.py | 104 ++- .../r.param.scale/tests/r_param_scale_test.py | 1 - .../testsuite/test_r_param_scale.py | 671 ------------------ 3 files changed, 83 insertions(+), 693 deletions(-) delete mode 100644 raster/r.param.scale/testsuite/test_r_param_scale.py diff --git a/raster/r.param.scale/benchmark/benchmark_r_param_scale_nprocs.py b/raster/r.param.scale/benchmark/benchmark_r_param_scale_nprocs.py index 48d833907c5..2b4f5f05369 100644 --- a/raster/r.param.scale/benchmark/benchmark_r_param_scale_nprocs.py +++ b/raster/r.param.scale/benchmark/benchmark_r_param_scale_nprocs.py @@ -1,56 +1,118 @@ """Benchmarking of r.param.scale - -@author Kaushik Raja +raster (2D) """ -import os -from subprocess import DEVNULL +# import os -import grass.benchmark as bm from grass.exceptions import CalledModuleError, GrassError from grass.pygrass.modules import Module +import grass.benchmark as bm + +# Baselines held fixed while one dimension is swept. +BASE_MAPSIZE = 50e6 # cells +BASE_WINDOW = 9 # window side, must be odd +BASE_MEMORY = 300 # MB +MAPSIZES = [10e6, 50e6, 100e6] +WINDOWS = [3, 9, 19, 39] +MEMORIES = [10, 100, 300, 1000] +METRICS = ["time", "speedup", "efficiency"] def main(): - results = [] + # Pin one thread per physical core so threads do not share a core's + # hyperthreads. The OpenMP runtime reads these when it initializes in each + # r.param.scale subprocess, so they must be set before the first run; the + # subprocesses inherit them from this environment. + # os.environ["OMP_PROC_BIND"] = "close" + # os.environ["OMP_PLACES"] = "cores" - # Users can add more or modify existing reference maps - benchmark(10000, "r.param.scale_100M", results) + # Sweep raster size at the baseline window and memory. + results = [] + for mapsize in MAPSIZES: + benchmark( + size=int(mapsize**0.5), + window=BASE_WINDOW, + memory=BASE_MEMORY, + label=f"r.param.scale_{int(mapsize / 1e6)}M", + results=results, + ) + plot(results, "rastersize") - bm.nprocs_plot(results, filename="r_param_scale_benchmark_nprocs.svg") + # Sweep window size at the baseline raster size and memory. + results = [] + for window in WINDOWS: + benchmark( + size=int(BASE_MAPSIZE**0.5), + window=window, + memory=BASE_MEMORY, + label=f"r.param.scale_window_{window}", + results=results, + ) + plot(results, "window") + # Sweep memory at the baseline raster size and window. + results = [] + for memory in MEMORIES: + benchmark( + size=int(BASE_MAPSIZE**0.5), + window=BASE_WINDOW, + memory=memory, + label=f"r.param.scale_memory_{memory}MB", + results=results, + ) + plot(results, "memory") -def benchmark(size, label, results): - reference = "r_param_scale_reference_map" - output = "benchmark_r_param_scale_nprocs" +def benchmark(size, window, memory, label, results): + reference = "benchmark_r_param_scale_reference" + output = "benchmark_r_param_scale" generate_map(rows=size, cols=size, fname=reference) module = Module( "r.param.scale", input=reference, output=output, - method="elev", - size=31, - memory=300, + method="feature", + size=window, + memory=memory, run_=False, - stdout_=DEVNULL, overwrite=True, ) results.append( - bm.benchmark_nprocs(module, label=label, max_nprocs=os.cpu_count(), repeat=3) + bm.benchmark_nprocs( + module, + label=label, + max_nprocs=12, + repeat=3, + ) + ) + Module( + "g.remove", + quiet=True, + flags="f", + type="raster", + pattern="benchmark_r_param_scale*", ) - Module("g.remove", quiet=True, flags="f", type="raster", name=(reference, output)) + + +def plot(results, sweep): + for metric in METRICS: + bm.nprocs_plot( + results, + filename=f"r_param_scale_{sweep}_{metric}.svg", + title=f"r.param.scale {sweep} {metric}", + metric=metric, + ) def generate_map(rows, cols, fname): - Module("g.region", flags="p", s=0, n=rows, w=0, e=cols, res=1) + Module("g.region", flags="p", n=rows, e=cols, res=1, w=0, s=0) # Generate using r.random.surface if r.surf.fractal fails try: print("Generating reference map using r.surf.fractal...") - Module("r.surf.fractal", output=fname) + Module("r.surf.fractal", output=fname, overwrite=True) except (CalledModuleError, GrassError): print("r.surf.fractal fails, using r.random.surface instead...") - Module("r.random.surface", output=fname) + Module("r.random.surface", output=fname, overwrite=True) if __name__ == "__main__": diff --git a/raster/r.param.scale/tests/r_param_scale_test.py b/raster/r.param.scale/tests/r_param_scale_test.py index 26795b14576..06a4f17ca6f 100644 --- a/raster/r.param.scale/tests/r_param_scale_test.py +++ b/raster/r.param.scale/tests/r_param_scale_test.py @@ -168,7 +168,6 @@ def _combos(): yield method, size -@pytest.mark.skip(reason="Waiting for a fix to the libgmath race, see #7539") @pytest.mark.parametrize(("method", "size"), list(_combos())) def test_param_scale_matches_reference(param_scale_session, method, size): """r.param.scale stats must match the main reference values.""" diff --git a/raster/r.param.scale/testsuite/test_r_param_scale.py b/raster/r.param.scale/testsuite/test_r_param_scale.py deleted file mode 100644 index ba5bf948bfd..00000000000 --- a/raster/r.param.scale/testsuite/test_r_param_scale.py +++ /dev/null @@ -1,671 +0,0 @@ -"""Test of r.param.scale parallelization (disk-strip path). - -Verifies that r.param.scale produces the same output when run threaded -(nprocs=4) as it does serially (nprocs=1), and locks in reference -univariate statistics for a sweep of methods, window sizes and the -constrained (-c) flag. - -Self-contained: the input DEM is generated in setUpClass with r.mapcalc -from a deterministic polynomial in row()/col() (pure +, -, * with decimal -constants, no trig and no rand), so it is bit-identical on every platform -and GRASS version. A small cubic term makes local curvature vary across -the raster while keeping the fit well conditioned. - -The region (400 x 400) together with memory=1 forces the module past its -memory budget, so it uses the multi-band per-thread file-descriptor -disk-strip path (brows < proc_rows) rather than holding the map in RAM. -The same memory setting is used for the serial and threaded runs. - -The reference values below were generated from a genuine single-threaded -run (OMP_NUM_THREADS=1, nprocs=1): the one-time LU decomposition relies on -a library pivot search that is only deterministic on a single thread, so -the reference must come from a single-threaded run. - -Used dataset: nc_spm_08_grass7 (only as a host location; the DEM is -generated, the NC rasters are not used). - -Author: Kaushik Raja -""" - -from grass.gunittest.case import TestCase -import grass.script as gs - - -REFERENCE = { - "elev_s3": { - "n": 158404, - "null_cells": 1596, - "cells": 160000, - "min": -564.762512, - "max": 1491.000666, - "range": 2055.763178, - "mean": 426.151076, - "mean_of_abs": 480.263031516931, - "stddev": 442.861043442243, - "variance": 196125.903798752, - "coeff_var": 103.921136982473, - "sum": 67504035.042704, - }, - "elev_s3_c": { - "n": 158404, - "null_cells": 1596, - "cells": 160000, - "min": -564.762512, - "max": 1491.000666, - "range": 2055.763178, - "mean": 426.151076, - "mean_of_abs": 480.263031516931, - "stddev": 442.861043442243, - "variance": 196125.903798752, - "coeff_var": 103.921136982473, - "sum": 67504035.042704, - }, - "elev_s9": { - "n": 153664, - "null_cells": 6336, - "cells": 160000, - "min": -546.478838, - "max": 1467.508536, - "range": 2013.987374, - "mean": 425.280891, - "mean_of_abs": 476.390592098136, - "stddev": 434.633318407327, - "variance": 188906.121469765, - "coeff_var": 102.199117713786, - "sum": 65350362.834624, - }, - "elev_s9_c": { - "n": 153664, - "null_cells": 6336, - "cells": 160000, - "min": -546.478838, - "max": 1467.508536, - "range": 2013.987374, - "mean": 425.280891, - "mean_of_abs": 476.390592098136, - "stddev": 434.633318407327, - "variance": 188906.121469765, - "coeff_var": 102.199117713786, - "sum": 65350362.834624, - }, - "elev_s15": { - "n": 148996, - "null_cells": 11004, - "cells": 160000, - "min": -528.372536, - "max": 1444.287922, - "range": 1972.660458, - "mean": 424.423924, - "mean_of_abs": 472.603709161615, - "stddev": 426.473749373164, - "variance": 181879.858904405, - "coeff_var": 100.482966500532, - "sum": 63237466.980304, - }, - "elev_s15_c": { - "n": 148996, - "null_cells": 11004, - "cells": 160000, - "min": -528.372536, - "max": 1444.287922, - "range": 1972.660458, - "mean": 424.423924, - "mean_of_abs": 472.603709161615, - "stddev": 426.473749373164, - "variance": 181879.858904405, - "coeff_var": 100.482966500532, - "sum": 63237466.980304, - }, - "slope_s3": { - "n": 158404, - "null_cells": 1596, - "cells": 160000, - "min": 31.143096450782, - "max": 83.9114645026956, - "range": 52.7683680519136, - "mean": 72.4250151761828, - "mean_of_abs": 72.4250151761828, - "stddev": 10.9080822146547, - "variance": 118.986257601667, - "coeff_var": 15.0612080482406, - "sum": 11472412.1039681, - }, - "slope_s3_c": { - "n": 158404, - "null_cells": 1596, - "cells": 160000, - "min": 31.143096450782, - "max": 83.9114645026956, - "range": 52.7683680519136, - "mean": 72.4250151761828, - "mean_of_abs": 72.4250151761828, - "stddev": 10.9080822146547, - "variance": 118.986257601667, - "coeff_var": 15.0612080482406, - "sum": 11472412.1039681, - }, - "slope_s9": { - "n": 153664, - "null_cells": 6336, - "cells": 160000, - "min": 31.944182807223, - "max": 83.8375409801301, - "range": 51.8933581729071, - "mean": 72.457518347825, - "mean_of_abs": 72.457518347825, - "stddev": 10.7142059615324, - "variance": 114.794209386136, - "coeff_var": 14.7868795479579, - "sum": 11134112.0994002, - }, - "slope_s9_c": { - "n": 153664, - "null_cells": 6336, - "cells": 160000, - "min": 31.944182807223, - "max": 83.8375409801301, - "range": 51.8933581729071, - "mean": 72.457518347825, - "mean_of_abs": 72.457518347825, - "stddev": 10.7142059615324, - "variance": 114.794209386136, - "coeff_var": 14.7868795479579, - "sum": 11134112.0994002, - }, - "slope_s15": { - "n": 148996, - "null_cells": 11004, - "cells": 160000, - "min": 32.7563367059122, - "max": 83.7623712622949, - "range": 51.0060345563827, - "mean": 72.4906681285151, - "mean_of_abs": 72.4906681285151, - "stddev": 10.5190698109164, - "variance": 110.650829686932, - "coeff_var": 14.5109295892647, - "sum": 10800819.5884762, - }, - "slope_s15_c": { - "n": 148996, - "null_cells": 11004, - "cells": 160000, - "min": 32.7563367059122, - "max": 83.7623712622949, - "range": 51.0060345563827, - "mean": 72.4906681285151, - "mean_of_abs": 72.4906681285151, - "stddev": 10.5190698109164, - "variance": 110.650829686932, - "coeff_var": 14.5109295892647, - "sum": 10800819.5884762, - }, - "profc_s3": { - "n": 158404, - "null_cells": 1596, - "cells": 160000, - "min": -0.00393749267654076, - "max": 0.00164711282024063, - "range": 0.00558460549678139, - "mean": -0.000342552958886785, - "mean_of_abs": 0.000491699464986382, - "stddev": 0.000829231153728635, - "variance": 6.87624306314123e-07, - "coeff_var": -242.073855214516, - "sum": -54.2617588995023, - }, - "profc_s3_c": { - "n": 158404, - "null_cells": 1596, - "cells": 160000, - "min": -0.00393749267654189, - "max": 0.00164711282024106, - "range": 0.00558460549678296, - "mean": -0.000342552958886785, - "mean_of_abs": 0.000491699464986382, - "stddev": 0.000829231153728635, - "variance": 6.87624306314123e-07, - "coeff_var": -242.073855214516, - "sum": -54.2617588995023, - }, - "profc_s9": { - "n": 153664, - "null_cells": 6336, - "cells": 160000, - "min": -0.00393703319719625, - "max": 0.00156307905481015, - "range": 0.0055001122520064, - "mean": -0.000343203532196393, - "mean_of_abs": 0.000487504566970183, - "stddev": 0.000819692530679842, - "variance": 6.71895844852323e-07, - "coeff_var": -238.835691880579, - "sum": -52.7380275714266, - }, - "profc_s9_c": { - "n": 153664, - "null_cells": 6336, - "cells": 160000, - "min": -0.00393703319719643, - "max": 0.00156307905480931, - "range": 0.00550011225200574, - "mean": -0.000343203532196394, - "mean_of_abs": 0.000487504566970183, - "stddev": 0.000819692530679842, - "variance": 6.71895844852324e-07, - "coeff_var": -238.835691880579, - "sum": -52.7380275714266, - }, - "profc_s15": { - "n": 148996, - "null_cells": 11004, - "cells": 160000, - "min": -0.00393547758337053, - "max": 0.0014817471036016, - "range": 0.00541722468697213, - "mean": -0.000343511486125908, - "mean_of_abs": 0.000483002036366709, - "stddev": 0.000809295470184593, - "variance": 6.54959158061301e-07, - "coeff_var": -235.594879027701, - "sum": -51.1818373868158, - }, - "profc_s15_c": { - "n": 148996, - "null_cells": 11004, - "cells": 160000, - "min": -0.00393547758337018, - "max": 0.00148174710360166, - "range": 0.00541722468697183, - "mean": -0.000343511486125908, - "mean_of_abs": 0.000483002036366709, - "stddev": 0.000809295470184592, - "variance": 6.54959158061301e-07, - "coeff_var": -235.594879027701, - "sum": -51.1818373868158, - }, - "crosc_s3": { - "n": 158404, - "null_cells": 1596, - "cells": 160000, - "min": -0.00798326486875052, - "max": 0.0140198635034536, - "range": 0.0220031283722041, - "mean": 0.00523063372892867, - "mean_of_abs": 0.00700064364287339, - "stddev": 0.0057384315690963, - "variance": 3.2929596873201e-05, - "coeff_var": 109.708151372924, - "sum": 828.553305197217, - }, - "crosc_s3_c": { - "n": 158404, - "null_cells": 1596, - "cells": 160000, - "min": -0.00798326486872778, - "max": 0.0140198635036052, - "range": 0.0220031283723329, - "mean": 0.00523063372892867, - "mean_of_abs": 0.00700064364287339, - "stddev": 0.0057384315690963, - "variance": 3.2929596873201e-05, - "coeff_var": 109.708151372924, - "sum": 828.553305197217, - }, - "crosc_s9": { - "n": 153664, - "null_cells": 6336, - "cells": 160000, - "min": -0.00790623053867785, - "max": 0.0139600812849915, - "range": 0.0218663118236694, - "mean": 0.00530129926236255, - "mean_of_abs": 0.00703005215756205, - "stddev": 0.00570183973099957, - "variance": 3.25109763180052e-05, - "coeff_var": 107.555515144763, - "sum": 814.618849851679, - }, - "crosc_s9_c": { - "n": 153664, - "null_cells": 6336, - "cells": 160000, - "min": -0.00790623053867993, - "max": 0.0139600812849679, - "range": 0.0218663118236478, - "mean": 0.00530129926236254, - "mean_of_abs": 0.00703005215756204, - "stddev": 0.00570183973099955, - "variance": 3.25109763180051e-05, - "coeff_var": 107.555515144763, - "sum": 814.618849851678, - }, - "crosc_s15": { - "n": 148996, - "null_cells": 11004, - "cells": 160000, - "min": -0.00782884301345861, - "max": 0.0139001591796006, - "range": 0.0217290021930592, - "mean": 0.00537273496983023, - "mean_of_abs": 0.0070588784164276, - "stddev": 0.00566360869239518, - "variance": 3.20764634205742e-05, - "coeff_var": 105.413885557324, - "sum": 800.516019564825, - }, - "crosc_s15_c": { - "n": 148996, - "null_cells": 11004, - "cells": 160000, - "min": -0.00782884301345516, - "max": 0.0139001591796018, - "range": 0.021729002193057, - "mean": 0.00537273496983024, - "mean_of_abs": 0.00705887841642761, - "stddev": 0.00566360869239518, - "variance": 3.20764634205743e-05, - "coeff_var": 105.413885557324, - "sum": 800.516019564826, - }, - "feature_s3": { - "n": 158404, - "null_cells": 1596, - "cells": 160000, - "min": 1.0, - "max": 5.0, - "range": 4.0, - "mean": 4.56869776015757, - "mean_of_abs": 4.56869776015757, - "stddev": 0.841755158740511, - "variance": 0.708551747266264, - "coeff_var": 18.4244001886323, - "sum": 723700.0, - }, - "feature_s3_c": { - "n": 158404, - "null_cells": 1596, - "cells": 160000, - "min": 1.0, - "max": 5.0, - "range": 4.0, - "mean": 4.56869776015757, - "mean_of_abs": 4.56869776015757, - "stddev": 0.841755158740511, - "variance": 0.708551747266264, - "coeff_var": 18.4244001886323, - "sum": 723700.0, - }, - "feature_s9": { - "n": 153664, - "null_cells": 6336, - "cells": 160000, - "min": 1.0, - "max": 5.0, - "range": 4.0, - "mean": 4.57504685547688, - "mean_of_abs": 4.57504685547688, - "stddev": 0.837584186678913, - "variance": 0.701547269774577, - "coeff_var": 18.3076635745539, - "sum": 703020.0, - }, - "feature_s9_c": { - "n": 153664, - "null_cells": 6336, - "cells": 160000, - "min": 1.0, - "max": 5.0, - "range": 4.0, - "mean": 4.57504685547688, - "mean_of_abs": 4.57504685547688, - "stddev": 0.837584186678913, - "variance": 0.701547269774577, - "coeff_var": 18.3076635745539, - "sum": 703020.0, - }, - "feature_s15": { - "n": 148996, - "null_cells": 11004, - "cells": 160000, - "min": 1.0, - "max": 5.0, - "range": 4.0, - "mean": 4.58155923648957, - "mean_of_abs": 4.58155923648957, - "stddev": 0.833230429104926, - "variance": 0.69427294798638, - "coeff_var": 18.1866125939988, - "sum": 682634.0, - }, - "feature_s15_c": { - "n": 148996, - "null_cells": 11004, - "cells": 160000, - "min": 1.0, - "max": 5.0, - "range": 4.0, - "mean": 4.58155923648957, - "mean_of_abs": 4.58155923648957, - "stddev": 0.833230429104926, - "variance": 0.69427294798638, - "coeff_var": 18.1866125939988, - "sum": 682634.0, - }, -} - - -class TestParamScale(TestCase): - """Parallel-vs-serial and reference-value regression for r.param.scale - on a generated DEM, exercising the disk-strip (multi-band) path. - """ - - # Generated input DEM and the fixed, reproducible region. - input_dem = "rps_test_dem" - region = {"n": 400, "s": 0, "e": 400, "w": 0, "res": 1} - dem_expr = ( - "double(100.0 + 0.5*col() + 0.4*row() " - "+ 0.003*col()*col() - 0.002*row()*row() + 0.0015*col()*row() " - "+ 0.00001*col()*col()*col() - 0.000008*row()*row()*row())" - ) - # Low memory budget forces the multi-band disk-strip path; identical for - # the serial and threaded runs so the comparison stays valid. - memory = 1 - - test_options = REFERENCE - to_remove = [] - - @classmethod - def setUpClass(cls): - cls.use_temp_region() - cls.runModule("g.region", **cls.region) - cls.runModule( - "r.mapcalc", - expression="{} = {}".format(cls.input_dem, cls.dem_expr), - overwrite=True, - ) - cls.to_remove.append(cls.input_dem) - - @classmethod - def tearDownClass(cls): - cls.del_temp_region() - if cls.to_remove: - cls.runModule( - "g.remove", flags="f", type="raster", name=",".join(cls.to_remove) - ) - - @staticmethod - def _key(method, size, constrained): - return "{}_s{}{}".format(method, size, "_c" if constrained else "") - - def _assert_identical(self, a, b, label): - """Null-aware bit-exact equality of two rasters. - - The diff is 1 wherever the null pattern differs OR (both non-null) - the values differ, and 0 otherwise. It is never null, so r.univar - does not skip cells. The differing-cell sum must be exactly 0, and - the two rasters must report the same non-null count. - """ - diff = "rps_diff_{}".format(label) - self.to_remove.append(diff) - expr = ( - "{d} = if(isnull({a}) != isnull({b}), 1, " - "if(isnull({a}), 0, if({a} != {b}, 1, 0)))" - ).format(d=diff, a=a, b=b) - self.runModule("r.mapcalc", expression=expr, overwrite=True) - dstats = gs.parse_command("r.univar", map=diff, flags="g") - differing = int(float(dstats.get("sum", 0) or 0)) - self.assertEqual( - differing, - 0, - msg="{} and {} differ in {} cells ({})".format(a, b, differing, label), - ) - a_n = int(gs.parse_command("r.univar", map=a, flags="g")["n"]) - b_n = int(gs.parse_command("r.univar", map=b, flags="g")["n"]) - self.assertEqual( - a_n, - b_n, - msg="non-null count mismatch ({} vs {}) for {}".format(a_n, b_n, label), - ) - - def _check_combo(self, method, size, constrained): - key = self._key(method, size, constrained) - serial = "rps_{}_serial".format(key) - threaded = "rps_{}_threaded".format(key) - self.to_remove.extend([serial, threaded]) - - common = { - "input": self.input_dem, - "method": method, - "size": size, - "memory": self.memory, - "overwrite": True, - } - serial_kw = dict(common, output=serial, nprocs=1) - threaded_kw = dict(common, output=threaded, nprocs=4) - if constrained: - serial_kw["flags"] = "c" - threaded_kw["flags"] = "c" - - self.assertModule("r.param.scale", **serial_kw) - self.assertModule("r.param.scale", **threaded_kw) - - # (a) reference-value regression guard (includes cell/null counts). - reference = self.test_options[key] - self.assertRasterFitsUnivar(raster=serial, reference=reference, precision=1e-5) - self.assertRasterFitsUnivar( - raster=threaded, reference=reference, precision=1e-5 - ) - # (b) null-aware bit-exact serial vs threaded. - self._assert_identical(serial, threaded, key) - - def _check_method(self, method): - for size in (3, 9, 15): - for constrained in (False, True): - self._check_combo(method, size, constrained) - - def test_elev(self): - """Method elev at sizes 3, 9, 15, with and without -c.""" - self._check_method("elev") - - def test_slope(self): - """Method slope at sizes 3, 9, 15, with and without -c.""" - self._check_method("slope") - - def test_profc(self): - """Method profc (profile curvature) at sizes 3, 9, 15, +/- -c.""" - self._check_method("profc") - - def test_crosc(self): - """Method crosc (cross-sectional curvature) at sizes 3, 9, 15, +/- -c.""" - self._check_method("crosc") - - def test_feature(self): - """Method feature (morphometric features) at sizes 3, 9, 15, +/- -c.""" - self._check_method("feature") - - def test_single_band_and_memory_invariance(self): - """Single-band path (memory=300) and invariance to band count. - - At size 15 the whole 400x400 map fits in memory=300 (one band), - whereas memory=1 forces multiple bands. For every method (default - and -c): the single-band serial and threaded outputs must be - identical, and the single-band serial output must be identical to - the multi-band (memory=1) serial output, i.e. the result is - invariant to how many bands the module uses. - """ - for method in ("elev", "slope", "profc", "crosc", "feature"): - for constrained in (False, True): - tag = "{}{}".format(method, "_c" if constrained else "") - s_m1 = "rps_inv_{}_s_m1".format(tag) - s_m300 = "rps_inv_{}_s_m300".format(tag) - t_m300 = "rps_inv_{}_t_m300".format(tag) - self.to_remove.extend([s_m1, s_m300, t_m300]) - - base = { - "input": self.input_dem, - "method": method, - "size": 15, - "overwrite": True, - } - runs = [ - dict(base, output=s_m1, nprocs=1, memory=1), - dict(base, output=s_m300, nprocs=1, memory=300), - dict(base, output=t_m300, nprocs=4, memory=300), - ] - if constrained: - for r in runs: - r["flags"] = "c" - for r in runs: - self.assertModule("r.param.scale", **r) - - # single-band serial == single-band threaded - self._assert_identical(s_m300, t_m300, "inv_{}_par".format(tag)) - # band-count invariance: one band (mem300) == multi band (mem1) - self._assert_identical(s_m300, s_m1, "inv_{}_mem".format(tag)) - # lock single-band output to the stored size-15 reference - reference = self.test_options[self._key(method, 15, constrained)] - self.assertRasterFitsUnivar( - raster=s_m300, reference=reference, precision=1e-5 - ) - self.assertRasterFitsUnivar( - raster=t_m300, reference=reference, precision=1e-5 - ) - - def test_mask(self): - """Masked sub-region: serial and threaded must match with a mask - active (exercises Rast_disable_omp_on_mask); masked cells are null - in both. Method elev, size 9, memory=1. - """ - maskmap = "rps_maskmap" - serial = "rps_mask_serial" - threaded = "rps_mask_threaded" - self.to_remove.extend([maskmap, serial, threaded]) - self.runModule( - "r.mapcalc", - expression=( - "{} = if(row() > 100 && row() < 300 " - "&& col() > 100 && col() < 300, 1, null())".format(maskmap) - ), - overwrite=True, - ) - self.runModule("r.mask", raster=maskmap) - try: - base = { - "input": self.input_dem, - "method": "elev", - "size": 9, - "overwrite": True, - "memory": 1, - } - self.assertModule("r.param.scale", output=serial, nprocs=1, **base) - self.assertModule("r.param.scale", output=threaded, nprocs=4, **base) - self._assert_identical(serial, threaded, "mask") - finally: - self.runModule("r.mask", flags="r") - - -if __name__ == "__main__": - from grass.gunittest.main import test - - test() From 2f32242981f314cc4558e255538c2294d3a7c038 Mon Sep 17 00:00:00 2001 From: Kaushik Raja Date: Mon, 6 Jul 2026 20:17:25 -0700 Subject: [PATCH 12/13] lib/gmath: Remove unsafe OpenMP pragma from G_ludcmp pivot search The pivot search in G_ludcmp parallelized its inner loop with an OpenMP pragma whose updates to the shared big and imax variables were unsynchronized, a data race. ThreadSanitizer reports the race directly, and it was reproducible as intermittent failures of the r.param.scale reference test, roughly one run in ten, whenever G_ludcmp ran with more than one thread from sequential code. Following maintainer direction on issue #7539, the pragma is removed rather than synchronized: it parallelized only a small one-time pivot loop, and having sequential callers spawn threads inside it is undesirable, especially on HPC systems. With the library routine now thread-safe by construction, the temporary workaround in r.param.scale that capped the one-time LU decomposition to a single thread is no longer needed and is removed as well. The output is verified bit-identical to the prior behavior and deterministic across thread counts: the r.param.scale reference test passes, the race-sensitive minic size 3 case produces a single checksum over repeated runs, and parallel and serial results match cell for cell. Related to #7539. --- lib/gmath/lu.c | 5 ----- raster/r.param.scale/process.c | 13 +------------ 2 files changed, 1 insertion(+), 17 deletions(-) diff --git a/lib/gmath/lu.c b/lib/gmath/lu.c index 84ab54d0050..bf8b8f90c6b 100644 --- a/lib/gmath/lu.c +++ b/lib/gmath/lu.c @@ -23,9 +23,6 @@ int G_ludcmp(double **a, int n, int *indx, double *d) vv = G_alloc_vector(n); *d = 1.0; - /* this pragma works, but doesn't really help speed things up */ - /* #pragma omp parallel for private(i, j, big, temp) shared(n, a, vv, - * is_singular) */ for (i = 0; i < n; i++) { big = 0.0; for (j = 0; j < n; j++) @@ -53,8 +50,6 @@ int G_ludcmp(double **a, int n, int *indx, double *d) } big = 0.0; - /* not very efficient, but this pragma helps speed things up a bit */ -#pragma omp parallel for private(i, k, sum, dum) shared(j, n, a, vv, big, imax) for (i = j; i < n; i++) { sum = a[i][j]; for (k = 0; k < j; k++) diff --git a/raster/r.param.scale/process.c b/raster/r.param.scale/process.c index b466426a40c..c4663c87ebc 100644 --- a/raster/r.param.scale/process.c +++ b/raster/r.param.scale/process.c @@ -63,22 +63,11 @@ void process(void) find_weight(weight_ptr); find_normal(normal_ptr, weight_ptr); - /* The one-time LU decomposition is run single-threaded because G_ludcmp's - * internal pivot search is not thread-safe (it has an unsynchronized shared - * reduction in lib/gmath/lu.c). Parallelizing this tiny one-time 6x6 solve - * gains nothing, and at nprocs>1 the race shifts the factorization by about - * 1 ULP, which then perturbs every per-cell solve. */ -#ifdef _OPENMP - int saved_threads = omp_get_max_threads(); - omp_set_num_threads(1); -#endif + /* One-time LU decomposition of the 6x6 normal-equations matrix. */ if (constrained) G_ludcmp(normal_ptr, 5, index_ptr, &temp); else G_ludcmp(normal_ptr, 6, index_ptr, &temp); -#ifdef _OPENMP - omp_set_num_threads(saved_threads); -#endif /*-----------------------------------------------------------------------*/ /* CHANGE 1: Per-thread input file descriptors. */ From 7fc3b28f95d2390e8522623fad8f99d6bab65c1f Mon Sep 17 00:00:00 2001 From: Anna Petrasova Date: Thu, 16 Jul 2026 04:02:33 -0400 Subject: [PATCH 13/13] r.param.scale: Add progress reporting and parallel tests, clean up comments --- raster/r.param.scale/Makefile | 4 + .../benchmark_r_param_scale_nprocs.py | 9 -- raster/r.param.scale/find_normal.c | 4 +- raster/r.param.scale/interface.c | 4 +- raster/r.param.scale/process.c | 127 +++++++----------- .../r.param.scale/tests/r_param_scale_test.py | 51 ++++++- 6 files changed, 105 insertions(+), 94 deletions(-) diff --git a/raster/r.param.scale/Makefile b/raster/r.param.scale/Makefile index fc5779d5413..e305a414dba 100644 --- a/raster/r.param.scale/Makefile +++ b/raster/r.param.scale/Makefile @@ -1,9 +1,13 @@ MODULE_TOPDIR = ../.. + PGM = r.param.scale + LIBES = $(GMATHLIB) $(RASTERLIB) $(GISLIB) $(MATHLIB) EXTRA_LIBS = $(OPENMP_LIBPATH) $(OPENMP_LIB) DEPENDENCIES = $(GMATHDEP) $(RASTERDEP) $(GISDEP) EXTRA_CFLAGS = $(OPENMP_CFLAGS) EXTRA_INC = $(OPENMP_INCPATH) + include $(MODULE_TOPDIR)/include/Make/Module.make + default: cmd diff --git a/raster/r.param.scale/benchmark/benchmark_r_param_scale_nprocs.py b/raster/r.param.scale/benchmark/benchmark_r_param_scale_nprocs.py index 2b4f5f05369..180a75ff404 100644 --- a/raster/r.param.scale/benchmark/benchmark_r_param_scale_nprocs.py +++ b/raster/r.param.scale/benchmark/benchmark_r_param_scale_nprocs.py @@ -2,8 +2,6 @@ raster (2D) """ -# import os - from grass.exceptions import CalledModuleError, GrassError from grass.pygrass.modules import Module import grass.benchmark as bm @@ -19,13 +17,6 @@ def main(): - # Pin one thread per physical core so threads do not share a core's - # hyperthreads. The OpenMP runtime reads these when it initializes in each - # r.param.scale subprocess, so they must be set before the first run; the - # subprocesses inherit them from this environment. - # os.environ["OMP_PROC_BIND"] = "close" - # os.environ["OMP_PLACES"] = "cores" - # Sweep raster size at the baseline window and memory. results = [] for mapsize in MAPSIZES: diff --git a/raster/r.param.scale/find_normal.c b/raster/r.param.scale/find_normal.c index 812383d7dbb..b447ccc6fa8 100644 --- a/raster/r.param.scale/find_normal.c +++ b/raster/r.param.scale/find_normal.c @@ -116,8 +116,8 @@ void find_obs(DCELL *z, /* Local window of elevs. */ for (row = 0; row < wsize; row++) for (col = 0; col < wsize; col++) { - /* w*z reused across all six sums; multiply order is kept */ - /* identical to the original so the result is unchanged. */ + /* w*z is shared by all six sums; the multiplication order + * affects rounding of the results. */ double wz; offset = row * wsize + col; diff --git a/raster/r.param.scale/interface.c b/raster/r.param.scale/interface.c index 5f33831aacd..1867cbcc245 100644 --- a/raster/r.param.scale/interface.c +++ b/raster/r.param.scale/interface.c @@ -32,7 +32,8 @@ void interface(int argc, char **argv) *parameter, /* Morphometric parameter to calculate. */ *expon, /* Inverse distance exponent for weight. */ *vert_sc, /* Vertical scaling factor. */ - *nprocs_opt, *mem_opt; + *nprocs_opt, /* Number of threads. */ + *mem_opt; /* Memory limit in MB. */ struct Flag *constr; /* Forces quadratic through the central */ @@ -55,6 +56,7 @@ void interface(int argc, char **argv) G_add_keyword(_("terrain")); G_add_keyword(_("elevation")); G_add_keyword(_("landform")); + G_add_keyword(_("parallel")); module->label = _("Extracts terrain parameters from a DEM."); module->description = _("Uses a multi-scale approach" " by taking fitting quadratic parameters to any " diff --git a/raster/r.param.scale/process.c b/raster/r.param.scale/process.c index c4663c87ebc..43c8c9425d0 100644 --- a/raster/r.param.scale/process.c +++ b/raster/r.param.scale/process.c @@ -1,10 +1,8 @@ /*****************************************************************************/ /*** process() ***/ -/*** Parallel (per-thread file descriptors): each thread opens its own ***/ -/*** file descriptor and reads its strip + halo directly from disk ***/ -/*** into a private buffer. No full map in RAM. Memory-light path for ***/ -/*** maps that do not fit in RAM. ***/ -/*** Based on the serial version by Jo Wood, Project ASSIST, 1993. ***/ +/*** Reads the input raster in bands and fits a quadratic surface ***/ +/*** in a moving window around each cell, in parallel with OpenMP. ***/ +/*** Based on the serial version by Jo Wood, Project ASSIST, 1993. ***/ /*****************************************************************************/ #include @@ -19,9 +17,9 @@ #include #endif -/* Rotate the per-thread ring: the oldest slot (0, smallest row index) is */ -/* recycled to become the newest slot (n-1). Mirrors r.neighbors rotate_bufs */ -/* (bufs.c:32-45). Pointer shuffle only, no data copy. */ +/* Rotate a ring of row buffers: the oldest slot (0, smallest row index) is + * recycled to become the newest slot (n-1). Pointer shuffle only, no data + * copy. */ static void rotate_ring(DCELL **ring, int n) { DCELL *temp = ring[0]; @@ -48,14 +46,10 @@ void process(void) else resoln = region.ns_res; - /*-----------------------------------------------------------------------*/ - /* CALCULATE LEAST SQUARES COEFFICIENTS (ONCE) */ - /* Same as serial. Weights, normal matrix, and LU decomposition depend */ - /* only on window size, not on data. Computed once, read-only by all */ - /* threads. */ - /*-----------------------------------------------------------------------*/ - - double *weight_ptr = (double *)G_malloc(DSQR(wsize) * sizeof(double)); + /* Weights, normal matrix, and its LU decomposition depend only on the + * window size, not on data. Computed once, read-only for all threads. */ + double *weight_ptr = + (double *)G_malloc((size_t)wsize * wsize * sizeof(double)); double **normal_ptr = dmatrix(0, 5, 0, 5); int *index_ptr = ivector(0, 5); double temp; @@ -69,22 +63,13 @@ void process(void) else G_ludcmp(normal_ptr, 6, index_ptr, &temp); - /*-----------------------------------------------------------------------*/ - /* CHANGE 1: Per-thread input file descriptors. */ - /* Each thread opens its own descriptor to the input raster. Each thread */ - /* will read its strip+halo directly from disk into its private buffer. */ - /* This is the r.neighbors pattern. No full_map array exists. */ - /*-----------------------------------------------------------------------*/ - + /* Each thread reads through its own descriptor to the input raster. */ int *fd_thread = G_malloc(sizeof(*fd_thread) * nprocs); for (int i = 0; i < nprocs; i++) fd_thread[i] = Rast_open_old(rast_in_name, ""); - /*-----------------------------------------------------------------------*/ - /* Band height from the memory cap. Reserve wsize input rows per thread */ - /* (the eventual ring-buffer cost), spend the rest on the band. */ - /*-----------------------------------------------------------------------*/ - + /* Band height from the memory cap: reserve wsize input rows per thread + * (the ring-buffer cost), spend the rest on the output band. */ size_t in_buf_size = sizeof(DCELL) * ncols * wsize * nprocs; /* memory= is in MiB; convert to bytes */ size_t out_buf_size = (size_t)memory * (1 << 20); @@ -98,15 +83,10 @@ void process(void) brows = proc_rows; if (brows < nprocs) brows = nprocs; - if (brows < 1) - brows = 1; - - /*-----------------------------------------------------------------------*/ - /* CHANGE 2: One shared band buffer (brows rows), reused each band, plus */ - /* a single NULL row for the EDGE borders. Threads write disjoint */ - /* band-relative rows; a serial loop writes each band in row order. */ - /*-----------------------------------------------------------------------*/ + /* One shared output buffer of brows rows, reused for each band, plus a + * single NULL row for the EDGE borders. Threads write disjoint rows of + * the band; a serial loop writes each band out in row order. */ DCELL *band_d = NULL, *null_d = NULL; CELL *band_c = NULL, *null_c = NULL; if (mparam != FEATURE) { @@ -120,10 +100,8 @@ void process(void) Rast_set_c_null_value(null_c, ncols); } - /* Top EDGE border rows are written as NULL. top_border is clamped to nrows - */ - /* so a region with fewer than 2*EDGE+1 rows cannot overflow (see bottom). - */ + /* Top EDGE border rows are written as NULL. top_border is clamped to + * nrows so a region with fewer than 2*EDGE+1 rows cannot overflow. */ int top_border = (EDGE < nrows) ? EDGE : nrows; for (row = 0; row < top_border; row++) { if (mparam != FEATURE) @@ -132,35 +110,30 @@ void process(void) Rast_put_row(fd_out, null_c, CELL_TYPE); } - /*-----------------------------------------------------------------------*/ - /* Per-thread scratch, allocated ONCE before the band loop and indexed by */ - /* tid (like r.neighbors allocate_bufs), instead of churning the allocator - */ - /* on every band: a wsize-row ring, the window matrix, and the obs vector. - */ - /*-----------------------------------------------------------------------*/ + /* Rows finished so far, shared across threads for progress reporting. */ + int computed = top_border; + /* Per-thread scratch, allocated once before the band loop and indexed + * by thread id: a wsize-row ring, the window matrix, and the observed + * vector. */ DCELL ***rings = G_malloc(sizeof(*rings) * nprocs); DCELL **window_ptrs = G_malloc(sizeof(*window_ptrs) * nprocs); double **obs_ptrs = G_malloc(sizeof(*obs_ptrs) * nprocs); - /* All ring rows live in one contiguous block; each thread's pointer */ - /* table (rings[t][k]) indexes into it, so rotate_ring only shuffles */ - /* pointers and there are no per-row allocations. */ + /* All ring rows live in one contiguous block; each thread's pointer + * table (rings[t][k]) indexes into it, so rotate_ring only shuffles + * pointers and there are no per-row allocations. */ DCELL *ring_block = G_malloc(sizeof(*ring_block) * nprocs * wsize * ncols); for (int t = 0; t < nprocs; t++) { rings[t] = G_malloc(sizeof(*rings[t]) * wsize); for (int k = 0; k < wsize; k++) rings[t][k] = ring_block + (size_t)(t * wsize + k) * ncols; - window_ptrs[t] = (DCELL *)G_malloc(DSQR(wsize) * sizeof(DCELL)); + window_ptrs[t] = + (DCELL *)G_malloc((size_t)wsize * wsize * sizeof(DCELL)); obs_ptrs[t] = dvector(0, 5); } - /*-----------------------------------------------------------------------*/ - /* CHANGE 3: Two-level structure. Outer serial loop over bands of brows */ - /* rows; inner parallel region splits the band across threads. Each */ - /* thread reads its sub-band + halo from disk into a private strip. */ - /*-----------------------------------------------------------------------*/ - + /* Outer serial loop over bands of brows rows; inner parallel region + * splits each band across threads. */ for (int band_start = EDGE; band_start < nrows - EDGE; band_start += brows) { int band_end = band_start + brows; @@ -179,9 +152,9 @@ void process(void) /* Split this band's rows [band_start .. band_end) across threads */ int my_start = - band_start + (int)((long)band_count * tid / nthreads); + band_start + (int)((size_t)band_count * tid / nthreads); int my_end = - band_start + (int)((long)band_count * (tid + 1) / nthreads); + band_start + (int)((size_t)band_count * (tid + 1) / nthreads); if (my_end > my_start) { /* Per-thread scratch, allocated once before the band loop. */ @@ -189,9 +162,9 @@ void process(void) DCELL *window_ptr = window_ptrs[tid]; double *obs_ptr = obs_ptrs[tid]; - /* Seed the ring with the first wsize-1 window rows; one new */ - /* row is read per output row below (read-once, mirrors */ - /* r.neighbors main.c:522-528). readrow = next disk row. */ + /* Seed the ring with the first wsize-1 window rows; one new + * row per output row is read below. readrow is the next row + * to read. */ int readrow = my_start - EDGE; for (int k = 0; k < wsize - 1; k++) { rotate_ring(ring, wsize); @@ -199,6 +172,8 @@ void process(void) } for (int orow = my_start; orow < my_end; orow++) { + G_percent(computed, nrows, 2); + size_t boff = (size_t)(orow - band_start) * ncols; DCELL *out_d = NULL; CELL *out_c = NULL; @@ -211,10 +186,9 @@ void process(void) Rast_set_c_null_value(out_c, ncols); } - /* Advance the ring by one row: rotate out the oldest, */ - /* read row orow+EDGE into the recycled newest slot; - * ring[k]*/ - /* then holds map row orow-EDGE+k (same layout as 3.1). */ + /* Advance the ring by one row: rotate out the oldest and + * read row orow+EDGE into the recycled newest slot; + * ring[k] then holds map row orow-EDGE+k. */ rotate_ring(ring, wsize); Rast_get_d_row(fd_thread[tid], ring[wsize - 1], readrow++); @@ -251,7 +225,6 @@ void process(void) if (found_null) continue; - /* Math block: unchanged from serial */ find_obs(window_ptr, obs_ptr, weight_ptr); if (constrained) G_lubksb(normal_ptr, 5, index_ptr, obs_ptr); @@ -266,6 +239,8 @@ void process(void) if (mparam == ELEV) *(out_d + col) += centre; } +#pragma omp atomic update + computed++; } } } /* end parallel */ @@ -283,23 +258,14 @@ void process(void) CELL_TYPE); } /* end band loop */ - /*-----------------------------------------------------------------------*/ - /* CHANGE 4: Close per-thread descriptors. */ - /*-----------------------------------------------------------------------*/ - for (int i = 0; i < nprocs; i++) Rast_close(fd_thread[i]); G_free(fd_thread); - /*-----------------------------------------------------------------------*/ - /* CHANGE 5: Bottom EDGE border rows are written as NULL. The computed */ - /* rows were written per band inside the loop above. */ - /*-----------------------------------------------------------------------*/ - - /* Bottom border: the remaining rows are NULL. With band_rows = max(0, */ - /* nrows-2*EDGE) computed rows, top_border + band_rows + bottom always sums - */ - /* to exactly nrows, even when the window is taller than the region. */ + /* Bottom border: the remaining rows are NULL. With band_rows = max(0, + * nrows - 2 * EDGE) computed rows, top_border + band_rows + bottom + * border always sums to exactly nrows, even when the window is taller + * than the region. */ int band_rows = (proc_rows > 0) ? proc_rows : 0; for (row = top_border + band_rows; row < nrows; row++) { if (mparam != FEATURE) @@ -307,6 +273,7 @@ void process(void) else Rast_put_row(fd_out, null_c, CELL_TYPE); } + G_percent(nrows, nrows, 2); if (mparam != FEATURE) { G_free(band_d); diff --git a/raster/r.param.scale/tests/r_param_scale_test.py b/raster/r.param.scale/tests/r_param_scale_test.py index 06a4f17ca6f..8d3a9403d6d 100644 --- a/raster/r.param.scale/tests/r_param_scale_test.py +++ b/raster/r.param.scale/tests/r_param_scale_test.py @@ -5,8 +5,9 @@ reference values captured from the current (main) state. This verifies that changes such as the OpenMP parallelization do not change the output. -Plain calls only (no nprocs, no memory), so the test runs identically on -main, where those options do not exist, and on the parallel branch. +The reference test runs with nprocs=1 to match the single-threaded +reference values. The parallel tests compare threaded and banded runs +against the serial run of the same binary. Reference provenance: the values below were generated with the main-branch r.param.scale binary at OMP_NUM_THREADS=1 (single threaded, so the internal @@ -180,6 +181,7 @@ def test_param_scale_matches_reference(param_scale_session, method, size): "output": output, "method": method, "size": size, + "nprocs": 1, "overwrite": True, } tools.r_param_scale(**kwargs) @@ -196,3 +198,48 @@ def test_param_scale_matches_reference(param_scale_session, method, size): # near-zero curvature outputs (profc/planc/minic) and zero stddev. for field in ("min", "max", "mean", "stddev", "sum"): assert stats[field] == pytest.approx(reference[field], rel=1e-6, abs=5e-8) + + +@pytest.mark.parametrize("method", ["slope", "feature"]) +@pytest.mark.parametrize("memory", [None, 0]) +def test_parallel_matches_serial(param_scale_session, method, memory): + """Threaded (and, with memory=0, banded) runs must match serial exactly.""" + tools = Tools(session=param_scale_session) + serial = f"serial_{method}" + tools.r_param_scale( + input="dem", output=serial, method=method, size=9, nprocs=1, overwrite=True + ) + + parallel = f"parallel_{method}_m{memory}" + kwargs = {} if memory is None else {"memory": memory} + tools.r_param_scale( + input="dem", + output=parallel, + method=method, + size=9, + nprocs=2, + overwrite=True, + **kwargs, + ) + + serial_stats = tools.r_univar(map=serial, format="json") + parallel_stats = tools.r_univar(map=parallel, format="json") + assert parallel_stats["null_cells"] == serial_stats["null_cells"] + assert parallel_stats["n"] == serial_stats["n"] + + diff = f"diff_{parallel}" + tools.r_mapcalc(expression=f"{diff} = abs({serial} - {parallel})") + diff_stats = tools.r_univar(map=diff, format="json") + assert diff_stats["n"] == serial_stats["n"] + assert diff_stats["max"] == 0 + + +def test_window_taller_than_region(param_scale_session): + """A window taller than the region must yield an all-null output.""" + tools = Tools(session=param_scale_session) + tools.r_param_scale( + input="dem", output="all_null", method="slope", size=51, overwrite=True + ) + stats = tools.r_univar(map="all_null", format="json") + assert stats["n"] == 0 + assert stats["null_cells"] == 2500