Add and use parallelisation abstraction

src/Distancer.cpp

@@ -1,10 +1,8 @@
 #include <Rcpp.h>
 
+#include "Parallel.h"
 #include "Distancer.h"
 
-#ifdef _OPENMP
-#include <omp.h>
-#endif
 
 double initialTransform (const double &x) { return x == 0.0 ? R_PosInf : 0.0; }
 
@@ -32,12 +30,8 @@ Array<double> * Distancer::run ()
     {
         const double sqPixdim = usePixdim ? pixdims[i] * pixdims[i] : 1.0;
 
-#ifdef _OPENMP
-#pragma omp parallel for
-#endif
         // Lines are independent, so can be processed in parallel
-        for (size_t j=0; j<result->countLines(i); j++)
-        {
+        PARALLEL_LOOP_START(j, result->countLines(i))
             // The vertices are the minima of a series of parabolas. The
             // intersections are the locations where they cross
             std::vector<int> vertices;
@@ -100,7 +94,7 @@ Array<double> * Distancer::run ()
                     dx *= pixdims[i];
                 it[l] = line[vertices[k]] + dx * dx;
             }
-        }
+        PARALLEL_LOOP_END
     }
 
     // Take the square-root of each value to get Euclidean distance. The

src/Parallel.h

@@ -0,0 +1,29 @@
+#ifndef PARALLEL_H_
+#define PARALLEL_H_
+
+#if defined(HAVE_LIBDISPATCH) && defined(HAVE_BLOCKS)
+
+#include <dispatch/dispatch.h>
+
+#define PARALLEL_LOOP_START(i,n) \
+    dispatch_apply(size_t(n), DISPATCH_APPLY_AUTO, ^(size_t i) {
+#define PARALLEL_LOOP_END });
+
+#elif _OPENMP
+
+#include <omp.h>
+
+#define PARALLEL_LOOP_START(i,n) \
+    _Pragma("omp parallel for")  \
+    for (size_t i=0; i<size_t(n); i++) {
+#define PARALLEL_LOOP_END }
+
+#else
+
+#define PARALLEL_LOOP_START(i,n) \
+    for (size_t i=0; i<size_t(n); i++) {
+#define PARALLEL_LOOP_END }
+
+#endif
+
+#endif

src/Resampler.cpp

@@ -1,11 +1,8 @@
 #include <Rcpp.h>
 
+#include "Parallel.h"
 #include "Resampler.h"
 
-#ifdef _OPENMP
-#include <omp.h>
-#endif
-
 // Presharpen data (i.e., calculate spline coefficients) along a single line
 // This function is slightly inscrutable but aims to be fast and general
 template <class InputIterator, class OutputIterator>
@@ -117,11 +114,7 @@ const std::vector<double> & Resampler::run (const Rcpp::NumericMatrix &locations
 
     samples.resize(nSamples);
 
-#ifdef _OPENMP
-#pragma omp parallel for
-#endif
-    for (int k=0; k<nSamples; k++)
-    {
+    PARALLEL_LOOP_START(k, nSamples)
         int_vector base(nDims);
         dbl_vector offset(nDims);
 
@@ -141,7 +134,7 @@ const std::vector<double> & Resampler::run (const Rcpp::NumericMatrix &locations
             }
         }
         samples[k] = samplePoint(base, offset, nDims-1);
-    }
+    PARALLEL_LOOP_END
 
     return samples;
 }
@@ -159,15 +152,11 @@ const std::vector<double> & Resampler::run (const std::vector<dbl_vector> &locat
         dims[i] = locations[i].size();
         Array<double> *result = new Array<double>(dims, NA_REAL);
 
-#ifdef _OPENMP
-#pragma omp parallel for
-#endif
-        for (int j=0; j<int(working->countLines(i)); j++)
-        {
+        PARALLEL_LOOP_START(j, working->countLines(i))
             CachedInterpolant interpolant(working->beginLine(j,i), working->endLine(j,i));
             interpolate(interpolant, locations[i], result->beginLine(j,i));
-        }
-
+        PARALLEL_LOOP_END
+        
         delete working;
         working = result;
     }