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