NFFT  3.3.0
taylor_nfft.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 
30 #include "config.h"
31 
32 #include <stdio.h>
33 #include <math.h>
34 #include <string.h>
35 #include <stdlib.h>
36 #ifdef HAVE_COMPLEX_H
37 #include <complex.h>
38 #endif
39 
40 #include "nfft3.h"
41 #include "infft.h"
42 
43 typedef struct
44 {
45  NFFT(plan) p; /* used for fftw and data */
46  INT *idx0; /* index of next neighbour of x_j on the oversampled regular grid */
47  R *deltax0; /* distance to the grid point */
48 } taylor_plan;
49 
61 static void taylor_init(taylor_plan *ths, int N, int M, int n, int m)
62 {
63  /* Note: no nfft precomputation! */
64  NFFT(init_guru)((NFFT(plan)*) ths, 1, &N, M, &n, m,
65  MALLOC_X | MALLOC_F_HAT | MALLOC_F | FFTW_INIT | FFT_OUT_OF_PLACE,
66  FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);
67 
68  ths->idx0 = (INT*) NFFT(malloc)((size_t)(M) * sizeof(INT));
69  ths->deltax0 = (R*) NFFT(malloc)((size_t)(M) * sizeof(R));
70 }
71 
79 static void taylor_precompute(taylor_plan *ths)
80 {
81  INT j;
82 
83  NFFT(plan)* cths = (NFFT(plan)*) ths;
84 
85  for (j = 0; j < cths->M_total; j++)
86  {
87  ths->idx0[j] = (LRINT(ROUND((cths->x[j] + K(0.5)) * (R)(cths->n[0])))
88  + cths->n[0] / 2) % cths->n[0];
89  ths->deltax0[j] = cths->x[j]
90  - (ROUND((cths->x[j] + K(0.5)) * (R)(cths->n[0])) / (R)(cths->n[0]) - K(0.5));
91  }
92 }
93 
101 static void taylor_finalize(taylor_plan *ths)
102 {
103  NFFT(free)(ths->deltax0);
104  NFFT(free)(ths->idx0);
105 
106  NFFT(finalize)((NFFT(plan)*) ths);
107 }
108 
119 static void taylor_trafo(taylor_plan *ths)
120 {
121  INT j, k, l, ll;
122  C *f, *f_hat, *g1;
123  R *deltax;
124  INT *idx;
125 
126  NFFT(plan) *cths = (NFFT(plan)*) ths;
127 
128  for (j = 0, f = cths->f; j < cths->M_total; j++)
129  *f++ = K(0.0);
130 
131  for (k = 0; k < cths->n_total; k++)
132  cths->g1[k] = K(0.0);
133 
134  for (k = -cths->N_total / 2, g1 = cths->g1 + cths->n_total
135  - cths->N_total / 2, f_hat = cths->f_hat; k < 0; k++)
136  (*g1++) = CPOW(-K2PI * II * (R)(k), (R)(cths->m)) * (*f_hat++);
137 
138  cths->g1[0] = cths->f_hat[cths->N_total / 2];
139 
140  for (k = 1, g1 = cths->g1 + 1, f_hat = cths->f_hat + cths->N_total / 2 + 1;
141  k < cths->N_total / 2; k++)
142  (*g1++) = CPOW(-K2PI * II * (R)(k), (R)(cths->m)) * (*f_hat++);
143 
144  for (l = cths->m - 1; l >= 0; l--)
145  {
146  for (k = -cths->N_total / 2, g1 = cths->g1 + cths->n_total
147  - cths->N_total / 2; k < 0; k++)
148  (*g1++) /= (-K2PI * II * (R)(k));
149 
150  for (k = 1, g1 = cths->g1 + 1; k < cths->N_total / 2; k++)
151  (*g1++) /= (-K2PI * II * (R)(k));
152 
153  FFTW(execute)(cths->my_fftw_plan1);
154 
155  ll = (l == 0 ? 1 : l);
156 
157  for (j = 0, f = cths->f, deltax = ths->deltax0, idx = ths->idx0;
158  j < cths->M_total; j++, f++)
159  (*f) = ((*f) * (*deltax++) + cths->g2[*idx++]) / (R)(ll);
160  }
161 }
162 
176 static void taylor_time_accuracy(int N, int M, int n, int m, int n_taylor,
177  int m_taylor, unsigned test_accuracy)
178 {
179  int r;
180  R t_ndft, t_nfft, t_taylor, t;
181  C *swapndft = NULL;
182  ticks t0, t1;
183 
184  taylor_plan tp;
185  NFFT(plan) np;
186 
187  printf("%d\t%d\t", N, M);
188 
189  taylor_init(&tp, N, M, n_taylor, m_taylor);
190 
191  NFFT(init_guru)(&np, 1, &N, M, &n, m,
192  PRE_PHI_HUT | PRE_FG_PSI | FFTW_INIT | FFT_OUT_OF_PLACE,
193  FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
194 
195  /* share nodes, input, and output vectors */
196  np.x = tp.p.x;
197  np.f_hat = tp.p.f_hat;
198  np.f = tp.p.f;
199 
200  /* output vector ndft */
201  if (test_accuracy)
202  swapndft = (C*) NFFT(malloc)((size_t)(M) * sizeof(C));
203 
204  /* init pseudo random nodes */
205  NFFT(vrand_shifted_unit_double)(np.x, np.M_total);
206 
207  /* nfft precomputation */
208  taylor_precompute(&tp);
209 
210  /* nfft precomputation */
211  if (np.flags & PRE_ONE_PSI)
212  NFFT(precompute_one_psi)(&np);
213 
214  /* init pseudo random Fourier coefficients */
215  NFFT(vrand_unit_complex)(np.f_hat, np.N_total);
216 
217  /* NDFT */
218  if (test_accuracy)
219  {
220  CSWAP(np.f, swapndft);
221 
222  t_ndft = K(0.0);
223  r = 0;
224  while (t_ndft < K(0.01))
225  {
226  r++;
227  t0 = getticks();
228  NFFT(trafo_direct)(&np);
229  t1 = getticks();
230  t = NFFT(elapsed_seconds)(t1, t0);
231  t_ndft += t;
232  }
233  t_ndft /= (R)(r);
234 
235  CSWAP(np.f, swapndft);
236  printf("%.2" __FES__ "\t", t_ndft);
237  }
238  else
239  printf("N/A\t\t");
240 
241  /* NFFT */
242  t_nfft = K(0.0);
243  r = 0;
244  while (t_nfft < K(0.01))
245  {
246  r++;
247  t0 = getticks();
248  NFFT(trafo)(&np);
249  t1 = getticks();
250  t = NFFT(elapsed_seconds)(t1, t0);
251  t_nfft += t;
252  }
253  t_nfft /= (R)(r);
254 
255  printf("%.2" __FES__ "\t%d\t%.2" __FES__ "\t", ((R)(n)) / ((R)(N)), m, t_nfft);
256 
257  if (test_accuracy)
258  printf("%.2" __FES__ "\t", NFFT(error_l_infty_complex)(swapndft, np.f, np.M_total));
259  else
260  printf("N/A\t\t");
261 
263  t_taylor = K(0.0);
264  r = 0;
265  while (t_taylor < K(0.01))
266  {
267  r++;
268  t0 = getticks();
269  taylor_trafo(&tp);
270  t1 = getticks();
271  t = NFFT(elapsed_seconds)(t1, t0);
272  t_taylor += t;
273  }
274  t_taylor /= (R)(r);
275 
276  printf("%.2" __FES__ "\t%d\t%.2" __FES__ "\t", ((R)(n_taylor)) / ((R)(N)), m_taylor, t_taylor);
277 
278  if (test_accuracy)
279  printf("%.2" __FES__ "\n", NFFT(error_l_infty_complex)(swapndft, np.f, np.M_total));
280  else
281  printf("N/A\t\n");
282 
283  fflush(stdout);
284 
285  /* finalise */
286  if (test_accuracy)
287  NFFT(free)(swapndft);
288 
289  NFFT(finalize)(&np);
290  taylor_finalize(&tp);
291 }
292 
293 int main(int argc, char **argv)
294 {
295  int l, m, trial;
296 
297  if (argc <= 2)
298  {
299  fprintf(stderr,
300  "taylor_nfft type first last trials sigma_nfft sigma_taylor.\n");
301  return EXIT_FAILURE;
302  }
303 
304  fprintf(stderr, "Testing the Nfft & a Taylor expansion based version.\n\n");
305  fprintf(stderr, "Columns: N, M, t_ndft, sigma_nfft, m_nfft, t_nfft, e_nfft");
306  fprintf(stderr, ", sigma_taylor, m_taylor, t_taylor, e_taylor\n");
307 
308  /* time vs. N = M */
309  if (atoi(argv[1]) == 0)
310  {
311  fprintf(stderr, "Fixed target accuracy, timings.\n\n");
312  int arg2 = atoi(argv[2]);
313  int arg3 = atoi(argv[3]);
314  int arg4 = atoi(argv[4]);
315  for (l = arg2; l <= arg3; l++)
316  {
317  int N = (int)(1U << l);
318  int M = (int)(1U << l);
319  int arg5 = (int)(atof(argv[5]) * N);
320  int arg6 = (int)(atof(argv[6]) * N);
321  for (trial = 0; trial < arg4; trial++)
322  {
323  taylor_time_accuracy(N, M, arg5, 6, arg6, 6, l <= 10 ? 1 : 0);
324  }
325  }
326  }
327 
328  /* error vs. m */
329  if (atoi(argv[1]) == 1)
330  {
331  int arg2 = atoi(argv[2]);
332  int arg3 = atoi(argv[3]);
333  int arg4 = atoi(argv[4]);
334  int N = atoi(argv[7]);
335  int arg5 = (int) (atof(argv[5]) * N);
336  int arg6 = (int) (atof(argv[6]) * N);
337  fprintf(stderr, "Fixed N=M=%d, error vs. m.\n\n", N);
338  for (m = arg2; m <= arg3; m++)
339  {
340  for (trial = 0; trial < arg4; trial++)
341  {
342  taylor_time_accuracy(N, N, arg5, m, arg6, m, 1);
343  }
344  }
345  }
346 
347  return EXIT_SUCCESS;
348 }
static void taylor_finalize(taylor_plan *ths)
Destroys a transform plan.
Definition: taylor_nfft.c:101
static void taylor_precompute(taylor_plan *ths)
Precomputation of weights and indices in Taylor expansion.
Definition: taylor_nfft.c:79
static void taylor_trafo(taylor_plan *ths)
Executes a Taylor-NFFT, see equation (1.1) in [Guide], computes fast and approximate by means of a Ta...
Definition: taylor_nfft.c:119
static void taylor_init(taylor_plan *ths, int N, int M, int n, int m)
Initialisation of a transform plan.
Definition: taylor_nfft.c:61
static void taylor_time_accuracy(int N, int M, int n, int m, int n_taylor, int m_taylor, unsigned test_accuracy)
Compares NDFT, NFFT, and Taylor-NFFT.
Definition: taylor_nfft.c:176
#define CSWAP(x, y)
Swap two vectors.
Definition: infft.h:152