00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include "config.h"
00021
00022 #include <stdio.h>
00023 #include <math.h>
00024 #include <string.h>
00025 #include <stdlib.h>
00026 #ifdef HAVE_COMPLEX_H
00027 #include <complex.h>
00028 #endif
00029
00030 #include "nfft3util.h"
00031 #include "nfft3.h"
00032 #include "infft.h"
00033
00034 int global_n;
00035 int global_d;
00036
00037 static int comp1(const void *x,const void *y)
00038 {
00039 return ((* (double*) x)<(* (double*) y)?-1:1);
00040 }
00041
00042 static int comp2(const void *x,const void *y)
00043 {
00044 int nx0,nx1,ny0,ny1;
00045 nx0=global_n*(* ((double*)x+0));
00046 nx1=global_n*(* ((double*)x+1));
00047 ny0=global_n*(* ((double*)y+0));
00048 ny1=global_n*(* ((double*)y+1));
00049
00050 if(nx0<ny0)
00051 return -1;
00052 else
00053 if(nx0==ny0)
00054 if(nx1<ny1)
00055 return -1;
00056 else
00057 return 1;
00058 else
00059 return 1;
00060 }
00061
00062 static int comp3(const void *x,const void *y)
00063 {
00064 int nx0,nx1,nx2,ny0,ny1,ny2;
00065 nx0=global_n*(* ((double*)x+0));
00066 nx1=global_n*(* ((double*)x+1));
00067 nx2=global_n*(* ((double*)x+2));
00068 ny0=global_n*(* ((double*)y+0));
00069 ny1=global_n*(* ((double*)y+1));
00070 ny2=global_n*(* ((double*)y+2));
00071
00072 if(nx0<ny0)
00073 return -1;
00074 else
00075 if(nx0==ny0)
00076 if(nx1<ny1)
00077 return -1;
00078 else
00079 if(nx1==ny1)
00080 if(nx2<ny2)
00081 return -1;
00082 else
00083 return 1;
00084 else
00085 return 1;
00086 else
00087 return 1;
00088 }
00089
00090 void measure_time_nfft(int d, int N, unsigned test_ndft)
00091 {
00092 int r, M, NN[d], nn[d];
00093 double t, t_fft, t_ndft, t_nfft;
00094 ticks t0, t1;
00095
00096 nfft_plan p;
00097 fftw_plan p_fft;
00098
00099 printf("\\verb+%d+&\t",(int)(log(N)/log(2)*d+0.5));
00100
00101 for(r=0,M=1;r<d;r++)
00102 {
00103 M=N*M;
00104 NN[r]=N;
00105 nn[r]=2*N;
00106 }
00107
00108 nfft_init_guru(&p, d, NN, M, nn, 2,
00109 PRE_PHI_HUT|
00110 PRE_FULL_PSI| MALLOC_F_HAT| MALLOC_X| MALLOC_F|
00111 FFTW_INIT| FFT_OUT_OF_PLACE,
00112 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00113
00114 p_fft=fftw_plan_dft(d, NN, p.f_hat, p.f, FFTW_FORWARD, FFTW_MEASURE);
00115
00117 nfft_vrand_shifted_unit_double(p.x, p.d*p.M_total);
00118
00119 global_n=nn[0];
00120 global_d=d;
00121 switch(d)
00122 {
00123 case 1: { qsort(p.x,p.M_total,d*sizeof(double),comp1); break; }
00124 case 2: { qsort(p.x,p.M_total,d*sizeof(double),comp2); break; }
00125 case 3: { qsort(p.x,p.M_total,d*sizeof(double),comp3); break; }
00126 }
00127
00128 nfft_precompute_one_psi(&p);
00129
00131 nfft_vrand_unit_complex(p.f_hat, p.N_total);
00132
00134 t_fft=0;
00135 r=0;
00136 while(t_fft<1)
00137 {
00138 r++;
00139 t0 = getticks();
00140 fftw_execute(p_fft);
00141 t1 = getticks();
00142 t = nfft_elapsed_seconds(t1,t0);
00143 t_fft+=t;
00144 }
00145 t_fft/=r;
00146
00147
00148 printf("\\verb+%.1e+ &\t",t_fft);
00149
00151 if(test_ndft)
00152 {
00153 t_ndft=0;
00154 r=0;
00155 while(t_ndft<1)
00156 {
00157 r++;
00158 t0 = getticks();
00159 nfft_trafo_direct(&p);
00160 t1 = getticks();
00161 t = nfft_elapsed_seconds(t1,t0);
00162 t_ndft+=t;
00163 }
00164 t_ndft/=r;
00165
00166 printf("\\verb+%.1e+ &\t",t_ndft);
00167 }
00168 else
00169
00170 printf("\\verb+*+\t&\t");
00171
00173 t_nfft=0;
00174 r=0;
00175 while(t_nfft<1)
00176 {
00177 r++;
00178 t0 = getticks();
00179 switch(d)
00180 {
00181 case 1: { nfft_trafo_1d(&p); break; }
00182 case 2: { nfft_trafo_2d(&p); break; }
00183 case 3: { nfft_trafo_3d(&p); break; }
00184 default: nfft_trafo(&p);
00185 }
00186 t1 = getticks();
00187 t = nfft_elapsed_seconds(t1,t0);
00188 t_nfft+=t;
00189 }
00190 t_nfft/=r;
00191
00192
00193 printf("\\verb+%.1e+ & \\verb+(%3.1f)+\\\\\n",t_nfft,t_nfft/t_fft);
00194
00195 fftw_destroy_plan(p_fft);
00196 nfft_finalize(&p);
00197 }
00198
00199 void measure_time_nfft_XXX2(int d, int N, unsigned test_ndft)
00200 {
00201 int r, M, NN[d], nn[d];
00202 double t, t_fft, t_ndft, t_nfft;
00203 ticks t0, t1;
00204
00205 nfft_plan p;
00206 fftw_plan p_fft;
00207
00208 printf("%d\t",(int)(log(N)/log(2)*d+0.5)); fflush(stdout);
00209
00210 for(r=0,M=1;r<d;r++)
00211 {
00212 M=N*M;
00213 NN[r]=N;
00214 nn[r]=2*N;
00215 }
00216
00217 nfft_init_guru(&p, d, NN, M, nn, 2,
00218 PRE_PHI_HUT|
00219 PRE_PSI|
00220 MALLOC_F_HAT| MALLOC_X| MALLOC_F|
00221 FFTW_INIT| FFT_OUT_OF_PLACE,
00222 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00223
00224 p_fft=fftw_plan_dft(d, NN, p.f_hat, p.f, FFTW_FORWARD, FFTW_MEASURE);
00225
00226 double _Complex *swapndft=(double _Complex*)nfft_malloc(p.M_total*sizeof(double _Complex));
00227
00229 nfft_vrand_shifted_unit_double(p.x, p.d*p.M_total);
00230
00231 qsort(p.x,p.M_total,d*sizeof(double),comp1);
00232
00233
00234 nfft_precompute_one_psi(&p);
00235
00237 nfft_vrand_unit_complex(p.f_hat, p.N_total);
00238
00240 t_fft=0;
00241 r=0;
00242 while(t_fft<0.1)
00243 {
00244 r++;
00245 t0 = getticks();
00246 fftw_execute(p_fft);
00247 t1 = getticks();
00248 t = nfft_elapsed_seconds(t1,t0);
00249 t_fft+=t;
00250 }
00251 t_fft/=r;
00252
00253 printf("%.1e\t",t_fft);
00254
00256 if(test_ndft)
00257 {
00258 NFFT_SWAP_complex(p.f,swapndft);
00259 t_ndft=0;
00260 r=0;
00261 while(t_ndft<0.1)
00262 {
00263 r++;
00264 t0 = getticks();
00265 nfft_trafo_direct(&p);
00266 t1 = getticks();
00267 t = nfft_elapsed_seconds(t1,t0);
00268 t_ndft+=t;
00269 }
00270 t_ndft/=r;
00271 printf("%.1e\t",t_ndft);
00272 NFFT_SWAP_complex(p.f,swapndft);
00273 }
00274 else
00275 printf("\t");
00276
00278 t_nfft=0;
00279 r=0;
00280 while(t_nfft<0.1)
00281 {
00282 r++;
00283 t0 = getticks();
00284 nfft_trafo(&p);
00285 t1 = getticks();
00286 t = nfft_elapsed_seconds(t1,t0);
00287 t_nfft+=t;
00288 }
00289 t_nfft/=r;
00290 printf("%.1e\t",t_nfft);
00291 if(test_ndft)
00292 printf("(%.1e)\t",X(error_l_2_complex)(swapndft, p.f, p.M_total));
00293
00295 t_nfft=0;
00296 r=0;
00297 while(t_nfft<0.1)
00298 {
00299 r++;
00300 t0 = getticks();
00301 nfft_trafo_1d(&p);
00302 t1 = getticks();
00303 t = nfft_elapsed_seconds(t1,t0);
00304 t_nfft+=t;
00305 }
00306 t_nfft/=r;
00307 printf("%.1e\t",t_nfft);
00308 if(test_ndft)
00309 printf("(%.1e)\t",X(error_l_2_complex)(swapndft, p.f, p.M_total));
00310
00311 printf("\n");
00312
00313 nfft_free(swapndft);
00314 fftw_destroy_plan(p_fft);
00315 nfft_finalize(&p);
00316 }
00317
00318 void measure_time_nfft_XXX3(int d, int N, unsigned test_ndft)
00319 {
00320 int r, M, NN[d], nn[d];
00321 double t, t_fft, t_ndft, t_nfft;
00322 ticks t0, t1;
00323
00324 nfft_plan p;
00325 fftw_plan p_fft;
00326
00327 printf("%d\t",(int)(log(N)/log(2)*d+0.5)); fflush(stdout);
00328
00329 for(r=0,M=1;r<d;r++)
00330 {
00331 M=N*M;
00332 NN[r]=N;
00333 nn[r]=2*N;
00334 }
00335
00336 nfft_init_guru(&p, d, NN, M, nn, 2,
00337 PRE_PHI_HUT|
00338 PRE_PSI|
00339 MALLOC_F_HAT| MALLOC_X| MALLOC_F|
00340 FFTW_INIT| FFT_OUT_OF_PLACE,
00341 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00342
00343 p_fft=fftw_plan_dft(d, NN, p.f, p.f_hat, FFTW_BACKWARD, FFTW_MEASURE);
00344
00345 double _Complex *swapndft=(double _Complex*)nfft_malloc(p.N_total*sizeof(double _Complex));
00346
00348 nfft_vrand_shifted_unit_double(p.x, p.d*p.M_total);
00349
00350 qsort(p.x,p.M_total,d*sizeof(double),comp1);
00351
00352
00353 nfft_precompute_one_psi(&p);
00354
00356 nfft_vrand_unit_complex(p.f, p.N_total);
00357
00359 t_fft=0;
00360 r=0;
00361 while(t_fft<0.1)
00362 {
00363 r++;
00364 t0 = getticks();
00365 fftw_execute(p_fft);
00366 t1 = getticks();
00367 t = nfft_elapsed_seconds(t1,t0);
00368 t_fft+=t;
00369 }
00370 t_fft/=r;
00371
00372 printf("%.1e\t",t_fft);
00373
00375 if(test_ndft)
00376 {
00377 NFFT_SWAP_complex(p.f_hat,swapndft);
00378 t_ndft=0;
00379 r=0;
00380 while(t_ndft<0.1)
00381 {
00382 r++;
00383 t0 = getticks();
00384 nfft_adjoint_direct(&p);
00385 t1 = getticks();
00386 t = nfft_elapsed_seconds(t1,t0);
00387 t_ndft+=t;
00388 }
00389 t_ndft/=r;
00390 printf("%.1e\t",t_ndft);
00391 NFFT_SWAP_complex(p.f_hat,swapndft);
00392 }
00393 else
00394 printf("\t");
00395
00397 t_nfft=0;
00398 r=0;
00399 while(t_nfft<0.1)
00400 {
00401 r++;
00402 t0 = getticks();
00403 nfft_adjoint(&p);
00404 t1 = getticks();
00405 t = nfft_elapsed_seconds(t1,t0);
00406 t_nfft+=t;
00407 }
00408 t_nfft/=r;
00409 printf("%.1e\t",t_nfft);
00410 if(test_ndft)
00411 printf("(%.1e)\t",X(error_l_2_complex)(swapndft, p.f_hat, p.N_total));
00412
00414 t_nfft=0;
00415 r=0;
00416 while(t_nfft<0.1)
00417 {
00418 r++;
00419 t0 = getticks();
00420 nfft_adjoint_1d(&p);
00421 t1 = getticks();
00422 t = nfft_elapsed_seconds(t1,t0);
00423 t_nfft+=t;
00424 }
00425 t_nfft/=r;
00426 printf("%.1e\t",t_nfft);
00427 if(test_ndft)
00428 printf("(%.1e)\t",X(error_l_2_complex)(swapndft, p.f_hat, p.N_total));
00429
00430 printf("\n");
00431
00432 nfft_free(swapndft);
00433 fftw_destroy_plan(p_fft);
00434 nfft_finalize(&p);
00435 }
00436
00437
00438
00439
00440 void measure_time_nfft_XXX4(int d, int N, unsigned test_ndft)
00441 {
00442 int r, M, NN[d], nn[d];
00443 double t, t_fft, t_ndft, t_nfft;
00444 ticks t0, t1;
00445
00446 nfft_plan p;
00447 fftw_plan p_fft;
00448
00449 printf("%d\t",(int)(log(N)/log(2)*d+0.5)); fflush(stdout);
00450
00451 for(r=0,M=1;r<d;r++)
00452 {
00453 M=N*M;
00454 NN[r]=N;
00455 nn[r]=2*N;
00456 }
00457
00458 nfft_init_guru(&p, d, NN, M, nn, 4,
00459 PRE_PHI_HUT|
00460 PRE_PSI|
00461 MALLOC_F_HAT| MALLOC_X| MALLOC_F|
00462 FFTW_INIT| FFT_OUT_OF_PLACE,
00463 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00464
00465 p_fft=fftw_plan_dft(d, NN, p.f_hat, p.f, FFTW_FORWARD, FFTW_MEASURE);
00466
00467 double _Complex *swapndft=(double _Complex*)nfft_malloc(p.M_total*sizeof(double _Complex));
00468
00470 nfft_vrand_shifted_unit_double(p.x, p.d*p.M_total);
00471
00472
00473
00474
00475
00476
00477
00478 nfft_precompute_one_psi(&p);
00479
00481 nfft_vrand_unit_complex(p.f_hat, p.N_total);
00482
00484 t_fft=0;
00485 r=0;
00486 while(t_fft<0.1)
00487 {
00488 r++;
00489 t0 = getticks();
00490 fftw_execute(p_fft);
00491 t1 = getticks();
00492 t = nfft_elapsed_seconds(t1,t0);
00493 t_fft+=t;
00494 }
00495 t_fft/=r;
00496
00497 printf("%.1e\t",t_fft);
00498
00500 nfft_vrand_unit_complex(p.f_hat, p.N_total);
00501
00503 if(test_ndft)
00504 {
00505 NFFT_SWAP_complex(p.f,swapndft);
00506 t_ndft=0;
00507 r=0;
00508 while(t_ndft<0.1)
00509 {
00510 r++;
00511 t0 = getticks();
00512 nfft_trafo_direct(&p);
00513 t1 = getticks();
00514 t = nfft_elapsed_seconds(t1,t0);
00515 t_ndft+=t;
00516 }
00517 t_ndft/=r;
00518 printf("%.1e\t",t_ndft);
00519
00520
00521
00522 NFFT_SWAP_complex(p.f,swapndft);
00523 }
00524 else
00525 printf("\t");
00526
00528 t_nfft=0;
00529 r=0;
00530 while(t_nfft<0.1)
00531 {
00532 r++;
00533 t0 = getticks();
00534 nfft_trafo(&p);
00535 t1 = getticks();
00536 t = nfft_elapsed_seconds(t1,t0);
00537 t_nfft+=t;
00538 }
00539 t_nfft/=r;
00540 printf("%.1e\t",t_nfft);
00541 if(test_ndft)
00542 printf("(%.1e)\t",X(error_l_2_complex)(swapndft, p.f, p.M_total));
00543
00544
00545
00547 t_nfft=0;
00548 r=0;
00549 while(t_nfft<0.1)
00550 {
00551 r++;
00552 t0 = getticks();
00553 nfft_trafo_2d(&p);
00554 t1 = getticks();
00555 t = nfft_elapsed_seconds(t1,t0);
00556 t_nfft+=t;
00557 }
00558 t_nfft/=r;
00559 printf("%.1e\t",t_nfft);
00560 if(test_ndft)
00561 printf("(%.1e)\t",X(error_l_2_complex)(swapndft, p.f, p.M_total));
00562
00563
00564
00565 printf("\n");
00566
00567 nfft_free(swapndft);
00568 fftw_destroy_plan(p_fft);
00569 nfft_finalize(&p);
00570 }
00571
00572
00573 void measure_time_nfft_XXX5(int d, int N, unsigned test_ndft)
00574 {
00575 int r, M, NN[d], nn[d];
00576 double t, t_fft, t_ndft, t_nfft;
00577 ticks t0, t1;
00578
00579 nfft_plan p;
00580 fftw_plan p_fft;
00581
00582 printf("%d\t",(int)(log(N)/log(2)*d+0.5)); fflush(stdout);
00583
00584 for(r=0,M=1;r<d;r++)
00585 {
00586 M=N*M;
00587 NN[r]=N;
00588 nn[r]=2*N;
00589 }
00590
00591 nfft_init_guru(&p, d, NN, M, nn, 4,
00592 PRE_PHI_HUT|
00593 PRE_PSI|
00594 MALLOC_F_HAT| MALLOC_X| MALLOC_F|
00595 FFTW_INIT| FFT_OUT_OF_PLACE,
00596 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00597
00598 p_fft=fftw_plan_dft(d, NN, p.f, p.f_hat, FFTW_FORWARD, FFTW_MEASURE);
00599
00600 double _Complex *swapndft=(double _Complex*)nfft_malloc(p.N_total*sizeof(double _Complex));
00601
00603 nfft_vrand_shifted_unit_double(p.x, p.d*p.M_total);
00604
00605
00606
00607 nfft_precompute_one_psi(&p);
00608
00610 nfft_vrand_unit_complex(p.f, p.M_total);
00611
00613 t_fft=0;
00614 r=0;
00615 while(t_fft<0.1)
00616 {
00617 r++;
00618 t0 = getticks();
00619 fftw_execute(p_fft);
00620 t1 = getticks();
00621 t = nfft_elapsed_seconds(t1,t0);
00622 t_fft+=t;
00623 }
00624 t_fft/=r;
00625
00626 printf("%.1e\t",t_fft);
00627
00629 nfft_vrand_unit_complex(p.f, p.M_total);
00630
00632 if(test_ndft)
00633 {
00634 NFFT_SWAP_complex(p.f_hat,swapndft);
00635 t_ndft=0;
00636 r=0;
00637 while(t_ndft<0.1)
00638 {
00639 r++;
00640 t0 = getticks();
00641 nfft_adjoint_direct(&p);
00642 t1 = getticks();
00643 t = nfft_elapsed_seconds(t1,t0);
00644 t_ndft+=t;
00645 }
00646 t_ndft/=r;
00647 printf("%.1e\t",t_ndft);
00648
00649
00650
00651 NFFT_SWAP_complex(p.f_hat,swapndft);
00652 }
00653 else
00654 printf("\t");
00655
00657 t_nfft=0;
00658 r=0;
00659 while(t_nfft<0.1)
00660 {
00661 r++;
00662 t0 = getticks();
00663 nfft_adjoint(&p);
00664 t1 = getticks();
00665 t = nfft_elapsed_seconds(t1,t0);
00666 t_nfft+=t;
00667 }
00668 t_nfft/=r;
00669 printf("%.1e\t",t_nfft);
00670 if(test_ndft)
00671 printf("(%.1e)\t",X(error_l_2_complex)(swapndft, p.f_hat, p.N_total));
00672
00673
00674
00676 t_nfft=0;
00677 r=0;
00678 while(t_nfft<0.1)
00679 {
00680 r++;
00681 t0 = getticks();
00682 nfft_adjoint_2d(&p);
00683 t1 = getticks();
00684 t = nfft_elapsed_seconds(t1,t0);
00685 t_nfft+=t;
00686 }
00687 t_nfft/=r;
00688 printf("%.1e\t",t_nfft);
00689 if(test_ndft)
00690 printf("(%.1e)\t",X(error_l_2_complex)(swapndft, p.f_hat, p.N_total));
00691
00692
00693
00694 printf("\n");
00695
00696 nfft_free(swapndft);
00697 fftw_destroy_plan(p_fft);
00698 nfft_finalize(&p);
00699 }
00700
00701
00702 void measure_time_nfft_XXX6(int d, int N, unsigned test_ndft)
00703 {
00704 int r, M, NN[d], nn[d];
00705 double t, t_fft, t_ndft, t_nfft;
00706 ticks t0, t1;
00707
00708 nfft_plan p;
00709 fftw_plan p_fft;
00710
00711 printf("%d\t",(int)(log(N)/log(2)*d+0.5)); fflush(stdout);
00712
00713 for(r=0,M=1;r<d;r++)
00714 {
00715 M=N*M;
00716 NN[r]=N;
00717 nn[r]=2*N;
00718 }
00719
00720 nfft_init_guru(&p, d, NN, M, nn, 2,
00721 PRE_PHI_HUT|
00722 FG_PSI|
00723 MALLOC_F_HAT| MALLOC_X| MALLOC_F|
00724 FFTW_INIT| FFT_OUT_OF_PLACE,
00725 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00726
00727 p_fft=fftw_plan_dft(d, NN, p.f_hat, p.f, FFTW_FORWARD, FFTW_MEASURE);
00728
00729 double _Complex *swapndft=(double _Complex*)nfft_malloc(p.M_total*sizeof(double _Complex));
00730
00732 nfft_vrand_shifted_unit_double(p.x, p.d*p.M_total);
00733
00734
00735
00736
00737 nfft_precompute_one_psi(&p);
00738
00740 nfft_vrand_unit_complex(p.f_hat, p.N_total);
00741
00743 t_fft=0;
00744 r=0;
00745 while(t_fft<0.1)
00746 {
00747 r++;
00748 t0 = getticks();
00749 fftw_execute(p_fft);
00750 t1 = getticks();
00751 t = nfft_elapsed_seconds(t1,t0);
00752 t_fft+=t;
00753 }
00754 t_fft/=r;
00755
00756 printf("%.1e\t",t_fft);
00757
00759 nfft_vrand_unit_complex(p.f_hat, p.N_total);
00760
00762 if(test_ndft)
00763 {
00764 NFFT_SWAP_complex(p.f,swapndft);
00765 t_ndft=0;
00766 r=0;
00767 while(t_ndft<0.1)
00768 {
00769 r++;
00770 t0 = getticks();
00771 nfft_trafo_direct(&p);
00772 t1 = getticks();
00773 t = nfft_elapsed_seconds(t1,t0);
00774 t_ndft+=t;
00775 }
00776 t_ndft/=r;
00777 printf("%.1e\t",t_ndft);
00778
00779
00780
00781 NFFT_SWAP_complex(p.f,swapndft);
00782 }
00783 else
00784 printf("\t");
00785
00787 t_nfft=0;
00788 r=0;
00789 while(t_nfft<0.1)
00790 {
00791 r++;
00792 t0 = getticks();
00793 nfft_trafo(&p);
00794 t1 = getticks();
00795 t = nfft_elapsed_seconds(t1,t0);
00796 t_nfft+=t;
00797 }
00798 t_nfft/=r;
00799 printf("%.1e\t",t_nfft);
00800 if(test_ndft)
00801 printf("(%.1e)\t",X(error_l_2_complex)(swapndft, p.f, p.M_total));
00802
00803
00804
00806 t_nfft=0;
00807 r=0;
00808 while(t_nfft<0.1)
00809 {
00810 r++;
00811 t0 = getticks();
00812 nfft_trafo_3d(&p);
00813 t1 = getticks();
00814 t = nfft_elapsed_seconds(t1,t0);
00815 t_nfft+=t;
00816 }
00817 t_nfft/=r;
00818 printf("%.1e\t",t_nfft);
00819 if(test_ndft)
00820 printf("(%.1e)\t",X(error_l_2_complex)(swapndft, p.f, p.M_total));
00821
00822
00823
00824 printf("\n");
00825
00826 nfft_free(swapndft);
00827 fftw_destroy_plan(p_fft);
00828 nfft_finalize(&p);
00829 }
00830
00831
00832 void measure_time_nfft_XXX7(int d, int N, unsigned test_ndft)
00833 {
00834 int r, M, NN[d], nn[d];
00835 double t, t_fft, t_ndft, t_nfft;
00836 ticks t0, t1;
00837
00838 nfft_plan p;
00839 fftw_plan p_fft;
00840
00841 printf("%d\t",(int)(log(N)/log(2)*d+0.5)); fflush(stdout);
00842
00843 for(r=0,M=1;r<d;r++)
00844 {
00845 M=N*M;
00846 NN[r]=N;
00847 nn[r]=2*N;
00848 }
00849
00850 nfft_init_guru(&p, d, NN, M, nn, 2,
00851 PRE_PHI_HUT|
00852 FG_PSI|
00853 MALLOC_F_HAT| MALLOC_X| MALLOC_F|
00854 FFTW_INIT| FFT_OUT_OF_PLACE,
00855 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00856
00857 p_fft=fftw_plan_dft(d, NN, p.f, p.f_hat, FFTW_FORWARD, FFTW_MEASURE);
00858
00859 double _Complex *swapndft=(double _Complex*)nfft_malloc(p.N_total*sizeof(double _Complex));
00860
00862 nfft_vrand_shifted_unit_double(p.x, p.d*p.M_total);
00863
00864
00865
00866 nfft_precompute_one_psi(&p);
00867
00869 nfft_vrand_unit_complex(p.f, p.M_total);
00870
00872 t_fft=0;
00873 r=0;
00874 while(t_fft<0.1)
00875 {
00876 r++;
00877 t0 = getticks();
00878 fftw_execute(p_fft);
00879 t1 = getticks();
00880 t = nfft_elapsed_seconds(t1,t0);
00881 t_fft+=t;
00882 }
00883 t_fft/=r;
00884
00885 printf("%.1e\t",t_fft);
00886
00888 nfft_vrand_unit_complex(p.f, p.M_total);
00889
00891 if(test_ndft)
00892 {
00893 NFFT_SWAP_complex(p.f_hat,swapndft);
00894 t_ndft=0;
00895 r=0;
00896 while(t_ndft<0.1)
00897 {
00898 r++;
00899 t0 = getticks();
00900 nfft_adjoint_direct(&p);
00901 t1 = getticks();
00902 t = nfft_elapsed_seconds(t1,t0);
00903 t_ndft+=t;
00904 }
00905 t_ndft/=r;
00906 printf("%.1e\t",t_ndft);
00907
00908
00909
00910 NFFT_SWAP_complex(p.f_hat,swapndft);
00911 }
00912 else
00913 printf("\t");
00914
00916 t_nfft=0;
00917 r=0;
00918 while(t_nfft<0.1)
00919 {
00920 r++;
00921 t0 = getticks();
00922 nfft_adjoint(&p);
00923 t1 = getticks();
00924 t = nfft_elapsed_seconds(t1,t0);
00925 t_nfft+=t;
00926 }
00927 t_nfft/=r;
00928 printf("%.1e\t",t_nfft);
00929 if(test_ndft)
00930 printf("(%.1e)\t",X(error_l_2_complex)(swapndft, p.f_hat, p.N_total));
00931
00932
00933
00935 t_nfft=0;
00936 r=0;
00937 while(t_nfft<0.1)
00938 {
00939 r++;
00940 t0 = getticks();
00941 nfft_adjoint_3d(&p);
00942 t1 = getticks();
00943 t = nfft_elapsed_seconds(t1,t0);
00944 t_nfft+=t;
00945 }
00946 t_nfft/=r;
00947 printf("%.1e\t",t_nfft);
00948 if(test_ndft)
00949 printf("(%.1e)\t",X(error_l_2_complex)(swapndft, p.f_hat, p.N_total));
00950
00951
00952
00953 printf("\n");
00954
00955 nfft_free(swapndft);
00956 fftw_destroy_plan(p_fft);
00957 nfft_finalize(&p);
00958 }
00959
00960 int main2(void)
00961 {
00962 int l,d,logIN;
00963
00964 for(l=3;l<=6;l++)
00965 {
00966 d=3;
00967 logIN=d*l;
00968 if(logIN<=15)
00969 {
00970 measure_time_nfft_XXX6(d,(1U<< (logIN/d)),1);
00971 measure_time_nfft_XXX7(d,(1U<< (logIN/d)),1);
00972 }
00973 else
00974 {
00975 measure_time_nfft_XXX6(d,(1U<< (logIN/d)),0);
00976 measure_time_nfft_XXX7(d,(1U<< (logIN/d)),0);
00977 }
00978 }
00979
00980 exit(-1);
00981
00982
00983 for(l=7;l<=12;l++)
00984 {
00985 d=2;
00986 logIN=d*l;
00987 if(logIN<=15)
00988 {
00989 measure_time_nfft_XXX4(d,(1U<< (logIN/d)),1);
00990 measure_time_nfft_XXX5(d,(1U<< (logIN/d)),1);
00991 }
00992 else
00993 {
00994 measure_time_nfft_XXX4(d,(1U<< (logIN/d)),0);
00995 measure_time_nfft_XXX5(d,(1U<< (logIN/d)),0);
00996 }
00997 }
00998
00999 exit(-1);
01000
01001 for(l=3;l<=12;l++)
01002 {
01003 logIN=l;
01004 if(logIN<=15)
01005 {
01006 measure_time_nfft_XXX2(1,(1U<< (logIN)),1);
01007 measure_time_nfft_XXX3(1,(1U<< (logIN)),1);
01008 }
01009 else
01010 {
01011 measure_time_nfft_XXX2(1,(1U<< (logIN)),0);
01012 measure_time_nfft_XXX3(1,(1U<< (logIN)),0);
01013 }
01014 }
01015
01016 exit(-1);
01017 }
01018
01019 int main(void)
01020 {
01021 int l,d,logIN;
01022
01023 printf("\\hline $l_N$ & FFT & NDFT & NFFT & NFFT/FFT\\\\\n");
01024 printf("\\hline \\hline \\multicolumn{5}{|c|}{$d=1$} \\\\ \\hline\n");
01025 for(l=3;l<=22;l++)
01026 {
01027 d=1;
01028 logIN=l;
01029 if(logIN<=15)
01030 measure_time_nfft(d,(1U<< (logIN/d)),1);
01031 else
01032 measure_time_nfft(d,(1U<< (logIN/d)),0);
01033
01034 fflush(stdout);
01035 }
01036
01037 printf("\\hline $l_N$ & FFT & NDFT & NFFT & NFFT/FFT\\\\\n");
01038 printf("\\hline \\hline \\multicolumn{5}{|c|}{$d=2$} \\\\ \\hline\n");
01039 for(l=3;l<=11;l++)
01040 {
01041 d=2;
01042 logIN=d*l;
01043 if(logIN<=15)
01044 measure_time_nfft(d,(1U<< (logIN/d)),1);
01045 else
01046 measure_time_nfft(d,(1U<< (logIN/d)),0);
01047
01048 fflush(stdout);
01049 }
01050
01051 printf("\\hline \\hline \\multicolumn{5}{|c|}{$d=3$} \\\\ \\hline\n");
01052 for(l=3;l<=7;l++)
01053 {
01054 d=3;
01055 logIN=d*l;
01056 if(logIN<=15)
01057 measure_time_nfft(d,(1U<< (logIN/d)),1);
01058 else
01059 measure_time_nfft(d,(1U<< (logIN/d)),0);
01060
01061 fflush(stdout);
01062 }
01063
01064 return 1;
01065 }