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);
68 ths->idx0 = (INT*) NFFT(malloc)((size_t)(M) *
sizeof(INT));
69 ths->deltax0 = (R*) NFFT(malloc)((size_t)(M) *
sizeof(R));
83 NFFT(plan)* cths = (NFFT(plan)*) ths;
85 for (j = 0; j < cths->M_total; j++)
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));
103 NFFT(free)(ths->deltax0);
104 NFFT(free)(ths->idx0);
106 NFFT(finalize)((NFFT(plan)*) ths);
126 NFFT(plan) *cths = (NFFT(plan)*) ths;
128 for (j = 0, f = cths->f; j < cths->M_total; j++)
131 for (k = 0; k < cths->n_total; k++)
132 cths->g1[k] = K(0.0);
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++);
138 cths->g1[0] = cths->f_hat[cths->N_total / 2];
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++);
144 for (l = cths->m - 1; l >= 0; l--)
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));
150 for (k = 1, g1 = cths->g1 + 1; k < cths->N_total / 2; k++)
151 (*g1++) /= (-K2PI * II * (R)(k));
153 FFTW(execute)(cths->my_fftw_plan1);
155 ll = (l == 0 ? 1 : l);
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);
177 int m_taylor,
unsigned test_accuracy)
180 R t_ndft, t_nfft, t_taylor, t;
187 printf(
"%d\t%d\t", N, M);
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);
197 np.f_hat = tp.p.f_hat;
202 swapndft = (C*) NFFT(malloc)((size_t)(M) *
sizeof(C));
205 NFFT(vrand_shifted_unit_double)(np.x, np.M_total);
211 if (np.flags & PRE_ONE_PSI)
212 NFFT(precompute_one_psi)(&np);
215 NFFT(vrand_unit_complex)(np.f_hat, np.N_total);
220 CSWAP(np.f, swapndft);
224 while (t_ndft < K(0.01))
228 NFFT(trafo_direct)(&np);
230 t = NFFT(elapsed_seconds)(t1, t0);
235 CSWAP(np.f, swapndft);
236 printf(
"%.2" __FES__
"\t", t_ndft);
244 while (t_nfft < K(0.01))
250 t = NFFT(elapsed_seconds)(t1, t0);
255 printf(
"%.2" __FES__
"\t%d\t%.2" __FES__
"\t", ((R)(n)) / ((R)(N)), m, t_nfft);
258 printf(
"%.2" __FES__
"\t", NFFT(error_l_infty_complex)(swapndft, np.f, np.M_total));
265 while (t_taylor < K(0.01))
271 t = NFFT(elapsed_seconds)(t1, t0);
276 printf(
"%.2" __FES__
"\t%d\t%.2" __FES__
"\t", ((R)(n_taylor)) / ((R)(N)), m_taylor, t_taylor);
279 printf(
"%.2" __FES__
"\n", NFFT(error_l_infty_complex)(swapndft, np.f, np.M_total));
287 NFFT(free)(swapndft);
293 int main(
int argc,
char **argv)
300 "taylor_nfft type first last trials sigma_nfft sigma_taylor.\n");
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");
309 if (atoi(argv[1]) == 0)
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++)
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++)
329 if (atoi(argv[1]) == 1)
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++)
340 for (trial = 0; trial < arg4; trial++)
static void taylor_finalize(taylor_plan *ths)
Destroys a transform plan.
static void taylor_precompute(taylor_plan *ths)
Precomputation of weights and indices in Taylor expansion.
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...
static void taylor_init(taylor_plan *ths, int N, int M, int n, int m)
Initialisation of a transform plan.
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.
#define CSWAP(x, y)
Swap two vectors.