Remapstat: OpenMP version.

src/Remapstat.cc

@@ -437,7 +437,7 @@ Varray2D<size_t> gen_mapdata(int gridID1, int gridID2)
   const auto gridtype2 = gridInqType(gridID2);
   auto grid2_is_reg2d = (gridtype2 == GRID_GAUSSIAN || gridtype2 == GRID_LONLAT);
 
-  const auto nlon2 = grid2_is_reg2d ? gridInqXsize(gridID2) : 0;
+  auto nlon2 = grid2_is_reg2d ? gridInqXsize(gridID2) : 0;
 
   int gridID2x = -1;
 
@@ -477,15 +477,23 @@ Varray2D<size_t> gen_mapdata(int gridID1, int gridID2)
 
   start = Options::cdoVerbose ? cdo_get_wtime() : 0;
 
-  size_t ndist = gridsize1;
+  size_t ndist_max = gridsize1;
+  if (gridsize1 > 1000000) ndist_max /= 4;
   std::vector<char> vmask(gridsize1, 0);
-  std::vector<size_t> adds(gridsize1);
-  std::vector<double> dist(gridsize1);
+  Varray2D<size_t> adds_2D(Threading::ompNumThreads, Varray<size_t>(ndist_max));
+  Varray2D<double> dist_2D(Threading::ompNumThreads, Varray<double>(ndist_max));
 
   StatInfo statInfo;
 
+#ifdef _OPENMP
+#pragma omp parallel for default(none) schedule(dynamic) shared(Options::cdoVerbose, gridsize2, gps, statInfo, mapdata, xvals1, yvals1, xvals2, yvals2, grid2_is_reg2d, ndist_max, nv, nlon2, xbounds2d, ybounds2d, xbounds, ybounds, adds_2D, dist_2D, vmask)
+#endif
   for (size_t i = 0; i < gridsize2; ++i)
     {
+      const auto ompthID = cdo_omp_get_thread_num();
+      auto &adds = adds_2D[ompthID];
+      auto &dist = dist_2D[ompthID];
+
       const auto lon2 = xvals2[i];
       const auto lat2 = yvals2[i];
 
@@ -493,7 +501,7 @@ Varray2D<size_t> gen_mapdata(int gridID1, int gridID2)
         ? calc_maxdist_rec2d(i, nlon2, lon2, lat2, xbounds2d, ybounds2d)
         : calc_maxdist(i, nv, lon2, lat2, xbounds, ybounds);
 
-      const auto nadds = gridPointSearchDistanceQnearest(gps, maxdist, lon2, lat2, ndist, adds.data(), dist.data());
+      const auto nadds = gridPointSearchDistanceQnearest(gps, maxdist, lon2, lat2, ndist_max, adds.data(), dist.data());
       // printf("%zu nadds %zu\n", i+1, nadds);
 
       const auto nvalues = grid2_is_reg2d