NFFT Logo 3.2.3
mpolar_fft_test.c
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2002, 2012 Jens Keiner, Stefan Kunis, Daniel Potts
3  *
4  * This program is free software; you can redistribute it and/or modify it under
5  * the terms of the GNU General Public License as published by the Free Software
6  * Foundation; either version 2 of the License, or (at your option) any later
7  * version.
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12  * details.
13  *
14  * You should have received a copy of the GNU General Public License along with
15  * this program; if not, write to the Free Software Foundation, Inc., 51
16  * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  */
18 
19 /* $Id: mpolar_fft_test.c 3896 2012-10-10 12:19:26Z tovo $ */
20 
30 #include "config.h"
31 
32 #include <math.h>
33 #include <stdlib.h>
34 #ifdef HAVE_COMPLEX_H
35 #include <complex.h>
36 #endif
37 
38 #include "nfft3util.h"
39 #include "nfft3.h"
40 #include "infft.h"
41 
48 double GLOBAL_elapsed_time;
49 
63 static int mpolar_grid(int T, int R, double *x, double *w)
64 {
65  int t, r;
66  double W;
67  int R2=2*ceil(sqrt(2.0)*R/2);
68  double xx, yy;
69  int M=0;
70 
71  for(t=-T/2; t<T/2; t++)
72  {
73  for(r=-R2/2; r<R2/2; r++)
74  {
75  xx = (double)r/R*cos(PI*t/T);
76  yy = (double)r/R*sin(PI*t/T);
77 
78  if ( ((-0.5-1.0/(double)R)<=xx) & (xx<=(0.5+1.0/(double)R)) &
79  ((-0.5-1.0/(double)R)<=yy) & (yy<=(0.5+1.0/(double)R)) )
80  {
81  x[2*M+0] = xx;
82  x[2*M+1] = yy;
83 
84  if (r==0)
85  w[M] = 1.0/4.0;
86  else
87  w[M] = fabs((double)r);
88 
89  M++;
90  }
91  }
92  }
93 
95  W=0.0;
96  for (t=0; t<M; t++)
97  W+=w[t];
98 
99  for (t=0; t<M; t++)
100  w[t]/=W;
101 
102  return M;
103 }
104 
106 static int mpolar_dft(fftw_complex *f_hat, int NN, fftw_complex *f, int T, int R, int m)
107 {
108  ticks t0, t1;
109  int j,k;
110  nfft_plan my_nfft_plan;
112  double *x, *w;
114  int N[2],n[2];
115  int M;
117  N[0]=NN; n[0]=2*N[0];
118  N[1]=NN; n[1]=2*N[1];
120  x = (double *)nfft_malloc(5*(T/2)*R*(sizeof(double)));
121  if (x==NULL)
122  return -1;
123 
124  w = (double *)nfft_malloc(5*(T*R)/4*(sizeof(double)));
125  if (w==NULL)
126  return -1;
127 
129  M=mpolar_grid(T,R,x,w);
130  nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
131  PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
132  FFTW_MEASURE| FFTW_DESTROY_INPUT);
133 
135  for(j=0;j<my_nfft_plan.M_total;j++)
136  {
137  my_nfft_plan.x[2*j+0] = x[2*j+0];
138  my_nfft_plan.x[2*j+1] = x[2*j+1];
139  }
140 
142  for(k=0;k<my_nfft_plan.N_total;k++)
143  my_nfft_plan.f_hat[k] = f_hat[k];
144 
145  t0 = getticks();
146 
148  nfft_trafo_direct(&my_nfft_plan);
149 
150  t1 = getticks();
151  GLOBAL_elapsed_time = nfft_elapsed_seconds(t1,t0);
152 
154  for(j=0;j<my_nfft_plan.M_total;j++)
155  f[j] = my_nfft_plan.f[j];
156 
158  nfft_finalize(&my_nfft_plan);
159  nfft_free(x);
160  nfft_free(w);
161 
162  return EXIT_SUCCESS;
163 }
164 
166 static int mpolar_fft(fftw_complex *f_hat, int NN, fftw_complex *f, int T, int R, int m)
167 {
168  ticks t0, t1;
169  int j,k;
170  nfft_plan my_nfft_plan;
172  double *x, *w;
174  int N[2],n[2];
175  int M;
177  N[0]=NN; n[0]=2*N[0];
178  N[1]=NN; n[1]=2*N[1];
180  x = (double *)nfft_malloc(5*T*R/2*(sizeof(double)));
181  if (x==NULL)
182  return -1;
183 
184  w = (double *)nfft_malloc(5*T*R/4*(sizeof(double)));
185  if (w==NULL)
186  return -1;
187 
189  M=mpolar_grid(T,R,x,w);
190  nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
191  PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
192  FFTW_MEASURE| FFTW_DESTROY_INPUT);
193 
196  for(j=0;j<my_nfft_plan.M_total;j++)
197  {
198  my_nfft_plan.x[2*j+0] = x[2*j+0];
199  my_nfft_plan.x[2*j+1] = x[2*j+1];
200  }
201 
203  if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
204  nfft_precompute_lin_psi(&my_nfft_plan);
205 
206  if(my_nfft_plan.nfft_flags & PRE_PSI)
207  nfft_precompute_psi(&my_nfft_plan);
208 
209  if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
210  nfft_precompute_full_psi(&my_nfft_plan);
211 
213  for(k=0;k<my_nfft_plan.N_total;k++)
214  my_nfft_plan.f_hat[k] = f_hat[k];
215 
216  t0 = getticks();
217 
219  nfft_trafo(&my_nfft_plan);
220 
221  t1 = getticks();
222  GLOBAL_elapsed_time = nfft_elapsed_seconds(t1,t0);
223 
225  for(j=0;j<my_nfft_plan.M_total;j++)
226  f[j] = my_nfft_plan.f[j];
227 
229  nfft_finalize(&my_nfft_plan);
230  nfft_free(x);
231  nfft_free(w);
232 
233  return EXIT_SUCCESS;
234 }
235 
237 static int inverse_mpolar_fft(fftw_complex *f, int T, int R, fftw_complex *f_hat, int NN, int max_i, int m)
238 {
239  ticks t0, t1;
240  int j,k;
241  nfft_plan my_nfft_plan;
242  solver_plan_complex my_infft_plan;
244  double *x, *w;
245  int l;
247  int N[2],n[2];
248  int M;
250  N[0]=NN; n[0]=2*N[0];
251  N[1]=NN; n[1]=2*N[1];
253  x = (double *)nfft_malloc(5*T*R/2*(sizeof(double)));
254  if (x==NULL)
255  return -1;
256 
257  w = (double *)nfft_malloc(5*T*R/4*(sizeof(double)));
258  if (w==NULL)
259  return -1;
260 
262  M=mpolar_grid(T,R,x,w);
263  nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
264  PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
265  FFTW_MEASURE| FFTW_DESTROY_INPUT);
266 
268  solver_init_advanced_complex(&my_infft_plan,(nfft_mv_plan_complex*)(&my_nfft_plan), CGNR | PRECOMPUTE_WEIGHT );
269 
271  for(j=0;j<my_nfft_plan.M_total;j++)
272  {
273  my_nfft_plan.x[2*j+0] = x[2*j+0];
274  my_nfft_plan.x[2*j+1] = x[2*j+1];
275  my_infft_plan.y[j] = f[j];
276  my_infft_plan.w[j] = w[j];
277  }
278 
280  if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
281  nfft_precompute_lin_psi(&my_nfft_plan);
282 
283  if(my_nfft_plan.nfft_flags & PRE_PSI)
284  nfft_precompute_psi(&my_nfft_plan);
285 
286  if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
287  nfft_precompute_full_psi(&my_nfft_plan);
288 
289 
291  if(my_infft_plan.flags & PRECOMPUTE_DAMP)
292  for(j=0;j<my_nfft_plan.N[0];j++)
293  for(k=0;k<my_nfft_plan.N[1];k++)
294  {
295  my_infft_plan.w_hat[j*my_nfft_plan.N[1]+k]=
296  (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);
297  }
298 
300  for(k=0;k<my_nfft_plan.N_total;k++)
301  my_infft_plan.f_hat_iter[k] = 0.0 + _Complex_I*0.0;
302 
303  t0 = getticks();
304 
306  solver_before_loop_complex(&my_infft_plan);
307 
308  if (max_i<1)
309  {
310  l=1;
311  for(k=0;k<my_nfft_plan.N_total;k++)
312  my_infft_plan.f_hat_iter[k] = my_infft_plan.p_hat_iter[k];
313  }
314  else
315  {
316  for(l=1;l<=max_i;l++)
317  {
318  solver_loop_one_step_complex(&my_infft_plan);
319  }
320  }
321 
322  t1 = getticks();
323  GLOBAL_elapsed_time = nfft_elapsed_seconds(t1,t0);
324 
326  for(k=0;k<my_nfft_plan.N_total;k++)
327  f_hat[k] = my_infft_plan.f_hat_iter[k];
328 
330  solver_finalize_complex(&my_infft_plan);
331  nfft_finalize(&my_nfft_plan);
332  nfft_free(x);
333  nfft_free(w);
334 
335  return EXIT_SUCCESS;
336 }
337 
339 static int comparison_fft(FILE *fp, int N, int T, int R)
340 {
341  ticks t0, t1;
342  fftw_plan my_fftw_plan;
343  fftw_complex *f_hat,*f;
344  int m,k;
345  double t_fft, t_dft_mpolar;
346 
347  f_hat = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*N*N);
348  f = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*(T*R/4)*5);
349 
350  my_fftw_plan = fftw_plan_dft_2d(N,N,f_hat,f,FFTW_BACKWARD,FFTW_MEASURE);
351 
352  for(k=0; k<N*N; k++)
353  f_hat[k] = (((double)rand())/RAND_MAX) + _Complex_I* (((double)rand())/RAND_MAX);
354 
355  t0 = getticks();
356  for(m=0;m<65536/N;m++)
357  {
358  fftw_execute(my_fftw_plan);
359  /* touch */
360  f_hat[2]=2*f_hat[0];
361  }
362  t1 = getticks();
363  GLOBAL_elapsed_time = nfft_elapsed_seconds(t1,t0);
364  t_fft=N*GLOBAL_elapsed_time/65536;
365 
366  if(N<256)
367  {
368  mpolar_dft(f_hat,N,f,T,R,1);
369  t_dft_mpolar=GLOBAL_elapsed_time;
370  }
371 
372  for (m=3; m<=9; m+=3)
373  {
374  if((m==3)&&(N<256))
375  fprintf(fp,"%d\t&\t&\t%1.1e&\t%1.1e&\t%d\t",N,t_fft,t_dft_mpolar,m);
376  else
377  if(m==3)
378  fprintf(fp,"%d\t&\t&\t%1.1e&\t &\t%d\t",N,t_fft,m);
379  else
380  fprintf(fp," \t&\t&\t &\t &\t%d\t",m);
381 
382  printf("N=%d\tt_fft=%1.1e\tt_dft_mpolar=%1.1e\tm=%d\t",N,t_fft,t_dft_mpolar,m);
383 
384  mpolar_fft(f_hat,N,f,T,R,m);
385  fprintf(fp,"%1.1e&\t",GLOBAL_elapsed_time);
386  printf("t_mpolar=%1.1e\t",GLOBAL_elapsed_time);
387  inverse_mpolar_fft(f,T,R,f_hat,N,2*m,m);
388  if(m==9)
389  fprintf(fp,"%1.1e\\\\\\hline\n",GLOBAL_elapsed_time);
390  else
391  fprintf(fp,"%1.1e\\\\\n",GLOBAL_elapsed_time);
392  printf("t_impolar=%1.1e\n",GLOBAL_elapsed_time);
393  }
394 
395  fflush(fp);
396 
397  nfft_free(f);
398  nfft_free(f_hat);
399 
400  return EXIT_SUCCESS;
401 }
402 
404 int main(int argc,char **argv)
405 {
406  int N;
407  int T, R;
408  int M;
409  double *x, *w;
410  fftw_complex *f_hat, *f, *f_direct, *f_tilde;
411  int k;
412  int max_i;
413  int m;
414  double temp1, temp2, E_max=0.0;
415  FILE *fp1, *fp2;
416  char filename[30];
417  int logN;
418 
419  if( argc!=4 )
420  {
421  printf("mpolar_fft_test N T R \n");
422  printf("\n");
423  printf("N mpolar FFT of size NxN \n");
424  printf("T number of slopes \n");
425  printf("R number of offsets \n");
426 
428  printf("\nHence, comparison FFTW, mpolar FFT and inverse mpolar FFT\n");
429  fp1=fopen("mpolar_comparison_fft.dat","w");
430  if (fp1==NULL)
431  return(-1);
432  for (logN=4; logN<=8; logN++)
433  comparison_fft(fp1,(1U<< logN), 3*(1U<< logN), 3*(1U<< (logN-1)));
434  fclose(fp1);
435 
436  exit(-1);
437  }
438 
439  N = atoi(argv[1]);
440  T = atoi(argv[2]);
441  R = atoi(argv[3]);
442  printf("N=%d, modified polar grid with T=%d, R=%d => ",N,T,R);
443 
444  x = (double *)nfft_malloc(5*T*R/2*(sizeof(double)));
445  w = (double *)nfft_malloc(5*T*R/4*(sizeof(double)));
446 
447  f_hat = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*N*N);
448  f = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*1.25*T*R); /* 4/pi*log(1+sqrt(2)) = 1.122... < 1.25 */
449  f_direct = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*1.25*T*R);
450  f_tilde = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*N*N);
451 
453  M=mpolar_grid(T,R,x,w); printf("M=%d.\n",M);
454 
456  fp1=fopen("input_data_r.dat","r");
457  fp2=fopen("input_data_i.dat","r");
458  if ((fp1==NULL) || (fp2==NULL))
459  return(-1);
460  for(k=0;k<N*N;k++)
461  {
462  fscanf(fp1,"%le ",&temp1);
463  fscanf(fp2,"%le ",&temp2);
464  f_hat[k]=temp1+ _Complex_I*temp2;
465  }
466  fclose(fp1);
467  fclose(fp2);
468 
470  mpolar_dft(f_hat,N,f_direct,T,R,1);
471  // mpolar_fft(f_hat,N,f_direct,T,R,12);
472 
474  printf("\nTest of the mpolar FFT: \n");
475  fp1=fopen("mpolar_fft_error.dat","w+");
476  for (m=1; m<=12; m++)
477  {
479  mpolar_fft(f_hat,N,f,T,R,m);
480 
482  E_max=X(error_l_infty_complex)(f_direct,f,M);
483  printf("m=%2d: E_max = %e\n",m,E_max);
484  fprintf(fp1,"%e\n",E_max);
485  }
486  fclose(fp1);
487 
489  for (m=3; m<=9; m+=3)
490  {
491  printf("\nTest of the inverse mpolar FFT for m=%d: \n",m);
492  sprintf(filename,"mpolar_ifft_error%d.dat",m);
493  fp1=fopen(filename,"w+");
494  for (max_i=0; max_i<=20; max_i+=2)
495  {
497  inverse_mpolar_fft(f_direct,T,R,f_tilde,N,max_i,m);
498 
500  E_max=X(error_l_infty_complex)(f_hat,f_tilde,N*N);
501  printf("%3d iterations: E_max = %e\n",max_i,E_max);
502  fprintf(fp1,"%e\n",E_max);
503  }
504  fclose(fp1);
505  }
506 
508  nfft_free(x);
509  nfft_free(w);
510  nfft_free(f_hat);
511  nfft_free(f);
512  nfft_free(f_direct);
513  nfft_free(f_tilde);
514 
515  return 0;
516 }
517 /* \} */
Generated on Tue Apr 30 2013 by Doxygen 1.8.1