NFFT  3.3.0
linogram_fft_test.c
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2002, 2015 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$ */
20 
30 #include <math.h>
31 #include <stdlib.h>
32 #include <complex.h>
33 
34 #define NFFT_PRECISION_DOUBLE
35 
36 #include "nfft3mp.h"
37 
44 NFFT_R GLOBAL_elapsed_time;
45 
49 static int linogram_grid(int T, int rr, NFFT_R *x, NFFT_R *w)
50 {
51  int t, r;
52  NFFT_R W = (NFFT_R) T * (((NFFT_R) rr / NFFT_K(2.0)) * ((NFFT_R) rr / NFFT_K(2.0)) + NFFT_K(1.0) / NFFT_K(4.0));
53 
54  for (t = -T / 2; t < T / 2; t++)
55  {
56  for (r = -rr / 2; r < rr / 2; r++)
57  {
58  if (t < 0)
59  {
60  x[2 * ((t + T / 2) * rr + (r + rr / 2)) + 0] = (NFFT_R) (r) / (NFFT_R)(rr);
61  x[2 * ((t + T / 2) * rr + (r + rr / 2)) + 1] = NFFT_K(4.0) * ((NFFT_R)(t) + (NFFT_R)(T) / NFFT_K(4.0)) / (NFFT_R)(T)
62  * (NFFT_R)(r) / (NFFT_R)(rr);
63  }
64  else
65  {
66  x[2 * ((t + T / 2) * rr + (r + rr / 2)) + 0] = -NFFT_K(4.0) * ((NFFT_R)(t) - (NFFT_R)(T) / NFFT_K(4.0)) / (NFFT_R)(T)
67  * (NFFT_R)(r) / (NFFT_R)(rr);
68  x[2 * ((t + T / 2) * rr + (r + rr / 2)) + 1] = (NFFT_R) r / (NFFT_R)(rr);
69  }
70  if (r == 0)
71  w[(t + T / 2) * rr + (r + rr / 2)] = NFFT_K(1.0) / NFFT_K(4.0) / W;
72  else
73  w[(t + T / 2) * rr + (r + rr / 2)] = NFFT_M(fabs)((NFFT_R) r) / W;
74  }
75  }
76 
77  return T * rr;
78 }
79 
81 static int linogram_dft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int rr, int m)
82 {
83  double t0, t1;
84  int j, k;
85  NFFT(plan) my_nfft_plan;
87  NFFT_R *x, *w;
89  int N[2], n[2];
90  int M = T * rr;
92  N[0] = NN;
93  n[0] = 2 * N[0];
94  N[1] = NN;
95  n[1] = 2 * N[1];
97  x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * rr) * (sizeof(NFFT_R)));
98  if (x == NULL)
99  return EXIT_FAILURE;
100 
101  w = (NFFT_R *) NFFT(malloc)((size_t)(T * rr) * (sizeof(NFFT_R)));
102  if (w == NULL)
103  return EXIT_FAILURE;
104 
106  NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, m,
107  PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | MALLOC_F | FFTW_INIT
108  | FFT_OUT_OF_PLACE,
109  FFTW_MEASURE | FFTW_DESTROY_INPUT);
110 
112  linogram_grid(T, rr, x, w);
113  for (j = 0; j < my_nfft_plan.M_total; j++)
114  {
115  my_nfft_plan.x[2 * j + 0] = x[2 * j + 0];
116  my_nfft_plan.x[2 * j + 1] = x[2 * j + 1];
117  }
118 
120  for (k = 0; k < my_nfft_plan.N_total; k++)
121  my_nfft_plan.f_hat[k] = f_hat[k];
122 
124  t0 = NFFT(clock_gettime_seconds)();
125  NFFT(trafo_direct)(&my_nfft_plan);
126  t1 = NFFT(clock_gettime_seconds)();
127  GLOBAL_elapsed_time = (t1 - t0);
128 
130  for (j = 0; j < my_nfft_plan.M_total; j++)
131  f[j] = my_nfft_plan.f[j];
132 
134  NFFT(finalize)(&my_nfft_plan);
135  NFFT(free)(x);
136  NFFT(free)(w);
137 
138  return EXIT_SUCCESS;
139 }
140 
142 static int linogram_fft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int rr, int m)
143 {
144  double t0, t1;
145  int j, k;
146  NFFT(plan) my_nfft_plan;
148  NFFT_R *x, *w;
150  int N[2], n[2];
151  int M = T * rr;
153  N[0] = NN;
154  n[0] = 2 * N[0];
155  N[1] = NN;
156  n[1] = 2 * N[1];
158  x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * rr) * (sizeof(NFFT_R)));
159  if (x == NULL)
160  return EXIT_FAILURE;
161 
162  w = (NFFT_R *) NFFT(malloc)((size_t)(T * rr) * (sizeof(NFFT_R)));
163  if (w == NULL)
164  return EXIT_FAILURE;
165 
167  NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, m,
168  PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | MALLOC_F | FFTW_INIT
169  | FFT_OUT_OF_PLACE,
170  FFTW_MEASURE | FFTW_DESTROY_INPUT);
171 
173  linogram_grid(T, rr, x, w);
174  for (j = 0; j < my_nfft_plan.M_total; j++)
175  {
176  my_nfft_plan.x[2 * j + 0] = x[2 * j + 0];
177  my_nfft_plan.x[2 * j + 1] = x[2 * j + 1];
178  }
179 
181  if (my_nfft_plan.flags & PRE_LIN_PSI)
182  NFFT(precompute_lin_psi)(&my_nfft_plan);
183 
184  if (my_nfft_plan.flags & PRE_PSI)
185  NFFT(precompute_psi)(&my_nfft_plan);
186 
187  if (my_nfft_plan.flags & PRE_FULL_PSI)
188  NFFT(precompute_full_psi)(&my_nfft_plan);
189 
191  for (k = 0; k < my_nfft_plan.N_total; k++)
192  my_nfft_plan.f_hat[k] = f_hat[k];
193 
195  t0 = NFFT(clock_gettime_seconds)();
196  NFFT(trafo)(&my_nfft_plan);
197  t1 = NFFT(clock_gettime_seconds)();
198  GLOBAL_elapsed_time = (t1 - t0);
199 
201  for (j = 0; j < my_nfft_plan.M_total; j++)
202  f[j] = my_nfft_plan.f[j];
203 
205  NFFT(finalize)(&my_nfft_plan);
206  NFFT(free)(x);
207  NFFT(free)(w);
208 
209  return EXIT_SUCCESS;
210 }
211 
213 static int inverse_linogram_fft(NFFT_C *f, int T, int rr, NFFT_C *f_hat, int NN,
214  int max_i, int m)
215 {
216  double t0, t1;
217  int j, k;
218  NFFT(plan) my_nfft_plan;
219  SOLVER(plan_complex) my_infft_plan;
221  NFFT_R *x, *w;
222  int l;
224  int N[2], n[2];
225  int M = T * rr;
227  N[0] = NN;
228  n[0] = 2 * N[0];
229  N[1] = NN;
230  n[1] = 2 * N[1];
232  x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * rr) * (sizeof(NFFT_R)));
233  if (x == NULL)
234  return EXIT_FAILURE;
235 
236  w = (NFFT_R *) NFFT(malloc)((size_t)(T * rr) * (sizeof(NFFT_R)));
237  if (w == NULL)
238  return EXIT_FAILURE;
239 
241  NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, m,
242  PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | MALLOC_F | FFTW_INIT
243  | FFT_OUT_OF_PLACE,
244  FFTW_MEASURE | FFTW_DESTROY_INPUT);
245 
247  SOLVER(init_advanced_complex)(&my_infft_plan,
248  (NFFT(mv_plan_complex)*) (&my_nfft_plan), CGNR | PRECOMPUTE_WEIGHT);
249 
251  linogram_grid(T, rr, x, w);
252  for (j = 0; j < my_nfft_plan.M_total; j++)
253  {
254  my_nfft_plan.x[2 * j + 0] = x[2 * j + 0];
255  my_nfft_plan.x[2 * j + 1] = x[2 * j + 1];
256  my_infft_plan.y[j] = f[j];
257  my_infft_plan.w[j] = w[j];
258  }
259 
261  if (my_nfft_plan.flags & PRE_LIN_PSI)
262  NFFT(precompute_lin_psi)(&my_nfft_plan);
263 
264  if (my_nfft_plan.flags & PRE_PSI)
265  NFFT(precompute_psi)(&my_nfft_plan);
266 
267  if (my_nfft_plan.flags & PRE_FULL_PSI)
268  NFFT(precompute_full_psi)(&my_nfft_plan);
269 
271  if (my_infft_plan.flags & PRECOMPUTE_DAMP)
272  for (j = 0; j < my_nfft_plan.N[0]; j++)
273  for (k = 0; k < my_nfft_plan.N[1]; k++)
274  {
275  my_infft_plan.w_hat[j * my_nfft_plan.N[1] + k] = (
276  NFFT_M(sqrt)(
277  NFFT_M(pow)((NFFT_R)(j - my_nfft_plan.N[0] / 2), NFFT_K(2.0))
278  + NFFT_M(pow)((NFFT_R)(k - my_nfft_plan.N[1] / 2), NFFT_K(2.0)))
279  > (NFFT_R)(my_nfft_plan.N[0] / 2) ? 0 : 1);
280  }
281 
283  for (k = 0; k < my_nfft_plan.N_total; k++)
284  my_infft_plan.f_hat_iter[k] = NFFT_K(0.0) + _Complex_I * NFFT_K(0.0);
285 
286  t0 = NFFT(clock_gettime_seconds)();
288  SOLVER(before_loop_complex)(&my_infft_plan);
289 
290  if (max_i < 1)
291  {
292  l = 1;
293  for (k = 0; k < my_nfft_plan.N_total; k++)
294  my_infft_plan.f_hat_iter[k] = my_infft_plan.p_hat_iter[k];
295  }
296  else
297  {
298  for (l = 1; l <= max_i; l++)
299  {
300  SOLVER(loop_one_step_complex)(&my_infft_plan);
301  }
302  }
303 
304  t1 = NFFT(clock_gettime_seconds)();
305  GLOBAL_elapsed_time = (t1 - t0);
306 
308  for (k = 0; k < my_nfft_plan.N_total; k++)
309  f_hat[k] = my_infft_plan.f_hat_iter[k];
310 
312  SOLVER(finalize_complex)(&my_infft_plan);
313  NFFT(finalize)(&my_nfft_plan);
314  NFFT(free)(x);
315  NFFT(free)(w);
316 
317  return EXIT_SUCCESS;
318 }
319 
321 static int comparison_fft(FILE *fp, int N, int T, int rr)
322 {
323  double t0, t1;
324  FFTW(plan) my_fftw_plan;
325  NFFT_C *f_hat, *f;
326  int m, k;
327  NFFT_R t_fft, t_dft_linogram;
328 
329  f_hat = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(N * N));
330  f = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)((T * rr / 4) * 5));
331 
332  my_fftw_plan = FFTW(plan_dft_2d)(N, N, f_hat, f, FFTW_BACKWARD, FFTW_MEASURE);
333 
334  for (k = 0; k < N * N; k++)
335  f_hat[k] = NFFT(drand48)() + _Complex_I * NFFT(drand48)();
336 
337  t0 = NFFT(clock_gettime_seconds)();
338  for (m = 0; m < 65536 / N; m++)
339  {
340  FFTW(execute)(my_fftw_plan);
341  /* touch */
342  f_hat[2] = NFFT_K(2.0) * f_hat[0];
343  }
344  t1 = NFFT(clock_gettime_seconds)();
345  GLOBAL_elapsed_time = (t1 - t0);
346  t_fft = (NFFT_R)(N) * GLOBAL_elapsed_time / NFFT_K(65536.0);
347 
348  if (N < 256)
349  {
350  linogram_dft(f_hat, N, f, T, rr, m);
351  t_dft_linogram = GLOBAL_elapsed_time;
352  }
353 
354  for (m = 3; m <= 9; m += 3)
355  {
356  if ((m == 3) && (N < 256))
357  fprintf(fp, "%d\t&\t&\t%1.1" NFFT__FES__ "&\t%1.1" NFFT__FES__ "&\t%d\t", N, t_fft, t_dft_linogram,
358  m);
359  else if (m == 3)
360  fprintf(fp, "%d\t&\t&\t%1.1" NFFT__FES__ "&\t &\t%d\t", N, t_fft, m);
361  else
362  fprintf(fp, " \t&\t&\t &\t &\t%d\t", m);
363 
364  printf("N=%d\tt_fft=%1.1" NFFT__FES__ "\tt_dft_linogram=%1.1" NFFT__FES__ "\tm=%d\t", N, t_fft,
365  t_dft_linogram, m);
366 
367  linogram_fft(f_hat, N, f, T, rr, m);
368  fprintf(fp, "%1.1" NFFT__FES__ "&\t", GLOBAL_elapsed_time);
369  printf("t_linogram=%1.1" NFFT__FES__ "\t", GLOBAL_elapsed_time);
370  inverse_linogram_fft(f, T, rr, f_hat, N, m + 3, m);
371  if (m == 9)
372  fprintf(fp, "%1.1" NFFT__FES__ "\\\\\\hline\n", GLOBAL_elapsed_time);
373  else
374  fprintf(fp, "%1.1" NFFT__FES__ "\\\\\n", GLOBAL_elapsed_time);
375  printf("t_ilinogram=%1.1" NFFT__FES__ "\n", GLOBAL_elapsed_time);
376  }
377 
378  fflush(fp);
379 
380  NFFT(free)(f);
381  NFFT(free)(f_hat);
382 
383  return EXIT_SUCCESS;
384 }
385 
387 int main(int argc, char **argv)
388 {
389  int N;
390  int T, rr;
391  int M;
392  NFFT_R *x, *w;
393  NFFT_C *f_hat, *f, *f_direct, *f_tilde;
394  int k;
395  int max_i;
396  int m;
397  NFFT_R temp1, temp2, E_max = NFFT_K(0.0);
398  FILE *fp1, *fp2;
399  char filename[30];
400  int logN;
401 
402  if (argc != 4)
403  {
404  printf("linogram_fft_test N T R \n");
405  printf("\n");
406  printf("N linogram FFT of size NxN \n");
407  printf("T number of slopes \n");
408  printf("R number of offsets \n");
409 
411  printf("\nHence, comparison FFTW, linogram FFT and inverse linogram FFT\n");
412  fp1 = fopen("linogram_comparison_fft.dat", "w");
413  if (fp1 == NULL)
414  return (-1);
415  for (logN = 4; logN <= 8; logN++)
416  comparison_fft(fp1, (int)(1U << logN), 3 * (int)(1U << logN),
417  3 * (int)(1U << (logN - 1)));
418  fclose(fp1);
419 
420  return EXIT_FAILURE;
421  }
422 
423  N = atoi(argv[1]);
424  T = atoi(argv[2]);
425  rr = atoi(argv[3]);
426  printf("N=%d, linogram grid with T=%d, R=%d => ", N, T, rr);
427 
428  x = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * rr / 2) * (sizeof(NFFT_R)));
429  w = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * rr / 4) * (sizeof(NFFT_R)));
430 
431  f_hat = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(N * N));
432  f = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(5 * T * rr / 4)); /* 4/pi*log(1+sqrt(2)) = 1.122... < 1.25 */
433  f_direct = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(5 * T * rr / 4));
434  f_tilde = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(N * N));
435 
437  M = linogram_grid(T, rr, x, w);
438  printf("M=%d.\n", M);
439 
441  fp1 = fopen("input_data_r.dat", "r");
442  fp2 = fopen("input_data_i.dat", "r");
443  if ((fp1 == NULL) || (fp2 == NULL))
444  return EXIT_FAILURE;
445  for (k = 0; k < N * N; k++)
446  {
447  fscanf(fp1, NFFT__FR__ " ", &temp1);
448  fscanf(fp2, NFFT__FR__ " ", &temp2);
449  f_hat[k] = temp1 + _Complex_I * temp2;
450  }
451  fclose(fp1);
452  fclose(fp2);
453 
455  linogram_dft(f_hat, N, f_direct, T, rr, 1);
456  // linogram_fft(f_hat,N,f_direct,T,R,12);
457 
459  printf("\nTest of the linogram FFT: \n");
460  fp1 = fopen("linogram_fft_error.dat", "w+");
461  for (m = 1; m <= 12; m++)
462  {
464  linogram_fft(f_hat, N, f, T, rr, m);
465 
467  E_max = NFFT(error_l_infty_complex)(f_direct, f, M);
468  //E_max=NFFT(error_l_2_complex)(f_direct,f,M);
469 
470  printf("m=%2d: E_max = %" NFFT__FES__ "\n", m, E_max);
471  fprintf(fp1, "%" NFFT__FES__ "\n", E_max);
472  }
473  fclose(fp1);
474 
476  for (m = 3; m <= 9; m += 3)
477  {
478  printf("\nTest of the inverse linogram FFT for m=%d: \n", m);
479  sprintf(filename, "linogram_ifft_error%d.dat", m);
480  fp1 = fopen(filename, "w+");
481  for (max_i = 0; max_i <= 20; max_i += 2)
482  {
484  inverse_linogram_fft(f_direct, T, rr, f_tilde, N, max_i, m);
485 
487  E_max = NFFT(error_l_infty_complex)(f_hat, f_tilde, N * N);
488  printf("%3d iterations: E_max = %" NFFT__FES__ "\n", max_i, E_max);
489  fprintf(fp1, "%" NFFT__FES__ "\n", E_max);
490  }
491  fclose(fp1);
492  }
493 
495  NFFT(free)(x);
496  NFFT(free)(w);
497  NFFT(free)(f_hat);
498  NFFT(free)(f);
499  NFFT(free)(f_direct);
500  NFFT(free)(f_tilde);
501 
502  return EXIT_SUCCESS;
503 }
504 /* \} */
static int max_i(int a, int b)
max
Definition: fastsum.c:48
static int inverse_linogram_fft(NFFT_C *f, int T, int rr, NFFT_C *f_hat, int NN, int max_i, int m)
NFFT-based inverse pseudo-polar FFT.
static int linogram_dft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int rr, int m)
discrete pseudo-polar FFT
static int comparison_fft(FILE *fp, int N, int T, int rr)
Comparison of the FFTW, linogram FFT, and inverse linogram FFT.
static int linogram_fft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int rr, int m)
NFFT-based pseudo-polar FFT.
static int linogram_grid(int T, int rr, NFFT_R *x, NFFT_R *w)
Generates the points x with weights w for the linogram grid with T slopes and R offsets.