NFFT 3.6.0
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
43NFFT_R GLOBAL_elapsed_time;
44
58static 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
102static 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 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
164static 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,
191 PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | MALLOC_F | FFTW_INIT
192 | FFT_OUT_OF_PLACE,
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
238static 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,
268 PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | MALLOC_F | FFTW_INIT
269 | FFT_OUT_OF_PLACE,
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
347static 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
412int 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/* \} */
max_i
static int max_i(int a, int b)
max
Definition fastsum.c:46
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:347
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
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:238
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:164
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_HAT
#define MALLOC_F_HAT
Definition nfft3.h:194
MALLOC_X
#define MALLOC_X
Definition nfft3.h:193
PRE_FULL_PSI
#define PRE_FULL_PSI
Definition nfft3.h:192
FFT_OUT_OF_PLACE
#define FFT_OUT_OF_PLACE
Definition nfft3.h:196
PRE_PSI
#define PRE_PSI
Definition nfft3.h:191
MALLOC_F
#define MALLOC_F
Definition nfft3.h:195
PRE_LIN_PSI
#define PRE_LIN_PSI
Definition nfft3.h:189
FFTW_INIT
#define FFTW_INIT
Definition nfft3.h:197
PRE_PHI_HUT
#define PRE_PHI_HUT
Definition nfft3.h:187
CGNR
#define CGNR
Definition nfft3.h:808
PRECOMPUTE_DAMP
#define PRECOMPUTE_DAMP
Definition nfft3.h:812
PRECOMPUTE_WEIGHT
#define PRECOMPUTE_WEIGHT
Definition nfft3.h:811
Generated by doxygen 1.9.8