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