Spreadinterponly (CPU and GPU)

docs/matlabhelp.doc

@@ -20,6 +20,7 @@
      ms    number of Fourier modes computed, may be even or odd;
            in either case, mode range is integers lying in [-ms/2, (ms-1)/2]
      opts.modeord: 0 (CMCL increasing mode ordering, default), 1 (FFT ordering)
+     opts.spreadinterponly: 0 (perform NUFFT, default), 1 (only spread/interp)
      opts   optional struct with optional fields controlling the following:
      opts.debug:   0 (silent, default), 1 (timing breakdown), 2 (debug info).
      opts.spread_debug: spreader: 0 (no text, default), 1 (some), or 2 (lots)
@@ -66,6 +67,7 @@
      isign if >=0, uses + sign in exponential, otherwise - sign.
      eps   relative precision requested (generally between 1e-15 and 1e-1)
      opts.modeord: 0 (CMCL increasing mode ordering, default), 1 (FFT ordering)
+     opts.spreadinterponly: 0 (perform NUFFT, default), 1 (only spread/interp)
      opts   optional struct with optional fields controlling the following:
      opts.debug:   0 (silent, default), 1 (timing breakdown), 2 (debug info).
      opts.spread_debug: spreader: 0 (no text, default), 1 (some), or 2 (lots)
@@ -162,6 +164,7 @@
      ms,mt  number of Fourier modes requested in x & y; each may be even or odd.
             In either case the mode range is integers lying in [-m/2, (m-1)/2]
      opts.modeord: 0 (CMCL increasing mode ordering, default), 1 (FFT ordering)
+     opts.spreadinterponly: 0 (perform NUFFT, default), 1 (only spread/interp)
      opts   optional struct with optional fields controlling the following:
      opts.debug:   0 (silent, default), 1 (timing breakdown), 2 (debug info).
      opts.spread_debug: spreader: 0 (no text, default), 1 (some), or 2 (lots)
@@ -211,6 +214,7 @@
      isign if >=0, uses + sign in exponential, otherwise - sign.
      eps   relative precision requested (generally between 1e-15 and 1e-1)
      opts.modeord: 0 (CMCL increasing mode ordering, default), 1 (FFT ordering)
+     opts.spreadinterponly: 0 (perform NUFFT, default), 1 (only spread/interp)
      opts   optional struct with optional fields controlling the following:
      opts.debug:   0 (silent, default), 1 (timing breakdown), 2 (debug info).
      opts.spread_debug: spreader: 0 (no text, default), 1 (some), or 2 (lots)
@@ -310,6 +314,7 @@
            even or odd.
            In either case the mode range is integers lying in [-m/2, (m-1)/2]
      opts.modeord: 0 (CMCL increasing mode ordering, default), 1 (FFT ordering)
+     opts.spreadinterponly: 0 (perform NUFFT, default), 1 (only spread/interp)
      opts   optional struct with optional fields controlling the following:
      opts.debug:   0 (silent, default), 1 (timing breakdown), 2 (debug info).
      opts.spread_debug: spreader: 0 (no text, default), 1 (some), or 2 (lots)
@@ -361,6 +366,7 @@
      isign if >=0, uses + sign in exponential, otherwise - sign.
      eps   relative precision requested (generally between 1e-15 and 1e-1)
      opts.modeord: 0 (CMCL increasing mode ordering, default), 1 (FFT ordering)
+     opts.spreadinterponly: 0 (perform NUFFT, default), 1 (only spread/interp)
      opts   optional struct with optional fields controlling the following:
      opts.debug:   0 (silent, default), 1 (timing breakdown), 2 (debug info).
      opts.spread_debug: spreader: 0 (no text, default), 1 (some), or 2 (lots)
@@ -493,6 +499,7 @@
      opts.floatprec: library precision to use, 'double' (default) or 'single'.
      for type 1 and 2 only, the following opts fields are also relevant:
      opts.modeord: 0 (CMCL increasing mode ordering, default), 1 (FFT ordering)
+     opts.spreadinterponly: 0 (perform NUFFT, default), 1 (only spread/interp)
  Outputs:
      plan            finufft_plan object (opaque pointer)
 

docs/opts.rst

@@ -92,6 +92,21 @@ Data handling options
 
   .. note:: The index *sets* are the same in the two ``modeord`` choices; their ordering differs only by a cyclic shift. The FFT ordering cyclically shifts the CMCL indices $\mbox{floor}(N/2)$ to the left (often called an "fftshift").
 
+**spreadinterponly**: [only has effect for type 1 or 2.] If ``0`` do the NUFFT as intended.
+If ``1``, omit the FFT and deconvolution (diagonal division by kernel Fourier transform) steps, thus
+returning *garbage answers as a NUFFT*, but allowing experts to perform an isolated
+spreading (if type 1) or interpolation (if type 2) by hijacking the usual FINUFFT API.
+The kernel (width and shape parameter)
+is determined by ``tol`` and ``opts.upsampfac``, just as it would be in an actual NUFFT.
+However, the kernel is not accessible in any other way, leaving the user to figure out how to
+make use of this interface to extract the actual kernel function.
+This provides a convenient (if slightly hacky) interface to our ``spreadinterp`` module
+(looping over multiple vectors if ``ntransf>1``).
+The known use-case here is estimating so-called density compensation, conventionally used in
+MRI (see `MRI-NUFFT <https://mind-inria.github.io/mri-nufft/nufft.html>`_),
+although it might also be useful in spectral Ewald.
+
+
 
 Diagnostic options
 ~~~~~~~~~~~~~~~~~~~~~~~

examples/spreadinterponly1d.cpp

@@ -0,0 +1,89 @@
+// this is all you must include for the finufft lib...
+#include <finufft.h>
+
+// also used in this example...
+#include <cassert>
+#include <chrono>
+#include <complex>
+#include <cstdio>
+#include <stdlib.h>
+#include <vector>
+using namespace std;
+using namespace std::chrono;
+
+int main(int argc, char *argv[])
+/* Example of double-prec spread/interp only tasks, with basic math tests.
+   Complex I/O arrays, but kernel is real.  Barnett 1/8/25
+   The math tests are:
+   1) for spread, check sum of spread kernel masses is as expected from one
+   times the strengths (ie testing the zero-frequency component in NUFFT).
+   2) for interp, check each interp kernel mass is the same as from one.
+   Without knowing the kernel, this is about all that can be done.
+   (Better math tests would be, ironically, to wrap the spreader/interpolator
+   into a NUFFT and test that :) But we already have that in FINUFFT.)
+
+   Compile and run (static library case):
+
+   g++ spreadinterponly1d.cpp -I../include ../lib-static/libfinufft.a -o
+   spreadinterponly1d -lfftw3 -lfftw3_omp && ./spreadinterponly1d
+
+   See: spreadtestnd for usage of internal (non FINUFFT-API) spread/interp.
+*/
+{
+  int M = 1e7; // number of nonuniform points
+  int N = 1e7; // size of regular grid
+  finufft_opts opts;
+  finufft_default_opts(&opts);
+  opts.spreadinterponly = 1;    // task: the following two control kernel used...
+  double tol            = 1e-9; // tolerance for (real) kernel shape design only
+  opts.upsampfac        = 2.0;  // pretend upsampling factor (really no upsampling)
+  // opts.spread_kerevalmeth = 0;  // would be needed for any nonstd upsampfac
+
+  complex<double> I = complex<double>(0.0, 1.0); // the imaginary unit
+  vector<double> x(M);                           // input
+  vector<complex<double>> c(M);                  // input
+  vector<complex<double>> F(N);                  // output (spread to this array)
+
+  // first spread M=1 single unit-strength at the origin, to get its total mass...
+  x[0]       = 0.0;
+  c[0]       = 1.0;
+  int unused = 1;
+  int ier    = finufft1d1(1, &x[0], &c[0], unused, tol, N, &F[0], &opts);
+  if (ier > 1) return ier;
+  complex<double> kersum = 0.0;
+  for (auto Fk : F) kersum += Fk; // kernel mass
+
+  // Now generate random nonuniform points (x) and complex strengths (c)...
+  for (int j = 0; j < M; ++j) {
+    x[j] = M_PI * (2 * ((double)rand() / RAND_MAX) - 1); // uniform random in [-pi,pi)
+    c[j] =
+        2 * ((double)rand() / RAND_MAX) - 1 + I * (2 * ((double)rand() / RAND_MAX) - 1);
+  }
+
+  opts.debug = 1;
+  auto t0    = steady_clock::now(); // now spread with all M pts... (dir=1)
+  ier        = finufft1d1(M, &x[0], &c[0], unused, tol, N, &F[0], &opts);
+  double t   = (steady_clock::now() - t0) / 1.0s;
+  if (ier > 1) return ier;
+  complex<double> csum = 0.0; // tot input strength
+  for (auto cj : c) csum += cj;
+  complex<double> mass = 0.0; // tot output mass
+  for (auto Fk : F) mass += Fk;
+  double relerr = abs(mass - kersum * csum) / abs(mass);
+  printf("1D spread-only, double-prec, %.3g s (%.3g NU pt/sec), ier=%d, mass err %.3g\n",
+         t, M / t, ier, relerr);
+
+  for (auto &Fk : F) Fk = complex<double>{1.0, 0.0}; // unit grid input
+  opts.debug = 0;
+  t0         = steady_clock::now(); // now interp to all M pts...  (dir=2)
+  ier        = finufft1d2(M, &x[0], &c[0], unused, tol, N, &F[0], &opts);
+  t          = (steady_clock::now() - t0) / 1.0s;
+  if (ier > 1) return ier;
+  csum = 0.0; // tot output
+  for (auto cj : c) csum += cj;
+  double maxerr = 0.0;
+  for (auto cj : c) maxerr = max(maxerr, abs(cj - kersum));
+  printf("1D interp-only, double-prec, %.3g s (%.3g NU pt/sec), ier=%d, max err %.3g\n",
+         t, M / t, ier, maxerr / abs(kersum));
+  return 0;
+}

include/cufinufft/impl.h

@@ -157,7 +157,7 @@ int cufinufft_makeplan_impl(int type, int dim, int *nmodes, int iflag, int ntran
   }
   if (d_plan->opts.gpu_spreadinterponly) {
     // XNOR implementation below with boolean logic.
-    if ((d_plan->opts.upsampfac != 1) == (type != 3)) {
+    if ((d_plan->opts.upsampfac != 1.0) == (type != 3)) {
       ier = FINUFFT_ERR_SPREADONLY_UPSAMP_INVALID;
       goto finalize;
     }
@@ -197,7 +197,6 @@ int cufinufft_makeplan_impl(int type, int dim, int *nmodes, int iflag, int ntran
     printf("[cufinufft] shared memory required for the spreader: %ld\n", mem_required);
   }
 
-
   // dynamically request the maximum amount of shared memory available
   // for the spreader
 
@@ -765,8 +764,8 @@ int cufinufft_setpts_impl(int M, T *d_kx, T *d_ky, T *d_kz, int N, T *d_s, T *d_
           thrust::cuda::par.on(stream), phase_iterator, phase_iterator + N,
           d_plan->deconv, d_plan->deconv,
           [c1, c2, c3, d1, d2, d3, realsign] __host__ __device__(
-              const thrust::tuple<T, T, T> tuple, cuda_complex<T> deconv)
-              -> cuda_complex<T> {
+              const thrust::tuple<T, T, T> tuple,
+              cuda_complex<T> deconv) -> cuda_complex<T> {
             // d2 and d3 are 0 if dim < 2 and dim < 3
             const auto phase = c1 * (thrust::get<0>(tuple) + d1) +
                                c2 * (thrust::get<1>(tuple) + d2) +

include/finufft.fh

@@ -7,7 +7,7 @@ c     erase chkbnds 1/7/25.
       type finufft_opts
 
 c     data handling opts...
-      integer modeord
+      integer modeord, spreadinterponly
 
 c     diagnostic opts...
       integer debug, spread_debug, showwarn

include/finufft/spreadinterp.h

@@ -52,11 +52,12 @@ FINUFFT_EXPORT int FINUFFT_CDECL spreadinterpSorted(
 template<typename T>
 FINUFFT_EXPORT T FINUFFT_CDECL evaluate_kernel(T x, const finufft_spread_opts &opts);
 template<typename T>
-FINUFFT_EXPORT T FINUFFT_CDECL evaluate_kernel_horner(T x, const finufft_spread_opts &opts);
+FINUFFT_EXPORT T FINUFFT_CDECL evaluate_kernel_horner(T x,
+                                                      const finufft_spread_opts &opts);
 template<typename T>
-FINUFFT_EXPORT int FINUFFT_CDECL setup_spreader(finufft_spread_opts &opts, T eps,
-                                                double upsampfac, int kerevalmeth,
-                                                int debug, int showwarn, int dim);
+FINUFFT_EXPORT int FINUFFT_CDECL setup_spreader(
+    finufft_spread_opts &opts, T eps, double upsampfac, int kerevalmeth, int debug,
+    int showwarn, int dim, int spreadinterponly);
 
 } // namespace spreadinterp
 } // namespace finufft

include/finufft_mod.f90

@@ -8,7 +8,7 @@ module finufft_mod
 type finufft_opts
 
    ! data handling opts...
-   integer(kind=C_INT) :: modeord
+   integer(kind=C_INT) :: modeord, spreadinterponly
 
    ! diagnostic opts...
    integer(kind=C_INT) :: debug, spread_debug, showwarn

include/finufft_opts.h

@@ -9,8 +9,10 @@ typedef struct finufft_opts { // defaults see finufft_core.cpp:finufft_default_o
   // sphinx tag (don't remove): @opts_start
   // FINUFFT options:
   // data handling opts...
-  int modeord; // (type 1,2 only): 0 CMCL-style increasing mode order
-               //                  1 FFT-style mode order
+  int modeord;          // (type 1,2 only): 0 CMCL-style increasing mode order
+                        //                  1 FFT-style mode order
+  int spreadinterponly; // (type 1,2 only): 0 (default) do full NUFFT
+                        //                  1 only do spread or interp
 
   // diagnostic opts...
   int debug;        // 0 silent, 1 some timing/debug, or 2 more