Adding options to switch kernel at runtime

CMakeLists.txt

@@ -83,6 +83,8 @@ if(FINUFFT_USE_CPU)
     add_subdirectory(src)
 endif()
 
+add_subdirectory(src/common)
+
 if(FINUFFT_USE_CUDA)
     include(cuda_setup)
     add_subdirectory(src/cuda)

cmake/utils.cmake

@@ -103,7 +103,7 @@ function(enable_asan target)
 endfunction()
 
 function(finufft_link_test target)
-    target_link_libraries(${target} PRIVATE finufft::finufft)
+    target_link_libraries(${target} PRIVATE finufft::finufft finufft_common)
     if(FINUFFT_USE_DUCC0)
         target_compile_definitions(${target} PRIVATE FINUFFT_USE_DUCC0)
     endif()

examples/cuda/CMakeLists.txt

@@ -5,6 +5,6 @@ foreach(srcfile ${example_src})
     get_filename_component(executable ${executable} NAME)
     add_executable(${executable} ${srcfile})
     target_include_directories(${executable} PUBLIC ${CUFINUFFT_INCLUDE_DIRS})
-    target_link_libraries(${executable} cufinufft CUDA::cufft CUDA::cudart)
+    target_link_libraries(${executable} PRIVATE cufinufft CUDA::cufft CUDA::cudart $<BUILD_INTERFACE:finufft_common>)
     target_compile_features(${executable} PRIVATE cxx_std_17)
 endforeach()

include/cufinufft/utils.h

@@ -89,7 +89,7 @@ class WithCudaDevice {
 };
 
 // math helpers whose source is in src/utils.cpp
-long next235beven(long n, long b);
+FINUFFT_EXPORT long next235beven(long n, long b);
 
 /**
  * does a complex atomic add on a shared memory address

include/finufft.fh

@@ -20,6 +20,7 @@ c     DEPRECATED: spread_kerpad kept for ABI compatibility, ignored by library
       real*8 upsampfac
       integer spread_thread, maxbatchsize, spread_nthr_atomic
       integer spread_max_sp_size
+      integer spread_kernel
       integer fftw_lock_fun, fftw_unlock_fun, fftw_lock_data
 
       end type

include/finufft/spreadinterp.h

@@ -21,9 +21,9 @@ namespace finufft {
 namespace spreadinterp {
 
 template<typename T>
-FINUFFT_EXPORT_TEST int setup_spreader(finufft_spread_opts &opts, T eps, double upsampfac,
-                                       int kerevalmeth, int debug, int showwarn,
-                                       int spreadinterponly, int dim);
+FINUFFT_EXPORT_TEST int setup_spreader(
+    finufft_spread_opts &opts, T eps, double upsampfac, int kerevalmeth, int debug,
+    int showwarn, int spreadinterponly, int dim, int kernel_type = 0);
 
 int spreadcheck(UBIGINT N1, UBIGINT N2, UBIGINT N3, const finufft_spread_opts &opts);
 template<typename T>

include/finufft_common/kernel.h

@@ -6,6 +6,8 @@
 #include <vector>
 
 #include <finufft_common/constants.h>
+#include <finufft_common/utils.h>
+#include <finufft_spread_opts.h>
 
 namespace finufft::kernel {
 
@@ -56,15 +58,34 @@ template<class T, class F> std::vector<T> fit_monomials(F &&f, int n, T a, T b)
   return c;
 }
 
-template<typename T> T evaluate_kernel(T x, T beta, T c) {
-  /* ES ("exp sqrt") kernel evaluation at single real argument:
-      phi(x) = exp(beta.(sqrt(1 - (2x/n_s)^2) - 1)),    for |x| < nspread/2
-     related to an asymptotic approximation to the Kaiser--Bessel, itself an
-     approximation to prolate spheroidal wavefunction (PSWF) of order 0.
-     This is the "reference implementation", used by eg finufft/onedim_* 2/17/17.
-     Rescaled so max is 1, Barnett 7/21/24
+template<typename T> T evaluate_kernel(T x, T beta, T c, int kernel_type = 0) {
+  /* Kernel evaluation at single real argument.
+     kernel_type == 0 : ES ("exp sqrt") kernel (default)
+       phi_ES(x) = exp(beta*(sqrt(1 - c*x^2) - 1))
+     kernel_type == 1 : Kaiser--Bessel (KB) kernel
+       phi_KB(x) = I_0(beta*sqrt(1 - c*x^2)) / I_0(beta)
+     Note: `std::cyl_bessel_i` from <cmath> is used for I_0.
+     Rescaled so max is 1.
   */
+  if (kernel_type == 1) {
+    // Kaiser--Bessel (normalized by I0(beta)). Use std::cyl_bessel_i from <cmath>.
+    const T inner        = std::sqrt(T(1) - c * x * x);
+    const T arg          = beta * inner;
+    const double i0_arg  = ::finufft::common::cyl_bessel_i(0, static_cast<double>(arg));
+    const double i0_beta = ::finufft::common::cyl_bessel_i(0, static_cast<double>(beta));
+    return static_cast<T>(i0_arg / i0_beta);
+  }
+
+  // default to ES
   return std::exp(beta * (std::sqrt(T(1) - c * x * x) - T(1)));
 }
 
+FINUFFT_EXPORT int compute_kernel_ns(double upsampfac, double tol, int kernel_type,
+                                     const finufft_spread_opts &opts);
+
+FINUFFT_EXPORT void initialize_kernel_params(finufft_spread_opts &opts, double upsampfac,
+                                             double tol, int kernel_type);
+
+FINUFFT_EXPORT double sigma_max_tol(double upsampfac, int kernel_type, int max_ns);
+
 } // namespace finufft::kernel

include/finufft_common/utils.h

@@ -13,6 +13,9 @@
 FINUFFT_EXPORT_TEST void gaussquad(int n, double *xgl, double *wgl);
 std::tuple<double, double> leg_eval(int n, double x);
 
+// Series implementation of the modified Bessel function of the first kind I_nu(x)
+double cyl_bessel_i(double nu, double x) noexcept;
+
 // helper to generate the integer sequence in range [Start, End]
 template<int Offset, typename Seq> struct offset_seq;
 
@@ -100,7 +103,7 @@
 }
 
 template<typename Tuple, std::size_t... I>
 auto extract_seqs_impl(const Tuple &t, std::index_sequence<I...>) {
   using T = std::remove_reference_t<Tuple>;
   return std::make_tuple(typename std::tuple_element_t<I, T>::seq_type{}...);
 }

include/finufft_mod.f90

@@ -21,6 +21,7 @@ module finufft_mod
    real(kind=C_DOUBLE) :: upsampfac
    integer(kind=C_INT) :: spread_thread, maxbatchsize
    integer(kind=C_INT) :: spread_nthr_atomic, spread_max_sp_size
+   integer(kind=C_INT) :: spread_kernel
    integer(kind=C_SIZE_T) :: fftw_lock_fun, fftw_unlock_fun, fftw_lock_data
    !  really, last should be type(C_PTR) :: etc, but fails to print nicely
 

include/finufft_opts.h

@@ -32,6 +32,8 @@ typedef struct finufft_opts { // defaults see finufft_core.cpp:finufft_default_o
   int spread_nthr_atomic; // if >=0, threads above which spreader OMP critical goes
                           // atomic
   int spread_max_sp_size; // if >0, overrides spreader (dir=1) max subproblem size
+  int spread_function;    // (dev only) 0:DEFAULT, (do not change), there is no guarantee
+                          // what non-zero values do and behaviour can change anytime
                           // sphinx tag (don't remove): @opts_end
 
   // User can provide their own FFTW planner lock functions for thread safety

include/finufft_spread_opts.h

@@ -26,6 +26,8 @@ typedef struct finufft_spread_opts {
   double ES_beta;
   double ES_halfwidth;
   double ES_c;
+  // Kernel selector: 0 = ES (default), 1 = Kaiser--Bessel (KB)
+  int kernel_type; // default 0
 } finufft_spread_opts;
 
 #endif // FINUFFT_SPREAD_OPTS_H

makefile

@@ -152,7 +152,7 @@ STATICLIB = lib-static/$(LIBNAME).a
 ABSDYNLIB = $(FINUFFT)$(DYNLIB)
 
 # spreader objs
-SOBJS = src/finufft_utils.o src/utils.o src/spreadinterp.o
+SOBJS = src/finufft_utils.o src/spreadinterp.o src/common/utils.o src/common/kernel.o
 
 # all lib dual-precision objs (note DUCC_OBJS empty if unused)
 OBJS = $(SOBJS) src/fft.o src/finufft_core.o src/c_interface.o fortran/finufftfort.o $(DUCC_OBJS)
@@ -283,10 +283,10 @@ test/%: test/%.cpp $(DYNLIB)
 test/%f: test/%.cpp $(DYNLIB)
 	$(CXX) $(CXXFLAGS) ${LDFLAGS} -DSINGLE $< $(ABSDYNLIB) $(LIBSFFT) -o $@
 # low-level tests that are cleaner if depend on only specific objects...
-test/testutils: test/testutils.cpp src/finufft_utils.o src/utils.o
-	$(CXX) $(CXXFLAGS) ${LDFLAGS} test/testutils.cpp src/finufft_utils.o src/utils.o $(LIBS) -o test/testutils
-test/testutilsf: test/testutils.cpp src/finufft_utils.o src/utils.o
-	$(CXX) $(CXXFLAGS) ${LDFLAGS} -DSINGLE test/testutils.cpp src/finufft_utils.o src/utils.o $(LIBS) -o test/testutilsf
+test/testutils: test/testutils.cpp src/finufft_utils.o src/common/utils.o
+	$(CXX) $(CXXFLAGS) ${LDFLAGS} test/testutils.cpp src/finufft_utils.o src/common/utils.o $(LIBS) -o test/testutils
+test/testutilsf: test/testutils.cpp src/finufft_utils.o src/common/utils.o
+	$(CXX) $(CXXFLAGS) ${LDFLAGS} -DSINGLE test/testutils.cpp src/finufft_utils.o src/common/utils.o $(LIBS) -o test/testutilsf
 
 # make sure all double-prec test executables ready for testing
 TESTS := $(basename $(wildcard test/*.cpp))

matlab/finufft.mw

@@ -103,6 +103,9 @@ $     }
 $     else if (strcmp(fname[ifield],"spread_max_sp_size") == 0) {
 $       oc->spread_max_sp_size = (int)round(*mxGetPr(mxGetFieldByNumber(om,idx,ifield)));
 $     }
+$     else if (strcmp(fname[ifield],"spread_kernel") == 0) {
+$       oc->spread_kernel = (int)round(*mxGetPr(mxGetFieldByNumber(om,idx,ifield)));
+$     }
 $     else if (strcmp(fname[ifield],"spreadinterponly") == 0) {
 $       oc->spreadinterponly = (int)round(*mxGetPr(mxGetFieldByNumber(om,idx,ifield)));
 $     }

perftest/cuda/CMakeLists.txt

@@ -1,6 +1,6 @@
 add_executable(cuperftest cuperftest.cu)
 target_include_directories(cuperftest PUBLIC ${CUFINUFFT_INCLUDE_DIRS})
-target_link_libraries(cuperftest cufinufft CUDA::cufft CUDA::cudart)
+target_link_libraries(cuperftest PRIVATE cufinufft CUDA::cufft CUDA::cudart $<BUILD_INTERFACE:finufft_common>)
 target_compile_features(cuperftest PRIVATE cxx_std_17)
 set_target_properties(
     cuperftest

python/finufft/finufft/_finufft.py

@@ -86,6 +86,7 @@ class FinufftOpts(ctypes.Structure):
                       ('maxbatchsize', c_int),
                       ('spread_nthr_atomic', c_int),
                       ('spread_max_sp_size', c_int),
+                      ('spread_kernel', c_int),
                       ('fftw_lock_fun', c_void_p),
                       ('fftw_unlock_fun', c_void_p),
                       ('fftw_lock_data', c_void_p)]

src/CMakeLists.txt

@@ -23,7 +23,6 @@ set(FINUFFT_SOURCES
     finufft_core.cpp
     c_interface.cpp
     finufft_utils.cpp
-    utils.cpp
 )
 
 if(FINUFFT_BUILD_FORTRAN)
@@ -48,7 +47,7 @@ if(FINUFFT_USE_DUCC0)
     target_compile_definitions(finufft PRIVATE FINUFFT_USE_DUCC0)
 endif()
 
-target_link_libraries(finufft PRIVATE $<BUILD_INTERFACE:finufft_fftlibs xsimd>)
+target_link_libraries(finufft PRIVATE $<BUILD_INTERFACE:finufft_fftlibs xsimd finufft_common>)
 if(FINUFFT_USE_OPENMP)
     target_link_libraries(finufft PRIVATE OpenMP::OpenMP_CXX)
     if(NOT FINUFFT_STATIC_LINKING)

src/common/CMakeLists.txt

@@ -0,0 +1,29 @@
+cmake_minimum_required(VERSION 3.24)
+
+# Build a small static helper library for CPU-only kernel parameter helpers
+# Source is in ../kernel.cpp as requested
+set(FINUFFT_COMMON_SOURCES kernel.cpp utils.cpp)
+
+add_library(finufft_common STATIC ${FINUFFT_COMMON_SOURCES})
+
+# The public include directory is the top-level include/
+target_include_directories(
+    finufft_common
+    PUBLIC $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/include> $<INSTALL_INTERFACE:include>
+)
+set_target_properties(finufft_common PROPERTIES POSITION_INDEPENDENT_CODE ON)
+
+# Visibility and compile options consistent with finufft
+if(FINUFFT_SHARED_LINKING)
+    target_compile_definitions(finufft_common PRIVATE FINUFFT_DLL)
+endif()
+
+target_compile_features(finufft_common PRIVATE cxx_std_17)
+
+# If you need -lm etc.
+find_library(MATH_LIBRARY m)
+if(MATH_LIBRARY)
+    target_link_libraries(finufft_common PRIVATE ${MATH_LIBRARY})
+endif()
+
+set_target_properties(finufft_common PROPERTIES CXX_VISIBILITY_PRESET hidden VISIBILITY_INLINES_HIDDEN YES)

src/common/kernel.cpp

@@ -0,0 +1,79 @@
+#include <algorithm>
+#include <cmath>
+#include <cstdio>
+#include <finufft_common/common.h>
+#include <finufft_common/kernel.h>
+#include <finufft_spread_opts.h>
+
+namespace finufft::kernel {
+
+int compute_kernel_ns(double upsampfac, double tol, int kernel_type,
+                      const finufft_spread_opts &opts) {
+  // Note: opts is unused here but kept for API consistency if needed later
+  int ns;
+  if (upsampfac == 2.0) {
+    ns = (int)std::ceil(-std::log10(tol / 10.0));
+  } else {
+    ns = (int)std::ceil(
+        -std::log(tol) / (finufft::common::PI * std::sqrt(1.0 - 1.0 / upsampfac)));
+  }
+  ns = std::max(2, ns);
+  return ns;
+}
+
+void initialize_kernel_params(finufft_spread_opts &opts, double upsampfac, double tol,
+                              int kernel_type) {
+  int ns;
+  // Respect any pre-set opts.nspread (e.g., clipped by caller). If it's <=0 compute it.
+  if (opts.nspread > 0) {
+    ns = opts.nspread;
+  } else {
+    ns = compute_kernel_ns(upsampfac, tol, kernel_type, opts);
+  }
+
+  opts.kernel_type  = kernel_type;
+  opts.nspread      = ns; // ensure opts is populated with the (possibly clipped) ns
+  opts.ES_halfwidth = (double)ns / 2.0;
+  opts.ES_c         = 4.0 / (double)(ns * ns);
+
+  if (kernel_type == 1) {
+    // Kaiser-Bessel (KB)
+    // Formula from Beatty et al. 2005
+    double tmp   = (double)ns * (double)ns / (upsampfac * upsampfac);
+    double term2 = (upsampfac - 0.5) * (upsampfac - 0.5);
+    opts.ES_beta = common::PI * std::sqrt(tmp * term2 - 0.8);
+  } else {
+    // Exponential of Semicircle (ES)
+    double betaoverns = 2.30;
+    if (ns == 2)
+      betaoverns = 2.20;
+    else if (ns == 3)
+      betaoverns = 2.26;
+    else if (ns == 4)
+      betaoverns = 2.38;
+
+    if (upsampfac != 2.0) {
+      double gamma = 0.97;
+      betaoverns   = gamma * common::PI * (1.0 - 1.0 / (2.0 * upsampfac));
+    }
+    opts.ES_beta = betaoverns * (double)ns;
+  }
+
+  if (opts.debug) {
+    const char *kname = (kernel_type == 1) ? "KB" : "ES";
+    printf("setup_spreader: using spread kernel type %d (%s)\n", kernel_type, kname);
+    printf("setup_spreader eps=%.3g sigma=%.3g (%s): chose ns=%d beta=%.3g\n", tol,
+           upsampfac, kname, ns, opts.ES_beta);
+  }
+}
+
+double sigma_max_tol(double upsampfac, int kernel_type, int max_ns) {
+  if (upsampfac == 2.0) {
+    return 10.0 * std::pow(10.0, -(double)max_ns);
+  } else {
+    return std::exp(
+        -(double)max_ns * finufft::common::PI * std::sqrt(1.0 - 1.0 / upsampfac));
+  }
+}
+
+} // namespace finufft::kernel

src/utils.cpp → src/common/utils.cpp

@@ -1,5 +1,19 @@
+
 #include <cmath>
 #include <finufft_common/common.h>
+#include <limits>
+
+// Prefer the standard library's special-math `cyl_bessel_i` when available.
+#if defined(__has_include)
+#if __has_include(<version>)
+#include <version>
+#endif
+#endif
+// Feature-test macro for special math functions (if available in the standard
+// library implementation). Fall back to our series implementation otherwise.
+#if defined(__cpp_lib_math_special_functions)
+#define FINUFFT_HAVE_STD_CYL_BESSEL_I 1
+#endif
 
 #ifdef __CUDACC__
 #include <cufinufft/types.h>
@@ -24,7 +38,7 @@ void gaussquad(int n, double *xgl, double *wgl) {
   for (int i = 0; i < n / 2; i++) { // Loop through nodes
     convcount = 0;
     x         = std::cos((2 * i + 1) * PI / (2 * n)); // Initial guess: Chebyshev node
-    while (true) {                               // Newton iteration
+    while (true) {                                    // Newton iteration
       auto [p, dp] = leg_eval(n, x);
       dx           = -p / dp;
       x += dx; // Newton step
@@ -71,6 +85,34 @@ std::tuple<double, double> leg_eval(int n, double x) {
   return {p2, n * (x * p2 - p1) / (x * x - 1)};
 }
 
+double cyl_bessel_i(double nu, double x) noexcept {
+#if defined(FINUFFT_HAVE_STD_CYL_BESSEL_I)
+  return std::cyl_bessel_i(nu, x);
+#else
+  if (x == 0.0) {
+    if (nu == 0.0) return 1.0;
+    return 0.0;
+  }
+
+  const double halfx = x / 2.0;
+  double term        = std::pow(halfx, nu) / std::tgamma(nu + 1.0); // k = 0
+  double sum         = term;
+
+  static constexpr auto eps      = std::numeric_limits<double>::epsilon() * 10.0;
+  static constexpr auto max_iter = 100000;
+
+  for (int k = 1; k < max_iter; ++k) {
+    term *= (halfx * halfx) / (static_cast<double>(k) * (nu + static_cast<double>(k)));
+    sum += term;
+
+    if (std::abs(term) < eps * std::abs(sum)) {
+      break;
+    }
+  }
+  return sum;
+#endif
+}
+
 } // namespace common
 } // namespace finufft