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