Commit 69f14dc1 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

cleanups, ChangeLog

parent f6c28023
Pipeline #106732 passed with stages
in 26 minutes and 25 seconds
0.19.0:
- fft:
- new function `r2r_fftw`, which supports FFTW's halfcomplex storage scheme.
- C++ sources for FFT calculation now have their own subdirectory.
......
......@@ -333,10 +333,16 @@ template<typename T> class ConvolverPlan
// prefetching distance
static constexpr size_t pfdist=2;
template<size_t supp> void interpolx(const mav<T,3> &cube,
template<size_t supp> void interpolx(size_t supp_, const mav<T,3> &cube,
size_t itheta0, size_t iphi0, const mav<T,1> &theta, const mav<T,1> &phi,
const mav<T,1> &psi, mav<T,1> &signal) const
{
if constexpr (supp>=8)
if (supp_<=supp/2) return interpolx<supp/2>(supp_, cube, itheta0, iphi0, theta, phi, psi, signal);
if constexpr (supp>4)
if (supp_<supp) return interpolx<supp-1>(supp_, cube, itheta0, iphi0, theta, phi, psi, signal);
MR_assert(supp_==supp, "requested support ou of range");
MR_assert(cube.stride(2)==1, "last axis of cube must be contiguous");
MR_assert(phi.shape(0)==theta.shape(0), "array shape mismatch");
MR_assert(psi.shape(0)==theta.shape(0), "array shape mismatch");
......@@ -397,10 +403,16 @@ template<typename T> class ConvolverPlan
}
});
}
template<size_t supp> void deinterpolx(mav<T,3> &cube,
template<size_t supp> void deinterpolx(size_t supp_, mav<T,3> &cube,
size_t itheta0, size_t iphi0, const mav<T,1> &theta, const mav<T,1> &phi,
const mav<T,1> &psi, const mav<T,1> &signal) const
{
if constexpr (supp>=8)
if (supp_<=supp/2) return deinterpolx<supp/2>(supp_, cube, itheta0, iphi0, theta, phi, psi, signal);
if constexpr (supp>4)
if (supp_<supp) return deinterpolx<supp-1>(supp_, cube, itheta0, iphi0, theta, phi, psi, signal);
MR_assert(supp_==supp, "requested support ou of range");
MR_assert(cube.stride(2)==1, "last axis of cube must be contiguous");
MR_assert(phi.shape(0)==theta.shape(0), "array shape mismatch");
MR_assert(psi.shape(0)==theta.shape(0), "array shape mismatch");
......@@ -642,54 +654,16 @@ template<typename T> class ConvolverPlan
size_t iphi0, const mav<T,1> &theta, const mav<T,1> &phi,
const mav<T,1> &psi, mav<T,1> &signal) const
{
if constexpr(is_same<T,double>::value)
switch(kernel->support())
{
case 9: interpolx< 9>(cube, itheta0, iphi0, theta, phi, psi, signal); return;
case 10: interpolx<10>(cube, itheta0, iphi0, theta, phi, psi, signal); return;
case 11: interpolx<11>(cube, itheta0, iphi0, theta, phi, psi, signal); return;
case 12: interpolx<12>(cube, itheta0, iphi0, theta, phi, psi, signal); return;
case 13: interpolx<13>(cube, itheta0, iphi0, theta, phi, psi, signal); return;
case 14: interpolx<14>(cube, itheta0, iphi0, theta, phi, psi, signal); return;
case 15: interpolx<15>(cube, itheta0, iphi0, theta, phi, psi, signal); return;
case 16: interpolx<16>(cube, itheta0, iphi0, theta, phi, psi, signal); return;
}
switch(kernel->support())
{
case 4: interpolx<4>(cube, itheta0, iphi0, theta, phi, psi, signal); break;
case 5: interpolx<5>(cube, itheta0, iphi0, theta, phi, psi, signal); break;
case 6: interpolx<6>(cube, itheta0, iphi0, theta, phi, psi, signal); break;
case 7: interpolx<7>(cube, itheta0, iphi0, theta, phi, psi, signal); break;
case 8: interpolx<8>(cube, itheta0, iphi0, theta, phi, psi, signal); break;
default: MR_fail("must not happen");
}
constexpr size_t maxsupp = is_same<T, double>::value ? 16 : 8;
interpolx<maxsupp>(kernel->support(), cube, itheta0, iphi0, theta, phi, psi, signal);
}
void deinterpol(mav<T,3> &cube, size_t itheta0,
size_t iphi0, const mav<T,1> &theta, const mav<T,1> &phi,
const mav<T,1> &psi, const mav<T,1> &signal) const
{
if constexpr(is_same<T,double>::value)
switch(kernel->support())
{
case 9: deinterpolx< 9>(cube, itheta0, iphi0, theta, phi, psi, signal); return;
case 10: deinterpolx<10>(cube, itheta0, iphi0, theta, phi, psi, signal); return;
case 11: deinterpolx<11>(cube, itheta0, iphi0, theta, phi, psi, signal); return;
case 12: deinterpolx<12>(cube, itheta0, iphi0, theta, phi, psi, signal); return;
case 13: deinterpolx<13>(cube, itheta0, iphi0, theta, phi, psi, signal); return;
case 14: deinterpolx<14>(cube, itheta0, iphi0, theta, phi, psi, signal); return;
case 15: deinterpolx<15>(cube, itheta0, iphi0, theta, phi, psi, signal); return;
case 16: deinterpolx<16>(cube, itheta0, iphi0, theta, phi, psi, signal); return;
}
switch(kernel->support())
{
case 4: deinterpolx<4>(cube, itheta0, iphi0, theta, phi, psi, signal); break;
case 5: deinterpolx<5>(cube, itheta0, iphi0, theta, phi, psi, signal); break;
case 6: deinterpolx<6>(cube, itheta0, iphi0, theta, phi, psi, signal); break;
case 7: deinterpolx<7>(cube, itheta0, iphi0, theta, phi, psi, signal); break;
case 8: deinterpolx<8>(cube, itheta0, iphi0, theta, phi, psi, signal); break;
default: MR_fail("must not happen");
}
constexpr size_t maxsupp = is_same<T, double>::value ? 16 : 8;
deinterpolx<maxsupp>(kernel->support(), cube, itheta0, iphi0, theta, phi, psi, signal);
}
void updateSlm(mav<complex<T>,2> &vslm, const mav<complex<T>,2> &vblm,
......
......@@ -935,8 +935,14 @@ template<typename Tcalc, typename Tacc, typename Tms, typename Timg> class Param
}
template<size_t SUPP, bool wgrid> [[gnu::hot]] void x2grid_c_helper
(mav<complex<Tcalc>,2> &grid, size_t p0, double w0)
(size_t supp, mav<complex<Tcalc>,2> &grid, size_t p0, double w0)
{
if constexpr (SUPP>=8)
if (supp<=SUPP/2) return x2grid_c_helper<SUPP/2, wgrid>(supp, grid, p0, w0);
if constexpr (SUPP>4)
if (supp<SUPP) return x2grid_c_helper<SUPP-1, wgrid>(supp, grid, p0, w0);
MR_assert(supp==SUPP, "requested support ou of range");
vector<mutex> locks(nu);
execDynamic(ranges.size(), nthreads, wgrid ? SUPP : 1, [&](Scheduler &sched)
......@@ -1019,33 +1025,19 @@ auto ix = ix_+ranges.size()/2; if (ix>=ranges.size()) ix -=ranges.size();
size_t p0, double w0=-1)
{
checkShape(grid.shape(), {nu, nv});
if constexpr (is_same<Tacc, double>::value)
switch(supp)
{
case 9: x2grid_c_helper< 9, wgrid>(grid, p0, w0); return;
case 10: x2grid_c_helper<10, wgrid>(grid, p0, w0); return;
case 11: x2grid_c_helper<11, wgrid>(grid, p0, w0); return;
case 12: x2grid_c_helper<12, wgrid>(grid, p0, w0); return;
case 13: x2grid_c_helper<13, wgrid>(grid, p0, w0); return;
case 14: x2grid_c_helper<14, wgrid>(grid, p0, w0); return;
case 15: x2grid_c_helper<15, wgrid>(grid, p0, w0); return;
case 16: x2grid_c_helper<16, wgrid>(grid, p0, w0); return;
}
switch(supp)
{
case 4: x2grid_c_helper< 4, wgrid>(grid, p0, w0); return;
case 5: x2grid_c_helper< 5, wgrid>(grid, p0, w0); return;
case 6: x2grid_c_helper< 6, wgrid>(grid, p0, w0); return;
case 7: x2grid_c_helper< 7, wgrid>(grid, p0, w0); return;
case 8: x2grid_c_helper< 8, wgrid>(grid, p0, w0); return;
default: MR_fail("must not happen");
}
constexpr size_t maxsupp = is_same<Tacc, double>::value ? 16 : 8;
x2grid_c_helper<maxsupp, wgrid>(supp, grid, p0, w0);
}
template<size_t SUPP, bool wgrid> [[gnu::hot]] void grid2x_c_helper
(const mav<complex<Tcalc>,2> &grid, size_t p0, double w0)
(size_t supp, const mav<complex<Tcalc>,2> &grid, size_t p0, double w0)
{
if constexpr (SUPP>=8)
if (supp<=SUPP/2) return grid2x_c_helper<SUPP/2, wgrid>(supp, grid, p0, w0);
if constexpr (SUPP>4)
if (supp<SUPP) return grid2x_c_helper<SUPP-1, wgrid>(supp, grid, p0, w0);
MR_assert(supp==SUPP, "requested support ou of range");
// Loop over sampling points
execDynamic(ranges.size(), nthreads, wgrid ? SUPP : 1, [&](Scheduler &sched)
{
......@@ -1124,28 +1116,8 @@ auto ix = ix_+ranges.size()/2; if (ix>=ranges.size()) ix -=ranges.size();
size_t p0, double w0=-1)
{
checkShape(grid.shape(), {nu, nv});
if constexpr (is_same<Tcalc, double>::value)
switch(supp)
{
case 9: grid2x_c_helper< 9, wgrid>(grid, p0, w0); return;
case 10: grid2x_c_helper<10, wgrid>(grid, p0, w0); return;
case 11: grid2x_c_helper<11, wgrid>(grid, p0, w0); return;
case 12: grid2x_c_helper<12, wgrid>(grid, p0, w0); return;
case 13: grid2x_c_helper<13, wgrid>(grid, p0, w0); return;
case 14: grid2x_c_helper<14, wgrid>(grid, p0, w0); return;
case 15: grid2x_c_helper<15, wgrid>(grid, p0, w0); return;
case 16: grid2x_c_helper<16, wgrid>(grid, p0, w0); return;
}
switch(supp)
{
case 4: grid2x_c_helper< 4, wgrid>(grid, p0, w0); return;
case 5: grid2x_c_helper< 5, wgrid>(grid, p0, w0); return;
case 6: grid2x_c_helper< 6, wgrid>(grid, p0, w0); return;
case 7: grid2x_c_helper< 7, wgrid>(grid, p0, w0); return;
case 8: grid2x_c_helper< 8, wgrid>(grid, p0, w0); return;
default: MR_fail("must not happen");
}
constexpr size_t maxsupp = is_same<Tcalc, double>::value ? 16 : 8;
grid2x_c_helper<maxsupp, wgrid>(supp, grid, p0, w0);
}
void apply_global_corrections(mav<Timg,2> &dirty)
......
Supports Markdown
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