Skip to content

GitLab

Commit 9b34409f authored by Martin Reinecke's avatar Martin Reinecke
Browse files

MAVs can manage their own memory if requested

parent 1d0a0e1d
Hide whitespace changes
Inline Side-by-side
mr_util/fft.h
View file @ 9b34409f
......@@ -2898,8 +2898,7 @@ template<typename T> void c2r(const cfmav<std::complex<T>> &in,
return c2r(in, out, axes[0], forward, fct, nthreads);
util::sanity_check_cr(in, out, axes);
if (in.size()==0) return;
aligned_array<std::complex<T>> tmp(in.size());
fmav<std::complex<T>> atmp(tmp.data(), in);
fmav<std::complex<T>> atmp(in.shape());
auto newaxes = shape_t({axes.begin(), --axes.end()});
c2c(in, atmp, newaxes, forward, T(1), nthreads);
c2r(atmp, out, axes.back(), forward, fct, nthreads);
......@@ -2933,9 +2932,7 @@ template<typename T> void r2r_genuine_hartley(const cfmav<T> &in,
if (in.size()==0) return;
shape_t tshp(in.shape());
tshp[axes.back()] = tshp[axes.back()]/2+1;
auto tinfo = fmav_info(tshp);
aligned_array<std::complex<T>> tdata(tinfo.size());
fmav<std::complex<T>> atmp(tdata.data(), tinfo);
fmav<std::complex<T>> atmp(tshp);
r2c(in, atmp, axes, true, fct, nthreads);
simple_iter iin(atmp);
rev_iter iout(out, axes);
......
mr_util/mav.h
View file @ 9b34409f
......@@ -25,6 +25,7 @@
#include <cstdlib>
#include <array>
#include <vector>
#include <memory>
#include "mr_util/error_handling.h"
namespace mr {
......@@ -85,10 +86,14 @@ class fmav_info
{ return shp==other.shp; }
};
template<typename T, size_t ndim> class cmav;
// "mav" stands for "multidimensional array view"
template<typename T> class cfmav: public fmav_info
{
protected:
using Tsp = shared_ptr<vector<T>>;
Tsp ptr;
T *d;
public:
......@@ -96,16 +101,26 @@ template<typename T> class cfmav: public fmav_info
: fmav_info(shp_, str_), d(const_cast<T *>(d_)) {}
cfmav(const T *d_, const shape_t &shp_)
: fmav_info(shp_), d(const_cast<T *>(d_)) {}
cfmav(const shape_t &shp_)
: fmav_info(shp_), ptr(make_unique<vector<T>>(size())),
d(const_cast<T *>(ptr->data())) {}
cfmav(const T* d_, const fmav_info &info)
: fmav_info(info), d(const_cast<T *>(d_)) {}
cfmav(const cfmav &other) = default;
cfmav(cfmav &&other) = default;
// Not for public use!
cfmav(const T *d_, const Tsp &p, const shape_t &shp_, const stride_t &str_)
: fmav_info(shp_, str_), ptr(p), d(const_cast<T *>(d_)) {}
template<typename I> const T &operator[](I i) const
{ return d[i]; }
const T *data() const
{ return d; }
};
template<typename T> class fmav: public cfmav<T>
{
protected:
using Tsp = shared_ptr<vector<T>>;
using parent = cfmav<T>;
using parent::d;
using parent::shp;
......@@ -116,8 +131,16 @@ template<typename T> class fmav: public cfmav<T>
: parent(d_, shp_, str_) {}
fmav(T *d_, const shape_t &shp_)
: parent(d_, shp_) {}
fmav(const shape_t &shp_)
: parent(shp_) {}
fmav(T* d_, const fmav_info &info)
: parent(d_, info) {}
fmav(const fmav &other) = default;
fmav(fmav &&other) = default;
// Not for public use!
fmav(T *d_, const Tsp &p, const shape_t &shp_, const stride_t &str_)
: parent(d_, p, shp_, str_) {}
template<typename I> T &operator[](I i) const
{ return d[i]; }
using parent::shape;
......@@ -129,31 +152,40 @@ template<typename T> class fmav: public cfmav<T>
using parent::conformable;
};
template<typename T, size_t ndim> class cmav
{
static_assert((ndim>0) && (ndim<3), "only supports 1D and 2D arrays");
static_assert((ndim>0) && (ndim<4), "only supports 1D, 2D, and 3D arrays");
protected:
T *d;
using Tsp = shared_ptr<vector<T>>;
array<size_t, ndim> shp;
array<ptrdiff_t, ndim> str;
Tsp ptr;
T *d;
public:
cmav(const T *d_, const array<size_t,ndim> &shp_,
const array<ptrdiff_t,ndim> &str_)
: d(const_cast<T *>(d_)), shp(shp_), str(str_) {}
: shp(shp_), str(str_), d(const_cast<T *>(d_)) {}
cmav(const T *d_, const array<size_t,ndim> &shp_)
: d(const_cast<T *>(d_)), shp(shp_)
: shp(shp_), d(const_cast<T *>(d_))
{
str[ndim-1]=1;
for (size_t i=2; i<=ndim; ++i)
str[ndim-i] = str[ndim-i+1]*shp[ndim-i+1];
}
cmav(const array<size_t,ndim> &shp_)
: shp(shp_), ptr(make_unique<vector<T>>(size())), d(ptr->data())
{
str[ndim-1]=1;
for (size_t i=2; i<=ndim; ++i)
str[ndim-i] = str[ndim-i+1]*shp[ndim-i+1];
}
cmav(const cmav &other) = default;
cmav(cmav &&other) = default;
operator cfmav<T>() const
{
return cfmav<T>(d, {shp.begin(), shp.end()}, {str.begin(), str.end()});
return cfmav<T>(d, ptr, {shp.begin(), shp.end()}, {str.begin(), str.end()});
}
const T &operator[](size_t i) const
{ return operator()(i); }
......@@ -167,6 +199,11 @@ template<typename T, size_t ndim> class cmav
static_assert(ndim==2, "ndim must be 2");
return d[str[0]*i + str[1]*j];
}
const T &operator()(size_t i, size_t j, size_t k) const
{
static_assert(ndim==3, "ndim must be 3");
return d[str[0]*i + str[1]*j + str[2]*k];
}
size_t shape(size_t i) const { return shp[i]; }
const array<size_t,ndim> &shape() const { return shp; }
size_t size() const
......@@ -196,8 +233,10 @@ template<typename T, size_t ndim> class cmav
template<typename T, size_t ndim> class mav: public cmav<T, ndim>
{
protected:
using Tsp = shared_ptr<vector<T>>;
using parent = cmav<T, ndim>;
using parent::d;
using parent::ptr;
using parent::shp;
using parent::str;
......@@ -207,9 +246,13 @@ template<typename T, size_t ndim> class mav: public cmav<T, ndim>
: parent(d_, shp_, str_) {}
mav(T *d_, const array<size_t,ndim> &shp_)
: parent(d_, shp_) {}
mav(const array<size_t,ndim> &shp_)
: parent(shp_) {}
mav(const mav &other) = default;
mav(mav &&other) = default;
operator fmav<T>() const
{
return fmav<T>(d, {shp.begin(), shp.end()}, {str.begin(), str.end()});
return fmav<T>(d, ptr, {shp.begin(), shp.end()}, {str.begin(), str.end()});
}
T &operator[](size_t i) const
{ return operator()(i); }
......@@ -223,6 +266,11 @@ template<typename T, size_t ndim> class mav: public cmav<T, ndim>
static_assert(ndim==2, "ndim must be 2");
return d[str[0]*i + str[1]*j];
}
T &operator()(size_t i, size_t j, size_t k) const
{
static_assert(ndim==3, "ndim must be 3");
return d[str[0]*i + str[1]*j + str[2]*k];
}
using parent::shape;
using parent::stride;
using parent::size;
......@@ -240,6 +288,11 @@ template<typename T, size_t ndim> class mav: public cmav<T, ndim>
for (size_t i=0; i<shp[0]; ++i)
for (size_t j=0; j<shp[1]; ++j)
d[str[0]*i + str[1]*j] = val;
else if (ndim==3)
for (size_t i=0; i<shp[0]; ++i)
for (size_t j=0; j<shp[1]; ++j)
for (size_t k=0; k<shp[2]; ++k)
d[str[0]*i + str[1]*j + str[2]*k] = val;
}
using parent::conformable;
};
......
nifty_gridder/gridder_cxx.h
View file @ 9b34409f
......@@ -59,27 +59,6 @@ template<size_t ndim> void checkShape
template<typename T> inline T fmod1 (T v)
{ return v-floor(v); }
template<typename T, size_t ndim> class tmpStorage
{
private:
vector<T> d;
mav<T,ndim> mav_;
static size_t prod(const array<size_t,ndim> &shp)
{
size_t res=1;
for (auto v: shp) res*=v;
return res;
}
public:
tmpStorage(const array<size_t,ndim> &shp)
: d(prod(shp)), mav_(d.data(), shp) {}
mav<T,ndim> &getMav() { return mav_; }
void fill(const T & val)
{ std::fill(d.begin(), d.end(), val); }
};
//
// Start of real gridder functionality
//
......@@ -98,8 +77,8 @@ template<typename T> void complex2hartley
for (size_t v=0; v<nv; ++v)
{
size_t xv = (v==0) ? 0 : nv-v;
grid2(u,v) += T(0.5)*(grid( u, v).real()+grid( u, v).imag()+
grid(xu,xv).real()-grid(xu,xv).imag());
grid2(u,v) = T(0.5)*(grid( u, v).real()+grid( u, v).imag()+
grid(xu,xv).real()-grid(xu,xv).imag());
}
}
});
......@@ -456,8 +435,7 @@ class GridderConfig
const mav<T,2> &dirty) const
{
checkShape(grid.shape(), {nu,nv});
tmpStorage<T,2> tmpdat({nu,nv});
auto tmav = tmpdat.getMav();
mav<T,2> tmav({nu,nv});
hartley2_2D<T>(grid, tmav, nthreads);
grid2dirty_post(tmav, dirty);
}
......@@ -1080,8 +1058,7 @@ template<typename T, typename Serv> void x2dirty(
WgridHelper<Serv> hlp(gconf, srv, verbosity);
double dw = hlp.DW();
dirty.fill(0);
tmpStorage<complex<T>,2> grid_({gconf.Nu(),gconf.Nv()});
auto grid=grid_.getMav();
mav<complex<T>,2> grid({gconf.Nu(),gconf.Nv()});
while(hlp.advance()) // iterate over w planes
{
if (hlp.Nvis()==0) continue;
......@@ -1099,12 +1076,10 @@ template<typename T, typename Serv> void x2dirty(
<< " visibilities" << endl;
if (verbosity>0) cout << "Using " << gconf.Nthreads() << " threads" << endl;
tmpStorage<complex<T>,2> grid_({gconf.Nu(), gconf.Nv()});
auto grid=grid_.getMav();
grid_.fill(0.);
mav<complex<T>,2> grid({gconf.Nu(), gconf.Nv()});
grid.fill(0.);
x2grid_c(gconf, srv, grid);
tmpStorage<T,2> rgrid_(grid.shape());
auto rgrid=rgrid_.getMav();
mav<T,2> rgrid(grid.shape());
complex2hartley(grid, rgrid, gconf.Nthreads());
gconf.grid2dirty(rgrid, dirty);
}
......@@ -1121,15 +1096,13 @@ template<typename T, typename Serv> void dirty2x(
if (verbosity>0) cout << "Degridding using improved w-stacking" << endl;
WgridHelper<Serv> hlp(gconf, srv, verbosity);
double dw = hlp.DW();
tmpStorage<T,2> tdirty_({nx_dirty,ny_dirty});
auto tdirty=tdirty_.getMav();
mav<T,2> tdirty({nx_dirty,ny_dirty});
for (size_t i=0; i<nx_dirty; ++i)
for (size_t j=0; j<ny_dirty; ++j)
tdirty(i,j) = dirty(i,j);
// correct for w gridding etc.
apply_global_corrections(gconf, tdirty, ES_Kernel(gconf.Supp(), gconf.Ofactor(), nthreads), dw, true);
tmpStorage<complex<T>,2> grid_({gconf.Nu(),gconf.Nv()});
auto grid=grid_.getMav();
mav<complex<T>,2> grid({gconf.Nu(),gconf.Nv()});
while(hlp.advance()) // iterate over w planes
{
if (hlp.Nvis()==0) continue;
......@@ -1144,11 +1117,9 @@ template<typename T, typename Serv> void dirty2x(
<< " visibilities" << endl;
if (verbosity>0) cout << "Using " << gconf.Nthreads() << " threads" << endl;
tmpStorage<T,2> grid_({gconf.Nu(), gconf.Nv()});
auto grid=grid_.getMav();
mav<T,2> grid({gconf.Nu(), gconf.Nv()});
gconf.dirty2grid(dirty, grid);
tmpStorage<complex<T>,2> grid2_(grid.shape());
auto grid2=grid2_.getMav();
mav<complex<T>,2> grid2(grid.shape());
hartley2complex(grid, grid2, gconf.Nthreads());
grid2x_c(gconf, grid2, srv);
}
......
nifty_gridder/setup.py
View file @ 9b34409f
......@@ -19,12 +19,12 @@ python_module_link_args = []
if sys.platform == 'darwin':
import distutils.sysconfig
extra_compile_args += ['--std=c++11', '--stdlib=libc++', '-mmacosx-version-min=10.9']
extra_compile_args += ['--std=c++14', '--stdlib=libc++', '-mmacosx-version-min=10.9']
vars = distutils.sysconfig.get_config_vars()
vars['LDSHARED'] = vars['LDSHARED'].replace('-bundle', '')
python_module_link_args+=['-bundle']
else:
extra_compile_args += ['--std=c++11', '-march=native', '-O3', '-ffast-math']
extra_compile_args += ['--std=c++14', '-march=native', '-O3', '-ffast-math']
python_module_link_args += ['-march=native', '-Wl,-rpath,$ORIGIN', '-ffast-math']
# if you don't want debugging info, add "-s" to python_module_link_args
......
pypocketfft/setup.py
View file @ 9b34409f
......@@ -13,7 +13,7 @@ class _deferred_pybind11_include(object):
include_dirs = ['..', _deferred_pybind11_include(True),
_deferred_pybind11_include()]
extra_compile_args = ['--std=c++11', '-march=native', '-ffast-math', '-O3']
extra_compile_args = ['--std=c++14', '-march=native', '-ffast-math', '-O3']
python_module_link_args = []
define_macros = []
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment