NFFT  3.4.1
mpolar_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 
29 #include <math.h>
30 #include <stdlib.h>
31 #include <complex.h>
32 
33 #define NFFT_PRECISION_DOUBLE
34 
35 #include "nfft3mp.h"
36 
43 NFFT_R GLOBAL_elapsed_time;
44 
58 static int mpolar_grid(int T, int S, NFFT_R *x, NFFT_R *w)
59 {
60  int t, r;
61  NFFT_R W;
62  int R2 = 2 * (int)(NFFT_M(lrint)(NFFT_M(ceil)(NFFT_M(sqrt)(NFFT_K(2.0)) * (NFFT_R)(S) / NFFT_K(2.0))));
63  NFFT_R xx, yy;
64  int M = 0;
65 
66  for (t = -T / 2; t < T / 2; t++)
67  {
68  for (r = -R2 / 2; r < R2 / 2; r++)
69  {
70  xx = (NFFT_R) (r) / (NFFT_R)(S) * NFFT_M(cos)(NFFT_KPI * (NFFT_R)(t) / (NFFT_R)(T));
71  yy = (NFFT_R) (r) / (NFFT_R)(S) * NFFT_M(sin)(NFFT_KPI * (NFFT_R)(t) / (NFFT_R)(T));
72 
73  if (((-NFFT_K(0.5) - NFFT_K(1.0) / (NFFT_R) S) <= xx) & (xx <= (NFFT_K(0.5) + NFFT_K(1.0) / (NFFT_R) S))
74  & ((-NFFT_K(0.5) - NFFT_K(1.0) / (NFFT_R) S) <= yy)
75  & (yy <= (NFFT_K(0.5) + NFFT_K(1.0) / (NFFT_R) S)))
76  {
77  x[2 * M + 0] = xx;
78  x[2 * M + 1] = yy;
79 
80  if (r == 0)
81  w[M] = NFFT_K(1.0) / NFFT_K(4.0);
82  else
83  w[M] = NFFT_M(fabs)((NFFT_R) r);
84 
85  M++;
86  }
87  }
88  }
89 
91  W = NFFT_K(0.0);
92  for (t = 0; t < M; t++)
93  W += w[t];
94 
95  for (t = 0; t < M; t++)
96  w[t] /= W;
97 
98  return M;
99 }
100 
102 static int mpolar_dft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int S, int m)
103 {
104  double t0, t1;
105  int j, k;
106  NFFT(plan) my_nfft_plan;
108  NFFT_R *x, *w;
110  int N[2], n[2];
111  int M;
113  N[0] = NN;
114  n[0] = 2 * N[0];
115  N[1] = NN;
116  n[1] = 2 * N[1];
118  x = (NFFT_R *) NFFT(malloc)((size_t)(5 * (T / 2) * S) * (sizeof(NFFT_R)));
119  if (x == NULL)
120  return EXIT_FAILURE;
121 
122  w = (NFFT_R *) NFFT(malloc)((size_t)(5 * (T * S) / 4) * (sizeof(NFFT_R)));
123  if (w == NULL)
124  return EXIT_FAILURE;
125 
127  M = mpolar_grid(T, S, x, w);
128  NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, m,
129  PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | MALLOC_F | FFTW_INIT
130  | FFT_OUT_OF_PLACE,
131  FFTW_MEASURE | FFTW_DESTROY_INPUT);
132 
134  for (j = 0; j < my_nfft_plan.M_total; j++)
135  {
136  my_nfft_plan.x[2 * j + 0] = x[2 * j + 0];
137  my_nfft_plan.x[2 * j + 1] = x[2 * j + 1];
138  }
139 
141  for (k = 0; k < my_nfft_plan.N_total; k++)
142  my_nfft_plan.f_hat[k] = f_hat[k];
143 
144  t0 = NFFT(clock_gettime_seconds)();
145 
147  NFFT(trafo_direct)(&my_nfft_plan);
148 
149  t1 = NFFT(clock_gettime_seconds)();
150  GLOBAL_elapsed_time = (t1 - t0);
151 
153  for (j = 0; j < my_nfft_plan.M_total; j++)
154  f[j] = my_nfft_plan.f[j];
155 
157  NFFT(finalize)(&my_nfft_plan);
158  NFFT(free)(x);
159  NFFT(free)(w);
160 
161  return EXIT_SUCCESS;
162 }
163 
165 static int mpolar_fft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int S, int m)
166 {
167  double t0, t1;
168  int j, k;
169  NFFT(plan) my_nfft_plan;
171  NFFT_R *x, *w;
173  int N[2], n[2];
174  int M;
176  N[0] = NN;
177  n[0] = 2 * N[0];
178  N[1] = NN;
179  n[1] = 2 * N[1];
181  x = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * S / 2) * (sizeof(NFFT_R)));
182  if (x == NULL)
183  return EXIT_FAILURE;
184 
185  w = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * S / 4) * (sizeof(NFFT_R)));
186  if (w == NULL)
187  return EXIT_FAILURE;
188 
190  M = mpolar_grid(T, S, x, w);
191  NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, m,
192  PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | MALLOC_F | FFTW_INIT
193  | FFT_OUT_OF_PLACE,
194  FFTW_MEASURE | FFTW_DESTROY_INPUT);
195 
198  for (j = 0; j < my_nfft_plan.M_total; j++)
199  {
200  my_nfft_plan.x[2 * j + 0] = x[2 * j + 0];
201  my_nfft_plan.x[2 * j + 1] = x[2 * j + 1];
202  }
203 
205  if (my_nfft_plan.flags & PRE_LIN_PSI)
206  NFFT(precompute_lin_psi)(&my_nfft_plan);
207 
208  if (my_nfft_plan.flags & PRE_PSI)
209  NFFT(precompute_psi)(&my_nfft_plan);
210 
211  if (my_nfft_plan.flags & PRE_FULL_PSI)
212  NFFT(precompute_full_psi)(&my_nfft_plan);
213 
215  for (k = 0; k < my_nfft_plan.N_total; k++)
216  my_nfft_plan.f_hat[k] = f_hat[k];
217 
218  t0 = NFFT(clock_gettime_seconds)();
219 
221  NFFT(trafo)(&my_nfft_plan);
222 
223  t1 = NFFT(clock_gettime_seconds)();
224  GLOBAL_elapsed_time = (t1 - t0);
225 
227  for (j = 0; j < my_nfft_plan.M_total; j++)
228  f[j] = my_nfft_plan.f[j];
229 
231  NFFT(finalize)(&my_nfft_plan);
232  NFFT(free)(x);
233  NFFT(free)(w);
234 
235  return EXIT_SUCCESS;
236 }
237 
239 static int inverse_mpolar_fft(NFFT_C *f, int T, int S, NFFT_C *f_hat, int NN, int max_i,
240  int m)
241 {
242  double t0, t1;
243  int j, k;
244  NFFT(plan) my_nfft_plan;
245  SOLVER(plan_complex) my_infft_plan;
247  NFFT_R *x, *w;
248  int l;
250  int N[2], n[2];
251  int M;
253  N[0] = NN;
254  n[0] = 2 * N[0];
255  N[1] = NN;
256  n[1] = 2 * N[1];
258  x = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * S / 2) * (sizeof(NFFT_R)));
259  if (x == NULL)
260  return EXIT_FAILURE;
261 
262  w = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * S / 4) * (sizeof(NFFT_R)));
263  if (w == NULL)
264  return EXIT_FAILURE;
265 
267  M = mpolar_grid(T, S, x, w);
268  NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, m,
269  PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | MALLOC_F | FFTW_INIT
270  | FFT_OUT_OF_PLACE,
271  FFTW_MEASURE | FFTW_DESTROY_INPUT);
272 
274  SOLVER(init_advanced_complex)(&my_infft_plan,
275  (NFFT(mv_plan_complex)*) (&my_nfft_plan), CGNR | PRECOMPUTE_WEIGHT);
276 
278  for (j = 0; j < my_nfft_plan.M_total; j++)
279  {
280  my_nfft_plan.x[2 * j + 0] = x[2 * j + 0];
281  my_nfft_plan.x[2 * j + 1] = x[2 * j + 1];
282  my_infft_plan.y[j] = f[j];
283  my_infft_plan.w[j] = w[j];
284  }
285 
287  if (my_nfft_plan.flags & PRE_LIN_PSI)
288  NFFT(precompute_lin_psi)(&my_nfft_plan);
289 
290  if (my_nfft_plan.flags & PRE_PSI)
291  NFFT(precompute_psi)(&my_nfft_plan);
292 
293  if (my_nfft_plan.flags & PRE_FULL_PSI)
294  NFFT(precompute_full_psi)(&my_nfft_plan);
295 
297  if (my_infft_plan.flags & PRECOMPUTE_DAMP)
298  for (j = 0; j < my_nfft_plan.N[0]; j++)
299  for (k = 0; k < my_nfft_plan.N[1]; k++)
300  {
301  my_infft_plan.w_hat[j * my_nfft_plan.N[1] + k] = (
302  NFFT_M(sqrt)(
303  NFFT_M(pow)((NFFT_R)(j - my_nfft_plan.N[0] / 2), NFFT_K(2.0))
304  + NFFT_M(pow)((NFFT_R)(k - my_nfft_plan.N[1] / 2), NFFT_K(2.0)))
305  > (NFFT_R)(my_nfft_plan.N[0] / 2) ? NFFT_K(0.0) : NFFT_K(1.0));
306  }
307 
309  for (k = 0; k < my_nfft_plan.N_total; k++)
310  my_infft_plan.f_hat_iter[k] = NFFT_K(0.0) + _Complex_I * NFFT_K(0.0);
311 
312  t0 = NFFT(clock_gettime_seconds)();
313 
315  SOLVER(before_loop_complex)(&my_infft_plan);
316 
317  if (max_i < 1)
318  {
319  l = 1;
320  for (k = 0; k < my_nfft_plan.N_total; k++)
321  my_infft_plan.f_hat_iter[k] = my_infft_plan.p_hat_iter[k];
322  }
323  else
324  {
325  for (l = 1; l <= max_i; l++)
326  {
327  SOLVER(loop_one_step_complex)(&my_infft_plan);
328  }
329  }
330 
331  t1 = NFFT(clock_gettime_seconds)();
332  GLOBAL_elapsed_time = (t1 - t0);
333 
335  for (k = 0; k < my_nfft_plan.N_total; k++)
336  f_hat[k] = my_infft_plan.f_hat_iter[k];
337 
339  SOLVER(finalize_complex)(&my_infft_plan);
340  NFFT(finalize)(&my_nfft_plan);
341  NFFT(free)(x);
342  NFFT(free)(w);
343 
344  return EXIT_SUCCESS;
345 }
346 
348 static int comparison_fft(FILE *fp, int N, int T, int S)
349 {
350  double t0, t1;
351  FFTW(plan) my_fftw_plan;
352  NFFT_C *f_hat, *f;
353  int m, k;
354  NFFT_R t_fft, t_dft_mpolar;
355 
356  f_hat = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(N * N));
357  f = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)((T * S / 4) * 5));
358 
359  my_fftw_plan = FFTW(plan_dft_2d)(N, N, f_hat, f, FFTW_BACKWARD, FFTW_MEASURE);
360 
361  for (k = 0; k < N * N; k++)
362  f_hat[k] = NFFT(drand48)() + _Complex_I * NFFT(drand48)();
363 
364  t0 = NFFT(clock_gettime_seconds)();
365  for (m = 0; m < 65536 / N; m++)
366  {
367  FFTW(execute)(my_fftw_plan);
368  /* touch */
369  f_hat[2] = NFFT_K(2.0) * f_hat[0];
370  }
371  t1 = NFFT(clock_gettime_seconds)();
372  GLOBAL_elapsed_time = (t1 - t0);
373  t_fft = (NFFT_R)(N) * GLOBAL_elapsed_time / NFFT_K(65536.0);
374 
375  if (N < 256)
376  {
377  mpolar_dft(f_hat, N, f, T, S, 1);
378  t_dft_mpolar = GLOBAL_elapsed_time;
379  }
380 
381  for (m = 3; m <= 9; m += 3)
382  {
383  if ((m == 3) && (N < 256))
384  fprintf(fp, "%d\t&\t&\t%1.1" NFFT__FES__ "&\t%1.1" NFFT__FES__ "&\t%d\t", N, t_fft, t_dft_mpolar, m);
385  else if (m == 3)
386  fprintf(fp, "%d\t&\t&\t%1.1" NFFT__FES__ "&\t &\t%d\t", N, t_fft, m);
387  else
388  fprintf(fp, " \t&\t&\t &\t &\t%d\t", m);
389 
390  printf("N=%d\tt_fft=%1.1" NFFT__FES__ "\tt_dft_mpolar=%1.1" NFFT__FES__ "\tm=%d\t", N, t_fft,
391  t_dft_mpolar, m);
392 
393  mpolar_fft(f_hat, N, f, T, S, m);
394  fprintf(fp, "%1.1" NFFT__FES__ "&\t", GLOBAL_elapsed_time);
395  printf("t_mpolar=%1.1" NFFT__FES__ "\t", GLOBAL_elapsed_time);
396  inverse_mpolar_fft(f, T, S, f_hat, N, 2 * m, m);
397  if (m == 9)
398  fprintf(fp, "%1.1" NFFT__FES__ "\\\\\\hline\n", GLOBAL_elapsed_time);
399  else
400  fprintf(fp, "%1.1" NFFT__FES__ "\\\\\n", GLOBAL_elapsed_time);
401  printf("t_impolar=%1.1" NFFT__FES__ "\n", GLOBAL_elapsed_time);
402  }
403 
404  fflush(fp);
405 
406  NFFT(free)(f);
407  NFFT(free)(f_hat);
408 
409  return EXIT_SUCCESS;
410 }
411 
413 int main(int argc, char **argv)
414 {
415  int N;
416  int T, S;
417  int M;
418  NFFT_R *x, *w;
419  NFFT_C *f_hat, *f, *f_direct, *f_tilde;
420  int k;
421  int max_i;
422  int m;
423  NFFT_R temp1, temp2, E_max = NFFT_K(0.0);
424  FILE *fp1, *fp2;
425  char filename[30];
426  int logN;
427 
428  if (argc != 4)
429  {
430  printf("mpolar_fft_test N T R \n");
431  printf("\n");
432  printf("N mpolar FFT of size NxN \n");
433  printf("T number of slopes \n");
434  printf("R number of offsets \n");
435 
437  printf("\nHence, comparison FFTW, mpolar FFT and inverse mpolar FFT\n");
438  fp1 = fopen("mpolar_comparison_fft.dat", "w");
439  if (fp1 == NULL)
440  return (-1);
441  for (logN = 4; logN <= 8; logN++)
442  comparison_fft(fp1, (int)(1U << logN), 3 * (int)(1U << logN),
443  3 * (int)(1U << (logN - 1)));
444  fclose(fp1);
445 
446  exit(EXIT_FAILURE);
447  }
448 
449  N = atoi(argv[1]);
450  T = atoi(argv[2]);
451  S = atoi(argv[3]);
452  printf("N=%d, modified polar grid with T=%d, R=%d => ", N, T, S);
453 
454  x = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * S / 2) * (sizeof(NFFT_R)));
455  w = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * S / 4) * (sizeof(NFFT_R)));
456 
457  f_hat = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(N * N));
458  f = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(1.25 * T * S)); /* 4/pi*log(1+sqrt(2)) = 1.122... < 1.25 */
459  f_direct = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(1.25 * T * S));
460  f_tilde = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(N * N));
461 
463  M = mpolar_grid(T, S, x, w);
464  printf("M=%d.\n", M);
465 
467  fp1 = fopen("input_data_r.dat", "r");
468  fp2 = fopen("input_data_i.dat", "r");
469  if ((fp1 == NULL) || (fp2 == NULL))
470  return (-1);
471  for (k = 0; k < N * N; k++)
472  {
473  fscanf(fp1, NFFT__FR__ " ", &temp1);
474  fscanf(fp2, NFFT__FR__ " ", &temp2);
475  f_hat[k] = temp1 + _Complex_I * temp2;
476  }
477  fclose(fp1);
478  fclose(fp2);
479 
481  mpolar_dft(f_hat, N, f_direct, T, S, 1);
482  // mpolar_fft(f_hat,N,f_direct,T,R,12);
483 
485  printf("\nTest of the mpolar FFT: \n");
486  fp1 = fopen("mpolar_fft_error.dat", "w+");
487  for (m = 1; m <= 12; m++)
488  {
490  mpolar_fft(f_hat, N, f, T, S, m);
491 
493  E_max = NFFT(error_l_infty_complex)(f_direct, f, M);
494  printf("m=%2d: E_max = %" NFFT__FES__ "\n", m, E_max);
495  fprintf(fp1, "%" NFFT__FES__ "\n", E_max);
496  }
497  fclose(fp1);
498 
500  for (m = 3; m <= 9; m += 3)
501  {
502  printf("\nTest of the inverse mpolar FFT for m=%d: \n", m);
503  sprintf(filename, "mpolar_ifft_error%d.dat", m);
504  fp1 = fopen(filename, "w+");
505  for (max_i = 0; max_i <= 20; max_i += 2)
506  {
508  inverse_mpolar_fft(f_direct, T, S, f_tilde, N, max_i, m);
509 
511  E_max = NFFT(error_l_infty_complex)(f_hat, f_tilde, N * N);
512  printf("%3d iterations: E_max = %" NFFT__FES__ "\n", max_i, E_max);
513  fprintf(fp1, "%" NFFT__FES__ "\n", E_max);
514  }
515  fclose(fp1);
516  }
517 
519  NFFT(free)(x);
520  NFFT(free)(w);
521  NFFT(free)(f_hat);
522  NFFT(free)(f);
523  NFFT(free)(f_direct);
524  NFFT(free)(f_tilde);
525 
526  return EXIT_SUCCESS;
527 }
528 /* \} */
max_i
static int max_i(int a, int b)
max
Definition: fastsum.c:46
PRECOMPUTE_DAMP
#define PRECOMPUTE_DAMP
Definition: nfft3.h:792
MALLOC_X
#define MALLOC_X
Definition: nfft3.h:199
MALLOC_F_HAT
#define MALLOC_F_HAT
Definition: nfft3.h:200
PRECOMPUTE_WEIGHT
#define PRECOMPUTE_WEIGHT
Definition: nfft3.h:791
inverse_mpolar_fft
static int inverse_mpolar_fft(NFFT_C *f, int T, int S, NFFT_C *f_hat, int NN, int max_i, int m)
inverse NFFT-based mpolar FFT
Definition: mpolar_fft_test.c:239
mpolar_fft
static int mpolar_fft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int S, int m)
NFFT-based mpolar FFT.
Definition: mpolar_fft_test.c:165
mpolar_dft
static int mpolar_dft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int S, int m)
discrete mpolar FFT
Definition: mpolar_fft_test.c:102
FFTW_INIT
#define FFTW_INIT
Definition: nfft3.h:203
main
int main(int argc, char **argv)
test program for various parameters
Definition: mpolar_fft_test.c:413
mpolar_grid
static int mpolar_grid(int T, int S, NFFT_R *x, NFFT_R *w)
Generates the points with weights for the modified polar grid with angles and offsets...
Definition: mpolar_fft_test.c:58
MALLOC_F
#define MALLOC_F
Definition: nfft3.h:201
FFT_OUT_OF_PLACE
#define FFT_OUT_OF_PLACE
Definition: nfft3.h:202
PRE_LIN_PSI
#define PRE_LIN_PSI
Definition: nfft3.h:195
PRE_PSI
#define PRE_PSI
Definition: nfft3.h:197
CGNR
#define CGNR
Definition: nfft3.h:788
comparison_fft
static int comparison_fft(FILE *fp, int N, int T, int S)
Comparison of the FFTW, mpolar FFT, and inverse mpolar FFT.
Definition: mpolar_fft_test.c:348
PRE_FULL_PSI
#define PRE_FULL_PSI
Definition: nfft3.h:198
PRE_PHI_HUT
#define PRE_PHI_HUT
Definition: nfft3.h:193
Generated on Thu Apr 12 2018 21:00:49 for NFFT by   doxygen 1.8.6