merge master

Author: lu1and10

CHANGELOG

@@ -3,9 +3,11 @@ If not stated, FINUFFT is assumed (cuFINUFFT <=1.3 is listed separately).
 
 * cufinufft now supports modeord(type 1,2 only): 0 CMCL-style increasing mode
   order, 1 FFT-style mode order.
+* New doc page: migration guide from NFFT3 (2d1 case only).
 * New foldrescale, removes [-3pi,3pi) restriction on NU points, and slight
   speedup at large tols. Deprecates both opts.chkbnds and error code
-  FINUFFT_ERR_SPREAD_PTS_OUT_RANGE.  PR #440 (Marco Barbone + Martin Reinecke)
+  FINUFFT_ERR_SPREAD_PTS_OUT_RANGE. Also inlined kernel eval code, increases
+  compile of spreadinterp.cpp to 10s. PR #440 (Marco Barbone + Martin Reinecke)
 * CPU plan stage allows any # threads, warns if > omp_get_max_threads(); or
   if single-threaded fixes nthr=1 and warns opts.nthreads>1 attempt.
   Sort now respects spread_opts.sort_threads not nthreads. Supercedes PR 431.

docs/index.rst

@@ -35,6 +35,7 @@ Documentation contents
    python_gpu
    julia
    changelog
+   nfft_migr
    cufinufft_migration
    devnotes
    related

docs/nfft_migr.rst

@@ -0,0 +1,78 @@
+.. _nfft_migr:
+
+Migration guide from NFFT3
+==========================
+
+Here we outline how the C/C++ user can replace NUFFT calls to the popular library
+`Chemnitz NFFT3 library <https://www-user.tu-chemnitz.de/~potts/nfft/>`_ with
+FINUFFT CPU calls to achieve the same effect, possibly with more performance
+or less RAM usage.
+See [KKP] in the :ref:`references<refs>` for more about NFFT3, and [FIN] for
+some performance and RAM comparisons performed using the codes available in 2018.
+We use the `NFFT source on GitHub <https://github.com/NFFT/nfft>`_, version 3.5.4alpha.
+So far we only discuss:
+
+ * the adjoint NFFT (a.k.a. type 1)
+
+Also of interest (but not yet demonstrated below) is:
+
+ * the forward NFFT transform (a.k.a. type 2)
+ * the nonuniform to nonuniform NNFFT (a.k.a. type 3)
+   
+ .. note:: The NFFT3 library can do more things---real-valued data, sphere, rotation group, hyperbolic cross, inverse transforms---none of which FINUFFT can yet do directly (although our three transforms can be used as components in such tasks). We do not address those here.
+
+Migrating a 2D adjoint transform (type 1) in C from NFFT3 to FINUFFT
+--------------------------------------------------------------------
+
+We need to start with the simplest example of using NFFT3 on "user data" generated
+using plain, transparent, C commands (rather than relying on NFFT3-supplied
+data-generation, direct transform, and printing utilities as in the
+NFFT example :file:`examples/nfft/simple_test.c`, or even its simplest
+version
+at https://www-user.tu-chemnitz.de/~potts/nfft/download/nfft_simple_test.tar.gz ).
+We choose 2D since it is the simplest rectangular
+case that illustrates how to get the transposed
+coordinates (or mode array) ordering correct.
+After installing NFFT3 one should be able to compile and run the following:
+
+.. literalinclude:: ../tutorial/nfft2d1_test.c
+   :language: c
+
+This is a basic example, running single-threaded, at the highest precision
+(using ``nfft_init_guru`` would allow more control.)
+It demonstrates: i) the NFFT3 user must write their data into arrays allocated
+by ``nfft_plan``,
+ii) the single nonuniform point coordinate array
+is interleaved ($x_1, y_1, x_2, y_2, \dots, x_M, y_M$),
+iii) the output mode ordering is C (row-major) rather than Fortran (column-major;
+this affects how to convert frequency indices into the output array index),
+and iv) there is an extra factor of $2\pi$ in the exponent relative
+to the FINUFFT definition, because NFFT3 assumes a 1-periodic input domain.
+The code is found in our :file:`tutorial/nfft2d1_test.c`. Running the executable gives:
+
+::
+
+ 2D type 1 (NFFT3) done in 0.589 s: f_hat[-17,33]=86.0632804289+-350.023846367i, rel err 9.93e-14
+
+To show how to migrate this, we write a self-contained code that generates exactly
+the same "user data" (same random seed), then uses FINUFFT to do the transform
+to achieve exactly the same ``f_hat`` output array (in row-major C ordering).
+This entails scaling and swapping the nonequispaced coordinates just before sending
+to FINUFFT. Here is the corresponding C code (compare to the above):
+
+.. literalinclude:: ../tutorial/migrate2d1_test.c
+   :language: c
+
+The fact that NFFT3 uses row-major mode arrays whereas FINUFFT uses column-major has
+been handled here by swapping the input $x$ and $y$ coordinates and array sizes in the
+FINUFFT call. (Equivalently, this could have been achieved by transposing the ``f_hat``
+output array. We recommend the former route since it saves memory.) Running the
+executable gives:
+
+::
+
+ 2D type 1 (FINUFFT) in 0.0787 s: f_hat[-17,33]=86.0632804289+-350.023846367i, rel err 9.58e-14
+
+Comparing to the above, we see the same answer to all shown digits, a similar error for this tested output entry, plus a 7.5$\times$ speed-up. (Both use a single thread, tested on the same AMD 5700U laptop.) The user may of course now set a coarser (larger) value for ``tol`` and see a further speed-up.
+
+We believe that the above gives the essentials of how to convert your code from using NFFT3 to FINUFFT. Please read our documentation, especially the guru interface if multiple related transforms are required, then post a GitHub Issue if you are still stuck.

src/spreadinterp.cpp

@@ -468,6 +468,8 @@ int interpSorted(BIGINT* sort_indices,BIGINT N1, BIGINT N2, BIGINT N3,
 	  }
 
 	  else{
+            // due to ns being padded up to next multiple of 4 in the eval_kernel_vec_Horner
+            // and writing zeros out to this padded size, these must occur in the order x,y,z...
 	    eval_kernel_vec_Horner(ker1,x1,ns,opts);
 	    if (ndims > 1) eval_kernel_vec_Horner(ker2,x2,ns,opts);  
 	    if (ndims > 2) eval_kernel_vec_Horner(ker3,x3,ns,opts);
@@ -979,6 +981,8 @@ void spread_subproblem_2d(BIGINT off1,BIGINT off2,BIGINT size1,BIGINT size2,
       set_kernel_args(kernel_args+ns, x2, opts);
       evaluate_kernel_vector(kernel_values, kernel_args, opts, 2*ns);
     } else {
+      // due to ns being padded up to next multiple of 4 in the eval_kernel_vec_Horner
+      // and writing zeros out to this padded size, these must occur in the order x,y...
       eval_kernel_vec_Horner(ker1,x1,ns,opts);
       eval_kernel_vec_Horner(ker2,x2,ns,opts);
     }
@@ -1038,6 +1042,8 @@ void spread_subproblem_3d(BIGINT off1,BIGINT off2,BIGINT off3,BIGINT size1,
       set_kernel_args(kernel_args+2*ns, x3, opts);
       evaluate_kernel_vector(kernel_values, kernel_args, opts, 3*ns);
     } else {
+      // due to ns being padded up to next multiple of 4 in the eval_kernel_vec_Horner
+      // and writing zeros out to this padded size, these must occur in the order x,y,z...
       eval_kernel_vec_Horner(ker1,x1,ns,opts);
       eval_kernel_vec_Horner(ker2,x2,ns,opts);
       eval_kernel_vec_Horner(ker3,x3,ns,opts);

tutorial/migrate2d1_test.c

@@ -0,0 +1,60 @@
+/* Minimal example of 2D "adjoint" (aka type 1) NUFFT using FINUFFT
+   to match NFFT3's conventions, and the same user-given data as
+   for nfft2d1_test.c.   Single-threaded, timer.   Barnett 5/18/24
+   To compile (assuming FINUFFT include and lib in path):
+   gcc migrate2d1_test.c -o migrate2d1_test -lfinufft -lfftw3 -lm
+ */
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <complex.h>
+#include <time.h>
+#include <finufft.h>
+
+int main(void) {
+  int N[2] = {300,200};         // N0, N1 output shape in nfft3 sense
+  int M = 500000;               // num. nonuniform input points
+  double tol = 1e-13;           // user must choose (unlike nfft3's simple call)
+  
+  // user allocates all external arrays (and no internal ones)
+  double* x = (double *)malloc(sizeof(double)*M);  // x (0th) coords only here
+  double* y = (double *)malloc(sizeof(double)*M);  // y (1st) coords need separate ptr
+  double complex* f = (double complex*)malloc(sizeof(double complex)*M);
+  double complex* f_hat = (double complex*)malloc(sizeof(double complex)*N[0]*N[1]);  // output
+  
+  // start with exactly the same "user data" as in nfft2d1_test.c...
+  srand(0);                     // fix seed
+  for (int j=0; j<M; ++j) {     // nonequispaced pts, and strengths f...
+    x[j] = (double)rand()/RAND_MAX - 0.5;    // x unif rand in [-1/2,1/2)
+    y[j] = (double)rand()/RAND_MAX - 0.5;    // y "
+    f[j] = 2*((double)rand()/RAND_MAX)-1 + I*(2*((double)rand()/RAND_MAX)-1);
+  }
+      
+  clock_t before = clock();
+
+  // do transform, includes precompute, writing to f_hat...
+  for (int j=0; j<M; ++j) {          // change user coords so finufft same as nfft3
+    x[j] *= 2*M_PI; y[j] *= 2*M_PI;  // scales from 1-periodic to 2pi-periodic
+  }
+  finufft_opts opts;                 // opts struct
+  finufft_default_opts(&opts);       // set default opts (must start with this)
+  opts.nthreads = 1;                 // enforce single-thread
+  int ier = finufft2d1(M,y,x,f,+1,tol,N[1],N[0],f_hat,&opts);  // both x,y and N0,N1 swapped!
+
+  double secs = (clock()-before)/(double)CLOCKS_PER_SEC;
+
+  // now test that f_hat is as it would have been if original data were sent to nfft3...
+  int kx=-17, ky=33;            // check one output f_hat(kx,ky) vs direct computation
+  int kxout = kx + N[0]/2;
+  int kyout = ky + N[1]/2;
+  int i = kyout + kxout*N[1];   // the output index (nfft3 convention, not FINUFFT's)
+  double complex f_hat_test = 0.0 + 0.0*I;
+  for (int j=0; j<M; ++j)       // since x,y were mult by 2pi, no such factor here...
+    f_hat_test += f[j] * cexp(I*((double)kx*x[j]+(double)ky*y[j]));
+  double err = cabs(f_hat[i] - f_hat_test) / cabs(f_hat_test);
+  printf("2D type 1 (FINUFFT) in %.3g s: f_hat[%d,%d]=%.12g+%.12gi, rel err %.3g\n",
+         secs,kx,ky,creal(f_hat[i]),cimag(f_hat[i]),err);
+  
+  free(x); free(y); free(f); free(f_hat);   // user deallocates own I/O arrays
+  return ier;
+}

tutorial/nfft2d1_test.c

@@ -0,0 +1,48 @@
+/* Minimal example of 2D "adjoint" (aka type 1) NUFFT using NFFT3 library,
+   with user-given data, single-threaded, plus timer.  Barnett 5/17/24
+   To compile (assuming nfft3 installed):
+   gcc nfft2d1_test.c -o nfft2d1_test -lnfft3 -lfftw3 -lm
+ */
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <complex.h>
+#include <time.h>
+#include "nfft3.h"
+
+int main(void) {
+  int N[2] = {300,200};         // N1, N2 output mode numbers
+  int M = 500000;               // num. nonuniform input points
+  int dim=2;
+  nfft_plan p;
+  nfft_init(&p, dim, N, M);     // allocates user I/O arrays too
+  // make some "user data" (we must use arrays that nfft allocated)...
+  srand(0);                     // fix seed
+  for (int j=0; j<M; ++j) {     // nonequispaced pts, and strengths f...
+    p.x[2*j] = (double)rand()/RAND_MAX - 0.5;    // x unif rand in [-1/2,1/2)
+    p.x[2*j+1] = (double)rand()/RAND_MAX - 0.5;  // y "
+    p.f[j] = 2*((double)rand()/RAND_MAX)-1 + I*(2*((double)rand()/RAND_MAX)-1);
+  }
+
+  clock_t before = clock();
+  
+  if(p.flags & PRE_ONE_PSI)     // precompute psi, the entries of the matrix B
+    nfft_precompute_one_psi(&p);
+  nfft_adjoint(&p);             // do transform, write to p.f_hat
+
+  double secs = (clock()-before)/(double)CLOCKS_PER_SEC;
+  
+  int kx=-17, ky=33;            // check one output f_hat(kx,ky) vs direct computation
+  int kxout = kx + N[0]/2;
+  int kyout = ky + N[1]/2;
+  int i = kyout + kxout*N[1];   // output index: array ordered x slow, y fast
+  double complex f_hat_test = 0.0 + 0.0*I;
+  for (int j=0; j<M; ++j)       // 2pi fac; p.x array is x interleaved with y...
+    f_hat_test += p.f[j] * cexp(2*M_PI*I*((double)kx*p.x[2*j]+(double)ky*p.x[2*j+1]));
+  double err = cabs(p.f_hat[i] - f_hat_test) / cabs(f_hat_test);
+  printf("2D type 1 (NFFT3) done in %.3g s: f_hat[%d,%d]=%.12g+%.12gi, rel err %.3g\n",
+         secs,kx,ky,creal(p.f_hat[i]),cimag(p.f_hat[i]),err);
+  
+  nfft_finalize(&p);            // free the plan
+  return 0;
+}