separate out 1D FFTs

Makefile.am

@@ -20,6 +20,7 @@ libmrutil_la_SOURCES = \
   mr_util/error_handling.h \
   mr_util/space_filling.cc \
   mr_util/space_filling.h \
+  mr_util/fft1d.h \
   mr_util/fft.h \
   mr_util/gl_integrator.h \
   mr_util/mav.h \

libsharp2/sharp.cc

@@ -29,7 +29,7 @@
 #include <atomic>
 #include <memory>
 #include "mr_util/math_utils.h"
-#include "mr_util/fft.h"
+#include "mr_util/fft1d.h"
 #include "libsharp2/sharp_internal.h"
 #include "libsharp2/sharp_almhelpers.h"
 #include "libsharp2/sharp_geomhelpers.h"
@@ -80,7 +80,7 @@ struct ringhelper
   double phi0_;
   vector<dcmplx> shiftarr;
   size_t s_shift;
-  unique_ptr<mr::detail_fft::rfftp<double>> plan;
+  unique_ptr<pocketfft_r<double>> plan;
   size_t length;
   bool norot;
   ringhelper() : length(0) {}
@@ -101,7 +101,7 @@ struct ringhelper
       }
     if (nph!=length)
       {
-      plan.reset(new mr::detail_fft::rfftp<double>(nph));
+      plan.reset(new pocketfft_r<double>(nph));
       length=nph;
       }
     }

libsharp2/sharp_geomhelpers.cc

@@ -29,7 +29,7 @@
 #include <vector>
 #include "libsharp2/sharp_geomhelpers.h"
 #include "mr_util/gl_integrator.h"
-#include "mr_util/fft.h"
+#include "mr_util/fft1d.h"
 #include "mr_util/error_handling.h"
 #include "mr_util/math_utils.h"
 
@@ -234,8 +234,10 @@ unique_ptr<sharp_geom_info> sharp_make_fejer1_geom_info (size_t nrings, size_t p
     weight[2*k  ]=2./(1.-4.*k*k)*sin((k*pi)/nrings);
     }
   if ((nrings&1)==0) weight[nrings-1]=0.;
-  auto tmp = fmav(weight.data(),{nrings});
-  mr::r2r_fftpack(tmp, tmp, {0}, false, false, 1.);
+  {
+  pocketfft_r<double> plan(nrings);
+  plan.exec(weight.data(), 1., false);
+  }
 
   for (size_t m=0; m<(nrings+1)/2; ++m)
     {
@@ -265,8 +267,10 @@ unique_ptr<sharp_geom_info> sharp_make_cc_geom_info (size_t nrings, size_t pprin
   for (size_t k=1; k<=(n/2-1); ++k)
     weight[2*k-1]=2./(1.-4.*k*k) + dw;
   weight[2*(n/2)-1]=(n-3.)/(2*(n/2)-1) -1. -dw*((2-(n&1))*n-1);
-  auto tmp = fmav(weight.data(),{n});
-  mr::r2r_fftpack(tmp, tmp, {0}, false, false, 1.);
+  {
+  pocketfft_r<double> plan(n);
+  plan.exec(weight.data(), 1., false);
+  }
   weight[n]=weight[0];
 
   for (size_t m=0; m<(nrings+1)/2; ++m)
@@ -290,8 +294,10 @@ static vector<double> get_dh_weights(size_t nrings)
   for (size_t k=1; k<=(nrings/2-1); ++k)
     weight[2*k-1]=2./(1.-4.*k*k);
   weight[2*(nrings/2)-1]=(nrings-3.)/(2*(nrings/2)-1) -1.;
-  auto tmp = fmav(weight.data(),{nrings});
-  mr::r2r_fftpack(tmp, tmp, {0}, false, false, 1.);
+  {
+  pocketfft_r<double> plan(nrings);
+  plan.exec(weight.data(), 1., false);
+  }
   return weight;
   }
 

pyHealpix/setup.py

@@ -40,7 +40,7 @@ def get_extension_modules():
                                '../libsharp2/sharp.cc', '../libsharp2/sharp_core.cc', '../libsharp2/sharp_geomhelpers.cc',
                                '../libsharp2/sharp_almhelpers.cc','../libsharp2/sharp_ylmgen.cc','../Healpix_cxx/healpix_base.cc',
                                '../Healpix_cxx/healpix_tables.cc'],
-                      depends=['../mr_util/fft.h', '../mr_util/mav.h', '../mr_util/threading.h',
+                      depends=['../mr_util/fft1d.h', '../mr_util/mav.h', '../mr_util/threading.h',
                                '../mr_util/aligned_array.h', '../mr_util/simd.h',
                                '../mr_util/cmplx.h', '../mr_util/unity_roots.h', '../mr_util/error_handling.h',
                                'setup.py'],

pypocketfft/pypocketfft.cc

@@ -380,7 +380,7 @@ PyObject * good_size(PyObject * /*self*/, PyObject * args)
   const auto n = static_cast<size_t>(n_);
   using namespace mr::detail_fft;
   return PyLong_FromSize_t(
-    real ? util::good_size_real(n) : util::good_size_cmplx(n));
+    real ? util1d::good_size_real(n) : util1d::good_size_cmplx(n));
   }
 
 const char *pypocketfft_DS = R"""(Fast Fourier and Hartley transforms.

pypocketfft/setup.py

@@ -37,7 +37,7 @@ def get_extension_modules():
     return [Extension('pypocketfft',
                       language='c++',
                       sources=['pypocketfft.cc','../mr_util/threading.cc'],
-                      depends=['../mr_util/fft.h', '../mr_util/mav.h', '../mr_util/threading.h',
+                      depends=['../mr_util/fft.h', '../mr_util/fft1d.h', '../mr_util/mav.h', '../mr_util/threading.h',
                                '../mr_util/aligned_array.h', '../mr_util/simd.h',
                                '../mr_util/cmplx.h', '../mr_util/unity_roots.h', '../mr_util/error_handling.h',
                                'setup.py'],

pysharp/setup.py

@@ -37,7 +37,7 @@ def get_extension_modules():
                       sources=['pysharp.cc','../mr_util/threading.cc','../mr_util/string_utils.cc',
                                '../libsharp2/sharp.cc', '../libsharp2/sharp_core.cc', '../libsharp2/sharp_geomhelpers.cc',
                                '../libsharp2/sharp_almhelpers.cc','../libsharp2/sharp_ylmgen.cc'],
-                      depends=['../mr_util/string_utils.h', '../mr_util/fft.h', '../mr_util/mav.h', '../mr_util/threading.h',
+                      depends=['../mr_util/string_utils.h', '../mr_util/fft1d.h', '../mr_util/mav.h', '../mr_util/threading.h',
                                '../mr_util/aligned_array.h', '../mr_util/simd.h',
                                '../mr_util/cmplx.h', '../mr_util/unity_roots.h', '../mr_util/error_handling.h',
                                'setup.py'],