diff --git a/splinterpolator.h b/splinterpolator.h
index 87a9c957168498ed49255f837aca1a261f927ff9..53fe8e04b5c866b0ba7b229b722156a2a5482a86 100644
--- a/splinterpolator.h
+++ b/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
}
}
}