00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
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
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);
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
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
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