NFFT Logo 3.2.3
linogram_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: linogram_fft_test.c 3896 2012-10-10 12:19:26Z tovo $ */
20 
29 #include "config.h"
30 
31 #include <math.h>
32 #include <stdlib.h>
33 #ifdef HAVE_COMPLEX_H
34 #include <complex.h>
35 #endif
36 
37 #include "nfft3util.h"
38 #include "nfft3.h"
39 #include "infft.h"
40 
47 double GLOBAL_elapsed_time;
48 
52 static int linogram_grid(int T, int R, double *x, double *w)
53 {
54  int t, r;
55  double W=(double)T*(((double)R/2.0)*((double)R/2.0)+1.0/4.0);
56 
57  for(t=-T/2; t<T/2; t++)
58  {
59  for(r=-R/2; r<R/2; r++)
60  {
61  if(t<0)
62  {
63  x[2*((t+T/2)*R+(r+R/2))+0] = (double)r/R;
64  x[2*((t+T/2)*R+(r+R/2))+1] = (double)4*(t+T/4)/T*r/R;
65  }
66  else
67  {
68  x[2*((t+T/2)*R+(r+R/2))+0] = -(double)4*(t-T/4)/T*r/R;
69  x[2*((t+T/2)*R+(r+R/2))+1] = (double)r/R;
70  }
71  if (r==0)
72  w[(t+T/2)*R+(r+R/2)] = 1.0/4.0/W;
73  else
74  w[(t+T/2)*R+(r+R/2)] = fabs((double)r)/W;
75  }
76  }
77 
78  return T*R;
79 }
80 
82 static int linogram_dft(fftw_complex *f_hat, int NN, fftw_complex *f, int T, int R, int m)
83 {
84  ticks t0, t1;
85  int j,k;
86  nfft_plan my_nfft_plan;
88  double *x, *w;
90  int N[2],n[2];
91  int M=T*R;
93  N[0]=NN; n[0]=2*N[0];
94  N[1]=NN; n[1]=2*N[1];
96  x = (double *)nfft_malloc(2*T*R*(sizeof(double)));
97  if (x==NULL)
98  return -1;
99 
100  w = (double *)nfft_malloc(T*R*(sizeof(double)));
101  if (w==NULL)
102  return -1;
103 
105  nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
106  PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
107  FFTW_MEASURE| FFTW_DESTROY_INPUT);
108 
110  linogram_grid(T,R,x,w);
111  for(j=0;j<my_nfft_plan.M_total;j++)
112  {
113  my_nfft_plan.x[2*j+0] = x[2*j+0];
114  my_nfft_plan.x[2*j+1] = x[2*j+1];
115  }
116 
117 
119  for(k=0;k<my_nfft_plan.N_total;k++)
120  my_nfft_plan.f_hat[k] = f_hat[k];
121 
123  t0 = getticks();
124  nfft_trafo_direct(&my_nfft_plan);
125  t1 = getticks();
126  GLOBAL_elapsed_time = nfft_elapsed_seconds(t1,t0);
127 
129  for(j=0;j<my_nfft_plan.M_total;j++)
130  f[j] = my_nfft_plan.f[j];
131 
133  nfft_finalize(&my_nfft_plan);
134  nfft_free(x);
135  nfft_free(w);
136 
137  return EXIT_SUCCESS;
138 }
139 
141 static int linogram_fft(fftw_complex *f_hat, int NN, fftw_complex *f, int T, int R, int m)
142 {
143  ticks t0, t1;
144  int j,k;
145  nfft_plan my_nfft_plan;
147  double *x, *w;
149  int N[2],n[2];
150  int M=T*R;
152  N[0]=NN; n[0]=2*N[0];
153  N[1]=NN; n[1]=2*N[1];
155  x = (double *)nfft_malloc(2*T*R*(sizeof(double)));
156  if (x==NULL)
157  return -1;
158 
159  w = (double *)nfft_malloc(T*R*(sizeof(double)));
160  if (w==NULL)
161  return -1;
162 
164  nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
165  PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
166  FFTW_MEASURE| FFTW_DESTROY_INPUT);
167 
169  linogram_grid(T,R,x,w);
170  for(j=0;j<my_nfft_plan.M_total;j++)
171  {
172  my_nfft_plan.x[2*j+0] = x[2*j+0];
173  my_nfft_plan.x[2*j+1] = x[2*j+1];
174  }
175 
177  if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
178  nfft_precompute_lin_psi(&my_nfft_plan);
179 
180  if(my_nfft_plan.nfft_flags & PRE_PSI)
181  nfft_precompute_psi(&my_nfft_plan);
182 
183  if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
184  nfft_precompute_full_psi(&my_nfft_plan);
185 
187  for(k=0;k<my_nfft_plan.N_total;k++)
188  my_nfft_plan.f_hat[k] = f_hat[k];
189 
191  t0 = getticks();
192  nfft_trafo(&my_nfft_plan);
193  t1 = getticks();
194  GLOBAL_elapsed_time = nfft_elapsed_seconds(t1,t0);
195 
197  for(j=0;j<my_nfft_plan.M_total;j++)
198  f[j] = my_nfft_plan.f[j];
199 
201  nfft_finalize(&my_nfft_plan);
202  nfft_free(x);
203  nfft_free(w);
204 
205  return EXIT_SUCCESS;
206 }
207 
209 static int inverse_linogram_fft(fftw_complex *f, int T, int R, fftw_complex *f_hat, int NN, int max_i, int m)
210 {
211  ticks t0, t1;
212  int j,k;
213  nfft_plan my_nfft_plan;
214  solver_plan_complex my_infft_plan;
216  double *x, *w;
217  int l;
219  int N[2],n[2];
220  int M=T*R;
222  N[0]=NN; n[0]=2*N[0];
223  N[1]=NN; n[1]=2*N[1];
225  x = (double *)nfft_malloc(2*T*R*(sizeof(double)));
226  if (x==NULL)
227  return -1;
228 
229  w = (double *)nfft_malloc(T*R*(sizeof(double)));
230  if (w==NULL)
231  return -1;
232 
234  nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
235  PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
236  FFTW_MEASURE| FFTW_DESTROY_INPUT);
237 
239  solver_init_advanced_complex(&my_infft_plan,(nfft_mv_plan_complex*)(&my_nfft_plan), CGNR | PRECOMPUTE_WEIGHT );
240 
242  linogram_grid(T,R,x,w);
243  for(j=0;j<my_nfft_plan.M_total;j++)
244  {
245  my_nfft_plan.x[2*j+0] = x[2*j+0];
246  my_nfft_plan.x[2*j+1] = x[2*j+1];
247  my_infft_plan.y[j] = f[j];
248  my_infft_plan.w[j] = w[j];
249  }
250 
252  if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
253  nfft_precompute_lin_psi(&my_nfft_plan);
254 
255  if(my_nfft_plan.nfft_flags & PRE_PSI)
256  nfft_precompute_psi(&my_nfft_plan);
257 
258  if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
259  nfft_precompute_full_psi(&my_nfft_plan);
260 
262  if(my_infft_plan.flags & PRECOMPUTE_DAMP)
263  for(j=0;j<my_nfft_plan.N[0];j++)
264  for(k=0;k<my_nfft_plan.N[1];k++)
265  {
266  my_infft_plan.w_hat[j*my_nfft_plan.N[1]+k]=
267  (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);
268  }
269 
271  for(k=0;k<my_nfft_plan.N_total;k++)
272  my_infft_plan.f_hat_iter[k] = 0.0 + _Complex_I*0.0;
273 
274  t0 = getticks();
276  solver_before_loop_complex(&my_infft_plan);
277 
278  if (max_i<1)
279  {
280  l=1;
281  for(k=0;k<my_nfft_plan.N_total;k++)
282  my_infft_plan.f_hat_iter[k] = my_infft_plan.p_hat_iter[k];
283  }
284  else
285  {
286  for(l=1;l<=max_i;l++)
287  {
288  solver_loop_one_step_complex(&my_infft_plan);
289  }
290  }
291 
292  t1 = getticks();
293  GLOBAL_elapsed_time = nfft_elapsed_seconds(t1,t0);
294 
296  for(k=0;k<my_nfft_plan.N_total;k++)
297  f_hat[k] = my_infft_plan.f_hat_iter[k];
298 
300  solver_finalize_complex(&my_infft_plan);
301  nfft_finalize(&my_nfft_plan);
302  nfft_free(x);
303  nfft_free(w);
304 
305  return EXIT_SUCCESS;
306 }
307 
309 int comparison_fft(FILE *fp, int N, int T, int R)
310 {
311  ticks t0, t1;
312  fftw_plan my_fftw_plan;
313  fftw_complex *f_hat,*f;
314  int m,k;
315  double t_fft, t_dft_linogram;
316 
317  f_hat = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*N*N);
318  f = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*(T*R/4)*5);
319 
320  my_fftw_plan = fftw_plan_dft_2d(N,N,f_hat,f,FFTW_BACKWARD,FFTW_MEASURE);
321 
322  for(k=0; k<N*N; k++)
323  f_hat[k] = (((double)rand())/RAND_MAX) + _Complex_I* (((double)rand())/RAND_MAX);
324 
325  t0 = getticks();
326  for(m=0;m<65536/N;m++)
327  {
328  fftw_execute(my_fftw_plan);
329  /* touch */
330  f_hat[2]=2*f_hat[0];
331  }
332  t1 = getticks();
333  GLOBAL_elapsed_time = nfft_elapsed_seconds(t1,t0);
334  t_fft=N*GLOBAL_elapsed_time/65536;
335 
336  if(N<256)
337  {
338  linogram_dft(f_hat,N,f,T,R,m);
339  t_dft_linogram=GLOBAL_elapsed_time;
340  }
341 
342  for (m=3; m<=9; m+=3)
343  {
344  if((m==3)&&(N<256))
345  fprintf(fp,"%d\t&\t&\t%1.1e&\t%1.1e&\t%d\t",N,t_fft,t_dft_linogram,m);
346  else
347  if(m==3)
348  fprintf(fp,"%d\t&\t&\t%1.1e&\t &\t%d\t",N,t_fft,m);
349  else
350  fprintf(fp," \t&\t&\t &\t &\t%d\t",m);
351 
352  printf("N=%d\tt_fft=%1.1e\tt_dft_linogram=%1.1e\tm=%d\t",N,t_fft,t_dft_linogram,m);
353 
354  linogram_fft(f_hat,N,f,T,R,m);
355  fprintf(fp,"%1.1e&\t",GLOBAL_elapsed_time);
356  printf("t_linogram=%1.1e\t",GLOBAL_elapsed_time);
357  inverse_linogram_fft(f,T,R,f_hat,N,m+3,m);
358  if(m==9)
359  fprintf(fp,"%1.1e\\\\\\hline\n",GLOBAL_elapsed_time);
360  else
361  fprintf(fp,"%1.1e\\\\\n",GLOBAL_elapsed_time);
362  printf("t_ilinogram=%1.1e\n",GLOBAL_elapsed_time);
363  }
364 
365  fflush(fp);
366 
367  nfft_free(f);
368  nfft_free(f_hat);
369 
370  return EXIT_SUCCESS;
371 }
372 
374 int main(int argc,char **argv)
375 {
376  int N;
377  int T, R;
378  int M;
379  double *x, *w;
380  fftw_complex *f_hat, *f, *f_direct, *f_tilde;
381  int k;
382  int max_i;
383  int m;
384  double temp1, temp2, E_max=0.0;
385  FILE *fp1, *fp2;
386  char filename[30];
387  int logN;
388 
389  if( argc!=4 )
390  {
391  printf("linogram_fft_test N T R \n");
392  printf("\n");
393  printf("N linogram FFT of size NxN \n");
394  printf("T number of slopes \n");
395  printf("R number of offsets \n");
396 
398  printf("\nHence, comparison FFTW, linogram FFT and inverse linogram FFT\n");
399  fp1=fopen("linogram_comparison_fft.dat","w");
400  if (fp1==NULL)
401  return(-1);
402  for (logN=4; logN<=8; logN++)
403  comparison_fft(fp1,(1U<< logN), 3*(1U<< logN), 3*(1U<< (logN-1)));
404  fclose(fp1);
405 
406  exit(-1);
407  }
408 
409  N = atoi(argv[1]);
410  T = atoi(argv[2]);
411  R = atoi(argv[3]);
412  printf("N=%d, linogram grid with T=%d, R=%d => ",N,T,R);
413 
414  x = (double *)nfft_malloc(5*T*R/2*(sizeof(double)));
415  w = (double *)nfft_malloc(5*T*R/4*(sizeof(double)));
416 
417  f_hat = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*N*N);
418  f = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*5*T*R/4); /* 4/pi*log(1+sqrt(2)) = 1.122... < 1.25 */
419  f_direct = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*5*T*R/4);
420  f_tilde = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*N*N);
421 
423  M=linogram_grid(T,R,x,w); printf("M=%d.\n",M);
424 
426  fp1=fopen("input_data_r.dat","r");
427  fp2=fopen("input_data_i.dat","r");
428  if ((fp1==NULL) || (fp2==NULL))
429  return(-1);
430  for(k=0;k<N*N;k++)
431  {
432  fscanf(fp1,"%le ",&temp1);
433  fscanf(fp2,"%le ",&temp2);
434  f_hat[k]=temp1+_Complex_I*temp2;
435  }
436  fclose(fp1);
437  fclose(fp2);
438 
440  linogram_dft(f_hat,N,f_direct,T,R,1);
441  // linogram_fft(f_hat,N,f_direct,T,R,12);
442 
444  printf("\nTest of the linogram FFT: \n");
445  fp1=fopen("linogram_fft_error.dat","w+");
446  for (m=1; m<=12; m++)
447  {
449  linogram_fft(f_hat,N,f,T,R,m);
450 
452  E_max=X(error_l_infty_complex)(f_direct,f,M);
453  //E_max=X(error_l_2_complex)(f_direct,f,M);
454 
455  printf("m=%2d: E_max = %e\n",m,E_max);
456  fprintf(fp1,"%e\n",E_max);
457  }
458  fclose(fp1);
459 
461  for (m=3; m<=9; m+=3)
462  {
463  printf("\nTest of the inverse linogram FFT for m=%d: \n",m);
464  sprintf(filename,"linogram_ifft_error%d.dat",m);
465  fp1=fopen(filename,"w+");
466  for (max_i=0; max_i<=20; max_i+=2)
467  {
469  inverse_linogram_fft(f_direct,T,R,f_tilde,N,max_i,m);
470 
472  E_max=X(error_l_infty_complex)(f_hat,f_tilde,N*N);
473  printf("%3d iterations: E_max = %e\n",max_i,E_max);
474  fprintf(fp1,"%e\n",E_max);
475  }
476  fclose(fp1);
477  }
478 
480  nfft_free(x);
481  nfft_free(w);
482  nfft_free(f_hat);
483  nfft_free(f);
484  nfft_free(f_direct);
485  nfft_free(f_tilde);
486 
487  return 0;
488 }
489 /* \} */
Generated on Tue Apr 30 2013 by Doxygen 1.8.1