-
Notifications
You must be signed in to change notification settings - Fork 73
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
docs: migration guide from nfft3 (two new tutorial C codes)
- Loading branch information
Showing
5 changed files
with
190 additions
and
1 deletion.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -35,6 +35,7 @@ Documentation contents | |
python_gpu | ||
julia | ||
changelog | ||
nfft_migr | ||
cufinufft_migration | ||
devnotes | ||
related | ||
|
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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. |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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; | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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; | ||
} |