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/CMakeLists.txt b/raster/CMakeLists.txt index ff92138b765..a3af02b2a71 100644 --- a/raster/CMakeLists.txt +++ b/raster/CMakeLists.txt @@ -411,7 +411,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) diff --git a/raster/r.param.scale/Makefile b/raster/r.param.scale/Makefile index 6cebd8a5d93..e305a414dba 100644 --- a/raster/r.param.scale/Makefile +++ b/raster/r.param.scale/Makefile @@ -3,7 +3,10 @@ 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 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..180a75ff404 --- /dev/null +++ b/raster/r.param.scale/benchmark/benchmark_r_param_scale_nprocs.py @@ -0,0 +1,110 @@ +"""Benchmarking of r.param.scale +raster (2D) +""" + +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(): + # 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") + + # 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, 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="feature", + size=window, + memory=memory, + run_=False, + overwrite=True, + ) + results.append( + 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*", + ) + + +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", 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, overwrite=True) + except (CalledModuleError, GrassError): + print("r.surf.fractal fails, using r.random.surface instead...") + Module("r.random.surface", output=fname, overwrite=True) + + +if __name__ == "__main__": + main() 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..b447ccc6fa8 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 is shared by all six sums; the multiplication order + * affects rounding of the results. */ + 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 2dac926caad..1867cbcc245 100644 --- a/raster/r.param.scale/interface.c +++ b/raster/r.param.scale/interface.c @@ -31,7 +31,9 @@ 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, /* Number of threads. */ + *mem_opt; /* Memory limit in MB. */ struct Flag *constr; /* Forces quadratic through the central */ @@ -54,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 " @@ -67,6 +70,8 @@ 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); + mem_opt = G_define_standard_option(G_OPT_MEMORYMB); constr = G_define_flag(); @@ -130,6 +135,11 @@ 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); + 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 8bb76b3ec67..29b8b6ca977 100644 --- a/raster/r.param.scale/main.c +++ b/raster/r.param.scale/main.c @@ -19,10 +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. */ + mparam, /* Morphometric parameter to calculate. */ + 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 f9f6d5e80f5..4816e830c45 100644 --- a/raster/r.param.scale/param.h +++ b/raster/r.param.scale/param.h @@ -85,10 +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. */ +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 6a0da1dbd87..43c8c9425d0 100644 --- a/raster/r.param.scale/process.c +++ b/raster/r.param.scale/process.c @@ -1,15 +1,8 @@ /*****************************************************************************/ - -/*** ***/ - /*** process() ***/ - -/*** Reads in a raster map row by row for processing. ***/ - -/*** Jo Wood, Project ASSIST, V1.0 7th February 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 @@ -20,56 +13,32 @@ #include "param.h" #include "nrutil.h" -void process(void) -{ - - /*-----------------------------------------------------------------------*/ - /* INITIALISE */ - - /*-----------------------------------------------------------------------*/ +#if defined(_OPENMP) +#include +#endif - 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 */ +/* 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]; + 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 */ + 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); + nrows = Rast_window_rows(); + ncols = Rast_window_cols(); - G_get_window(®ion); /* Fill out the region structure (the */ - /* geographical limits etc.) */ - - 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; @@ -77,210 +46,253 @@ void process(void) else resoln = region.ns_res; - /*-----------------------------------------------------------------------*/ - /* RESERVE MEMORY TO HOLD Z VALUES AND MATRICES */ + /* 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; - /*-----------------------------------------------------------------------*/ + find_weight(weight_ptr); + find_normal(normal_ptr, weight_ptr); - row_in = (DCELL *)G_malloc(ncols * sizeof(DCELL) * wsize); - /* Reserve `wsize' rows of memory. */ + /* 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); + /* 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 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); + 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; + + /* 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) { - row_out = - Rast_allocate_buf(DCELL_TYPE); /* Initialise output row buffer. */ - Rast_set_d_null_value(row_out, ncols); + 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 { - featrow_out = - Rast_allocate_buf(CELL_TYPE); /* Initialise output row buffer. */ - Rast_set_c_null_value(featrow_out, ncols); + band_c = G_malloc(sizeof(*band_c) * brows * ncols); + null_c = Rast_allocate_buf(CELL_TYPE); + Rast_set_c_null_value(null_c, 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. ---*/ - - find_weight(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. ---*/ - - find_normal(normal_ptr, weight_ptr); - - /*--- Apply LU decomposition to normal equations. ---*/ - - 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); - */ + /* 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) + Rast_put_row(fd_out, null_d, DCELL_TYPE); + else + Rast_put_row(fd_out, null_c, CELL_TYPE); } - else { - G_ludcmp(normal_ptr, 6, index_ptr, &temp); - /* disp_matrix(normal_ptr,obs_ptr,obs_ptr,6); - */ + /* 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. */ + 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((size_t)wsize * wsize * sizeof(DCELL)); + obs_ptrs[t] = dvector(0, 5); } - /*-----------------------------------------------------------------------*/ - /* PROCESS INPUT RASTER AND WRITE OUT RASTER LINE BY LINE */ + /* 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; + 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; +#if defined(_OPENMP) + nthreads = omp_get_num_threads(); + tid = omp_get_thread_num(); +#endif + + /* Split this band's rows [band_start .. band_end) across threads */ + int my_start = + band_start + (int)((size_t)band_count * tid / nthreads); + int my_end = + band_start + (int)((size_t)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 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); + Rast_get_d_row(fd_thread[tid], ring[wsize - 1], readrow++); + } - /*-----------------------------------------------------------------------*/ + for (int orow = my_start; orow < my_end; orow++) { + G_percent(computed, nrows, 2); - 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); - else { - Rast_set_d_null_value(row_out + col, 1); - } - 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)) { - if (mparam == FEATURE) - Rast_set_c_null_value(featrow_out + col, 1); - else { - Rast_set_d_null_value(row_out + col, 1); - } - found_null = TRUE; - break; + 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); } - *(window_ptr + (wind_row * wsize) + wind_col) = - *(window_cell)-centre; - } - } - if (found_null) - continue; + /* 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++); + + 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; + } - /*--- Use LU back substitution to solve normal equations. ---*/ - find_obs(window_ptr, obs_ptr, weight_ptr); + 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; + } + } + if (found_null) + continue; - /* disp_wind(window_ptr); - disp_matrix(normal_ptr,obs_ptr,obs_ptr,6); - */ + 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 (constrained) { - G_lubksb(normal_ptr, 5, index_ptr, obs_ptr); - /* - disp_matrix(normal_ptr,obs_ptr,obs_ptr,5); - */ - } + if (mparam == FEATURE) + *(out_c + col) = (CELL)feature(obs_ptr); + else + *(out_d + col) = param(mparam, obs_ptr); - else { - G_lubksb(normal_ptr, 6, index_ptr, obs_ptr); - /* - disp_matrix(normal_ptr,obs_ptr,obs_ptr,6); - */ + if (mparam == ELEV) + *(out_d + col) += centre; + } +#pragma omp atomic update + computed++; + } } + } /* end parallel */ - /*--- Calculate terrain parameter based on quad. coefficients. ---*/ - if (mparam == FEATURE) - *(featrow_out + col) = (CELL)feature(obs_ptr); - else - *(row_out + col) = param(mparam, obs_ptr); - - if (mparam == ELEV) - *(row_out + col) += centre; /* Add central elevation back */ - } - + /* Serial in-order write of this band's computed rows. */ 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++) { + 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 */ + + for (int i = 0; i < nprocs; i++) + Rast_close(fd_thread[i]); + G_free(fd_thread); + + /* 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) - Rast_put_row(fd_out, row_out, - DCELL_TYPE); /* Write out the edge cells as NULL. */ + Rast_put_row(fd_out, null_d, DCELL_TYPE); else - Rast_put_row(fd_out, featrow_out, - CELL_TYPE); /* Write out the edge cells as NULL. */ + Rast_put_row(fd_out, null_c, CELL_TYPE); } + G_percent(nrows, nrows, 2); - /*-----------------------------------------------------------------------*/ - /* FREE MEMORY USED TO STORE RASTER ROWS, LOCAL WINDOW AND MATRICES */ - - /*-----------------------------------------------------------------------*/ - - G_free(row_in); - if (mparam != FEATURE) - G_free(row_out); - else - G_free(featrow_out); - - G_free(window_ptr); + if (mparam != FEATURE) { + G_free(band_d); + G_free(null_d); + } + else { + G_free(band_c); + G_free(null_c); + } + for (int t = 0; t < nprocs; t++) { + 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); + G_free(weight_ptr); free_dmatrix(normal_ptr, 0, 5, 0, 5); - free_dvector(obs_ptr, 0, 5); free_ivector(index_ptr, 0, 5); } 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 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..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 @@ -168,7 +169,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.""" @@ -181,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) @@ -197,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