Follow-up to the big reorg PR

CMakeLists.txt

@@ -169,7 +169,7 @@ if(FINUFFT_USE_CPU)
   set(XTL_VERSION 0.7.7)
   set(XSIMD_VERSION 13.0.0)
   # using latest ducc0 version for now as it fixes MacOS GCC build
-  set(DUCC0_VERSION b0beb85e03982344a31ebb119758d7aa3ef5d362)
+  set(DUCC0_VERSION ducc0_0_35_0)
   set(FINUFFT_FFTW_LIBRARIES)
   include(cmake/setupXSIMD.cmake)
   if(FINUFFT_USE_DUCC0)
@@ -258,7 +258,7 @@ if(FINUFFT_USE_CPU)
       contrib/legendre_rule_fast.cpp
       src/fft.cpp
       src/finufft_core.cpp
-      src/simpleinterfaces.cpp
+      src/c_interface.cpp
       fortran/finufftfort.cpp)
   else()
     add_library(
@@ -268,7 +268,7 @@ if(FINUFFT_USE_CPU)
       contrib/legendre_rule_fast.cpp
       src/fft.cpp
       src/finufft_core.cpp
-      src/simpleinterfaces.cpp
+      src/c_interface.cpp
       fortran/finufftfort.cpp)
   endif()
   set_finufft_options(finufft)

include/finufft.h

@@ -2,8 +2,7 @@
 
 // This contains both single and double precision user-facing commands.
 // "macro-safe" rewrite, including the plan object, Barnett 5/21/22-6/7/22.
-// They will clobber any prior macros starting FINUFFT*, so in the lib/test
-// sources finufft.h must be included before defs.h
+// They will clobber any prior macros starting FINUFFT*.
 
 /* Devnotes.
    A) Two precisions done by including the "either precision" headers twice.

include/finufft/dirft.h

@@ -1,22 +1,36 @@
 #ifndef DIRFT_H
 #define DIRFT_H
 
-#include <finufft/defs.h>
+#include <finufft/finufft_core.h>
 
-void dirft1d1(BIGINT nj, FLT *x, CPX *c, int isign, BIGINT ms, CPX *f);
-void dirft1d2(BIGINT nj, FLT *x, CPX *c, int iflag, BIGINT ms, CPX *f);
-void dirft1d3(BIGINT nj, FLT *x, CPX *c, int iflag, BIGINT nk, FLT *s, CPX *f);
+template<typename T>
+void dirft1d1(BIGINT nj, T *x, std::complex<T> *c, int isign, BIGINT ms,
+              std::complex<T> *f);
+template<typename T>
+void dirft1d2(BIGINT nj, T *x, std::complex<T> *c, int iflag, BIGINT ms,
+              std::complex<T> *f);
+template<typename T>
+void dirft1d3(BIGINT nj, T *x, std::complex<T> *c, int iflag, BIGINT nk, T *s,
+              std::complex<T> *f);
 
-void dirft2d1(BIGINT nj, FLT *x, FLT *y, CPX *c, int iflag, BIGINT ms, BIGINT mt, CPX *f);
-void dirft2d2(BIGINT nj, FLT *x, FLT *y, CPX *c, int iflag, BIGINT ms, BIGINT mt, CPX *f);
-void dirft2d3(BIGINT nj, FLT *x, FLT *y, CPX *c, int iflag, BIGINT nk, FLT *s, FLT *t,
-              CPX *f);
+template<typename T>
+void dirft2d1(BIGINT nj, T *x, T *y, std::complex<T> *c, int iflag, BIGINT ms, BIGINT mt,
+              std::complex<T> *f);
+template<typename T>
+void dirft2d2(BIGINT nj, T *x, T *y, std::complex<T> *c, int iflag, BIGINT ms, BIGINT mt,
+              std::complex<T> *f);
+template<typename T>
+void dirft2d3(BIGINT nj, T *x, T *y, std::complex<T> *c, int iflag, BIGINT nk, T *s, T *t,
+              std::complex<T> *f);
 
-void dirft3d1(BIGINT nj, FLT *x, FLT *y, FLT *z, CPX *c, int iflag, BIGINT ms, BIGINT mt,
-              BIGINT mu, CPX *f);
-void dirft3d2(BIGINT nj, FLT *x, FLT *y, FLT *z, CPX *c, int iflag, BIGINT ms, BIGINT mt,
-              BIGINT mu, CPX *f);
-void dirft3d3(BIGINT nj, FLT *x, FLT *y, FLT *z, CPX *c, int iflag, BIGINT nk, FLT *s,
-              FLT *t, FLT *u, CPX *f);
+template<typename T>
+void dirft3d1(BIGINT nj, T *x, T *y, T *z, std::complex<T> *c, int iflag, BIGINT ms,
+              BIGINT mt, BIGINT mu, std::complex<T> *f);
+template<typename T>
+void dirft3d2(BIGINT nj, T *x, T *y, T *z, std::complex<T> *c, int iflag, BIGINT ms,
+              BIGINT mt, BIGINT mu, std::complex<T> *f);
+template<typename T>
+void dirft3d3(BIGINT nj, T *x, T *y, T *z, std::complex<T> *c, int iflag, BIGINT nk, T *s,
+              T *t, T *u, std::complex<T> *f);
 
 #endif

include/finufft/fft.h

@@ -10,13 +10,12 @@ template<typename T> class Finufft_FFT_plan {
 public:
   [[maybe_unused]] Finufft_FFT_plan(void (*)(void *) = nullptr,
                                     void (*)(void *) = nullptr, void * = nullptr) {}
+  // deleting these operations to be consistent with the FFTW plans (seel below)
+  Finufft_FFT_plan(const Finufft_FFT_plan &)            = delete;
+  Finufft_FFT_plan &operator=(const Finufft_FFT_plan &) = delete;
   [[maybe_unused]] void plan(const std::vector<int> & /*dims*/, size_t /*batchSize*/,
                              std::complex<T> * /*ptr*/, int /*sign*/, int /*options*/,
                              int /*nthreads*/) {}
-  [[maybe_unused]] static std::complex<T> *alloc_complex(size_t N) {
-    return new std::complex<T>[N];
-  }
-  [[maybe_unused]] static void free(std::complex<T> *ptr) { delete[] ptr; }
 
   [[maybe_unused]] static void forget_wisdom() {}
   [[maybe_unused]] static void cleanup() {}
@@ -63,11 +62,15 @@ template<> struct Finufft_FFT_plan<float> {
 #endif
     unlock();
   }
+  // we have raw pointers in the object (the FFTW plan).
+  // If we allow copying those, we end up destroying the plans multiple times.
+  Finufft_FFT_plan(const Finufft_FFT_plan &) = delete;
   [[maybe_unused]] ~Finufft_FFT_plan() {
     lock();
     fftwf_destroy_plan(plan_);
     unlock();
   }
+  Finufft_FFT_plan &operator=(const Finufft_FFT_plan &) = delete;
 
   void plan
       [[maybe_unused]] (const std::vector<int> &dims, size_t batchSize,
@@ -84,12 +87,6 @@ template<> struct Finufft_FFT_plan<float> {
                                 1, int(nf), sign, unsigned(options));
     unlock();
   }
-  static std::complex<float> *alloc_complex [[maybe_unused]] (size_t N) {
-    return reinterpret_cast<std::complex<float> *>(fftwf_alloc_complex(N));
-  }
-  static void free [[maybe_unused]] (std::complex<float> *ptr) {
-    if (ptr) fftwf_free(reinterpret_cast<fftwf_complex *>(ptr));
-  }
   void execute [[maybe_unused]] () { fftwf_execute(plan_); }
 
   static void forget_wisdom [[maybe_unused]] () { fftwf_forget_wisdom(); }
@@ -131,11 +128,13 @@ template<> struct Finufft_FFT_plan<double> {
 #endif
     unlock();
   }
+  Finufft_FFT_plan(const Finufft_FFT_plan &) = delete;
   [[maybe_unused]] ~Finufft_FFT_plan() {
     lock();
     fftw_destroy_plan(plan_);
     unlock();
   }
+  Finufft_FFT_plan &operator=(const Finufft_FFT_plan &) = delete;
 
   void plan
       [[maybe_unused]] (const std::vector<int> &dims, size_t batchSize,
@@ -152,12 +151,6 @@ template<> struct Finufft_FFT_plan<double> {
                                sign, unsigned(options));
     unlock();
   }
-  static std::complex<double> *alloc_complex [[maybe_unused]] (size_t N) {
-    return reinterpret_cast<std::complex<double> *>(fftw_alloc_complex(N));
-  }
-  static void free [[maybe_unused]] (std::complex<double> *ptr) {
-    fftw_free(reinterpret_cast<fftw_complex *>(ptr));
-  }
   void execute [[maybe_unused]] () { fftw_execute(plan_); }
 
   static void forget_wisdom [[maybe_unused]] () { fftw_forget_wisdom(); }

include/finufft/finufft_core.h

@@ -61,6 +61,7 @@
 
 #include <finufft_errors.h>
 #include <memory>
+#include <xsimd/xsimd.hpp>
 
 // All indexing in library that potentially can exceed 2^31 uses 64-bit signed.
 // This includes all calling arguments (eg M,N) that could be huge someday.
@@ -95,6 +96,10 @@
 // values for M = nj (also nk in type 3)...
 inline constexpr BIGINT MAX_NU_PTS = BIGINT(1e14);
 
+// We define our own PI here because M_PI is not actually part of standard C++
+inline constexpr double PI      = 3.14159265358979329;
+inline constexpr double INV_2PI = 0.159154943091895336;
+
 // ----- OpenMP macros which also work when omp not present -----
 // Allows compile-time switch off of openmp, so compilation without any openmp
 // is done (Note: _OPENMP is automatically set by -fopenmp compile flag)
@@ -138,7 +143,8 @@
 
   // These default and delete specifications just state the obvious,
   // but are here to silence compiler warnings.
-  FINUFFT_PLAN_T() = default;
+  FINUFFT_PLAN_T(int type, int dim, const BIGINT *n_modes, int iflag, int ntrans, TF tol,
+                 finufft_opts *opts, int &ier);
   // Copy construction and assignent are already deleted implicitly
   // because of the unique_ptr member.
   FINUFFT_PLAN_T(const FINUFFT_PLAN_T &)            = delete;
@@ -159,37 +165,40 @@
   BIGINT mu;               // number of modes in z (3) direction = N3
   BIGINT N;                // total # modes (prod of above three)
 
-  BIGINT nf1;              // size of internal fine grid in x (1) direction
-  BIGINT nf2;              // " y (2)
-  BIGINT nf3;              // " z (3)
-  BIGINT nf;               // total # fine grid points (product of the above three)
+  BIGINT nf1 = 1;          // size of internal fine grid in x (1) direction
+  BIGINT nf2 = 1;          // " y (2)
+  BIGINT nf3 = 1;          // " z (3)
+  BIGINT nf  = 1;          // total # fine grid points (product of the above three)
 
   int fftSign;             // sign in exponential for NUFFT defn, guaranteed to be +-1
 
   std::vector<TF> phiHat1; // FT of kernel in t1,2, on x-axis mode grid
   std::vector<TF> phiHat2; // " y-axis.
   std::vector<TF> phiHat3; // " z-axis.
 
-  TC *fwBatch = nullptr;   // (batches of) fine grid(s) for FFTW to plan
-                           // & act on. Usually the largest working array
+  std::vector<TC, xsimd::aligned_allocator<TC, 64>> fwBatch; // (batches of) fine grid(s)
+                                                             // for the FFT to plan & act
+                                                             // on. Usually the largest
+                                                             // working array.
 
   std::vector<BIGINT> sortIndices; // precomputed NU pt permutation, speeds spread/interp
   bool didSort;                    // whether binsorting used (false: identity perm used)
 
   TF *X = nullptr, *Y = nullptr, *Z = nullptr; // for t1,2: ptr to user-supplied NU pts
-                                               // (no new allocs). for t3: allocated as
-                                               // "primed" (scaled) src pts x'_j, etc
+                                               // (no new allocs). for t3: points to
+                                               // "primed" (scaled) Xp, Yp, Zp
 
   // type 3 specific
   TF *S = nullptr, *T = nullptr, *U = nullptr; // pointers to user's target NU pts arrays
                                                // (no new allocs)
   std::vector<TC> prephase;                    // pre-phase, for all input NU pts
   std::vector<TC> deconv;     // reciprocal of kernel FT, phase, all output NU pts
   std::vector<TC> CpBatch;    // working array of prephased strengths
-  std::vector<TF> Sp, Tp, Up; // internal primed targs (s'_k, etc),
-                              // allocated
+  std::vector<TF> Xp, Yp, Zp; // internal primed NU points (x'_j, etc)
+  std::vector<TF> Sp, Tp, Up; // internal primed targs (s'_k, etc)
   type3params<TF> t3P;        // groups together type 3 shift, scale, phase, parameters
-  FINUFFT_PLAN_T<TF> *innerT2plan = nullptr; // ptr used for type 2 in step 2 of type 3
+  std::unique_ptr<FINUFFT_PLAN_T<TF>> innerT2plan; // ptr used for type 2 in step 2 of
+                                                   // type 3
 
   // other internal structs
   std::unique_ptr<Finufft_FFT_plan<TF>> fftPlan;
@@ -204,10 +213,5 @@
 template<typename TF>
 int finufft_makeplan_t(int type, int dim, const BIGINT *n_modes, int iflag, int ntrans,
                        TF tol, FINUFFT_PLAN_T<TF> **pp, finufft_opts *opts);
-template<typename TF>
-int finufft_setpts_t(FINUFFT_PLAN_T<TF> *p, BIGINT nj, TF *xj, TF *yj, TF *zj, BIGINT nk,
-                     TF *s, TF *t, TF *u);
-template<typename TF>
-int finufft_execute_t(FINUFFT_PLAN_T<TF> *p, std::complex<TF> *cj, std::complex<TF> *fk);
 
 #endif // FINUFFT_CORE_H