NFFT Logo 3.2.2
nfft_times.c
00001 /*
00002  * Copyright (c) 2002, 2012 Jens Keiner, Stefan Kunis, Daniel Potts
00003  *
00004  * This program is free software; you can redistribute it and/or modify it under
00005  * the terms of the GNU General Public License as published by the Free Software
00006  * Foundation; either version 2 of the License, or (at your option) any later
00007  * version.
00008  *
00009  * This program is distributed in the hope that it will be useful, but WITHOUT
00010  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00011  * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
00012  * details.
00013  *
00014  * You should have received a copy of the GNU General Public License along with
00015  * this program; if not, write to the Free Software Foundation, Inc., 51
00016  * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
00017  */
00018 
00019 /* $Id: nfft_times.c 3896 2012-10-10 12:19:26Z tovo $ */
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   //  printf("\\verb+%.1e+ & \\verb+(%.1f)+ &\t",t_fft,t_fft/(p.N_total*(log(N)/log(2)*d))*auxC);
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       //printf("\\verb+%.1e+ & \\verb+(%d)+&\t",t_ndft,(int)round(t_ndft/(p.N_total*p.N_total)*auxC));
00166       printf("\\verb+%.1e+ &\t",t_ndft);
00167     }
00168   else
00169     //    printf("\\verb+*+\t&\t&\t");
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   //  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);
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   //nfft_vpr_double(p.x,p.M_total,"nodes x");
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   //nfft_vpr_double(p.x,p.M_total,"nodes x");
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   //for(j=0;j<2*M;j+=2)
00473   //   p.x[j]=0.5*p.x[j]+0.25*(p.x[j]>=0?1:-1);
00474 
00475   //qsort(p.x,p.M_total,d*sizeof(double),comp1);
00476   //nfft_vpr_double(p.x,p.M_total,"nodes x");
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       //printf("f=%e+i%e\t",creal(p.f[0]),cimag(p.f[0]));
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   //printf("f=%e+i%e\t",creal(p.f[0]),cimag(p.f[0]));
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   //printf("f=%e+i%e\t",creal(p.f[0]),cimag(p.f[0]));
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   //sort_nodes(p.x,p.d,p.M_total,
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       //printf("\nf_hat=%e+i%e\t",creal(p.f_hat[0]),cimag(p.f_hat[0]));
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   //printf("\nf_hat=%e+i%e\t",creal(p.f_hat[0]),cimag(p.f_hat[0]));
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   //printf("\nf_hat=%e+i%e\t",creal(p.f_hat[0]),cimag(p.f_hat[0]));
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   //qsort(p.x,p.M_total,d*sizeof(double),comp1);
00735   //nfft_vpr_double(p.x,p.M_total,"nodes x");
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       //printf("f=%e+i%e\t",creal(p.f[0]),cimag(p.f[0]));
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   //printf("f=%e+i%e\t",creal(p.f[0]),cimag(p.f[0]));
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   //printf("f=%e+i%e\t",creal(p.f[0]),cimag(p.f[0]));
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   //sort_nodes(p.x,p.d,p.M_total,
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       //printf("\nf_hat=%e+i%e\t",creal(p.f_hat[0]),cimag(p.f_hat[0]));
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   //printf("\nf_hat=%e+i%e\t",creal(p.f_hat[0]),cimag(p.f_hat[0]));
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   //printf("\nf_hat=%e+i%e\t",creal(p.f_hat[0]),cimag(p.f_hat[0]));
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 }

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409