fix: avoid unnecessary 2GB allocation in makeplan for FFTW_ESTIMATE

Since d4d08a4 (PR #633, between v2.3.0 and v2.4.1), makeplan allocates
a temporary std::vector for FFTW planning that is never read. For large
3D grids (e.g. 250^3 with upsamp=2.0 → 500^3 = 2GB), this causes
~250ms of page faults for a buffer that neither FFTW_ESTIMATE nor
DUCC0 touches.

Fix: remove the ptr parameter from plan() and let each backend manage
its own buffer internally. For FFTW_ESTIMATE (the default), a small
dummy is used since FFTW never touches the buffer. For FFTW_MEASURE
etc., a full buffer is allocated for trial FFTs. The DUCC0 stub
remains a no-op.

Author: DiamonDinoia

include/finufft/fft.h

@@ -14,8 +14,7 @@ template<typename T> class Finufft_FFT_plan {
   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*/) {}
+                             int /*sign*/, int /*options*/, int /*nthreads*/) {}
 
   [[maybe_unused]] static void forget_wisdom() {}
   [[maybe_unused]] static void cleanup() {}
@@ -28,6 +27,7 @@ template<typename T> class Finufft_FFT_plan {
 #include <complex>
 #include <fftw3.h> // (after complex) needed so can typedef FFTW_CPX
 //clang-format on
+#include <array>
 #include <mutex>
 
 template<typename T> class Finufft_FFT_plan {};
@@ -73,11 +73,17 @@ template<> struct Finufft_FFT_plan<float> {
   }
   Finufft_FFT_plan &operator=(const Finufft_FFT_plan &) = delete;
 
-  void plan
-      [[maybe_unused]] (const std::vector<int> &dims, size_t batchSize,
-                        std::complex<float> *ptr, int sign, int options, int nthreads) {
+  void plan [[maybe_unused]] (const std::vector<int> &dims, size_t batchSize, int sign,
+                              int options, int nthreads) {
     uint64_t nf = 1;
     for (auto i : dims) nf *= i;
+    // FFTW_ESTIMATE never touches the buffer; FFTW_MEASURE etc. run trial FFTs.
+    // Use a 1-element dummy for ESTIMATE, full aligned buffer otherwise.
+    using cpxf = std::complex<float>;
+    std::array<cpxf, 1> dummy{};
+    std::vector<cpxf> buf(options & FFTW_ESTIMATE ? 0 : nf * batchSize);
+    auto *ptr =
+        reinterpret_cast<fftwf_complex *>(buf.empty() ? dummy.data() : buf.data());
     lock();
     // Destroy existing plans before creating new ones (handles re-planning)
     if (plan_) {
@@ -91,14 +97,12 @@ template<> struct Finufft_FFT_plan<float> {
 #ifdef _OPENMP
     fftwf_plan_with_nthreads(nthreads);
 #endif
-    plan_     = fftwf_plan_many_dft(int(dims.size()), dims.data(), int(batchSize),
-                                    reinterpret_cast<fftwf_complex *>(ptr), nullptr, 1,
-                                    int(nf), reinterpret_cast<fftwf_complex *>(ptr), nullptr,
-                                    1, int(nf), sign, unsigned(options));
-    plan_adj_ = fftwf_plan_many_dft(int(dims.size()), dims.data(), int(batchSize),
-                                    reinterpret_cast<fftwf_complex *>(ptr), nullptr, 1,
-                                    int(nf), reinterpret_cast<fftwf_complex *>(ptr),
-                                    nullptr, 1, int(nf), -sign, unsigned(options));
+    plan_     = fftwf_plan_many_dft(int(dims.size()), dims.data(), int(batchSize), ptr,
+                                    nullptr, 1, int(nf), ptr, nullptr, 1, int(nf), sign,
+                                    unsigned(options));
+    plan_adj_ = fftwf_plan_many_dft(int(dims.size()), dims.data(), int(batchSize), ptr,
+                                    nullptr, 1, int(nf), ptr, nullptr, 1, int(nf), -sign,
+                                    unsigned(options));
     unlock();
   }
   void execute [[maybe_unused]] (std::complex<float> *data) const {
@@ -158,11 +162,16 @@ template<> struct Finufft_FFT_plan<double> {
   }
   Finufft_FFT_plan &operator=(const Finufft_FFT_plan &) = delete;
 
-  void plan
-      [[maybe_unused]] (const std::vector<int> &dims, size_t batchSize,
-                        std::complex<double> *ptr, int sign, int options, int nthreads) {
+  void plan [[maybe_unused]] (const std::vector<int> &dims, size_t batchSize, int sign,
+                              int options, int nthreads) {
     uint64_t nf = 1;
     for (auto i : dims) nf *= i;
+    // FFTW_ESTIMATE never touches the buffer; FFTW_MEASURE etc. run trial FFTs.
+    // Use a 1-element dummy for ESTIMATE, full aligned buffer otherwise.
+    using cpxd = std::complex<double>;
+    std::array<cpxd, 1> dummy{};
+    std::vector<cpxd> buf(options & FFTW_ESTIMATE ? 0 : nf * batchSize);
+    auto *ptr = reinterpret_cast<fftw_complex *>(buf.empty() ? dummy.data() : buf.data());
     lock();
     // Destroy existing plans before creating new ones (handles re-planning)
     if (plan_) {
@@ -176,14 +185,12 @@ template<> struct Finufft_FFT_plan<double> {
 #ifdef _OPENMP
     fftw_plan_with_nthreads(nthreads);
 #endif
-    plan_     = fftw_plan_many_dft(int(dims.size()), dims.data(), int(batchSize),
-                                   reinterpret_cast<fftw_complex *>(ptr), nullptr, 1, int(nf),
-                                   reinterpret_cast<fftw_complex *>(ptr), nullptr, 1, int(nf),
-                                   sign, unsigned(options));
-    plan_adj_ = fftw_plan_many_dft(int(dims.size()), dims.data(), int(batchSize),
-                                   reinterpret_cast<fftw_complex *>(ptr), nullptr, 1,
-                                   int(nf), reinterpret_cast<fftw_complex *>(ptr),
-                                   nullptr, 1, int(nf), -sign, unsigned(options));
+    plan_ =
+        fftw_plan_many_dft(int(dims.size()), dims.data(), int(batchSize), ptr, nullptr, 1,
+                           int(nf), ptr, nullptr, 1, int(nf), sign, unsigned(options));
+    plan_adj_ =
+        fftw_plan_many_dft(int(dims.size()), dims.data(), int(batchSize), ptr, nullptr, 1,
+                           int(nf), ptr, nullptr, 1, int(nf), -sign, unsigned(options));
     unlock();
   }
   void execute [[maybe_unused]] (std::complex<double> *data) const {

src/finufft_core.cpp

@@ -764,8 +764,7 @@ template<typename TF> int FINUFFT_PLAN_T<TF>::init_grid_kerFT_FFT() {
     timer.restart(); // plan the FFTW (to act in-place on the workspace fwBatch)
     int nthr_fft  = opts.nthreads;
     const auto ns = gridsize_for_fft(*this);
-    std::vector<TC, xsimd::aligned_allocator<TC, 64>> fwBatch(nf() * batchSize);
-    fftPlan->plan(ns, batchSize, fwBatch.data(), fftSign, opts.fftw, nthr_fft);
+    fftPlan->plan(ns, batchSize, fftSign, opts.fftw, nthr_fft);
     if (opts.debug)
       printf("[%s] FFT plan (mode %d, nthr=%d):\t%.3g s\n", __func__, opts.fftw, nthr_fft,
              timer.elapsedsec());