r.param.scale: parallelize with OpenMP

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++)

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)

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
 

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()

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);
 

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;
         }
 }
 

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);

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.         */

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)

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.         */