Commit 6d591c2b authored by Martin Reinecke's avatar Martin Reinecke

seems to work

parent 1ef56637
......@@ -34,6 +34,7 @@
#include "ducc0/infra/threading.h"
#include "ducc0/infra/useful_macros.h"
#include "ducc0/infra/mav.h"
#include "ducc0/infra/simd.h"
#include "ducc0/math/es_kernel.h"
namespace ducc0 {
......@@ -536,8 +537,9 @@ template<typename T, typename T2=complex<T>> class Helper
public:
const T2 *p0r;
T2 *p0w;
T kernel[64] DUCC0_ALIGNED(64);
T kernel[64];
static constexpr size_t vlen=native_simd<T>::size();
native_simd<T> simd_kernel[64/vlen];
Helper(const GridderConfig &gconf_, const T2 *grid_r_, T2 *grid_w_,
vector<std::mutex> &locks_, double w0_=-1, double dw_=-1)
......@@ -551,7 +553,7 @@ template<typename T, typename T2=complex<T>> class Helper
w0(w0_),
xdw(T(1)/dw_),
nexp(2*supp + do_w_gridding),
nvecs(vlen*((nexp+vlen-1)/vlen)),
nvecs((nexp+vlen-1)/vlen),
locks(locks_)
{}
~Helper() { if (grid_w) dump(); }
......@@ -571,16 +573,13 @@ template<typename T, typename T2=complex<T>> class Helper
kernel[i+supp] = T(y0+i*xsupp);
}
if (do_w_gridding)
kernel[2*supp] = min(T(1), T(xdw*xsupp*abs(w0-in.w)));
for (size_t i=nexp; i<nvecs; ++i)
kernel[2*supp] = T(xdw*xsupp*abs(w0-in.w));
for (size_t i=nexp; i<nvecs*vlen; ++i)
kernel[i]=0;
memcpy(simd_kernel, kernel, nvecs*vlen*sizeof(T));
for (size_t i=0; i<nvecs; ++i)
{
kernel[i] = T(1) - kernel[i]*kernel[i];
kernel[i] = (kernel[i]<0) ? T(-200.) : beta*(sqrt(kernel[i])-T(1));
}
for (size_t i=0; i<nvecs; ++i)
kernel[i] = exp(kernel[i]);
simd_kernel[i] = esk(simd_kernel[i], beta);
memcpy(kernel, simd_kernel, nvecs*vlen*sizeof(T));
if ((iu0<bu0) || (iv0<bv0) || (iu0+supp>bu0+su) || (iv0+supp>bv0+sv))
{
if (grid_w) { dump(); fill(wbuf.begin(), wbuf.end(), T(0)); }
......
......@@ -25,6 +25,7 @@
#include <cmath>
#include <vector>
#include <array>
#include "ducc0/infra/simd.h"
#include "ducc0/infra/error_handling.h"
#include "ducc0/infra/threading.h"
#include "ducc0/math/constants.h"
......@@ -36,6 +37,28 @@ namespace detail_es_kernel {
using namespace std;
template<typename T> T esk (T v, T beta)
{
auto tmp = (1-v)*(1+v);
auto tmp2 = tmp>=0;
return tmp2*exp(T(beta)*(sqrt(tmp*tmp2)-1));
}
template<typename T> class exponator
{
public:
T operator()(T val) { return exp(val); }
};
template<typename T> native_simd<T> esk (native_simd<T> v, T beta)
{
auto tmp = (T(1)-v)*(T(1)+v);
auto mask = tmp<T(0);
where(mask,tmp)*=T(0);
auto res = (beta*(sqrt(tmp)-T(1))).apply(exponator<T>());
where (mask,res)*=T(0);
return res;
}
/* This class implements the "exponential of semicircle" gridding kernel
described in https://arxiv.org/abs/1808.06736 */
class ES_Kernel
......@@ -62,11 +85,11 @@ class ES_Kernel
: ES_Kernel(supp_, 2., nthreads){}
template<typename T> T operator()(T v) const
{
auto tmp = (1-v)*(1+v);
auto tmp2 = tmp>=0;
return tmp2*exp(T(beta)*(sqrt(tmp*tmp2)-1));
}
{ return esk(v, T(beta)); }
template<typename T> native_simd<T> operator()(native_simd<T> v)
{ return esk(v, T(beta)); }
/* Compute correction factors for the ES gridding kernel
This implementation follows eqs. (3.8) to (3.10) of Barnett et al. 2018 */
double corfac(double v) const
......@@ -138,6 +161,7 @@ class ES_Kernel
}
using detail_es_kernel::esk;
using detail_es_kernel::ES_Kernel;
}
......
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