36 static int comp1(
const void *x,
const void *y)
38 return ((*(
const R*) x) < (*(
const R*) y) ? -1 : 1);
41 static int comp2(
const void *x,
const void *y)
43 int nx0, nx1, ny0, ny1;
44 nx0 = global_n * (int)LRINT(*((
const R*) x + 0));
45 nx1 = global_n * (int)LRINT(*((
const R*) x + 1));
46 ny0 = global_n * (
int)LRINT(*((const R*) y + 0));
47 ny1 = global_n * (
int)LRINT(*((const R*) y + 1));
60 static
int comp3(const
void *x, const
void *y)
62 int nx0, nx1, nx2, ny0, ny1, ny2;
63 nx0 = global_n * (int)LRINT(*((
const R*) x + 0));
64 nx1 = global_n * (int)LRINT(*((
const R*) x + 1));
65 nx2 = global_n * (
int)LRINT(*((const R*) x + 2));
66 ny0 = global_n * (
int)LRINT(*((const R*) y + 0));
67 ny1 = global_n * (
int)LRINT(*((const R*) y + 1));
68 ny2 = global_n * (
int)LRINT(*((const R*) y + 2));
86 static
void measure_time_nfft(
int d,
int N,
unsigned test_ndft)
88 int r, M, NN[d], nn[d];
89 R t, t_fft, t_ndft, t_nfft;
95 printf("\\verb+%ld+&\t", LRINT(LOG((R)(N)) / LOG((R)(2)) * (R)(d) + K(0.5)));
97 for (r = 0, M = 1; r < d; r++)
104 NFFT(init_guru)(&p, d, NN, M, nn, 2,
106 PRE_FULL_PSI | MALLOC_F_HAT | MALLOC_X | MALLOC_F |
107 FFTW_INIT | FFT_OUT_OF_PLACE,
108 FFTW_MEASURE | FFTW_DESTROY_INPUT);
110 p_fft = FFTW(plan_dft)(d, NN, p.f_hat, p.f, FFTW_FORWARD, FFTW_MEASURE);
113 NFFT(vrand_shifted_unit_double)(p.x, p.d * p.M_total);
122 qsort(p.x, (
size_t)(p.M_total), (
size_t)(d) *
sizeof(R), comp1);
127 qsort(p.x, (
size_t)(p.M_total), (
size_t)(d) *
sizeof(R), comp2);
132 qsort(p.x, (
size_t)(p.M_total), (
size_t)(d) *
sizeof(R), comp3);
137 NFFT(precompute_one_psi)(&p);
140 NFFT(vrand_unit_complex)(p.f_hat, p.N_total);
146 while (t_fft < K(1.0))
150 FFTW(execute)(p_fft);
152 t = NFFT(elapsed_seconds)(t1, t0);
158 printf(
"\\verb+%.1" __FES__
"+ &\t", t_fft);
165 while (t_ndft < K(1.0))
169 NFFT(trafo_direct)(&p);
171 t = NFFT(elapsed_seconds)(t1, t0);
176 printf(
"\\verb+%.1" __FES__
"+ &\t", t_ndft);
180 printf(
"\\verb+*+\t&\t");
185 while (t_nfft < K(1.0))
210 t = NFFT(elapsed_seconds)(t1, t0);
216 printf(
"\\verb+%.1" __FES__
"+ & \\verb+(%3.1" __FIS__
")+\\\\\n", t_nfft, t_nfft / t_fft);
218 FFTW(destroy_plan)(p_fft);
222 static void measure_time_nfft_XXX2(
int d,
int N,
unsigned test_ndft)
224 int r, M, NN[d], nn[d];
225 R t, t_fft, t_ndft, t_nfft;
231 printf("%ld\t", LRINT(LOG((R)(N)) / LOG((R)(2)) * (R)(d) + K(0.5)));
234 for (r = 0, M = 1; r < d; r++)
241 NFFT(init_guru)(&p, d, NN, M, nn, 2,
244 MALLOC_F_HAT | MALLOC_X | MALLOC_F |
245 FFTW_INIT | FFT_OUT_OF_PLACE,
246 FFTW_MEASURE | FFTW_DESTROY_INPUT);
248 p_fft = FFTW(plan_dft)(d, NN, p.f_hat, p.f, FFTW_FORWARD, FFTW_MEASURE);
250 C *swapndft = (C*) NFFT(malloc)((size_t)(p.M_total) *
sizeof(C));
253 NFFT(vrand_shifted_unit_double)(p.x, p.d * p.M_total);
255 qsort(p.x, (
size_t)(p.M_total), (
size_t)(d) *
sizeof(R), comp1);
258 NFFT(precompute_one_psi)(&p);
261 NFFT(vrand_unit_complex)(p.f_hat, p.N_total);
266 while (t_fft < K(0.1))
270 FFTW(execute)(p_fft);
272 t = NFFT(elapsed_seconds)(t1, t0);
277 printf(
"%.1" __FES__
"\t", t_fft);
282 CSWAP(p.f, swapndft);
285 while (t_ndft < K(0.1))
289 NFFT(trafo_direct)(&p);
291 t = NFFT(elapsed_seconds)(t1, t0);
295 printf(
"%.1" __FES__
"\t", t_ndft);
296 CSWAP(p.f, swapndft);
304 while (t_nfft < K(0.1))
310 t = NFFT(elapsed_seconds)(t1, t0);
314 printf(
"%.1" __FES__
"\t", t_nfft);
316 printf(
"(%.1" __FES__
")\t", NFFT(error_l_2_complex)(swapndft, p.f, p.M_total));
321 while (t_nfft < K(0.1))
327 t = NFFT(elapsed_seconds)(t1, t0);
331 printf(
"%.1" __FES__
"\t", t_nfft);
333 printf(
"(%.1" __FES__
")\t", NFFT(error_l_2_complex)(swapndft, p.f, p.M_total));
337 NFFT(free)(swapndft);
338 FFTW(destroy_plan)(p_fft);
342 static void measure_time_nfft_XXX3(
int d,
int N,
unsigned test_ndft)
344 int r, M, NN[d], nn[d];
345 R t, t_fft, t_ndft, t_nfft;
351 printf("%ld\t", LRINT(LOG((R)(N)) / LOG((R)(2)) * (R)(d) + K(0.5)));
354 for (r = 0, M = 1; r < d; r++)
361 NFFT(init_guru)(&p, d, NN, M, nn, 2,
364 MALLOC_F_HAT | MALLOC_X | MALLOC_F |
365 FFTW_INIT | FFT_OUT_OF_PLACE,
366 FFTW_MEASURE | FFTW_DESTROY_INPUT);
368 p_fft = FFTW(plan_dft)(d, NN, p.f, p.f_hat, FFTW_BACKWARD, FFTW_MEASURE);
370 C *swapndft = (C*) NFFT(malloc)((size_t)(p.N_total) *
sizeof(C));
373 NFFT(vrand_shifted_unit_double)(p.x, p.d * p.M_total);
375 qsort(p.x, (
size_t)(p.M_total), (
size_t)(d) *
sizeof(R), comp1);
378 NFFT(precompute_one_psi)(&p);
381 NFFT(vrand_unit_complex)(p.f, p.N_total);
386 while (t_fft < K(0.1))
390 FFTW(execute)(p_fft);
392 t = NFFT(elapsed_seconds)(t1, t0);
397 printf(
"%.1" __FES__
"\t", t_fft);
402 CSWAP(p.f_hat, swapndft);
405 while (t_ndft < K(0.1))
409 NFFT(adjoint_direct)(&p);
411 t = NFFT(elapsed_seconds)(t1, t0);
415 printf(
"%.1" __FES__
"\t", t_ndft);
416 CSWAP(p.f_hat, swapndft);
424 while (t_nfft < K(0.1))
430 t = NFFT(elapsed_seconds)(t1, t0);
434 printf(
"%.1" __FES__
"\t", t_nfft);
436 printf(
"(%.1" __FES__
")\t", NFFT(error_l_2_complex)(swapndft, p.f_hat, p.N_total));
441 while (t_nfft < K(0.1))
445 NFFT(adjoint_1d)(&p);
447 t = NFFT(elapsed_seconds)(t1, t0);
451 printf(
"%.1" __FES__
"\t", t_nfft);
453 printf(
"(%.1" __FES__
")\t", NFFT(error_l_2_complex)(swapndft, p.f_hat, p.N_total));
457 NFFT(free)(swapndft);
458 FFTW(destroy_plan)(p_fft);
462 static void measure_time_nfft_XXX4(
int d,
int N,
unsigned test_ndft)
464 int r, M, NN[d], nn[d];
465 R t, t_fft, t_ndft, t_nfft;
471 printf("%ld\t", LRINT(LOG((R)(N)) / LOG((R)(2)) * (R)(d) + K(0.5)));
474 for (r = 0, M = 1; r < d; r++)
481 NFFT(init_guru)(&p, d, NN, M, nn, 4,
484 MALLOC_F_HAT | MALLOC_X | MALLOC_F |
485 FFTW_INIT | FFT_OUT_OF_PLACE,
486 FFTW_MEASURE | FFTW_DESTROY_INPUT);
488 p_fft = FFTW(plan_dft)(d, NN, p.f_hat, p.f, FFTW_FORWARD, FFTW_MEASURE);
490 C *swapndft = (C*) NFFT(malloc)((size_t)(p.M_total) *
sizeof(C));
493 NFFT(vrand_shifted_unit_double)(p.x, p.d * p.M_total);
501 NFFT(precompute_one_psi)(&p);
504 NFFT(vrand_unit_complex)(p.f_hat, p.N_total);
509 while (t_fft < K(0.1))
513 FFTW(execute)(p_fft);
515 t = NFFT(elapsed_seconds)(t1, t0);
520 printf(
"%.1" __FES__
"\t", t_fft);
523 NFFT(vrand_unit_complex)(p.f_hat, p.N_total);
528 CSWAP(p.f, swapndft);
531 while (t_ndft < K(0.1))
535 NFFT(trafo_direct)(&p);
537 t = NFFT(elapsed_seconds)(t1, t0);
541 printf(
"%.1" __FES__
"\t", t_ndft);
545 CSWAP(p.f, swapndft);
553 while (t_nfft < K(0.1))
559 t = NFFT(elapsed_seconds)(t1, t0);
563 printf(
"%.1" __FES__
"\t", t_nfft);
565 printf(
"(%.1" __FES__
")\t", NFFT(error_l_2_complex)(swapndft, p.f, p.M_total));
572 while (t_nfft < K(0.1))
578 t = NFFT(elapsed_seconds)(t1, t0);
582 printf(
"%.1" __FES__
"\t", t_nfft);
584 printf(
"(%.1" __FES__
")\t", NFFT(error_l_2_complex)(swapndft, p.f, p.M_total));
590 NFFT(free)(swapndft);
591 FFTW(destroy_plan)(p_fft);
595 static void measure_time_nfft_XXX5(
int d,
int N,
unsigned test_ndft)
597 int r, M, NN[d], nn[d];
598 R t, t_fft, t_ndft, t_nfft;
604 printf("%ld\t", LRINT(LOG((R)(N)) / LOG((R)(2)) * (R)(d) + K(0.5)));
607 for (r = 0, M = 1; r < d; r++)
614 NFFT(init_guru)(&p, d, NN, M, nn, 4,
617 MALLOC_F_HAT | MALLOC_X | MALLOC_F |
618 FFTW_INIT | FFT_OUT_OF_PLACE,
619 FFTW_MEASURE | FFTW_DESTROY_INPUT);
621 p_fft = FFTW(plan_dft)(d, NN, p.f, p.f_hat, FFTW_FORWARD, FFTW_MEASURE);
623 C *swapndft = (C*) NFFT(malloc)((size_t)(p.N_total) *
sizeof(C));
626 NFFT(vrand_shifted_unit_double)(p.x, p.d * p.M_total);
630 NFFT(precompute_one_psi)(&p);
633 NFFT(vrand_unit_complex)(p.f, p.M_total);
638 while (t_fft < K(0.1))
642 FFTW(execute)(p_fft);
644 t = NFFT(elapsed_seconds)(t1, t0);
649 printf(
"%.1" __FES__
"\t", t_fft);
652 NFFT(vrand_unit_complex)(p.f, p.M_total);
657 CSWAP(p.f_hat, swapndft);
660 while (t_ndft < K(0.1))
664 NFFT(adjoint_direct)(&p);
666 t = NFFT(elapsed_seconds)(t1, t0);
670 printf(
"%.1" __FES__
"\t", t_ndft);
674 CSWAP(p.f_hat, swapndft);
682 while (t_nfft < K(0.1))
688 t = NFFT(elapsed_seconds)(t1, t0);
692 printf(
"%.1" __FES__
"\t", t_nfft);
694 printf(
"(%.1" __FES__
")\t", NFFT(error_l_2_complex)(swapndft, p.f_hat, p.N_total));
701 while (t_nfft < K(0.1))
705 NFFT(adjoint_2d)(&p);
707 t = NFFT(elapsed_seconds)(t1, t0);
711 printf(
"%.1" __FES__
"\t", t_nfft);
713 printf(
"(%.1" __FES__
")\t", NFFT(error_l_2_complex)(swapndft, p.f_hat, p.N_total));
719 NFFT(free)(swapndft);
720 FFTW(destroy_plan)(p_fft);
724 static void measure_time_nfft_XXX6(
int d,
int N,
unsigned test_ndft)
726 int r, M, NN[d], nn[d];
727 R t, t_fft, t_ndft, t_nfft;
733 printf("%ld\t", LRINT(LOG((R)(N)) / LOG((R)(2)) * (R)(d) + K(0.5)));
736 for (r = 0, M = 1; r < d; r++)
743 NFFT(init_guru)(&p, d, NN, M, nn, 2,
746 MALLOC_F_HAT | MALLOC_X | MALLOC_F |
747 FFTW_INIT | FFT_OUT_OF_PLACE,
748 FFTW_MEASURE | FFTW_DESTROY_INPUT);
750 p_fft = FFTW(plan_dft)(d, NN, p.f_hat, p.f, FFTW_FORWARD, FFTW_MEASURE);
752 C *swapndft = (C*) NFFT(malloc)((size_t)(p.M_total) *
sizeof(C));
755 NFFT(vrand_shifted_unit_double)(p.x, p.d * p.M_total);
760 NFFT(precompute_one_psi)(&p);
763 NFFT(vrand_unit_complex)(p.f_hat, p.N_total);
768 while (t_fft < K(0.1))
772 FFTW(execute)(p_fft);
774 t = NFFT(elapsed_seconds)(t1, t0);
779 printf(
"%.1" __FES__
"\t", t_fft);
782 NFFT(vrand_unit_complex)(p.f_hat, p.N_total);
787 CSWAP(p.f, swapndft);
790 while (t_ndft < K(0.1))
794 NFFT(trafo_direct)(&p);
796 t = NFFT(elapsed_seconds)(t1, t0);
800 printf(
"%.1" __FES__
"\t", t_ndft);
804 CSWAP(p.f, swapndft);
812 while (t_nfft < K(0.1))
818 t = NFFT(elapsed_seconds)(t1, t0);
822 printf(
"%.1" __FES__
"\t", t_nfft);
824 printf(
"(%.1" __FES__
")\t", NFFT(error_l_2_complex)(swapndft, p.f, p.M_total));
831 while (t_nfft < K(0.1))
837 t = NFFT(elapsed_seconds)(t1, t0);
841 printf(
"%.1" __FES__
"\t", t_nfft);
843 printf(
"(%.1" __FES__
")\t", NFFT(error_l_2_complex)(swapndft, p.f, p.M_total));
849 NFFT(free)(swapndft);
850 FFTW(destroy_plan)(p_fft);
854 static void measure_time_nfft_XXX7(
int d,
int N,
unsigned test_ndft)
856 int r, M, NN[d], nn[d];
857 R t, t_fft, t_ndft, t_nfft;
863 printf("%ld\t", LRINT(LOG((R)(N)) / LOG((R)(2)) * (R)(d) + K(0.5)));
866 for (r = 0, M = 1; r < d; r++)
873 NFFT(init_guru)(&p, d, NN, M, nn, 2,
876 MALLOC_F_HAT | MALLOC_X | MALLOC_F |
877 FFTW_INIT | FFT_OUT_OF_PLACE,
878 FFTW_MEASURE | FFTW_DESTROY_INPUT);
880 p_fft = FFTW(plan_dft)(d, NN, p.f, p.f_hat, FFTW_FORWARD, FFTW_MEASURE);
882 C *swapndft = (C*) NFFT(malloc)((size_t)(p.N_total) *
sizeof(C));
885 NFFT(vrand_shifted_unit_double)(p.x, p.d * p.M_total);
889 NFFT(precompute_one_psi)(&p);
892 NFFT(vrand_unit_complex)(p.f, p.M_total);
897 while (t_fft < K(0.1))
901 FFTW(execute)(p_fft);
903 t = NFFT(elapsed_seconds)(t1, t0);
908 printf(
"%.1" __FES__
"\t", t_fft);
911 NFFT(vrand_unit_complex)(p.f, p.M_total);
916 CSWAP(p.f_hat, swapndft);
919 while (t_ndft < K(0.1))
923 NFFT(adjoint_direct)(&p);
925 t = NFFT(elapsed_seconds)(t1, t0);
929 printf(
"%.1" __FES__
"\t", t_ndft);
933 CSWAP(p.f_hat, swapndft);
941 while (t_nfft < K(0.1))
947 t = NFFT(elapsed_seconds)(t1, t0);
951 printf(
"%.1" __FES__
"\t", t_nfft);
953 printf(
"(%.1" __FES__
")\t", NFFT(error_l_2_complex)(swapndft, p.f_hat, p.N_total));
960 while (t_nfft < K(0.1))
964 NFFT(adjoint_3d)(&p);
966 t = NFFT(elapsed_seconds)(t1, t0);
970 printf(
"%.1" __FES__
"\t", t_nfft);
972 printf(
"(%.1" __FES__
")\t", NFFT(error_l_2_complex)(swapndft, p.f_hat, p.N_total));
978 NFFT(free)(swapndft);
979 FFTW(destroy_plan)(p_fft);
1020 UNUSED(measure_time_nfft_XXX2);
1021 UNUSED(measure_time_nfft_XXX3);
1022 UNUSED(measure_time_nfft_XXX4);
1023 UNUSED(measure_time_nfft_XXX5);
1024 UNUSED(measure_time_nfft_XXX6);
1025 UNUSED(measure_time_nfft_XXX7);
1027 printf(
"\\hline $l_N$ & FFT & NDFT & NFFT & NFFT/FFT\\\\\n");
1028 printf(
"\\hline \\hline \\multicolumn{5}{|c|}{$d=1$} \\\\ \\hline\n");
1029 for (l = 3; l <= 22; l++)
1033 int N = (int)(1U << (logIN / d));
1034 measure_time_nfft(d, N, logIN <= 15 ? 1 : 0);
1039 printf(
"\\hline $l_N$ & FFT & NDFT & NFFT & NFFT/FFT\\\\\n");
1040 printf(
"\\hline \\hline \\multicolumn{5}{|c|}{$d=2$} \\\\ \\hline\n");
1041 for (l = 3; l <= 11; l++)
1045 int N = (int)(1U << (logIN / d));
1046 measure_time_nfft(d, N, logIN <= 15 ? 1 : 0);
1051 printf(
"\\hline \\hline \\multicolumn{5}{|c|}{$d=3$} \\\\ \\hline\n");
1052 for (l = 3; l <= 7; l++)
1056 int N = (int)(1U << (logIN / d));
1057 measure_time_nfft(d, N, logIN <= 15 ? 1 : 0);
#define UNUSED(x)
Dummy use of unused parameters to silence compiler warnings.
#define CSWAP(x, y)
Swap two vectors.