NFFT 3.5.3alpha
ndft_fast.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
25#include "config.h"
26
27#include <stdio.h>
28#include <math.h>
29#include <string.h>
30#include <stdlib.h>
31#ifdef HAVE_COMPLEX_H
32#include <complex.h>
33#endif
34
35#include "nfft3.h"
36#include "infft.h"
37
38static void ndft_horner_trafo(NFFT(plan) *ths)
39{
40 INT j, k;
41 C *f_hat_k, *f_j;
42 C exp_omega_0;
43
44 for (j = 0, f_j = ths->f; j < ths->M_total; j++, f_j++)
45 (*f_j) = K(0.0);
46
47 for (j = 0, f_j = ths->f; j < ths->M_total; j++, f_j++)
48 {
49 exp_omega_0 = CEXP(+K2PI * II * ths->x[j]);
50 for (k = 0, f_hat_k = ths->f_hat; k < ths->N[0]; k++, f_hat_k++)
51 {
52 (*f_j) += (*f_hat_k);
53 (*f_j) *= exp_omega_0;
54 }
55 (*f_j) *= CEXP(-KPI * II * (R)(ths->N[0]) * ths->x[j]);
56 }
57} /* ndft_horner_trafo */
58
59static void ndft_pre_full_trafo(NFFT(plan) *ths, C *A)
60{
61 INT j, k;
62 C *f_hat_k, *f_j;
63 C *A_local;
64
65 for (j = 0, f_j = ths->f; j < ths->M_total; j++, f_j++)
66 (*f_j) = K(0.0);
67
68 for (j = 0, f_j = ths->f, A_local = A; j < ths->M_total; j++, f_j++)
69 for (k = 0, f_hat_k = ths->f_hat; k < ths->N[0]; k++, f_hat_k++, A_local++)
70 (*f_j) += (*f_hat_k) * (*A_local);
71} /* ndft_pre_full_trafo */
72
73static void ndft_pre_full_init(NFFT(plan) *ths, C *A)
74{
75 INT j, k;
76 C *f_hat_k, *f_j, *A_local;
77
78 for (j = 0, f_j = ths->f, A_local = A; j < ths->M_total; j++, f_j++)
79 for (k = 0, f_hat_k = ths->f_hat; k < ths->N[0]; k++, f_hat_k++, A_local++)
80 (*A_local) = CEXP(
81 -K2PI * II * ((R) (k) - (R) (ths->N[0]) / K(2.0)) * ths->x[j]);
82
83} /* ndft_pre_full_init */
84
85static void ndft_time(int N, int M, unsigned test_ndft, unsigned test_pre_full)
86{
87 int r;
88 R t, t_ndft, t_horner, t_pre_full, t_nfft;
89 C *A = NULL;
90 ticks t0, t1;
91
92 NFFT(plan) np;
93
94 printf("%d\t%d\t", N, M);
95
96 NFFT(init_1d)(&np, N, M);
97
98 /* Initialize pseudo random nodes. */
99 NFFT(vrand_shifted_unit_double)(np.x, np.M_total);
100
101 if (test_pre_full)
102 {
103 A = (C*) NFFT(malloc)((size_t)(N * M) * sizeof(R));
104 ndft_pre_full_init(&np, A);
105 }
106
107 /* Initialize pseudo random Fourier coefficients. */
108 NFFT(vrand_unit_complex)(np.f_hat, np.N_total);
109
110 /* NDFT */
111 if (test_ndft)
112 {
113 t_ndft = K(0.0);
114 r = 0;
115 while (t_ndft < K(0.1))
116 {
117 r++;
118 t0 = getticks();
119 NFFT(trafo_direct)(&np);
120 t1 = getticks();
121 t = NFFT(elapsed_seconds)(t1, t0);
122 t_ndft += t;
123 }
124 t_ndft /= (R) (r);
125
126 printf("%.2" __FES__ "\t", t_ndft);
127 }
128 else
129 printf("N/A\t\t");
130
131 /* Horner NDFT */
132 t_horner = K(0.0);
133 r = 0;
134 while (t_horner < K(0.1))
135 {
136 r++;
137 t0 = getticks();
138 ndft_horner_trafo(&np);
139 t1 = getticks();
140 t = NFFT(elapsed_seconds)(t1, t0);
141 t_horner += t;
142 }
143 t_horner /= (R)(r);
144
145 printf("%.2" __FES__ "\t", t_horner);
146
147 /* Fully precomputed NDFT */
148 if (test_pre_full)
149 {
150 t_pre_full = K(0.0);
151 r = 0;
152 while (t_pre_full < K(0.1))
153 {
154 r++;
155 t0 = getticks();
156 ndft_pre_full_trafo(&np, A);
157 t1 = getticks();
158 t = NFFT(elapsed_seconds)(t1, t0);
159 t_pre_full += t;
160 }
161 t_pre_full /= (R)(r);
162
163 printf("%.2" __FES__ "\t", t_pre_full);
164 }
165 else
166 printf("N/A\t\t");
167
168 t_nfft = K(0.0);
169 r = 0;
170 while (t_nfft < K(0.1))
171 {
172 r++;
173 t0 = getticks();
174 NFFT(trafo)(&np);
175 t1 = getticks();
176 t = NFFT(elapsed_seconds)(t1, t0);
177 t_nfft += t;
178 }
179 t_nfft /= (R)(r);
180
181 printf("%.2" __FES__ "\n", t_nfft);
182
183 fflush(stdout);
184
185 if (test_pre_full)
186 NFFT(free)(A);
187
188 NFFT(finalize)(&np);
189}
190
191int main(int argc, char **argv)
192{
193 int l, trial;
194
195 if (argc < 4)
196 {
197 fprintf(stderr, "ndft_fast type first last trials\n");
198 return EXIT_FAILURE;
199 }
200 else
201 {
202 int arg2 = (atoi(argv[2]));
203 int arg3 = (atoi(argv[3]));
204 int arg4 = (atoi(argv[4]));
205 fprintf(stderr, "Testing ndft, Horner-like ndft, fully precomputed ndft.\n");
206 fprintf(stderr, "Columns: N, M, t_ndft, t_horner, t_pre_full, t_nfft\n\n");
207
208 /* time vs. N=M */
209 if (atoi(argv[1]) == 0)
210 {
211 for (l = arg2; l <= arg3; l++)
212 {
213 for (trial = 0; trial < arg4; trial++)
214 {
215 int N = (int)(1U << l);
216 int M = (int)(1U << l);
217 ndft_time(N, M, l <= 15 ? 1 : 0, l <= 13 ? 1 : 0);
218 }
219 }
220 }
221 }
222
223 return EXIT_SUCCESS;
224}
Internal header file for auxiliary definitions and functions.
Header file for the nfft3 library.