NFFT  3.5.3
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,
130  FFTW_MEASURE);
131 
133  for (j = 0; j < my_nfft_plan.M_total; j++)
134  {
135  my_nfft_plan.x[2 * j + 0] = x[2 * j + 0];
136  my_nfft_plan.x[2 * j + 1] = x[2 * j + 1];
137  }
138 
140  for (k = 0; k < my_nfft_plan.N_total; k++)
141  my_nfft_plan.f_hat[k] = f_hat[k];
142 
143  t0 = NFFT(clock_gettime_seconds)();
144 
146  NFFT(trafo_direct)(&my_nfft_plan);
147 
148  t1 = NFFT(clock_gettime_seconds)();
149  GLOBAL_elapsed_time = (t1 - t0);
150 
152  for (j = 0; j < my_nfft_plan.M_total; j++)
153  f[j] = my_nfft_plan.f[j];
154 
156  NFFT(finalize)(&my_nfft_plan);
157  NFFT(free)(x);
158  NFFT(free)(w);
159 
160  return EXIT_SUCCESS;
161 }
162 
164 static int mpolar_fft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int S, int m)
165 {
166  double t0, t1;
167  int j, k;
168  NFFT(plan) my_nfft_plan;
170  NFFT_R *x, *w;
172  int N[2], n[2];
173  int M;
175  N[0] = NN;
176  n[0] = 2 * N[0];
177  N[1] = NN;
178  n[1] = 2 * N[1];
180  x = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * S / 2) * (sizeof(NFFT_R)));
181  if (x == NULL)
182  return EXIT_FAILURE;
183 
184  w = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * S / 4) * (sizeof(NFFT_R)));
185  if (w == NULL)
186  return EXIT_FAILURE;
187 
189  M = mpolar_grid(T, S, x, w);
190  NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, m,
193  FFTW_MEASURE | FFTW_DESTROY_INPUT);
194 
197  for (j = 0; j < my_nfft_plan.M_total; j++)
198  {
199  my_nfft_plan.x[2 * j + 0] = x[2 * j + 0];
200  my_nfft_plan.x[2 * j + 1] = x[2 * j + 1];
201  }
202 
204  if (my_nfft_plan.flags & PRE_LIN_PSI)
205  NFFT(precompute_lin_psi)(&my_nfft_plan);
206 
207  if (my_nfft_plan.flags & PRE_PSI)
208  NFFT(precompute_psi)(&my_nfft_plan);
209 
210  if (my_nfft_plan.flags & PRE_FULL_PSI)
211  NFFT(precompute_full_psi)(&my_nfft_plan);
212 
214  for (k = 0; k < my_nfft_plan.N_total; k++)
215  my_nfft_plan.f_hat[k] = f_hat[k];
216 
217  t0 = NFFT(clock_gettime_seconds)();
218 
220  NFFT(trafo)(&my_nfft_plan);
221 
222  t1 = NFFT(clock_gettime_seconds)();
223  GLOBAL_elapsed_time = (t1 - t0);
224 
226  for (j = 0; j < my_nfft_plan.M_total; j++)
227  f[j] = my_nfft_plan.f[j];
228 
230  NFFT(finalize)(&my_nfft_plan);
231  NFFT(free)(x);
232  NFFT(free)(w);
233 
234  return EXIT_SUCCESS;
235 }
236 
238 static int inverse_mpolar_fft(NFFT_C *f, int T, int S, NFFT_C *f_hat, int NN, int max_i,
239  int m)
240 {
241  double t0, t1;
242  int j, k;
243  NFFT(plan) my_nfft_plan;
244  SOLVER(plan_complex) my_infft_plan;
246  NFFT_R *x, *w;
247  int l;
249  int N[2], n[2];
250  int M;
252  N[0] = NN;
253  n[0] = 2 * N[0];
254  N[1] = NN;
255  n[1] = 2 * N[1];
257  x = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * S / 2) * (sizeof(NFFT_R)));
258  if (x == NULL)
259  return EXIT_FAILURE;
260 
261  w = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * S / 4) * (sizeof(NFFT_R)));
262  if (w == NULL)
263  return EXIT_FAILURE;
264 
266  M = mpolar_grid(T, S, x, w);
267  NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, m,
270  FFTW_MEASURE | FFTW_DESTROY_INPUT);
271 
273  SOLVER(init_advanced_complex)(&my_infft_plan,
274  (NFFT(mv_plan_complex)*) (&my_nfft_plan), CGNR | PRECOMPUTE_WEIGHT);
275 
277  for (j = 0; j < my_nfft_plan.M_total; j++)
278  {
279  my_nfft_plan.x[2 * j + 0] = x[2 * j + 0];
280  my_nfft_plan.x[2 * j + 1] = x[2 * j + 1];
281  my_infft_plan.y[j] = f[j];
282  my_infft_plan.w[j] = w[j];
283  }
284 
286  if (my_nfft_plan.flags & PRE_LIN_PSI)
287  NFFT(precompute_lin_psi)(&my_nfft_plan);
288 
289  if (my_nfft_plan.flags & PRE_PSI)
290  NFFT(precompute_psi)(&my_nfft_plan);
291 
292  if (my_nfft_plan.flags & PRE_FULL_PSI)
293  NFFT(precompute_full_psi)(&my_nfft_plan);
294 
296  if (my_infft_plan.flags & PRECOMPUTE_DAMP)
297  for (j = 0; j < my_nfft_plan.N[0]; j++)
298  for (k = 0; k < my_nfft_plan.N[1]; k++)
299  {
300  my_infft_plan.w_hat[j * my_nfft_plan.N[1] + k] = (
301  NFFT_M(sqrt)(
302  NFFT_M(pow)((NFFT_R)(j - my_nfft_plan.N[0] / 2), NFFT_K(2.0))
303  + NFFT_M(pow)((NFFT_R)(k - my_nfft_plan.N[1] / 2), NFFT_K(2.0)))
304  > (NFFT_R)(my_nfft_plan.N[0] / 2) ? NFFT_K(0.0) : NFFT_K(1.0));
305  }
306 
308  for (k = 0; k < my_nfft_plan.N_total; k++)
309  my_infft_plan.f_hat_iter[k] = NFFT_K(0.0) + _Complex_I * NFFT_K(0.0);
310 
311  t0 = NFFT(clock_gettime_seconds)();
312 
314  SOLVER(before_loop_complex)(&my_infft_plan);
315 
316  if (max_i < 1)
317  {
318  l = 1;
319  for (k = 0; k < my_nfft_plan.N_total; k++)
320  my_infft_plan.f_hat_iter[k] = my_infft_plan.p_hat_iter[k];
321  }
322  else
323  {
324  for (l = 1; l <= max_i; l++)
325  {
326  SOLVER(loop_one_step_complex)(&my_infft_plan);
327  }
328  }
329 
330  t1 = NFFT(clock_gettime_seconds)();
331  GLOBAL_elapsed_time = (t1 - t0);
332 
334  for (k = 0; k < my_nfft_plan.N_total; k++)
335  f_hat[k] = my_infft_plan.f_hat_iter[k];
336 
338  SOLVER(finalize_complex)(&my_infft_plan);
339  NFFT(finalize)(&my_nfft_plan);
340  NFFT(free)(x);
341  NFFT(free)(w);
342 
343  return EXIT_SUCCESS;
344 }
345 
347 static int comparison_fft(FILE *fp, int N, int T, int S)
348 {
349  double t0, t1;
350  FFTW(plan) my_fftw_plan;
351  NFFT_C *f_hat, *f;
352  int m, k;
353  NFFT_R t_fft, t_dft_mpolar;
354 
355  f_hat = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(N * N));
356  f = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)((T * S / 4) * 5));
357 
358  my_fftw_plan = FFTW(plan_dft_2d)(N, N, f_hat, f, FFTW_BACKWARD, FFTW_MEASURE);
359 
360  for (k = 0; k < N * N; k++)
361  f_hat[k] = NFFT(drand48)() + _Complex_I * NFFT(drand48)();
362 
363  t0 = NFFT(clock_gettime_seconds)();
364  for (m = 0; m < 65536 / N; m++)
365  {
366  FFTW(execute)(my_fftw_plan);
367  /* touch */
368  f_hat[2] = NFFT_K(2.0) * f_hat[0];
369  }
370  t1 = NFFT(clock_gettime_seconds)();
371  GLOBAL_elapsed_time = (t1 - t0);
372  t_fft = (NFFT_R)(N) * GLOBAL_elapsed_time / NFFT_K(65536.0);
373 
374  if (N < 256)
375  {
376  mpolar_dft(f_hat, N, f, T, S, 1);
377  t_dft_mpolar = GLOBAL_elapsed_time;
378  }
379 
380  for (m = 3; m <= 9; m += 3)
381  {
382  if ((m == 3) && (N < 256))
383  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);
384  else if (m == 3)
385  fprintf(fp, "%d\t&\t&\t%1.1" NFFT__FES__ "&\t &\t%d\t", N, t_fft, m);
386  else
387  fprintf(fp, " \t&\t&\t &\t &\t%d\t", m);
388 
389  printf("N=%d\tt_fft=%1.1" NFFT__FES__ "\tt_dft_mpolar=%1.1" NFFT__FES__ "\tm=%d\t", N, t_fft,
390  t_dft_mpolar, m);
391 
392  mpolar_fft(f_hat, N, f, T, S, m);
393  fprintf(fp, "%1.1" NFFT__FES__ "&\t", GLOBAL_elapsed_time);
394  printf("t_mpolar=%1.1" NFFT__FES__ "\t", GLOBAL_elapsed_time);
395  inverse_mpolar_fft(f, T, S, f_hat, N, 2 * m, m);
396  if (m == 9)
397  fprintf(fp, "%1.1" NFFT__FES__ "\\\\\\hline\n", GLOBAL_elapsed_time);
398  else
399  fprintf(fp, "%1.1" NFFT__FES__ "\\\\\n", GLOBAL_elapsed_time);
400  printf("t_impolar=%1.1" NFFT__FES__ "\n", GLOBAL_elapsed_time);
401  }
402 
403  fflush(fp);
404 
405  NFFT(free)(f);
406  NFFT(free)(f_hat);
407 
408  return EXIT_SUCCESS;
409 }
410 
412 int main(int argc, char **argv)
413 {
414  int N;
415  int T, S;
416  int M;
417  NFFT_R *x, *w;
418  NFFT_C *f_hat, *f, *f_direct, *f_tilde;
419  int k;
420  int max_i;
421  int m;
422  NFFT_R temp1, temp2, E_max = NFFT_K(0.0);
423  FILE *fp1, *fp2;
424  char filename[30];
425  int logN;
426 
427  if (argc != 4)
428  {
429  printf("mpolar_fft_test N T R \n");
430  printf("\n");
431  printf("N mpolar FFT of size NxN \n");
432  printf("T number of slopes \n");
433  printf("R number of offsets \n");
434 
436  printf("\nHence, comparison FFTW, mpolar FFT and inverse mpolar FFT\n");
437  fp1 = fopen("mpolar_comparison_fft.dat", "w");
438  if (fp1 == NULL)
439  return (-1);
440  for (logN = 4; logN <= 8; logN++)
441  comparison_fft(fp1, (int)(1U << logN), 3 * (int)(1U << logN),
442  3 * (int)(1U << (logN - 1)));
443  fclose(fp1);
444 
445  exit(EXIT_FAILURE);
446  }
447 
448  N = atoi(argv[1]);
449  T = atoi(argv[2]);
450  S = atoi(argv[3]);
451  printf("N=%d, modified polar grid with T=%d, R=%d => ", N, T, S);
452 
453  x = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * S / 2) * (sizeof(NFFT_R)));
454  w = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * S / 4) * (sizeof(NFFT_R)));
455 
456  f_hat = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(N * N));
457  f = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(1.25 * T * S)); /* 4/pi*log(1+sqrt(2)) = 1.122... < 1.25 */
458  f_direct = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(1.25 * T * S));
459  f_tilde = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(N * N));
460 
462  M = mpolar_grid(T, S, x, w);
463  printf("M=%d.\n", M);
464 
466  fp1 = fopen("input_data_r.dat", "r");
467  fp2 = fopen("input_data_i.dat", "r");
468  if ((fp1 == NULL) || (fp2 == NULL))
469  return (-1);
470  for (k = 0; k < N * N; k++)
471  {
472  fscanf(fp1, NFFT__FR__ " ", &temp1);
473  fscanf(fp2, NFFT__FR__ " ", &temp2);
474  f_hat[k] = temp1 + _Complex_I * temp2;
475  }
476  fclose(fp1);
477  fclose(fp2);
478 
480  mpolar_dft(f_hat, N, f_direct, T, S, 1);
481  // mpolar_fft(f_hat,N,f_direct,T,R,12);
482 
484  printf("\nTest of the mpolar FFT: \n");
485  fp1 = fopen("mpolar_fft_error.dat", "w+");
486  for (m = 1; m <= 12; m++)
487  {
489  mpolar_fft(f_hat, N, f, T, S, m);
490 
492  E_max = NFFT(error_l_infty_complex)(f_direct, f, M);
493  printf("m=%2d: E_max = %" NFFT__FES__ "\n", m, E_max);
494  fprintf(fp1, "%" NFFT__FES__ "\n", E_max);
495  }
496  fclose(fp1);
497 
499  for (m = 3; m <= 9; m += 3)
500  {
501  printf("\nTest of the inverse mpolar FFT for m=%d: \n", m);
502  sprintf(filename, "mpolar_ifft_error%d.dat", m);
503  fp1 = fopen(filename, "w+");
504  for (max_i = 0; max_i <= 20; max_i += 2)
505  {
507  inverse_mpolar_fft(f_direct, T, S, f_tilde, N, max_i, m);
508 
510  E_max = NFFT(error_l_infty_complex)(f_hat, f_tilde, N * N);
511  printf("%3d iterations: E_max = %" NFFT__FES__ "\n", max_i, E_max);
512  fprintf(fp1, "%" NFFT__FES__ "\n", E_max);
513  }
514  fclose(fp1);
515  }
516 
518  NFFT(free)(x);
519  NFFT(free)(w);
520  NFFT(free)(f_hat);
521  NFFT(free)(f);
522  NFFT(free)(f_direct);
523  NFFT(free)(f_tilde);
524 
525  return EXIT_SUCCESS;
526 }
527 /* \} */
static int max_i(int a, int b)
max
Definition: fastsum.c:46
static int comparison_fft(FILE *fp, int N, int T, int S)
Comparison of the FFTW, mpolar FFT, and inverse mpolar FFT.
int main(int argc, char **argv)
test program for various parameters
static int mpolar_dft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int S, int m)
discrete 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
static int mpolar_fft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int S, int m)
NFFT-based mpolar FFT.
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.
#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