Added optional nthr to splinterpolator constructors. Will be used for deconvolution.

splinterpolator.h

@@ -14,6 +14,7 @@
 #include <vector>
 #include <string>
 #include <cmath>
+#include <thread>
 #include <iomanip>
 #include "armawrap/newmat.h"
 #include "miscmaths/miscmaths.h"
@@ -43,12 +44,21 @@ class Splinterpolator
 {
 public:
   // Constructors
-  Splinterpolator() : _valid(false), _own_coef(false), _coef(0), _cptr(0), _ndim(0) {}
+  Splinterpolator() : _valid(false), _own_coef(false), _coef(0), _cptr(0), _ndim(0), _nthr(1) {}
   Splinterpolator(const T *data, const std::vector<unsigned int>& dim, const std::vector<ExtrapolationType>& et, unsigned int order=3, bool copy_low_order=true, double prec=1e-8) : _valid(false), _own_coef(false), _coef(0), _cptr(0), _ndim(0)
   {
     common_construction(data,dim,order,prec,et,copy_low_order);
   }
-  Splinterpolator(const T *data, const std::vector<unsigned int>& dim, ExtrapolationType et=Zeros, unsigned int order=3, bool copy_low_order=true, double prec=1e-8) : _valid(false), _own_coef(false), _coef(0), _cptr(0), _ndim(0)
+  Splinterpolator(const T *data, const std::vector<unsigned int>& dim, ExtrapolationType et=Zeros, unsigned int order=3, bool copy_low_order=true, double prec=1e-8) : _valid(false), _own_coef(false), _coef(0), _cptr(0), _ndim(0), _nthr(1)
+  {
+    std::vector<ExtrapolationType>   ett(dim.size(),et);
+    common_construction(data,dim,order,prec,ett,copy_low_order);
+  }
+  Splinterpolator(const T *data, const std::vector<unsigned int>& dim, const std::vector<ExtrapolationType>& et, unsigned int order=3, bool copy_low_order=true, unsigned int nthr=1, double prec=1e-8) : _valid(false), _own_coef(false), _coef(0), _cptr(0), _ndim(0), _nthr(nthr)
+  {
+    common_construction(data,dim,order,prec,et,copy_low_order);
+  }
+  Splinterpolator(const T *data, const std::vector<unsigned int>& dim, ExtrapolationType et=Zeros, unsigned int order=3, bool copy_low_order=true, unsigned int nthr=1, double prec=1e-8) : _valid(false), _own_coef(false), _coef(0), _cptr(0), _ndim(0), _nthr(nthr)
   {
     std::vector<ExtrapolationType>   ett(dim.size(),et);
     common_construction(data,dim,order,prec,ett,copy_low_order);
@@ -194,6 +204,7 @@ private:
   const T                         *_cptr;         // Pointer to constant data. Used instead of _coef when we don't copy the data
   unsigned int                    _order;         // Order of splines
   unsigned int                    _ndim;          // # of non-singleton dimensions
+  unsigned int                    _nthr;          // Number of threads used for the deconvolution
   double                          _prec;          // Precision when dealing with edges
   std::vector<unsigned int>       _dim;           // Dimensions of data
   std::vector<ExtrapolationType>  _et;            // How to do extrapolation
@@ -205,6 +216,8 @@ private:
   void assign(const Splinterpolator<T>& src);
   bool calc_coef(const T *data, bool copy);
   void deconv_along(unsigned int dim);
+  void deconv_along_mt_helper(unsigned int dim, unsigned int mdim, unsigned int mstep, unsigned int offset, unsigned int step, 
+			      const std::vector<unsigned int>& rdim, const std::vector<unsigned int>& rstep);						
   T coef(int *indx) const;
   const T* coef_ptr() const {if (_own_coef) return(_coef); else return(_cptr); }
   unsigned int indx2indx(int indx, unsigned int d) const;
@@ -1368,6 +1381,7 @@ void Splinterpolator<T>::assign(const Splinterpolator<T>& src)
   _cptr = src._cptr;
   _order = src._order;
   _ndim = src._ndim;
+  _nthr = src._nthr;
   _prec = src._prec;
   _dim = src._dim;
   _et = src._et;
@@ -1438,16 +1452,50 @@ void Splinterpolator<T>::deconv_along(unsigned int dim)
     ss *= _dim[i];
   }
 
-  SplineColumn  col(mdim,mstep);          // Column helps us do the job
+  if (_nthr==1) { // If we are to run single-threaded
+    SplineColumn  col(mdim,mstep);          // Column helps us do the job
+    for (unsigned int l=0; l<rdim[3]; l++) {
+      for (unsigned int k=0; k<rdim[2]; k++) {
+	for (unsigned int j=0; j<rdim[1]; j++) {
+	  T *dp = _coef + l*rstep[3] + k*rstep[2] + j*rstep[1];
+	  for (unsigned int i=0; i<rdim[0]; i++, dp+=rstep[0]) {
+	    col.Get(dp);                          // Extract a column from the volume
+	    col.Deconv(_order,_et[dim],_prec);    // Deconvolve it
+	    col.Set(dp);                          // Put back the deconvolved column
+	  }
+	}
+      }
+    }
+  }
+  else { // We are running multi-threaded
+    std::vector<std::thread> threads(_nthr-1); // + main thread makes _nthr
+    for (unsigned int t=0; t<_nthr-1; t++) {
+      threads[t] = std::thread(&Splinterpolator::deconv_along_mt_helper,this,dim,mdim,mstep,t,_nthr,std::ref(rdim),std::ref(rstep));
+    }
+    deconv_along_mt_helper(dim,mdim,mstep,_nthr-1,_nthr,rdim,rstep);
+    std::for_each(threads.begin(),threads.end(),std::mem_fn(&std::thread::join)); // Join the threads
+  }
+  return;
+}
 
+template<class T>
+void Splinterpolator<T>::deconv_along_mt_helper(unsigned int                     dim,
+						unsigned int                     mdim, 
+						unsigned int                     mstep, 
+						unsigned int                     offset, // Offset into parallel dimension
+						unsigned int                     step,   // Step size in parallel dimension
+						const std::vector<unsigned int>& rdim,
+						const std::vector<unsigned int>& rstep)
+{
+  SplineColumn  col(mdim,mstep);          // Column helps us do the job
   for (unsigned int l=0; l<rdim[3]; l++) {
     for (unsigned int k=0; k<rdim[2]; k++) {
       for (unsigned int j=0; j<rdim[1]; j++) {
-        T *dp = _coef + l*rstep[3] + k*rstep[2] + j*rstep[1];
-	for (unsigned int i=0; i<rdim[0]; i++, dp+=rstep[0]) {
-	  col.Get(dp);                          // Extract a column from the volume
-          col.Deconv(_order,_et[dim],_prec);    // Deconvolve it
-          col.Set(dp);                          // Put back the deconvolved column
+	T *dp = _coef + l*rstep[3] + k*rstep[2] + j*rstep[1] + offset*rstep[0];
+	for (unsigned int i=offset; i<rdim[0]; i+=step, dp+=step*rstep[0]) {
+ 	  col.Get(dp);                          // Extract a column from the volume
+	  col.Deconv(_order,_et[dim],_prec);    // Deconvolve it
+	  col.Set(dp);                          // Put back the deconvolved column
 	}
       }
     }