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