40 static void ndft_horner_trafo(NFFT(plan) *ths)
46 for (j = 0, f_j = ths->f; j < ths->M_total; j++, f_j++)
49 for (j = 0, f_j = ths->f; j < ths->M_total; j++, f_j++)
51 exp_omega_0 = CEXP(+K2PI * II * ths->x[j]);
52 for (k = 0, f_hat_k = ths->f_hat; k < ths->N[0]; k++, f_hat_k++)
55 (*f_j) *= exp_omega_0;
57 (*f_j) *= CEXP(-KPI * II * (R)(ths->N[0]) * ths->x[j]);
61 static void ndft_pre_full_trafo(NFFT(plan) *ths, C *A)
67 for (j = 0, f_j = ths->f; j < ths->M_total; j++, f_j++)
70 for (j = 0, f_j = ths->f, A_local = A; j < ths->M_total; j++, f_j++)
71 for (k = 0, f_hat_k = ths->f_hat; k < ths->N[0]; k++, f_hat_k++, A_local++)
72 (*f_j) += (*f_hat_k) * (*A_local);
75 static void ndft_pre_full_init(NFFT(plan) *ths, C *A)
78 C *f_hat_k, *f_j, *A_local;
80 for (j = 0, f_j = ths->f, A_local = A; j < ths->M_total; j++, f_j++)
81 for (k = 0, f_hat_k = ths->f_hat; k < ths->N[0]; k++, f_hat_k++, A_local++)
83 -K2PI * II * ((R) (k) - (R) (ths->N[0]) / K(2.0)) * ths->x[j]);
87 static void ndft_time(
int N,
int M,
unsigned test_ndft,
unsigned test_pre_full)
90 R t, t_ndft, t_horner, t_pre_full, t_nfft;
96 printf("%d\t%d\t", N, M);
98 NFFT(init_1d)(&np, N, M);
101 NFFT(vrand_shifted_unit_double)(np.x, np.M_total);
105 A = (C*) NFFT(malloc)((size_t)(N * M) *
sizeof(R));
106 ndft_pre_full_init(&np, A);
110 NFFT(vrand_unit_complex)(np.f_hat, np.N_total);
117 while (t_ndft < K(0.1))
121 NFFT(trafo_direct)(&np);
123 t = NFFT(elapsed_seconds)(t1, t0);
128 printf(
"%.2" __FES__
"\t", t_ndft);
136 while (t_horner < K(0.1))
140 ndft_horner_trafo(&np);
142 t = NFFT(elapsed_seconds)(t1, t0);
147 printf(
"%.2" __FES__
"\t", t_horner);
154 while (t_pre_full < K(0.1))
158 ndft_pre_full_trafo(&np, A);
160 t = NFFT(elapsed_seconds)(t1, t0);
163 t_pre_full /= (R)(r);
165 printf(
"%.2" __FES__
"\t", t_pre_full);
172 while (t_nfft < K(0.1))
178 t = NFFT(elapsed_seconds)(t1, t0);
183 printf(
"%.2" __FES__
"\n", t_nfft);
193 int main(
int argc,
char **argv)
199 fprintf(stderr,
"ndft_fast type first last trials\n");
204 int arg2 = (atoi(argv[2]));
205 int arg3 = (atoi(argv[3]));
206 int arg4 = (atoi(argv[4]));
207 fprintf(stderr,
"Testing ndft, Horner-like ndft, fully precomputed ndft.\n");
208 fprintf(stderr,
"Columns: N, M, t_ndft, t_horner, t_pre_full, t_nfft\n\n");
211 if (atoi(argv[1]) == 0)
213 for (l = arg2; l <= arg3; l++)
215 for (trial = 0; trial < arg4; trial++)
217 int N = (int)(1U << l);
218 int M = (int)(1U << l);
219 ndft_time(N, M, l <= 15 ? 1 : 0, l <= 13 ? 1 : 0);