Documentation & changes to kereval0 to make it faster

Author: DiamonDinoia

src/spreadinterp.cpp

@@ -59,10 +59,11 @@ static FINUFFT_ALWAYS_INLINE auto ker_eval(FLT *FINUFFT_RESTRICT ker,
                                            const finufft_spread_opts &opts,
                                            const V... elems) noexcept;
 static FINUFFT_ALWAYS_INLINE FLT fold_rescale(FLT x, BIGINT N) noexcept;
-static FINUFFT_ALWAYS_INLINE void set_kernel_args(
-    FLT *args, FLT x, const finufft_spread_opts &opts) noexcept;
+template<uint8_t ns>
+static FINUFFT_ALWAYS_INLINE void set_kernel_args(FLT *args, FLT x) noexcept;
+template<uint8_t N>
 static FINUFFT_ALWAYS_INLINE void evaluate_kernel_vector(
-    FLT *ker, FLT *args, const finufft_spread_opts &opts, int N) noexcept;
+    FLT *ker, FLT *args, const finufft_spread_opts &opts) noexcept;
 template<uint8_t w, class simd_type = xsimd::make_sized_batch_t<
                         FLT, find_optimal_simd_width<FLT, w>()>> // aka ns
 static FINUFFT_ALWAYS_INLINE void eval_kernel_vec_Horner(
@@ -703,16 +704,15 @@ FLT evaluate_kernel(FLT x, const finufft_spread_opts &opts)
     return exp((FLT)opts.ES_beta * sqrt((FLT)1.0 - (FLT)opts.ES_c * x * x));
 }
 
-void set_kernel_args(FLT *args, FLT x, const finufft_spread_opts &opts) noexcept
+template<uint8_t ns>
+void set_kernel_args(FLT *args, FLT x) noexcept
 // Fills vector args[] with kernel arguments x, x+1, ..., x+ns-1.
 // needed for the vectorized kernel eval of Ludvig af K.
 {
-  int ns = opts.nspread;
   for (int i = 0; i < ns; i++) args[i] = x + (FLT)i;
 }
-
-void evaluate_kernel_vector(FLT *ker, FLT *args, const finufft_spread_opts &opts,
-                            const int N) noexcept
+template<uint8_t N>
+void evaluate_kernel_vector(FLT *ker, FLT *args, const finufft_spread_opts &opts) noexcept
 /* Evaluate ES kernel for a vector of N arguments; by Ludvig af K.
    If opts.kerpad true, args and ker must be allocated for Npad, and args is
    written to (to pad to length Npad), only first N outputs are correct.
@@ -742,8 +742,7 @@ void evaluate_kernel_vector(FLT *ker, FLT *args, const finufft_spread_opts &opts
     if (opts.kerpad) {
       // padded part should be zero, in spread_subproblem_nd_kernels, there are
       // out of bound writes to trg arrays
-      for (int i = N; i < Npad; ++i)
-        ker[i] = 0.0;
+      for (int i = N; i < Npad; ++i) ker[i] = 0.0;
     }
   } else {
     for (int i = 0; i < N; i++) // dummy for timing only
@@ -1798,8 +1797,8 @@ auto ker_eval(FLT *FINUFFT_RESTRICT ker, const finufft_spread_opts &opts,
     }
     if constexpr (kerevalmeth == 0) {
       alignas(simd_type::arch_type::alignment()) std::array<T, MAX_NSPREAD> kernel_args{};
-      set_kernel_args(kernel_args.data(), inputs[i], opts);
-      evaluate_kernel_vector(ker + (i * MAX_NSPREAD), kernel_args.data(), opts, ns);
+      set_kernel_args<ns>(kernel_args.data(), inputs[i]);
+      evaluate_kernel_vector<ns>(ker + (i * MAX_NSPREAD), kernel_args.data(), opts);
     }
   }
   return ker;

test/README

@@ -22,6 +22,14 @@ finufft{1,2,3}dmany_test{f}: accuracy/speed tests for vectorized transforms,
                    in a given dimension. Types 1, 2, and 3 are tested.
                    (exit code 0 is a pass).
                    Call with no arguments for argument documentation.
+finufft3dkernel_test{f} : test of the kernel evaluation methods.
+                   Only 3D NUFFT are tested.
+                   Types 1, 2, and 3 are tested.
+                   It requires Nmodes1 Nmodes2 Nmodes3 Nsrc
+                   sizes of the 3d grid and the number of sources.
+                   Optional arguments are the tolerance, debug flags,
+                   sort flag and upsampfac.
+                   (exit code 0 is a pass).
 dumbinputs{f} :    test of edge cases, invalid inputs, and plan interface.
                    No arguments needed (exit code 0 is a pass).
 testutils{f} :     test of utils module.

test/finufft3dkernel_test.cpp

@@ -4,15 +4,25 @@
 using namespace std;
 using namespace finufft::utils;
 
-const char *help[] = {"Tester for FINUFFT in 3d, all 3 types, either precision.",
-                      "",
-                      "Usage: finufft3d_test Nmodes1 Nmodes2 Nmodes3 Nsrc [tol [debug "
-                      "[spread_sort [upsampfac]]]]",
-                      "\teg:\tfinufft3d_test 100 200 50 1e6 1e-12 0 2 0.0 1e-11",
-                      "\tnotes:\tif errfail present, exit code 1 if any error > errfail",
-                      NULL};
-// Barnett 2/2/17 onwards.
-
+const char *help[] = {
+    "Tester for FINUFFT in 3d, all 3 types, either precision.",
+    "",
+    "Usage: finufft3d_test Nmodes1 Nmodes2 Nmodes3 Nsrc",
+    "\t[tol] error tolerance (default 1e-6)",
+    "\t[debug] (default 0) 0: silent, 1: text, 2: as 1 but also spreader",
+    "\t[spread_sort] (default 2) 0: don't sort NU pts, 1: do, 2: auto",
+    "\t[upsampfac] (default 2.0)",
+    "\teg: finufft3d_test 100 200 50 1e6 1e-12 0 2 0.0",
+    "\tnotes: exit code 1 if any error > tol",
+    nullptr};
+/**
+ * @brief Test the 3D NUFFT of type 1, 2, and 3.
+ * It evaluates the error of the kernel evaluation methods.
+ * It uses err(a,b)=||a-b||_2 / ||a||_2 as the error metric.
+ * It return FINUFFT error code if it is not 0.
+ * It returns 1 if any error exceeds tol.
+ * It returns 0 if test passes.
+ */
 int main(int argc, char *argv[]) {
   BIGINT M, N1, N2, N3;      // M = # srcs, N1,N2,N3 = # modes
   double w, tol      = 1e-6; // default
@@ -142,9 +152,6 @@ int main(int argc, char *argv[]) {
   std::vector<FLT> s(N); // targ freqs (1-cmpt)
   std::vector<FLT> t(N); // targ freqs (2-cmpt)
   std::vector<FLT> u(N); // targ freqs (3-cmpt)
-  FLT S1 = (FLT)N1 / 2;  // choose freq range sim to type 1
-  FLT S2 = (FLT)N2 / 2;
-  FLT S3 = (FLT)N3 / 2;
 
   timer.restart();
   printf("kerevalmeth 0:\n");
@@ -171,5 +178,7 @@ int main(int argc, char *argv[]) {
   err    = relerrtwonorm(N, F0.data(), F1.data());
   errmax = max(err, errmax);
   printf("\ttype 3 rel l2-err in F is %.3g\n", err);
+  // return 1 if any error exceeds tol
+  // or return finufft error code if it is not 0
   return (errmax > tol);
 }