NFFT  3.3.0
nfft_times.c
1 /*
2  * Copyright (c) 2002, 2015 Jens Keiner, Stefan Kunis, Daniel Potts
3  *
4  * This program is free software; you can redistribute it and/or modify it under
5  * the terms of the GNU General Public License as published by the Free Software
6  * Foundation; either version 2 of the License, or (at your option) any later
7  * version.
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12  * details.
13  *
14  * You should have received a copy of the GNU General Public License along with
15  * this program; if not, write to the Free Software Foundation, Inc., 51
16  * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  */
18 
19 /* $Id$ */
20 #include "config.h"
21 
22 #include <stdio.h>
23 #include <math.h>
24 #include <string.h>
25 #include <stdlib.h>
26 #ifdef HAVE_COMPLEX_H
27 #include <complex.h>
28 #endif
29 
30 #include "nfft3.h"
31 #include "infft.h"
32 
33 int global_n;
34 int global_d;
35 
36 static int comp1(const void *x, const void *y)
37 {
38  return ((*(const R*) x) < (*(const R*) y) ? -1 : 1);
39 }
40 
41 static int comp2(const void *x, const void *y)
42 {
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));
48 
49  if (nx0 < ny0)
50  return -1;
51  else if (nx0 == ny0)
52  if (nx1 < ny1)
53  return -1;
54  else
55  return 1;
56  else
57  return 1;
58 }
59 
60 static int comp3(const void *x, const void *y)
61 {
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));
69 
70  if (nx0 < ny0)
71  return -1;
72  else if (nx0 == ny0)
73  if (nx1 < ny1)
74  return -1;
75  else if (nx1 == ny1)
76  if (nx2 < ny2)
77  return -1;
78  else
79  return 1;
80  else
81  return 1;
82  else
83  return 1;
84 }
85 
86 static void measure_time_nfft(int d, int N, unsigned test_ndft)
87 {
88  int r, M, NN[d], nn[d];
89  R t, t_fft, t_ndft, t_nfft;
90  ticks t0, t1;
91 
92  NFFT(plan) p;
93  FFTW(plan) p_fft;
94 
95  printf("\\verb+%ld+&\t", LRINT(LOG((R)(N)) / LOG((R)(2)) * (R)(d) + K(0.5)));
96 
97  for (r = 0, M = 1; r < d; r++)
98  {
99  M = N * M;
100  NN[r] = N;
101  nn[r] = 2 * N;
102  }
103 
104  NFFT(init_guru)(&p, d, NN, M, nn, 2,
105  PRE_PHI_HUT |
106  PRE_FULL_PSI | MALLOC_F_HAT | MALLOC_X | MALLOC_F |
107  FFTW_INIT | FFT_OUT_OF_PLACE,
108  FFTW_MEASURE | FFTW_DESTROY_INPUT);
109 
110  p_fft = FFTW(plan_dft)(d, NN, p.f_hat, p.f, FFTW_FORWARD, FFTW_MEASURE);
111 
113  NFFT(vrand_shifted_unit_double)(p.x, p.d * p.M_total);
114 
115  global_n = nn[0];
116  global_d = d;
117 
118  switch (d)
119  {
120  case 1:
121  {
122  qsort(p.x, (size_t)(p.M_total), (size_t)(d) * sizeof(R), comp1);
123  break;
124  }
125  case 2:
126  {
127  qsort(p.x, (size_t)(p.M_total), (size_t)(d) * sizeof(R), comp2);
128  break;
129  }
130  case 3:
131  {
132  qsort(p.x, (size_t)(p.M_total), (size_t)(d) * sizeof(R), comp3);
133  break;
134  }
135  }
136 
137  NFFT(precompute_one_psi)(&p);
138 
139  /* init pseudo random Fourier coefficients */
140  NFFT(vrand_unit_complex)(p.f_hat, p.N_total);
141 
143  t_fft = K(0.0);
144  r = 0;
145 
146  while (t_fft < K(1.0))
147  {
148  r++;
149  t0 = getticks();
150  FFTW(execute)(p_fft);
151  t1 = getticks();
152  t = NFFT(elapsed_seconds)(t1, t0);
153  t_fft += t;
154  }
155  t_fft /= (R)(r);
156 
157  // printf("\\verb+%.1" __FES__ "+ & \\verb+(%.1f)+ &\t",t_fft,t_fft/(p.N_total*(log(N)/log(2)*d))*auxC);
158  printf("\\verb+%.1" __FES__ "+ &\t", t_fft);
159 
161  if (test_ndft)
162  {
163  t_ndft = K(0.0);
164  r = 0;
165  while (t_ndft < K(1.0))
166  {
167  r++;
168  t0 = getticks();
169  NFFT(trafo_direct)(&p);
170  t1 = getticks();
171  t = NFFT(elapsed_seconds)(t1, t0);
172  t_ndft += t;
173  }
174  t_ndft /= (R)(r);
175  //printf("\\verb+%.1" __FES__ "+ & \\verb+(%d)+&\t",t_ndft,(int)round(t_ndft/(p.N_total*p.N_total)*auxC));
176  printf("\\verb+%.1" __FES__ "+ &\t", t_ndft);
177  }
178  else
179  // printf("\\verb+*+\t&\t&\t");
180  printf("\\verb+*+\t&\t");
181 
183  t_nfft = K(0.0);
184  r = 0;
185  while (t_nfft < K(1.0))
186  {
187  r++;
188  t0 = getticks();
189  switch (d)
190  {
191  case 1:
192  {
193  NFFT(trafo_1d)(&p);
194  break;
195  }
196  case 2:
197  {
198  NFFT(trafo_2d)(&p);
199  break;
200  }
201  case 3:
202  {
203  NFFT(trafo_3d)(&p);
204  break;
205  }
206  default:
207  NFFT(trafo)(&p);
208  }
209  t1 = getticks();
210  t = NFFT(elapsed_seconds)(t1, t0);
211  t_nfft += t;
212  }
213  t_nfft /= (R)(r);
214 
215  // printf("\\verb+%.1" __FES__ "+ & \\verb+(%d)+ & \\verb+(%.1" __FES__ ")+\\\\\n",t_nfft,(int)round(t_nfft/(p.N_total*(log(N)/log(2)*d))*auxC),t_nfft/t_fft);
216  printf("\\verb+%.1" __FES__ "+ & \\verb+(%3.1" __FIS__ ")+\\\\\n", t_nfft, t_nfft / t_fft);
217 
218  FFTW(destroy_plan)(p_fft);
219  NFFT(finalize)(&p);
220 }
221 
222 static void measure_time_nfft_XXX2(int d, int N, unsigned test_ndft)
223 {
224  int r, M, NN[d], nn[d];
225  R t, t_fft, t_ndft, t_nfft;
226  ticks t0, t1;
227 
228  NFFT(plan) p;
229  FFTW(plan) p_fft;
230 
231  printf("%ld\t", LRINT(LOG((R)(N)) / LOG((R)(2)) * (R)(d) + K(0.5)));
232  fflush(stdout);
233 
234  for (r = 0, M = 1; r < d; r++)
235  {
236  M = N * M;
237  NN[r] = N;
238  nn[r] = 2 * N;
239  }
240 
241  NFFT(init_guru)(&p, d, NN, M, nn, 2,
242  PRE_PHI_HUT |
243  PRE_PSI |
244  MALLOC_F_HAT | MALLOC_X | MALLOC_F |
245  FFTW_INIT | FFT_OUT_OF_PLACE,
246  FFTW_MEASURE | FFTW_DESTROY_INPUT);
247 
248  p_fft = FFTW(plan_dft)(d, NN, p.f_hat, p.f, FFTW_FORWARD, FFTW_MEASURE);
249 
250  C *swapndft = (C*) NFFT(malloc)((size_t)(p.M_total) * sizeof(C));
251 
253  NFFT(vrand_shifted_unit_double)(p.x, p.d * p.M_total);
254 
255  qsort(p.x, (size_t)(p.M_total), (size_t)(d) * sizeof(R), comp1);
256  //nfft_vpr_double(p.x,p.M_total,"nodes x");
257 
258  NFFT(precompute_one_psi)(&p);
259 
261  NFFT(vrand_unit_complex)(p.f_hat, p.N_total);
262 
264  t_fft = K(0.0);
265  r = 0;
266  while (t_fft < K(0.1))
267  {
268  r++;
269  t0 = getticks();
270  FFTW(execute)(p_fft);
271  t1 = getticks();
272  t = NFFT(elapsed_seconds)(t1, t0);
273  t_fft += t;
274  }
275  t_fft /= (R)(r);
276 
277  printf("%.1" __FES__ "\t", t_fft);
278 
280  if (test_ndft)
281  {
282  CSWAP(p.f, swapndft);
283  t_ndft = K(0.0);
284  r = 0;
285  while (t_ndft < K(0.1))
286  {
287  r++;
288  t0 = getticks();
289  NFFT(trafo_direct)(&p);
290  t1 = getticks();
291  t = NFFT(elapsed_seconds)(t1, t0);
292  t_ndft += t;
293  }
294  t_ndft /= (R)(r);
295  printf("%.1" __FES__ "\t", t_ndft);
296  CSWAP(p.f, swapndft);
297  }
298  else
299  printf("\t");
300 
302  t_nfft = K(0.0);
303  r = 0;
304  while (t_nfft < K(0.1))
305  {
306  r++;
307  t0 = getticks();
308  NFFT(trafo)(&p);
309  t1 = getticks();
310  t = NFFT(elapsed_seconds)(t1, t0);
311  t_nfft += t;
312  }
313  t_nfft /= (R)(r);
314  printf("%.1" __FES__ "\t", t_nfft);
315  if (test_ndft)
316  printf("(%.1" __FES__ ")\t", NFFT(error_l_2_complex)(swapndft, p.f, p.M_total));
317 
319  t_nfft = K(0.0);
320  r = 0;
321  while (t_nfft < K(0.1))
322  {
323  r++;
324  t0 = getticks();
325  NFFT(trafo_1d)(&p);
326  t1 = getticks();
327  t = NFFT(elapsed_seconds)(t1, t0);
328  t_nfft += t;
329  }
330  t_nfft /= (R)(r);
331  printf("%.1" __FES__ "\t", t_nfft);
332  if (test_ndft)
333  printf("(%.1" __FES__ ")\t", NFFT(error_l_2_complex)(swapndft, p.f, p.M_total));
334 
335  printf("\n");
336 
337  NFFT(free)(swapndft);
338  FFTW(destroy_plan)(p_fft);
339  NFFT(finalize)(&p);
340 }
341 
342 static void measure_time_nfft_XXX3(int d, int N, unsigned test_ndft)
343 {
344  int r, M, NN[d], nn[d];
345  R t, t_fft, t_ndft, t_nfft;
346  ticks t0, t1;
347 
348  NFFT(plan) p;
349  FFTW(plan) p_fft;
350 
351  printf("%ld\t", LRINT(LOG((R)(N)) / LOG((R)(2)) * (R)(d) + K(0.5)));
352  fflush(stdout);
353 
354  for (r = 0, M = 1; r < d; r++)
355  {
356  M = N * M;
357  NN[r] = N;
358  nn[r] = 2 * N;
359  }
360 
361  NFFT(init_guru)(&p, d, NN, M, nn, 2,
362  PRE_PHI_HUT |
363  PRE_PSI |
364  MALLOC_F_HAT | MALLOC_X | MALLOC_F |
365  FFTW_INIT | FFT_OUT_OF_PLACE,
366  FFTW_MEASURE | FFTW_DESTROY_INPUT);
367 
368  p_fft = FFTW(plan_dft)(d, NN, p.f, p.f_hat, FFTW_BACKWARD, FFTW_MEASURE);
369 
370  C *swapndft = (C*) NFFT(malloc)((size_t)(p.N_total) * sizeof(C));
371 
373  NFFT(vrand_shifted_unit_double)(p.x, p.d * p.M_total);
374 
375  qsort(p.x, (size_t)(p.M_total), (size_t)(d) * sizeof(R), comp1);
376  //nfft_vpr_double(p.x,p.M_total,"nodes x");
377 
378  NFFT(precompute_one_psi)(&p);
379 
381  NFFT(vrand_unit_complex)(p.f, p.N_total);
382 
384  t_fft = K(0.0);
385  r = 0;
386  while (t_fft < K(0.1))
387  {
388  r++;
389  t0 = getticks();
390  FFTW(execute)(p_fft);
391  t1 = getticks();
392  t = NFFT(elapsed_seconds)(t1, t0);
393  t_fft += t;
394  }
395  t_fft /= (R)(r);
396 
397  printf("%.1" __FES__ "\t", t_fft);
398 
400  if (test_ndft)
401  {
402  CSWAP(p.f_hat, swapndft);
403  t_ndft = K(0.0);
404  r = 0;
405  while (t_ndft < K(0.1))
406  {
407  r++;
408  t0 = getticks();
409  NFFT(adjoint_direct)(&p);
410  t1 = getticks();
411  t = NFFT(elapsed_seconds)(t1, t0);
412  t_ndft += t;
413  }
414  t_ndft /= (R)(r);
415  printf("%.1" __FES__ "\t", t_ndft);
416  CSWAP(p.f_hat, swapndft);
417  }
418  else
419  printf("\t");
420 
422  t_nfft = K(0.0);
423  r = 0;
424  while (t_nfft < K(0.1))
425  {
426  r++;
427  t0 = getticks();
428  NFFT(adjoint)(&p);
429  t1 = getticks();
430  t = NFFT(elapsed_seconds)(t1, t0);
431  t_nfft += t;
432  }
433  t_nfft /= (R)(r);
434  printf("%.1" __FES__ "\t", t_nfft);
435  if (test_ndft)
436  printf("(%.1" __FES__ ")\t", NFFT(error_l_2_complex)(swapndft, p.f_hat, p.N_total));
437 
439  t_nfft = K(0.0);
440  r = 0;
441  while (t_nfft < K(0.1))
442  {
443  r++;
444  t0 = getticks();
445  NFFT(adjoint_1d)(&p);
446  t1 = getticks();
447  t = NFFT(elapsed_seconds)(t1, t0);
448  t_nfft += t;
449  }
450  t_nfft /= (R)(r);
451  printf("%.1" __FES__ "\t", t_nfft);
452  if (test_ndft)
453  printf("(%.1" __FES__ ")\t", NFFT(error_l_2_complex)(swapndft, p.f_hat, p.N_total));
454 
455  printf("\n");
456 
457  NFFT(free)(swapndft);
458  FFTW(destroy_plan)(p_fft);
459  NFFT(finalize)(&p);
460 }
461 
462 static void measure_time_nfft_XXX4(int d, int N, unsigned test_ndft)
463 {
464  int r, M, NN[d], nn[d];
465  R t, t_fft, t_ndft, t_nfft;
466  ticks t0, t1;
467 
468  NFFT(plan) p;
469  FFTW(plan) p_fft;
470 
471  printf("%ld\t", LRINT(LOG((R)(N)) / LOG((R)(2)) * (R)(d) + K(0.5)));
472  fflush(stdout);
473 
474  for (r = 0, M = 1; r < d; r++)
475  {
476  M = N * M;
477  NN[r] = N;
478  nn[r] = 2 * N;
479  }
480 
481  NFFT(init_guru)(&p, d, NN, M, nn, 4,
482  PRE_PHI_HUT |
483  PRE_PSI |
484  MALLOC_F_HAT | MALLOC_X | MALLOC_F |
485  FFTW_INIT | FFT_OUT_OF_PLACE,
486  FFTW_MEASURE | FFTW_DESTROY_INPUT);
487 
488  p_fft = FFTW(plan_dft)(d, NN, p.f_hat, p.f, FFTW_FORWARD, FFTW_MEASURE);
489 
490  C *swapndft = (C*) NFFT(malloc)((size_t)(p.M_total) * sizeof(C));
491 
493  NFFT(vrand_shifted_unit_double)(p.x, p.d * p.M_total);
494 
495  //for(j=0;j<2*M;j+=2)
496  // p.x[j]=0.5*p.x[j]+0.25*(p.x[j]>=0?1:-1);
497 
498  //qsort(p.x,p.M_total,d*sizeof(double),comp1);
499  //nfft_vpr_double(p.x,p.M_total,"nodes x");
500 
501  NFFT(precompute_one_psi)(&p);
502 
504  NFFT(vrand_unit_complex)(p.f_hat, p.N_total);
505 
507  t_fft = K(0.0);
508  r = 0;
509  while (t_fft < K(0.1))
510  {
511  r++;
512  t0 = getticks();
513  FFTW(execute)(p_fft);
514  t1 = getticks();
515  t = NFFT(elapsed_seconds)(t1, t0);
516  t_fft += t;
517  }
518  t_fft /= (R)(r);
519 
520  printf("%.1" __FES__ "\t", t_fft);
521 
523  NFFT(vrand_unit_complex)(p.f_hat, p.N_total);
524 
526  if (test_ndft)
527  {
528  CSWAP(p.f, swapndft);
529  t_ndft = K(0.0);
530  r = 0;
531  while (t_ndft < K(0.1))
532  {
533  r++;
534  t0 = getticks();
535  NFFT(trafo_direct)(&p);
536  t1 = getticks();
537  t = NFFT(elapsed_seconds)(t1, t0);
538  t_ndft += t;
539  }
540  t_ndft /= (R)(r);
541  printf("%.1" __FES__ "\t", t_ndft);
542 
543  //printf("f=%e+i%e\t",creal(p.f[0]),cimag(p.f[0]));
544 
545  CSWAP(p.f, swapndft);
546  }
547  else
548  printf("\t");
549 
551  t_nfft = K(0.0);
552  r = 0;
553  while (t_nfft < K(0.1))
554  {
555  r++;
556  t0 = getticks();
557  NFFT(trafo)(&p);
558  t1 = getticks();
559  t = NFFT(elapsed_seconds)(t1, t0);
560  t_nfft += t;
561  }
562  t_nfft /= (R)(r);
563  printf("%.1" __FES__ "\t", t_nfft);
564  if (test_ndft)
565  printf("(%.1" __FES__ ")\t", NFFT(error_l_2_complex)(swapndft, p.f, p.M_total));
566 
567  //printf("f=%e+i%e\t",creal(p.f[0]),cimag(p.f[0]));
568 
570  t_nfft = K(0.0);
571  r = 0;
572  while (t_nfft < K(0.1))
573  {
574  r++;
575  t0 = getticks();
576  NFFT(trafo_2d)(&p);
577  t1 = getticks();
578  t = NFFT(elapsed_seconds)(t1, t0);
579  t_nfft += t;
580  }
581  t_nfft /= (R)(r);
582  printf("%.1" __FES__ "\t", t_nfft);
583  if (test_ndft)
584  printf("(%.1" __FES__ ")\t", NFFT(error_l_2_complex)(swapndft, p.f, p.M_total));
585 
586  //printf("f=%e+i%e\t",creal(p.f[0]),cimag(p.f[0]));
587 
588  printf("\n");
589 
590  NFFT(free)(swapndft);
591  FFTW(destroy_plan)(p_fft);
592  NFFT(finalize)(&p);
593 }
594 
595 static void measure_time_nfft_XXX5(int d, int N, unsigned test_ndft)
596 {
597  int r, M, NN[d], nn[d];
598  R t, t_fft, t_ndft, t_nfft;
599  ticks t0, t1;
600 
601  NFFT(plan) p;
602  FFTW(plan) p_fft;
603 
604  printf("%ld\t", LRINT(LOG((R)(N)) / LOG((R)(2)) * (R)(d) + K(0.5)));
605  fflush(stdout);
606 
607  for (r = 0, M = 1; r < d; r++)
608  {
609  M = N * M;
610  NN[r] = N;
611  nn[r] = 2 * N;
612  }
613 
614  NFFT(init_guru)(&p, d, NN, M, nn, 4,
615  PRE_PHI_HUT |
616  PRE_PSI |
617  MALLOC_F_HAT | MALLOC_X | MALLOC_F |
618  FFTW_INIT | FFT_OUT_OF_PLACE,
619  FFTW_MEASURE | FFTW_DESTROY_INPUT);
620 
621  p_fft = FFTW(plan_dft)(d, NN, p.f, p.f_hat, FFTW_FORWARD, FFTW_MEASURE);
622 
623  C *swapndft = (C*) NFFT(malloc)((size_t)(p.N_total) * sizeof(C));
624 
626  NFFT(vrand_shifted_unit_double)(p.x, p.d * p.M_total);
627 
628  //sort_nodes(p.x,p.d,p.M_total,
629 
630  NFFT(precompute_one_psi)(&p);
631 
633  NFFT(vrand_unit_complex)(p.f, p.M_total);
634 
636  t_fft = K(0.0);
637  r = 0;
638  while (t_fft < K(0.1))
639  {
640  r++;
641  t0 = getticks();
642  FFTW(execute)(p_fft);
643  t1 = getticks();
644  t = NFFT(elapsed_seconds)(t1, t0);
645  t_fft += t;
646  }
647  t_fft /= (R)(r);
648 
649  printf("%.1" __FES__ "\t", t_fft);
650 
652  NFFT(vrand_unit_complex)(p.f, p.M_total);
653 
655  if (test_ndft)
656  {
657  CSWAP(p.f_hat, swapndft);
658  t_ndft = K(0.0);
659  r = 0;
660  while (t_ndft < K(0.1))
661  {
662  r++;
663  t0 = getticks();
664  NFFT(adjoint_direct)(&p);
665  t1 = getticks();
666  t = NFFT(elapsed_seconds)(t1, t0);
667  t_ndft += t;
668  }
669  t_ndft /= (R)(r);
670  printf("%.1" __FES__ "\t", t_ndft);
671 
672  //printf("\nf_hat=%e+i%e\t",creal(p.f_hat[0]),cimag(p.f_hat[0]));
673 
674  CSWAP(p.f_hat, swapndft);
675  }
676  else
677  printf("\t");
678 
680  t_nfft = K(0.0);
681  r = 0;
682  while (t_nfft < K(0.1))
683  {
684  r++;
685  t0 = getticks();
686  NFFT(adjoint)(&p);
687  t1 = getticks();
688  t = NFFT(elapsed_seconds)(t1, t0);
689  t_nfft += t;
690  }
691  t_nfft /= (R)(r);
692  printf("%.1" __FES__ "\t", t_nfft);
693  if (test_ndft)
694  printf("(%.1" __FES__ ")\t", NFFT(error_l_2_complex)(swapndft, p.f_hat, p.N_total));
695 
696  //printf("\nf_hat=%e+i%e\t",creal(p.f_hat[0]),cimag(p.f_hat[0]));
697 
699  t_nfft = K(0.0);
700  r = 0;
701  while (t_nfft < K(0.1))
702  {
703  r++;
704  t0 = getticks();
705  NFFT(adjoint_2d)(&p);
706  t1 = getticks();
707  t = NFFT(elapsed_seconds)(t1, t0);
708  t_nfft += t;
709  }
710  t_nfft /= (R)(r);
711  printf("%.1" __FES__ "\t", t_nfft);
712  if (test_ndft)
713  printf("(%.1" __FES__ ")\t", NFFT(error_l_2_complex)(swapndft, p.f_hat, p.N_total));
714 
715  //printf("\nf_hat=%e+i%e\t",creal(p.f_hat[0]),cimag(p.f_hat[0]));
716 
717  printf("\n");
718 
719  NFFT(free)(swapndft);
720  FFTW(destroy_plan)(p_fft);
721  NFFT(finalize)(&p);
722 }
723 
724 static void measure_time_nfft_XXX6(int d, int N, unsigned test_ndft)
725 {
726  int r, M, NN[d], nn[d];
727  R t, t_fft, t_ndft, t_nfft;
728  ticks t0, t1;
729 
730  NFFT(plan) p;
731  FFTW(plan) p_fft;
732 
733  printf("%ld\t", LRINT(LOG((R)(N)) / LOG((R)(2)) * (R)(d) + K(0.5)));
734  fflush(stdout);
735 
736  for (r = 0, M = 1; r < d; r++)
737  {
738  M = N * M;
739  NN[r] = N;
740  nn[r] = 2 * N;
741  }
742 
743  NFFT(init_guru)(&p, d, NN, M, nn, 2,
744  PRE_PHI_HUT |
745  FG_PSI |
746  MALLOC_F_HAT | MALLOC_X | MALLOC_F |
747  FFTW_INIT | FFT_OUT_OF_PLACE,
748  FFTW_MEASURE | FFTW_DESTROY_INPUT);
749 
750  p_fft = FFTW(plan_dft)(d, NN, p.f_hat, p.f, FFTW_FORWARD, FFTW_MEASURE);
751 
752  C *swapndft = (C*) NFFT(malloc)((size_t)(p.M_total) * sizeof(C));
753 
755  NFFT(vrand_shifted_unit_double)(p.x, p.d * p.M_total);
756 
757  //qsort(p.x,p.M_total,d*sizeof(double),comp1);
758  //nfft_vpr_double(p.x,p.M_total,"nodes x");
759 
760  NFFT(precompute_one_psi)(&p);
761 
763  NFFT(vrand_unit_complex)(p.f_hat, p.N_total);
764 
766  t_fft = K(0.0);
767  r = 0;
768  while (t_fft < K(0.1))
769  {
770  r++;
771  t0 = getticks();
772  FFTW(execute)(p_fft);
773  t1 = getticks();
774  t = NFFT(elapsed_seconds)(t1, t0);
775  t_fft += t;
776  }
777  t_fft /= (R)(r);
778 
779  printf("%.1" __FES__ "\t", t_fft);
780 
782  NFFT(vrand_unit_complex)(p.f_hat, p.N_total);
783 
785  if (test_ndft)
786  {
787  CSWAP(p.f, swapndft);
788  t_ndft = K(0.0);
789  r = 0;
790  while (t_ndft < K(0.1))
791  {
792  r++;
793  t0 = getticks();
794  NFFT(trafo_direct)(&p);
795  t1 = getticks();
796  t = NFFT(elapsed_seconds)(t1, t0);
797  t_ndft += t;
798  }
799  t_ndft /= (R)(r);
800  printf("%.1" __FES__ "\t", t_ndft);
801 
802  //printf("f=%e+i%e\t",creal(p.f[0]),cimag(p.f[0]));
803 
804  CSWAP(p.f, swapndft);
805  }
806  else
807  printf("\t");
808 
810  t_nfft = K(0.0);
811  r = 0;
812  while (t_nfft < K(0.1))
813  {
814  r++;
815  t0 = getticks();
816  NFFT(trafo)(&p);
817  t1 = getticks();
818  t = NFFT(elapsed_seconds)(t1, t0);
819  t_nfft += t;
820  }
821  t_nfft /= (R)(r);
822  printf("%.1" __FES__ "\t", t_nfft);
823  if (test_ndft)
824  printf("(%.1" __FES__ ")\t", NFFT(error_l_2_complex)(swapndft, p.f, p.M_total));
825 
826  //printf("f=%e+i%e\t",creal(p.f[0]),cimag(p.f[0]));
827 
829  t_nfft = K(0.0);
830  r = 0;
831  while (t_nfft < K(0.1))
832  {
833  r++;
834  t0 = getticks();
835  NFFT(trafo_3d)(&p);
836  t1 = getticks();
837  t = NFFT(elapsed_seconds)(t1, t0);
838  t_nfft += t;
839  }
840  t_nfft /= (R)(r);
841  printf("%.1" __FES__ "\t", t_nfft);
842  if (test_ndft)
843  printf("(%.1" __FES__ ")\t", NFFT(error_l_2_complex)(swapndft, p.f, p.M_total));
844 
845  //printf("f=%e+i%e\t",creal(p.f[0]),cimag(p.f[0]));
846 
847  printf("\n");
848 
849  NFFT(free)(swapndft);
850  FFTW(destroy_plan)(p_fft);
851  NFFT(finalize)(&p);
852 }
853 
854 static void measure_time_nfft_XXX7(int d, int N, unsigned test_ndft)
855 {
856  int r, M, NN[d], nn[d];
857  R t, t_fft, t_ndft, t_nfft;
858  ticks t0, t1;
859 
860  NFFT(plan) p;
861  FFTW(plan) p_fft;
862 
863  printf("%ld\t", LRINT(LOG((R)(N)) / LOG((R)(2)) * (R)(d) + K(0.5)));
864  fflush(stdout);
865 
866  for (r = 0, M = 1; r < d; r++)
867  {
868  M = N * M;
869  NN[r] = N;
870  nn[r] = 2 * N;
871  }
872 
873  NFFT(init_guru)(&p, d, NN, M, nn, 2,
874  PRE_PHI_HUT |
875  FG_PSI |
876  MALLOC_F_HAT | MALLOC_X | MALLOC_F |
877  FFTW_INIT | FFT_OUT_OF_PLACE,
878  FFTW_MEASURE | FFTW_DESTROY_INPUT);
879 
880  p_fft = FFTW(plan_dft)(d, NN, p.f, p.f_hat, FFTW_FORWARD, FFTW_MEASURE);
881 
882  C *swapndft = (C*) NFFT(malloc)((size_t)(p.N_total) * sizeof(C));
883 
885  NFFT(vrand_shifted_unit_double)(p.x, p.d * p.M_total);
886 
887  //sort_nodes(p.x,p.d,p.M_total,
888 
889  NFFT(precompute_one_psi)(&p);
890 
892  NFFT(vrand_unit_complex)(p.f, p.M_total);
893 
895  t_fft = K(0.0);
896  r = 0;
897  while (t_fft < K(0.1))
898  {
899  r++;
900  t0 = getticks();
901  FFTW(execute)(p_fft);
902  t1 = getticks();
903  t = NFFT(elapsed_seconds)(t1, t0);
904  t_fft += t;
905  }
906  t_fft /= (R)(r);
907 
908  printf("%.1" __FES__ "\t", t_fft);
909 
911  NFFT(vrand_unit_complex)(p.f, p.M_total);
912 
914  if (test_ndft)
915  {
916  CSWAP(p.f_hat, swapndft);
917  t_ndft = K(0.0);
918  r = 0;
919  while (t_ndft < K(0.1))
920  {
921  r++;
922  t0 = getticks();
923  NFFT(adjoint_direct)(&p);
924  t1 = getticks();
925  t = NFFT(elapsed_seconds)(t1, t0);
926  t_ndft += t;
927  }
928  t_ndft /= (R)(r);
929  printf("%.1" __FES__ "\t", t_ndft);
930 
931  //printf("\nf_hat=%e+i%e\t",creal(p.f_hat[0]),cimag(p.f_hat[0]));
932 
933  CSWAP(p.f_hat, swapndft);
934  }
935  else
936  printf("\t");
937 
939  t_nfft = K(0.0);
940  r = 0;
941  while (t_nfft < K(0.1))
942  {
943  r++;
944  t0 = getticks();
945  NFFT(adjoint)(&p);
946  t1 = getticks();
947  t = NFFT(elapsed_seconds)(t1, t0);
948  t_nfft += t;
949  }
950  t_nfft /= (R)(r);
951  printf("%.1" __FES__ "\t", t_nfft);
952  if (test_ndft)
953  printf("(%.1" __FES__ ")\t", NFFT(error_l_2_complex)(swapndft, p.f_hat, p.N_total));
954 
955  //printf("\nf_hat=%e+i%e\t",creal(p.f_hat[0]),cimag(p.f_hat[0]));
956 
958  t_nfft = K(0.0);
959  r = 0;
960  while (t_nfft < K(0.1))
961  {
962  r++;
963  t0 = getticks();
964  NFFT(adjoint_3d)(&p);
965  t1 = getticks();
966  t = NFFT(elapsed_seconds)(t1, t0);
967  t_nfft += t;
968  }
969  t_nfft /= (R)(r);
970  printf("%.1" __FES__ "\t", t_nfft);
971  if (test_ndft)
972  printf("(%.1" __FES__ ")\t", NFFT(error_l_2_complex)(swapndft, p.f_hat, p.N_total));
973 
974  //printf("\nf_hat=%e+i%e\t",creal(p.f_hat[0]),cimag(p.f_hat[0]));
975 
976  printf("\n");
977 
978  NFFT(free)(swapndft);
979  FFTW(destroy_plan)(p_fft);
980  NFFT(finalize)(&p);
981 }
982 
983 //static int main(void)
984 //{
985 // int l, d, logIN;
986 //
987 // for (l = 3; l <= 6; l++)
988 // {
989 // d = 3;
990 // logIN = d * l;
991 // int N = (int)(1U << (logIN / d));
992 // measure_time_nfft_XXX6(d, N, logIN <= 15 ? 1 : 0);
993 // measure_time_nfft_XXX7(d, N, logIN <= 15 ? 1 : 0);
994 // }
995 //
996 // for (l = 7; l <= 12; l++)
997 // {
998 // d = 2;
999 // logIN = d * l;
1000 // int N = (int)(1U << (logIN / d));
1001 // measure_time_nfft_XXX4(d, N, logIN <= 15 ? 1 : 0);
1002 // measure_time_nfft_XXX5(d, N, logIN <= 15 ? 1 : 0);
1003 // }
1004 //
1005 // for (l = 3; l <= 12; l++)
1006 // {
1007 // logIN = l;
1008 // int N = (int)(1U << (logIN));
1009 // measure_time_nfft_XXX2(1, N, logIN <= 15 ? 1 : 0);
1010 // measure_time_nfft_XXX3(1, N, logIN <= 15 ? 1 : 0);
1011 // }
1012 //
1013 // return EXIT_SUCCESS;
1014 //}
1015 
1016 int main(void)
1017 {
1018  int l, d, logIN;
1019 
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);
1026 
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++)
1030  {
1031  d = 1;
1032  logIN = l;
1033  int N = (int)(1U << (logIN / d));
1034  measure_time_nfft(d, N, logIN <= 15 ? 1 : 0);
1035 
1036  fflush(stdout);
1037  }
1038 
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++)
1042  {
1043  d = 2;
1044  logIN = d * l;
1045  int N = (int)(1U << (logIN / d));
1046  measure_time_nfft(d, N, logIN <= 15 ? 1 : 0);
1047 
1048  fflush(stdout);
1049  }
1050 
1051  printf("\\hline \\hline \\multicolumn{5}{|c|}{$d=3$} \\\\ \\hline\n");
1052  for (l = 3; l <= 7; l++)
1053  {
1054  d = 3;
1055  logIN = d * l;
1056  int N = (int)(1U << (logIN / d));
1057  measure_time_nfft(d, N, logIN <= 15 ? 1 : 0);
1058 
1059  fflush(stdout);
1060  }
1061 
1062  return 1;
1063 }
#define UNUSED(x)
Dummy use of unused parameters to silence compiler warnings.
Definition: infft.h:1383
#define CSWAP(x, y)
Swap two vectors.
Definition: infft.h:152