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,6 +10,9 @@ 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*/) {}
@@ -63,11 +66,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,
@@ -131,11 +138,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,

include/finufft/finufft_core.h

@@ -95,6 +95,17 @@
 // values for M = nj (also nk in type 3)...
 inline constexpr BIGINT MAX_NU_PTS = BIGINT(1e14);
 
+// MR: In the longer term I suggest to move
+// away from M_PI, which was never part of the standard.
+// Perhaps a constexpr pi in the namespace finufft, or a constexpr finufft_pi
+// if no namespaces are used?
+// In C++20 these constants will be part of the language, and the problem will go away.
+#ifndef M_PI // Windows apparently doesn't have this const
+#define M_PI 3.14159265358979329
+#endif
+#define M_1_2PI 0.159154943091895336
+#define M_2PI   6.28318530717958648
+
 // ----- 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 +149,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;
@@ -147,49 +159,53 @@
 
   int type;                // transform type (Rokhlin naming): 1,2 or 3
   int dim;                 // overall dimension: 1,2 or 3
   int ntrans;              // how many transforms to do at once (vector or "many" mode)
   BIGINT nj;               // num of NU pts in type 1,2 (for type 3, num input x pts)
   BIGINT nk;               // number of NU freq pts (type 3 only)
   TF tol;                  // relative user tolerance
   int batchSize;           // # strength vectors to group together for FFTW, etc
   int nbatch;              // how many batches done to cover all ntrans vectors
 
   BIGINT ms;               // number of modes in x (1) dir (historical CMCL name) = N1
   BIGINT mt;               // number of modes in y (2) direction = N2
   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
+  TC *fwBatch = nullptr;   // (batches of) fine grid(s) for the FFT to plan
+                           // & act on. Usually the largest working array.
+                           // Unfortunately this can't be a vector, since
+                           // in the presence of FFTW it must be allocated
+                           // by fftw_malloc.
 
   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 +220,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

include/finufft/spreadinterp.h

@@ -1,12 +1,13 @@
 // Defines interface to spreading/interpolation code.
 
-// Devnotes: see defs.h for definition of MAX_NSPREAD (as of 9/24/18).
+// Devnotes: see finufft_core.h for definition of MAX_NSPREAD
 // RESCALE macro moved to spreadinterp.cpp, 7/15/20.
 // finufft_spread_opts renamed 6/7/22.
 
 #ifndef SPREADINTERP_H
 #define SPREADINTERP_H
 
+#include <finufft/finufft_core.h>
 #include <finufft_spread_opts.h>
 
 /* Bitwise debugging timing flag (TF) defs; see finufft_spread_opts.flags.