NFFT Logo 3.2.2
mpolar_fft_test.c
Go to the documentation of this file.
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: mpolar_fft_test.c 3896 2012-10-10 12:19:26Z tovo $ */
00020 
00030 #include "config.h"
00031 
00032 #include <math.h>
00033 #include <stdlib.h>
00034 #ifdef HAVE_COMPLEX_H
00035 #include <complex.h>
00036 #endif
00037 
00038 #include "nfft3util.h"
00039 #include "nfft3.h"
00040 #include "infft.h"
00041 
00048 double GLOBAL_elapsed_time;
00049 
00063 static int mpolar_grid(int T, int R, double *x, double *w)
00064 {
00065   int t, r;
00066   double W;
00067   int R2=2*ceil(sqrt(2.0)*R/2);
00068   double xx, yy;
00069   int M=0;
00070 
00071   for(t=-T/2; t<T/2; t++)
00072   {
00073     for(r=-R2/2; r<R2/2; r++)
00074     {
00075       xx = (double)r/R*cos(PI*t/T);
00076       yy = (double)r/R*sin(PI*t/T);
00077 
00078       if ( ((-0.5-1.0/(double)R)<=xx) & (xx<=(0.5+1.0/(double)R)) &
00079         ((-0.5-1.0/(double)R)<=yy) & (yy<=(0.5+1.0/(double)R)) )
00080       {
00081         x[2*M+0] = xx;
00082         x[2*M+1] = yy;
00083 
00084         if (r==0)
00085           w[M] = 1.0/4.0;
00086         else
00087           w[M] = fabs((double)r);
00088 
00089         M++; 
00090       }
00091     }
00092   }
00093 
00095    W=0.0;
00096    for (t=0; t<M; t++)
00097       W+=w[t];
00098 
00099    for (t=0; t<M; t++)
00100     w[t]/=W;
00101 
00102   return M;                             
00103 }
00104 
00106 static int mpolar_dft(fftw_complex *f_hat, int NN, fftw_complex *f, int T, int R, int m)
00107 {
00108   ticks t0, t1;
00109   int j,k;                              
00110   nfft_plan my_nfft_plan;               
00112   double *x, *w;                        
00114   int N[2],n[2];
00115   int M;                                
00117   N[0]=NN; n[0]=2*N[0];                 
00118   N[1]=NN; n[1]=2*N[1];                 
00120   x = (double *)nfft_malloc(5*(T/2)*R*(sizeof(double)));
00121   if (x==NULL)
00122     return -1;
00123 
00124   w = (double *)nfft_malloc(5*(T*R)/4*(sizeof(double)));
00125   if (w==NULL)
00126     return -1;
00127 
00129   M=mpolar_grid(T,R,x,w);
00130   nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
00131                   PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00132                   FFTW_MEASURE| FFTW_DESTROY_INPUT);
00133 
00135   for(j=0;j<my_nfft_plan.M_total;j++)
00136   {
00137     my_nfft_plan.x[2*j+0] = x[2*j+0];
00138     my_nfft_plan.x[2*j+1] = x[2*j+1];
00139   }
00140 
00142   for(k=0;k<my_nfft_plan.N_total;k++)
00143     my_nfft_plan.f_hat[k] = f_hat[k];
00144 
00145   t0 = getticks();
00146 
00148   nfft_trafo_direct(&my_nfft_plan);
00149 
00150   t1 = getticks();
00151   GLOBAL_elapsed_time = nfft_elapsed_seconds(t1,t0);
00152 
00154   for(j=0;j<my_nfft_plan.M_total;j++)
00155     f[j] = my_nfft_plan.f[j];
00156 
00158   nfft_finalize(&my_nfft_plan);
00159   nfft_free(x);
00160   nfft_free(w);
00161 
00162   return EXIT_SUCCESS;
00163 }
00164 
00166 static int mpolar_fft(fftw_complex *f_hat, int NN, fftw_complex *f, int T, int R, int m)
00167 {
00168   ticks t0, t1;
00169   int j,k;                              
00170   nfft_plan my_nfft_plan;               
00172   double *x, *w;                        
00174   int N[2],n[2];
00175   int M;                                
00177   N[0]=NN; n[0]=2*N[0];                 
00178   N[1]=NN; n[1]=2*N[1];                 
00180   x = (double *)nfft_malloc(5*T*R/2*(sizeof(double)));
00181   if (x==NULL)
00182     return -1;
00183 
00184   w = (double *)nfft_malloc(5*T*R/4*(sizeof(double)));
00185   if (w==NULL)
00186     return -1;
00187 
00189   M=mpolar_grid(T,R,x,w);
00190   nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
00191                   PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00192                   FFTW_MEASURE| FFTW_DESTROY_INPUT);
00193 
00196   for(j=0;j<my_nfft_plan.M_total;j++)
00197   {
00198     my_nfft_plan.x[2*j+0] = x[2*j+0];
00199     my_nfft_plan.x[2*j+1] = x[2*j+1];
00200   }
00201 
00203   if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
00204     nfft_precompute_lin_psi(&my_nfft_plan);
00205 
00206   if(my_nfft_plan.nfft_flags & PRE_PSI)
00207     nfft_precompute_psi(&my_nfft_plan);
00208 
00209   if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
00210     nfft_precompute_full_psi(&my_nfft_plan);
00211 
00213   for(k=0;k<my_nfft_plan.N_total;k++)
00214     my_nfft_plan.f_hat[k] = f_hat[k];
00215 
00216   t0 = getticks();
00217 
00219   nfft_trafo(&my_nfft_plan);
00220 
00221   t1 = getticks();
00222   GLOBAL_elapsed_time = nfft_elapsed_seconds(t1,t0);
00223 
00225   for(j=0;j<my_nfft_plan.M_total;j++)
00226     f[j] = my_nfft_plan.f[j];
00227 
00229   nfft_finalize(&my_nfft_plan);
00230   nfft_free(x);
00231   nfft_free(w);
00232 
00233   return EXIT_SUCCESS;
00234 }
00235 
00237 static int inverse_mpolar_fft(fftw_complex *f, int T, int R, fftw_complex *f_hat, int NN, int max_i, int m)
00238 {
00239   ticks t0, t1;
00240   int j,k;                              
00241   nfft_plan my_nfft_plan;               
00242   solver_plan_complex my_infft_plan;             
00244   double *x, *w;                        
00245   int l;                                
00247   int N[2],n[2];
00248   int M;                                
00250   N[0]=NN; n[0]=2*N[0];                 
00251   N[1]=NN; n[1]=2*N[1];                 
00253   x = (double *)nfft_malloc(5*T*R/2*(sizeof(double)));
00254   if (x==NULL)
00255     return -1;
00256 
00257   w = (double *)nfft_malloc(5*T*R/4*(sizeof(double)));
00258   if (w==NULL)
00259     return -1;
00260 
00262   M=mpolar_grid(T,R,x,w);
00263   nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
00264                   PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00265                   FFTW_MEASURE| FFTW_DESTROY_INPUT);
00266 
00268     solver_init_advanced_complex(&my_infft_plan,(nfft_mv_plan_complex*)(&my_nfft_plan), CGNR | PRECOMPUTE_WEIGHT );
00269 
00271   for(j=0;j<my_nfft_plan.M_total;j++)
00272   {
00273     my_nfft_plan.x[2*j+0] = x[2*j+0];
00274     my_nfft_plan.x[2*j+1] = x[2*j+1];
00275     my_infft_plan.y[j]    = f[j];
00276     my_infft_plan.w[j]    = w[j];
00277   }
00278 
00280   if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
00281     nfft_precompute_lin_psi(&my_nfft_plan);
00282 
00283   if(my_nfft_plan.nfft_flags & PRE_PSI)
00284     nfft_precompute_psi(&my_nfft_plan);
00285 
00286   if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
00287     nfft_precompute_full_psi(&my_nfft_plan);
00288 
00289 
00291  if(my_infft_plan.flags & PRECOMPUTE_DAMP)
00292    for(j=0;j<my_nfft_plan.N[0];j++)
00293      for(k=0;k<my_nfft_plan.N[1];k++)
00294      {
00295         my_infft_plan.w_hat[j*my_nfft_plan.N[1]+k]=
00296           (sqrt(pow(j-my_nfft_plan.N[0]/2,2)+pow(k-my_nfft_plan.N[1]/2,2))>(my_nfft_plan.N[0]/2)?0:1);
00297      }
00298 
00300   for(k=0;k<my_nfft_plan.N_total;k++)
00301       my_infft_plan.f_hat_iter[k] = 0.0 + _Complex_I*0.0;
00302 
00303   t0 = getticks();
00304 
00306   solver_before_loop_complex(&my_infft_plan);
00307 
00308   if (max_i<1)
00309   {
00310     l=1;
00311     for(k=0;k<my_nfft_plan.N_total;k++)
00312       my_infft_plan.f_hat_iter[k] = my_infft_plan.p_hat_iter[k];
00313   }
00314   else
00315   {
00316     for(l=1;l<=max_i;l++)
00317     {
00318       solver_loop_one_step_complex(&my_infft_plan);
00319     }
00320   }
00321 
00322   t1 = getticks();
00323   GLOBAL_elapsed_time = nfft_elapsed_seconds(t1,t0);
00324 
00326   for(k=0;k<my_nfft_plan.N_total;k++)
00327     f_hat[k] = my_infft_plan.f_hat_iter[k];
00328 
00330   solver_finalize_complex(&my_infft_plan);
00331   nfft_finalize(&my_nfft_plan);
00332   nfft_free(x);
00333   nfft_free(w);
00334 
00335   return EXIT_SUCCESS;
00336 }
00337 
00339 static int comparison_fft(FILE *fp, int N, int T, int R)
00340 {
00341   ticks t0, t1;
00342   fftw_plan my_fftw_plan;
00343   fftw_complex *f_hat,*f;
00344   int m,k;
00345   double t_fft, t_dft_mpolar;
00346 
00347   f_hat = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*N*N);
00348   f     = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*(T*R/4)*5);
00349 
00350   my_fftw_plan = fftw_plan_dft_2d(N,N,f_hat,f,FFTW_BACKWARD,FFTW_MEASURE);
00351 
00352   for(k=0; k<N*N; k++)
00353     f_hat[k] = (((double)rand())/RAND_MAX) + _Complex_I* (((double)rand())/RAND_MAX);
00354 
00355   t0 = getticks();
00356   for(m=0;m<65536/N;m++)
00357     {
00358       fftw_execute(my_fftw_plan);
00359       /* touch */
00360       f_hat[2]=2*f_hat[0];
00361     }
00362   t1 = getticks();
00363   GLOBAL_elapsed_time = nfft_elapsed_seconds(t1,t0);
00364   t_fft=N*GLOBAL_elapsed_time/65536;
00365 
00366   if(N<256)
00367     {
00368       mpolar_dft(f_hat,N,f,T,R,1);
00369       t_dft_mpolar=GLOBAL_elapsed_time;
00370     }
00371 
00372   for (m=3; m<=9; m+=3)
00373     {
00374       if((m==3)&&(N<256))
00375         fprintf(fp,"%d\t&\t&\t%1.1e&\t%1.1e&\t%d\t",N,t_fft,t_dft_mpolar,m);
00376       else
00377         if(m==3)
00378     fprintf(fp,"%d\t&\t&\t%1.1e&\t       &\t%d\t",N,t_fft,m);
00379   else
00380     fprintf(fp,"  \t&\t&\t       &\t       &\t%d\t",m);
00381 
00382       printf("N=%d\tt_fft=%1.1e\tt_dft_mpolar=%1.1e\tm=%d\t",N,t_fft,t_dft_mpolar,m);
00383 
00384       mpolar_fft(f_hat,N,f,T,R,m);
00385       fprintf(fp,"%1.1e&\t",GLOBAL_elapsed_time);
00386       printf("t_mpolar=%1.1e\t",GLOBAL_elapsed_time);
00387       inverse_mpolar_fft(f,T,R,f_hat,N,2*m,m);
00388       if(m==9)
00389   fprintf(fp,"%1.1e\\\\\\hline\n",GLOBAL_elapsed_time);
00390       else
00391   fprintf(fp,"%1.1e\\\\\n",GLOBAL_elapsed_time);
00392       printf("t_impolar=%1.1e\n",GLOBAL_elapsed_time);
00393     }
00394 
00395   fflush(fp);
00396 
00397   nfft_free(f);
00398   nfft_free(f_hat);
00399 
00400   return EXIT_SUCCESS;
00401 }
00402 
00404 int main(int argc,char **argv)
00405 {
00406   int N;                                
00407   int T, R;                             
00408   int M;                                
00409   double *x, *w;                        
00410   fftw_complex *f_hat, *f, *f_direct, *f_tilde;
00411   int k;
00412   int max_i;                            
00413   int m;
00414   double temp1, temp2, E_max=0.0;
00415   FILE *fp1, *fp2;
00416   char filename[30];
00417   int logN;
00418 
00419   if( argc!=4 )
00420   {
00421     printf("mpolar_fft_test N T R \n");
00422     printf("\n");
00423     printf("N          mpolar FFT of size NxN    \n");
00424     printf("T          number of slopes          \n");
00425     printf("R          number of offsets         \n");
00426 
00428     printf("\nHence, comparison FFTW, mpolar FFT and inverse mpolar FFT\n");
00429     fp1=fopen("mpolar_comparison_fft.dat","w");
00430     if (fp1==NULL)
00431   return(-1);
00432     for (logN=4; logN<=8; logN++)
00433   comparison_fft(fp1,(1U<< logN), 3*(1U<< logN), 3*(1U<< (logN-1)));
00434     fclose(fp1);
00435 
00436     exit(-1);
00437   }
00438 
00439   N = atoi(argv[1]);
00440   T = atoi(argv[2]);
00441   R = atoi(argv[3]);
00442   printf("N=%d, modified polar grid with T=%d, R=%d => ",N,T,R);
00443 
00444   x = (double *)nfft_malloc(5*T*R/2*(sizeof(double)));
00445   w = (double *)nfft_malloc(5*T*R/4*(sizeof(double)));
00446 
00447   f_hat    = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*N*N);
00448   f        = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*1.25*T*R);  /* 4/pi*log(1+sqrt(2)) = 1.122... < 1.25 */
00449   f_direct = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*1.25*T*R);
00450   f_tilde  = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*N*N);
00451 
00453   M=mpolar_grid(T,R,x,w); printf("M=%d.\n",M);
00454 
00456   fp1=fopen("input_data_r.dat","r");
00457   fp2=fopen("input_data_i.dat","r");
00458   if ((fp1==NULL) || (fp2==NULL))
00459     return(-1);
00460   for(k=0;k<N*N;k++)
00461   {
00462     fscanf(fp1,"%le ",&temp1);
00463     fscanf(fp2,"%le ",&temp2);
00464     f_hat[k]=temp1+ _Complex_I*temp2;
00465   }
00466   fclose(fp1);
00467   fclose(fp2);
00468 
00470       mpolar_dft(f_hat,N,f_direct,T,R,1);
00471   //  mpolar_fft(f_hat,N,f_direct,T,R,12);
00472 
00474   printf("\nTest of the mpolar FFT: \n");
00475   fp1=fopen("mpolar_fft_error.dat","w+");
00476   for (m=1; m<=12; m++)
00477   {
00479     mpolar_fft(f_hat,N,f,T,R,m);
00480 
00482     E_max=X(error_l_infty_complex)(f_direct,f,M);
00483     printf("m=%2d: E_max = %e\n",m,E_max);
00484     fprintf(fp1,"%e\n",E_max);
00485   }
00486   fclose(fp1);
00487 
00489   for (m=3; m<=9; m+=3)
00490   {
00491     printf("\nTest of the inverse mpolar FFT for m=%d: \n",m);
00492     sprintf(filename,"mpolar_ifft_error%d.dat",m);
00493     fp1=fopen(filename,"w+");
00494     for (max_i=0; max_i<=20; max_i+=2)
00495     {
00497       inverse_mpolar_fft(f_direct,T,R,f_tilde,N,max_i,m);
00498 
00500       E_max=X(error_l_infty_complex)(f_hat,f_tilde,N*N);
00501       printf("%3d iterations: E_max = %e\n",max_i,E_max);
00502       fprintf(fp1,"%e\n",E_max);
00503     }
00504     fclose(fp1);
00505   }
00506 
00508   nfft_free(x);
00509   nfft_free(w);
00510   nfft_free(f_hat);
00511   nfft_free(f);
00512   nfft_free(f_direct);
00513   nfft_free(f_tilde);
00514 
00515   return 0;
00516 }
00517 /* \} */

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409