NFFT  3.3.0
flags.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 
29 #include "config.h"
30 
31 #include <stdio.h>
32 #include <math.h>
33 #include <string.h>
34 #include <stdlib.h>
35 #ifdef HAVE_COMPLEX_H
36 #include <complex.h>
37 #endif
38 
39 #include "nfft3.h"
40 #include "infft.h"
41 
42 #ifdef GAUSSIAN
43  unsigned test_fg=1;
44 #else
45  unsigned test_fg=0;
46 #endif
47 
48 #ifdef MEASURE_TIME_FFTW
49  unsigned test_fftw=1;
50 #else
51  unsigned test_fftw=0;
52 #endif
53 
54 #ifdef MEASURE_TIME
55  unsigned test=1;
56 #else
57  unsigned test=0;
58 #endif
59 
60 static void flags_cp(NFFT(plan) *dst, NFFT(plan) *src)
61 {
62  dst->x = src->x;
63  dst->f_hat = src->f_hat;
64  dst->f = src->f;
65  dst->g1 = src->g1;
66  dst->g2 = src->g2;
67  dst->my_fftw_plan1 = src->my_fftw_plan1;
68  dst->my_fftw_plan2 = src->my_fftw_plan2;
69 }
70 
71 static void time_accuracy(int d, int N, int M, int n, int m, unsigned test_ndft,
72  unsigned test_pre_full_psi)
73 {
74  int r, NN[d], nn[d];
75  R t_ndft, t, e;
76  C *swapndft = NULL;
77  ticks t0, t1;
78 
79  NFFT(plan) p;
80  NFFT(plan) p_pre_phi_hut;
81  NFFT(plan) p_fg_psi;
82  NFFT(plan) p_pre_lin_psi;
83  NFFT(plan) p_pre_fg_psi;
84  NFFT(plan) p_pre_psi;
85  NFFT(plan) p_pre_full_psi;
86 
87  printf("%d\t%d\t", d, N);
88 
89  for (r = 0; r < d; r++)
90  {
91  NN[r] = N;
92  nn[r] = n;
93  }
94 
95  /* output vector ndft */
96  if (test_ndft)
97  swapndft = (C*) NFFT(malloc)((size_t)(M) * sizeof(C));
98 
99  NFFT(init_guru)(&p, d, NN, M, nn, m,
100  MALLOC_X | MALLOC_F_HAT | MALLOC_F |
101  FFTW_INIT | FFT_OUT_OF_PLACE,
102  FFTW_MEASURE | FFTW_DESTROY_INPUT);
103 
105  NFFT(vrand_shifted_unit_double)(p.x, p.d * p.M_total);
106 
107  NFFT(init_guru)(&p_pre_phi_hut, d, NN, M, nn, m, PRE_PHI_HUT, 0);
108  flags_cp(&p_pre_phi_hut, &p);
109  NFFT(precompute_one_psi)(&p_pre_phi_hut);
110 
111  if (test_fg)
112  {
113  NFFT(init_guru)(&p_fg_psi, d, NN, M, nn, m, FG_PSI, 0);
114  flags_cp(&p_fg_psi, &p);
115  NFFT(precompute_one_psi)(&p_fg_psi);
116  }
117 
118  NFFT(init_guru)(&p_pre_lin_psi, d, NN, M, nn, m, PRE_LIN_PSI, 0);
119  flags_cp(&p_pre_lin_psi, &p);
120  NFFT(precompute_one_psi)(&p_pre_lin_psi);
121 
122  if (test_fg)
123  {
124  NFFT(init_guru)(&p_pre_fg_psi, d, NN, M, nn, m, PRE_FG_PSI, 0);
125  flags_cp(&p_pre_fg_psi, &p);
126  NFFT(precompute_one_psi)(&p_pre_fg_psi);
127  }
128 
129  NFFT(init_guru)(&p_pre_psi, d, NN, M, nn, m, PRE_PSI, 0);
130  flags_cp(&p_pre_psi, &p);
131  NFFT(precompute_one_psi)(&p_pre_psi);
132 
133  if (test_pre_full_psi)
134  {
135  NFFT(init_guru)(&p_pre_full_psi, d, NN, M, nn, m, PRE_FULL_PSI, 0);
136  flags_cp(&p_pre_full_psi, &p);
137  NFFT(precompute_one_psi)(&p_pre_full_psi);
138  }
139 
140  /* init pseudo random Fourier coefficients */
141  NFFT(vrand_unit_complex)(p.f_hat, p.N_total);
142 
143  /* NDFT */
144  if (test_ndft)
145  {
146  CSWAP(p.f, swapndft);
147 
148  t_ndft = K(0.0);
149  r = 0;
150  while (t_ndft < K(0.01))
151  {
152  r++;
153  t0 = getticks();
154  NFFT(trafo_direct)(&p);
155  t1 = getticks();
156  t = NFFT(elapsed_seconds)(t1, t0);
157  t_ndft += t;
158  }
159  t_ndft /= (R)(r);
160 
161  CSWAP(p.f, swapndft);
162  }
163  else
164  t_ndft = MKNAN("");
165 
166  /* NFFTs */
167  NFFT(trafo)(&p);
168  NFFT(trafo)(&p_pre_phi_hut);
169  if (test_fg)
170  NFFT(trafo)(&p_fg_psi);
171  else
172  p_fg_psi.MEASURE_TIME_t[2] = MKNAN("");
173  NFFT(trafo)(&p_pre_lin_psi);
174  if (test_fg)
175  NFFT(trafo)(&p_pre_fg_psi);
176  else
177  p_pre_fg_psi.MEASURE_TIME_t[2] = MKNAN("");
178  NFFT(trafo)(&p_pre_psi);
179  if (test_pre_full_psi)
180  NFFT(trafo)(&p_pre_full_psi);
181  else
182  p_pre_full_psi.MEASURE_TIME_t[2] = MKNAN("");
183 
184  if (test_fftw == 0)
185  p.MEASURE_TIME_t[1] = MKNAN("");
186 
187  if (test_ndft)
188  e = NFFT(error_l_2_complex)(swapndft, p.f, p.M_total);
189  else
190  e = MKNAN("");
191 
192  printf(
193  "%.2" __FES__ "\t%d\t%.2" __FES__ "\t%.2" __FES__ "\t%.2" __FES__ "\t%.2" __FES__ "\t%.2" __FES__ "\t%.2" __FES__ "\t%.2" __FES__ "\t%.2" __FES__ "\t%.2" __FES__ "\t%.2" __FES__ "\n",
194  t_ndft, m, e, p.MEASURE_TIME_t[0], p_pre_phi_hut.MEASURE_TIME_t[0],
195  p.MEASURE_TIME_t[1], p.MEASURE_TIME_t[2], p_fg_psi.MEASURE_TIME_t[2],
196  p_pre_lin_psi.MEASURE_TIME_t[2], p_pre_fg_psi.MEASURE_TIME_t[2],
197  p_pre_psi.MEASURE_TIME_t[2], p_pre_full_psi.MEASURE_TIME_t[2]);
198 
199  fflush(stdout);
200 
202  if (test_pre_full_psi)
203  NFFT(finalize)(&p_pre_full_psi);
204  NFFT(finalize)(&p_pre_psi);
205  if (test_fg)
206  NFFT(finalize)(&p_pre_fg_psi);
207  NFFT(finalize)(&p_pre_lin_psi);
208  if (test_fg)
209  NFFT(finalize)(&p_fg_psi);
210  NFFT(finalize)(&p_pre_phi_hut);
211  NFFT(finalize)(&p);
212 
213  if (test_ndft)
214  NFFT(free)(swapndft);
215 }
216 
217 static void accuracy_pre_lin_psi(int d, int N, int M, int n, int m, int K)
218 {
219  int r, NN[d], nn[d];
220  R e;
221  C *swapndft;
222 
223  NFFT(plan) p;
224 
225  for (r = 0; r < d; r++)
226  {
227  NN[r] = N;
228  nn[r] = n;
229  }
230 
231  /* output vector ndft */
232  swapndft = (C*) NFFT(malloc)((size_t)(M) * sizeof(C));
233 
234  NFFT(init_guru)(&p, d, NN, M, nn, m,
235  MALLOC_X | MALLOC_F_HAT | MALLOC_F |
236  PRE_PHI_HUT | PRE_LIN_PSI |
237  FFTW_INIT | FFT_OUT_OF_PLACE,
238  FFTW_MEASURE | FFTW_DESTROY_INPUT);
239 
241  NFFT(free)(p.psi);
242  p.K = K;
243  p.psi = (R*) NFFT(malloc)((size_t)((p.K + 1) * p.d) * sizeof(R));
244 
246  NFFT(precompute_one_psi)(&p);
247 
249  NFFT(vrand_shifted_unit_double)(p.x, p.d * p.M_total);
250 
252  NFFT(vrand_unit_complex)(p.f_hat, p.N_total);
253 
255  CSWAP(p.f, swapndft);
256  NFFT(trafo_direct)(&p);
257  CSWAP(p.f, swapndft);
258 
260  NFFT(trafo)(&p);
261  e = NFFT(error_l_2_complex)(swapndft, p.f, p.M_total);
262 
263  // printf("%d\t%d\t%d\t%d\t%.2e\n",d,N,m,K,e);
264  printf("$%.1" __FES__ "$&\t", e);
265 
266  fflush(stdout);
267 
269  NFFT(finalize)(&p);
270  NFFT(free)(swapndft);
271 }
272 
273 int main(int argc, char **argv)
274 {
275  int l, trial;
276 
277  if (argc <= 2)
278  {
279  fprintf(stderr, "flags type first last trials d m\n");
280  return EXIT_FAILURE;
281  }
282 
283  if ((test == 0) && (atoi(argv[1]) < 2))
284  {
285  fprintf(stderr, "MEASURE_TIME in infft.h not set\n");
286  return EXIT_FAILURE;
287  }
288 
289  fprintf(stderr, "Testing different precomputation schemes for the nfft.\n");
290  fprintf(stderr, "Columns: d, N=M, t_ndft, e_nfft, t_D, t_pre_phi_hut, ");
291  fprintf(stderr, "t_fftw, t_B, t_fg_psi, t_pre_lin_psi, t_pre_fg_psi, ");
292  fprintf(stderr, "t_pre_psi, t_pre_full_psi\n\n");
293 
294  int arg2 = atoi(argv[2]);
295  int arg3 = atoi(argv[3]);
296  int arg4 = atoi(argv[4]);
297 
298  /* time vs. N=M */
299  if (atoi(argv[1]) == 0)
300  {
301  int d = atoi(argv[5]);
302  int m = atoi(argv[6]);
303 
304  for (l = arg2; l <= arg3; l++)
305  {
306  int N = (int)(1U << l);
307  int M = (int)(1U << (d * l));
308  for (trial = 0; trial < arg4; trial++)
309  {
310  time_accuracy(d, N, M, 2 * N, m, 0, 0);
311  }
312  }
313  }
314  else if (atoi(argv[1]) == 1) /* accuracy vs. time */
315  {
316  int d = atoi(argv[5]);
317  int N = atoi(argv[6]);
318  int m;
319 
320  for (m = arg2; m <= arg3; m++)
321  {
322  for (trial = 0; trial < arg4; trial++)
323  {
324  time_accuracy(d, N, (int)(LRINT(POW((R)(N), (R)(d)))), 2 * N, m, 1, 1);
325  }
326  }
327  }
328  else if (atoi(argv[1]) == 2) /* accuracy vs. K for linear interpolation, assumes (m+1)|K */
329  {
330  int d = atoi(argv[5]);
331  int N = atoi(argv[6]);
332  int m = atoi(argv[7]);
333 
334  printf("$\\log_2(K/(m+1))$&\t");
335 
336  for (l = arg2; l < arg3; l++)
337  printf("$%d$&\t", l);
338 
339  printf("$%d$\\\\\n", arg3);
340 
341  printf("$\\tilde E_2$&\t");
342  for (l = arg2; l <= arg3; l++)
343  {
344  int x = (m + 1) * (int)(1U << l);
345  accuracy_pre_lin_psi(d, N, (int)(LRINT(POW((R)(N), (R)(d)))), 2 * N, m, x);
346  }
347 
348  printf("\n");
349  }
350 
351 
352  return EXIT_SUCCESS;
353 }
#define CSWAP(x, y)
Swap two vectors.
Definition: infft.h:152