Migration guide from NFFT3¶
Here we outline how the C/C++ user can replace NUFFT calls to the popular library Chemnitz NFFT3 library with FINUFFT CPU calls to achieve the same effect, possibly with more performance or less RAM usage. See [KKP] in the references 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, 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 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:
/* 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 "nfft3.h"
#include <complex.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.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;
}
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 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):
/* 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 <complex.h>
#include <finufft.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.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;
}
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.
