NFFT Logo 3.2.2
linogram_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: linogram_fft_test.c 3896 2012-10-10 12:19:26Z tovo $ */
00020 
00029 #include "config.h"
00030 
00031 #include <math.h>
00032 #include <stdlib.h>
00033 #ifdef HAVE_COMPLEX_H
00034 #include <complex.h>
00035 #endif
00036 
00037 #include "nfft3util.h"
00038 #include "nfft3.h"
00039 #include "infft.h"
00040 
00047 double GLOBAL_elapsed_time;
00048 
00052 static int linogram_grid(int T, int R, double *x, double *w)
00053 {
00054   int t, r;
00055   double W=(double)T*(((double)R/2.0)*((double)R/2.0)+1.0/4.0);
00056 
00057   for(t=-T/2; t<T/2; t++)
00058   {
00059     for(r=-R/2; r<R/2; r++)
00060     {
00061       if(t<0)
00062       {
00063         x[2*((t+T/2)*R+(r+R/2))+0] = (double)r/R;
00064         x[2*((t+T/2)*R+(r+R/2))+1] = (double)4*(t+T/4)/T*r/R;
00065       }
00066       else
00067       {
00068         x[2*((t+T/2)*R+(r+R/2))+0] = -(double)4*(t-T/4)/T*r/R;
00069         x[2*((t+T/2)*R+(r+R/2))+1] = (double)r/R;
00070       }
00071       if (r==0)
00072         w[(t+T/2)*R+(r+R/2)] = 1.0/4.0/W;
00073       else
00074         w[(t+T/2)*R+(r+R/2)] = fabs((double)r)/W;
00075     }
00076   }
00077 
00078   return T*R;                           
00079 }
00080 
00082 static int linogram_dft(fftw_complex *f_hat, int NN, fftw_complex *f, int T, int R, int m)
00083 {
00084   ticks t0, t1;
00085   int j,k;                              
00086   nfft_plan my_nfft_plan;               
00088   double *x, *w;                        
00090   int N[2],n[2];
00091   int M=T*R;                            
00093   N[0]=NN; n[0]=2*N[0];                 
00094   N[1]=NN; n[1]=2*N[1];                 
00096   x = (double *)nfft_malloc(2*T*R*(sizeof(double)));
00097   if (x==NULL)
00098     return -1;
00099 
00100   w = (double *)nfft_malloc(T*R*(sizeof(double)));
00101   if (w==NULL)
00102     return -1;
00103 
00105   nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
00106                   PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00107                   FFTW_MEASURE| FFTW_DESTROY_INPUT);
00108 
00110   linogram_grid(T,R,x,w);
00111   for(j=0;j<my_nfft_plan.M_total;j++)
00112   {
00113     my_nfft_plan.x[2*j+0] = x[2*j+0];
00114     my_nfft_plan.x[2*j+1] = x[2*j+1];
00115   }
00116 
00117 
00119   for(k=0;k<my_nfft_plan.N_total;k++)
00120     my_nfft_plan.f_hat[k] = f_hat[k];
00121 
00123   t0 = getticks();
00124   nfft_trafo_direct(&my_nfft_plan);
00125   t1 = getticks();
00126   GLOBAL_elapsed_time = nfft_elapsed_seconds(t1,t0);
00127 
00129   for(j=0;j<my_nfft_plan.M_total;j++)
00130     f[j] = my_nfft_plan.f[j];
00131 
00133   nfft_finalize(&my_nfft_plan);
00134   nfft_free(x);
00135   nfft_free(w);
00136 
00137   return EXIT_SUCCESS;
00138 }
00139 
00141 static int linogram_fft(fftw_complex *f_hat, int NN, fftw_complex *f, int T, int R, int m)
00142 {
00143   ticks t0, t1;
00144   int j,k;                              
00145   nfft_plan my_nfft_plan;               
00147   double *x, *w;                        
00149   int N[2],n[2];
00150   int M=T*R;                            
00152   N[0]=NN; n[0]=2*N[0];                 
00153   N[1]=NN; n[1]=2*N[1];                 
00155   x = (double *)nfft_malloc(2*T*R*(sizeof(double)));
00156   if (x==NULL)
00157     return -1;
00158 
00159   w = (double *)nfft_malloc(T*R*(sizeof(double)));
00160   if (w==NULL)
00161     return -1;
00162 
00164   nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
00165                   PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00166                   FFTW_MEASURE| FFTW_DESTROY_INPUT);
00167 
00169   linogram_grid(T,R,x,w);
00170   for(j=0;j<my_nfft_plan.M_total;j++)
00171   {
00172     my_nfft_plan.x[2*j+0] = x[2*j+0];
00173     my_nfft_plan.x[2*j+1] = x[2*j+1];
00174   }
00175 
00177   if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
00178     nfft_precompute_lin_psi(&my_nfft_plan);
00179 
00180   if(my_nfft_plan.nfft_flags & PRE_PSI)
00181     nfft_precompute_psi(&my_nfft_plan);
00182 
00183   if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
00184     nfft_precompute_full_psi(&my_nfft_plan);
00185 
00187   for(k=0;k<my_nfft_plan.N_total;k++)
00188     my_nfft_plan.f_hat[k] = f_hat[k];
00189 
00191   t0 = getticks();
00192   nfft_trafo(&my_nfft_plan);
00193   t1 = getticks();
00194   GLOBAL_elapsed_time = nfft_elapsed_seconds(t1,t0);
00195 
00197   for(j=0;j<my_nfft_plan.M_total;j++)
00198     f[j] = my_nfft_plan.f[j];
00199 
00201   nfft_finalize(&my_nfft_plan);
00202   nfft_free(x);
00203   nfft_free(w);
00204 
00205   return EXIT_SUCCESS;
00206 }
00207 
00209 static int inverse_linogram_fft(fftw_complex *f, int T, int R, fftw_complex *f_hat, int NN, int max_i, int m)
00210 {
00211   ticks t0, t1;
00212   int j,k;                              
00213   nfft_plan my_nfft_plan;               
00214   solver_plan_complex my_infft_plan;             
00216   double *x, *w;                        
00217   int l;                                
00219   int N[2],n[2];
00220   int M=T*R;                            
00222   N[0]=NN; n[0]=2*N[0];                 
00223   N[1]=NN; n[1]=2*N[1];                 
00225   x = (double *)nfft_malloc(2*T*R*(sizeof(double)));
00226   if (x==NULL)
00227     return -1;
00228 
00229   w = (double *)nfft_malloc(T*R*(sizeof(double)));
00230   if (w==NULL)
00231     return -1;
00232 
00234   nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
00235                   PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00236                   FFTW_MEASURE| FFTW_DESTROY_INPUT);
00237 
00239   solver_init_advanced_complex(&my_infft_plan,(nfft_mv_plan_complex*)(&my_nfft_plan), CGNR | PRECOMPUTE_WEIGHT );
00240 
00242   linogram_grid(T,R,x,w);
00243   for(j=0;j<my_nfft_plan.M_total;j++)
00244   {
00245     my_nfft_plan.x[2*j+0] = x[2*j+0];
00246     my_nfft_plan.x[2*j+1] = x[2*j+1];
00247     my_infft_plan.y[j]    = f[j];
00248     my_infft_plan.w[j]    = w[j];
00249   }
00250 
00252   if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
00253     nfft_precompute_lin_psi(&my_nfft_plan);
00254 
00255   if(my_nfft_plan.nfft_flags & PRE_PSI)
00256     nfft_precompute_psi(&my_nfft_plan);
00257 
00258   if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
00259     nfft_precompute_full_psi(&my_nfft_plan);
00260 
00262   if(my_infft_plan.flags & PRECOMPUTE_DAMP)
00263     for(j=0;j<my_nfft_plan.N[0];j++)
00264       for(k=0;k<my_nfft_plan.N[1];k++)
00265   {
00266     my_infft_plan.w_hat[j*my_nfft_plan.N[1]+k]=
00267         (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);
00268   }
00269 
00271   for(k=0;k<my_nfft_plan.N_total;k++)
00272     my_infft_plan.f_hat_iter[k] = 0.0 + _Complex_I*0.0;
00273 
00274   t0 = getticks();
00276   solver_before_loop_complex(&my_infft_plan);
00277 
00278   if (max_i<1)
00279   {
00280     l=1;
00281     for(k=0;k<my_nfft_plan.N_total;k++)
00282       my_infft_plan.f_hat_iter[k] = my_infft_plan.p_hat_iter[k];
00283   }
00284   else
00285   {
00286     for(l=1;l<=max_i;l++)
00287     {
00288       solver_loop_one_step_complex(&my_infft_plan);
00289     }
00290   }
00291 
00292   t1 = getticks();
00293   GLOBAL_elapsed_time = nfft_elapsed_seconds(t1,t0);
00294 
00296   for(k=0;k<my_nfft_plan.N_total;k++)
00297     f_hat[k] = my_infft_plan.f_hat_iter[k];
00298 
00300   solver_finalize_complex(&my_infft_plan);
00301   nfft_finalize(&my_nfft_plan);
00302   nfft_free(x);
00303   nfft_free(w);
00304 
00305   return EXIT_SUCCESS;
00306 }
00307 
00309 int comparison_fft(FILE *fp, int N, int T, int R)
00310 {
00311   ticks t0, t1;
00312   fftw_plan my_fftw_plan;
00313   fftw_complex *f_hat,*f;
00314   int m,k;
00315   double t_fft, t_dft_linogram;
00316 
00317   f_hat = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*N*N);
00318   f     = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*(T*R/4)*5);
00319 
00320   my_fftw_plan = fftw_plan_dft_2d(N,N,f_hat,f,FFTW_BACKWARD,FFTW_MEASURE);
00321 
00322   for(k=0; k<N*N; k++)
00323     f_hat[k] = (((double)rand())/RAND_MAX) + _Complex_I* (((double)rand())/RAND_MAX);
00324 
00325   t0 = getticks();
00326   for(m=0;m<65536/N;m++)
00327     {
00328       fftw_execute(my_fftw_plan);
00329       /* touch */
00330       f_hat[2]=2*f_hat[0];
00331     }
00332   t1 = getticks();
00333   GLOBAL_elapsed_time = nfft_elapsed_seconds(t1,t0);
00334   t_fft=N*GLOBAL_elapsed_time/65536;
00335 
00336   if(N<256)
00337     {
00338       linogram_dft(f_hat,N,f,T,R,m);
00339       t_dft_linogram=GLOBAL_elapsed_time;
00340     }
00341 
00342   for (m=3; m<=9; m+=3)
00343     {
00344       if((m==3)&&(N<256))
00345         fprintf(fp,"%d\t&\t&\t%1.1e&\t%1.1e&\t%d\t",N,t_fft,t_dft_linogram,m);
00346       else
00347         if(m==3)
00348     fprintf(fp,"%d\t&\t&\t%1.1e&\t       &\t%d\t",N,t_fft,m);
00349   else
00350     fprintf(fp,"  \t&\t&\t       &\t       &\t%d\t",m);
00351 
00352       printf("N=%d\tt_fft=%1.1e\tt_dft_linogram=%1.1e\tm=%d\t",N,t_fft,t_dft_linogram,m);
00353 
00354       linogram_fft(f_hat,N,f,T,R,m);
00355       fprintf(fp,"%1.1e&\t",GLOBAL_elapsed_time);
00356       printf("t_linogram=%1.1e\t",GLOBAL_elapsed_time);
00357       inverse_linogram_fft(f,T,R,f_hat,N,m+3,m);
00358       if(m==9)
00359   fprintf(fp,"%1.1e\\\\\\hline\n",GLOBAL_elapsed_time);
00360       else
00361   fprintf(fp,"%1.1e\\\\\n",GLOBAL_elapsed_time);
00362       printf("t_ilinogram=%1.1e\n",GLOBAL_elapsed_time);
00363     }
00364 
00365   fflush(fp);
00366 
00367   nfft_free(f);
00368   nfft_free(f_hat);
00369 
00370   return EXIT_SUCCESS;
00371 }
00372 
00374 int main(int argc,char **argv)
00375 {
00376   int N;                                
00377   int T, R;                             
00378   int M;                                
00379   double *x, *w;                        
00380   fftw_complex *f_hat, *f, *f_direct, *f_tilde;
00381   int k;
00382   int max_i;                            
00383   int m;
00384   double temp1, temp2, E_max=0.0;
00385   FILE *fp1, *fp2;
00386   char filename[30];
00387   int logN;
00388 
00389   if( argc!=4 )
00390   {
00391     printf("linogram_fft_test N T R \n");
00392     printf("\n");
00393     printf("N          linogram FFT of size NxN    \n");
00394     printf("T          number of slopes          \n");
00395     printf("R          number of offsets         \n");
00396 
00398     printf("\nHence, comparison FFTW, linogram FFT and inverse linogram FFT\n");
00399     fp1=fopen("linogram_comparison_fft.dat","w");
00400     if (fp1==NULL)
00401   return(-1);
00402     for (logN=4; logN<=8; logN++)
00403   comparison_fft(fp1,(1U<< logN), 3*(1U<< logN), 3*(1U<< (logN-1)));
00404     fclose(fp1);
00405 
00406     exit(-1);
00407   }
00408 
00409   N = atoi(argv[1]);
00410   T = atoi(argv[2]);
00411   R = atoi(argv[3]);
00412   printf("N=%d, linogram grid with T=%d, R=%d => ",N,T,R);
00413 
00414   x = (double *)nfft_malloc(5*T*R/2*(sizeof(double)));
00415   w = (double *)nfft_malloc(5*T*R/4*(sizeof(double)));
00416 
00417   f_hat    = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*N*N);
00418   f        = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*5*T*R/4);  /* 4/pi*log(1+sqrt(2)) = 1.122... < 1.25 */
00419   f_direct = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*5*T*R/4);
00420   f_tilde  = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*N*N);
00421 
00423   M=linogram_grid(T,R,x,w); printf("M=%d.\n",M);
00424 
00426   fp1=fopen("input_data_r.dat","r");
00427   fp2=fopen("input_data_i.dat","r");
00428   if ((fp1==NULL) || (fp2==NULL))
00429     return(-1);
00430   for(k=0;k<N*N;k++)
00431   {
00432     fscanf(fp1,"%le ",&temp1);
00433     fscanf(fp2,"%le ",&temp2);
00434     f_hat[k]=temp1+_Complex_I*temp2;
00435   }
00436   fclose(fp1);
00437   fclose(fp2);
00438 
00440       linogram_dft(f_hat,N,f_direct,T,R,1);
00441   //  linogram_fft(f_hat,N,f_direct,T,R,12);
00442 
00444   printf("\nTest of the linogram FFT: \n");
00445   fp1=fopen("linogram_fft_error.dat","w+");
00446   for (m=1; m<=12; m++)
00447   {
00449     linogram_fft(f_hat,N,f,T,R,m);
00450 
00452     E_max=X(error_l_infty_complex)(f_direct,f,M);
00453     //E_max=X(error_l_2_complex)(f_direct,f,M);
00454 
00455     printf("m=%2d: E_max = %e\n",m,E_max);
00456     fprintf(fp1,"%e\n",E_max);
00457   }
00458   fclose(fp1);
00459 
00461   for (m=3; m<=9; m+=3)
00462   {
00463     printf("\nTest of the inverse linogram FFT for m=%d: \n",m);
00464     sprintf(filename,"linogram_ifft_error%d.dat",m);
00465     fp1=fopen(filename,"w+");
00466     for (max_i=0; max_i<=20; max_i+=2)
00467     {
00469       inverse_linogram_fft(f_direct,T,R,f_tilde,N,max_i,m);
00470 
00472       E_max=X(error_l_infty_complex)(f_hat,f_tilde,N*N);
00473       printf("%3d iterations: E_max = %e\n",max_i,E_max);
00474       fprintf(fp1,"%e\n",E_max);
00475     }
00476     fclose(fp1);
00477   }
00478 
00480   nfft_free(x);
00481   nfft_free(w);
00482   nfft_free(f_hat);
00483   nfft_free(f);
00484   nfft_free(f_direct);
00485   nfft_free(f_tilde);
00486 
00487   return 0;
00488 }
00489 /* \} */

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409