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