Testing branch

lib/calc/xnmedian.c

@@ -44,8 +44,7 @@ static int dcmp(const void *aa, const void *bb)
 
 int f_nmedian(int argc, const int *argt, void **args)
 {
-    static void *array;
-    static int alloc;
+    void *array;
     int size = argc * Rast_cell_size(argt[0]);
     int i, j;
 
@@ -56,10 +55,7 @@ int f_nmedian(int argc, const int *argt, void **args)
         if (argt[i] != argt[0])
             return E_ARG_TYPE;
 
-    if (size > alloc) {
-        alloc = size;
-        array = G_realloc(array, size);
-    }
+    array = G_malloc(size);
 
     switch (argt[0]) {
     case CELL_TYPE: {

raster/r.mapcalc/CMakeLists.txt

@@ -30,7 +30,8 @@ build_program(
   OPTIONAL_DEPENDS
   Readline::Readline
   Readline::History
-  Threads::Threads)
+  Threads::Threads
+  OPENMP)
 
 build_program(
   NAME

raster/r.mapcalc/Makefile

@@ -12,8 +12,10 @@ r3_mapcalc_OBJS := $(filter-out map.o xcoor.o xres.o, $(AUTO_OBJS))
 include $(MODULE_TOPDIR)/include/Make/Multi.make
 
 EXTRA_CFLAGS = $(READLINEINCPATH) $(PTHREADINCPATH)
-LIBES2 = $(CALCLIB) $(GISLIB) $(RASTERLIB) $(BTREELIB) $(READLINELIBPATH) $(READLINELIB) $(HISTORYLIB) $(PTHREADLIBPATH) $(PTHREADLIB)
-LIBES3 = $(CALCLIB) $(RASTER3DLIB) $(GISLIB) $(RASTERLIB) $(BTREELIB) $(READLINELIBPATH) $(READLINELIB) $(HISTORYLIB) $(PTHREADLIBPATH) $(PTHREADLIB)
+LIBES2 = $(CALCLIB) $(GISLIB) $(RASTERLIB) $(BTREELIB) $(READLINELIBPATH) $(READLINELIB) $(HISTORYLIB) $(PTHREADLIBPATH) $(PTHREADLIB) $(OPENMP_LIBPATH) $(OPENMP_LIB)
+LIBES3 = $(CALCLIB) $(RASTER3DLIB) $(GISLIB) $(RASTERLIB) $(BTREELIB) $(READLINELIBPATH) $(READLINELIB) $(HISTORYLIB) $(PTHREADLIBPATH) $(PTHREADLIB) $(OPENMP_LIBPATH) $(OPENMP_LIB)
+EXTRA_CFLAGS = $(OPENMP_CFLAGS)
+EXTRA_INC = $(OPENMP_INCPATH)
 
 default: multi
 

raster/r.mapcalc/benchmark/benchmark_rmapcalc.py

@@ -0,0 +1,104 @@
+"""Benchmarking of r.mapcalc
+raster (2D)
+
+@author Chung-Yuan Liang, 2025
+"""
+
+from grass.exceptions import CalledModuleError
+from grass.pygrass.modules import Module
+
+import grass.benchmark as bm
+
+
+def main():
+    results = []
+    metrics = ["time", "speedup", "efficiency"]
+    mapsizes = [10e6, 50e6, 100e6]
+
+    # run benchmarks
+
+    for mapsize in mapsizes:
+        benchmark(
+            size=int(mapsize**0.5),
+            step=0,
+            results=results,
+        )
+
+    # plot results
+
+    for metric in metrics:
+        bm.nprocs_plot(
+            results,
+            filename=f"r_mapcalc_{metric}.svg",
+            title=f"r.mapcalc {metric}",
+            metric=metric,
+        )
+
+
+def benchmark(size, step, results):
+    map1 = "benchmark_r_mapcalc_map1"
+    map2 = "benchmark_r_mapcalc_map2"
+    output = "benchmark_r_mapcalc"
+
+    generate_map(rows=size, cols=size, fname=map1)
+    generate_map(rows=size, cols=size, fname=map2)
+    module = Module(
+        "r.mapcalc",
+        expression=f"{output}=({map1} + {map2})",
+        overwrite=True,
+    )
+
+    results.append(
+        bm.benchmark_nprocs(
+            module,
+            label=f"r.mapcalc_simple_{int((size * size) / 1e6)}M",
+            max_nprocs=12,
+            repeat=5,
+        )
+    )
+
+    module = Module(
+        "r.mapcalc",
+        expression=f"{output}=({map1} + {map2} - 2*{map2} + {map1}*{map2} - {map1}/{map2})/2",
+        overwrite=True,
+    )
+    results.append(
+        bm.benchmark_nprocs(
+            module,
+            label=f"r.mapcalc_complex_{int((size * size) / 1e6)}M",
+            max_nprocs=12,
+            repeat=5,
+        )
+    )
+
+    module = Module(
+        "r.mapcalc",
+        expression=f"{output}= if (({map1}[5, 5] + {map2}[-5, -5]) > 0, 1.6, 0.1)",
+        overwrite=True,
+    )
+    results.append(
+        bm.benchmark_nprocs(
+            module,
+            label=f"r.mapcalc_neighbor_{int((size * size) / 1e6)}M",
+            max_nprocs=12,
+            repeat=5,
+        )
+    )
+    Module(
+        "g.remove", quiet=True, flags="f", type="raster", pattern="benchmark_r_mapcalc*"
+    )
+
+
+def generate_map(rows, cols, fname):
+    Module("g.region", flags="p", rows=rows, cols=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, overwrite=True)
+    except CalledModuleError:
+        print("r.surf.fractal fails, using r.random.surface instead...")
+        Module("r.random.surface", output=fname, overwrite=True)
+
+
+if __name__ == "__main__":
+    main()