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
19 /*! \file flags.c
20 *
21 * \brief Testing the nfft.
22 *
23 * \author Stefan Kunis
24 *
25 * References: Time and Memory Requirements of the Nonequispaced FFT
26 */
27 #include "config.h"
28
29 #include <stdio.h>
30 #include <math.h>
31 #include <string.h>
32 #include <stdlib.h>
33 #ifdef HAVE_COMPLEX_H
34 #include <complex.h>
35 #endif
36
37 #include "nfft3.h"
38 #include "infft.h"
39
40 #ifdef GAUSSIAN
41 unsigned test_fg=1;
42 #else
43 unsigned test_fg=0;
44 #endif
45
46 #ifdef MEASURE_TIME_FFTW
47 unsigned test_fftw=1;
48 #else
49 unsigned test_fftw=0;
50 #endif
51
52 #ifdef MEASURE_TIME
53 unsigned test=1;
54 #else
55 unsigned test=0;
56 #endif
57
flags_cp(NFFT (plan)* dst,NFFT (plan)* src)58 static void flags_cp(NFFT(plan) *dst, NFFT(plan) *src)
59 {
60 dst->x = src->x;
61 dst->f_hat = src->f_hat;
62 dst->f = src->f;
63 dst->g1 = src->g1;
64 dst->g2 = src->g2;
65 dst->my_fftw_plan1 = src->my_fftw_plan1;
66 dst->my_fftw_plan2 = src->my_fftw_plan2;
67 }
68
time_accuracy(int d,int N,int M,int n,int m,unsigned test_ndft,unsigned test_pre_full_psi)69 static void time_accuracy(int d, int N, int M, int n, int m, unsigned test_ndft,
70 unsigned test_pre_full_psi)
71 {
72 int r, NN[d], nn[d];
73 R t_ndft, t, e;
74 C *swapndft = NULL;
75 ticks t0, t1;
76
77 NFFT(plan) p;
78 NFFT(plan) p_pre_phi_hut;
79 NFFT(plan) p_fg_psi;
80 NFFT(plan) p_pre_lin_psi;
81 NFFT(plan) p_pre_fg_psi;
82 NFFT(plan) p_pre_psi;
83 NFFT(plan) p_pre_full_psi;
84
85 printf("%d\t%d\t", d, N);
86
87 for (r = 0; r < d; r++)
88 {
89 NN[r] = N;
90 nn[r] = n;
91 }
92
93 /* output vector ndft */
94 if (test_ndft)
95 swapndft = (C*) NFFT(malloc)((size_t)(M) * sizeof(C));
96
97 NFFT(init_guru)(&p, d, NN, M, nn, m,
98 MALLOC_X | MALLOC_F_HAT | MALLOC_F |
99 FFTW_INIT | FFT_OUT_OF_PLACE,
100 FFTW_MEASURE | FFTW_DESTROY_INPUT);
101
102 /** init pseudo random nodes */
103 NFFT(vrand_shifted_unit_double)(p.x, p.d * p.M_total);
104
105 NFFT(init_guru)(&p_pre_phi_hut, d, NN, M, nn, m, PRE_PHI_HUT, 0);
106 flags_cp(&p_pre_phi_hut, &p);
107 NFFT(precompute_one_psi)(&p_pre_phi_hut);
108
109 if (test_fg)
110 {
111 NFFT(init_guru)(&p_fg_psi, d, NN, M, nn, m, FG_PSI, 0);
112 flags_cp(&p_fg_psi, &p);
113 NFFT(precompute_one_psi)(&p_fg_psi);
114 }
115
116 NFFT(init_guru)(&p_pre_lin_psi, d, NN, M, nn, m, PRE_LIN_PSI, 0);
117 flags_cp(&p_pre_lin_psi, &p);
118 NFFT(precompute_one_psi)(&p_pre_lin_psi);
119
120 if (test_fg)
121 {
122 NFFT(init_guru)(&p_pre_fg_psi, d, NN, M, nn, m, PRE_FG_PSI, 0);
123 flags_cp(&p_pre_fg_psi, &p);
124 NFFT(precompute_one_psi)(&p_pre_fg_psi);
125 }
126
127 NFFT(init_guru)(&p_pre_psi, d, NN, M, nn, m, PRE_PSI, 0);
128 flags_cp(&p_pre_psi, &p);
129 NFFT(precompute_one_psi)(&p_pre_psi);
130
131 if (test_pre_full_psi)
132 {
133 NFFT(init_guru)(&p_pre_full_psi, d, NN, M, nn, m, PRE_FULL_PSI, 0);
134 flags_cp(&p_pre_full_psi, &p);
135 NFFT(precompute_one_psi)(&p_pre_full_psi);
136 }
137
138 /* init pseudo random Fourier coefficients */
139 NFFT(vrand_unit_complex)(p.f_hat, p.N_total);
140
141 /* NDFT */
142 if (test_ndft)
143 {
144 CSWAP(p.f, swapndft);
145
146 t_ndft = K(0.0);
147 r = 0;
148 while (t_ndft < K(0.01))
149 {
150 r++;
151 t0 = getticks();
152 NFFT(trafo_direct)(&p);
153 t1 = getticks();
154 t = NFFT(elapsed_seconds)(t1, t0);
155 t_ndft += t;
156 }
157 t_ndft /= (R)(r);
158
159 CSWAP(p.f, swapndft);
160 }
161 else
162 t_ndft = MKNAN("");
163
164 /* NFFTs */
165 NFFT(trafo)(&p);
166 NFFT(trafo)(&p_pre_phi_hut);
167 if (test_fg)
168 NFFT(trafo)(&p_fg_psi);
169 else
170 p_fg_psi.MEASURE_TIME_t[2] = MKNAN("");
171 NFFT(trafo)(&p_pre_lin_psi);
172 if (test_fg)
173 NFFT(trafo)(&p_pre_fg_psi);
174 else
175 p_pre_fg_psi.MEASURE_TIME_t[2] = MKNAN("");
176 NFFT(trafo)(&p_pre_psi);
177 if (test_pre_full_psi)
178 NFFT(trafo)(&p_pre_full_psi);
179 else
180 p_pre_full_psi.MEASURE_TIME_t[2] = MKNAN("");
181
182 if (test_fftw == 0)
183 p.MEASURE_TIME_t[1] = MKNAN("");
184
185 if (test_ndft)
186 e = NFFT(error_l_2_complex)(swapndft, p.f, p.M_total);
187 else
188 e = MKNAN("");
189
190 printf(
191 "%.2" __FES__ "\t%d\t%.2" __FES__ "\t%.2" __FES__ "\t%.2" __FES__ "\t%.2" __FES__ "\t%.2" __FES__ "\t%.2" __FES__ "\t%.2" __FES__ "\t%.2" __FES__ "\t%.2" __FES__ "\t%.2" __FES__ "\n",
192 t_ndft, m, e, p.MEASURE_TIME_t[0], p_pre_phi_hut.MEASURE_TIME_t[0],
193 p.MEASURE_TIME_t[1], p.MEASURE_TIME_t[2], p_fg_psi.MEASURE_TIME_t[2],
194 p_pre_lin_psi.MEASURE_TIME_t[2], p_pre_fg_psi.MEASURE_TIME_t[2],
195 p_pre_psi.MEASURE_TIME_t[2], p_pre_full_psi.MEASURE_TIME_t[2]);
196
197 fflush(stdout);
198
199 /** finalise */
200 if (test_pre_full_psi)
201 NFFT(finalize)(&p_pre_full_psi);
202 NFFT(finalize)(&p_pre_psi);
203 if (test_fg)
204 NFFT(finalize)(&p_pre_fg_psi);
205 NFFT(finalize)(&p_pre_lin_psi);
206 if (test_fg)
207 NFFT(finalize)(&p_fg_psi);
208 NFFT(finalize)(&p_pre_phi_hut);
209 NFFT(finalize)(&p);
210
211 if (test_ndft)
212 NFFT(free)(swapndft);
213 }
214
accuracy_pre_lin_psi(int d,int N,int M,int n,int m,int K)215 static void accuracy_pre_lin_psi(int d, int N, int M, int n, int m, int K)
216 {
217 int r, NN[d], nn[d];
218 R e;
219 C *swapndft;
220
221 NFFT(plan) p;
222
223 for (r = 0; r < d; r++)
224 {
225 NN[r] = N;
226 nn[r] = n;
227 }
228
229 /* output vector ndft */
230 swapndft = (C*) NFFT(malloc)((size_t)(M) * sizeof(C));
231
232 NFFT(init_guru)(&p, d, NN, M, nn, m,
233 MALLOC_X | MALLOC_F_HAT | MALLOC_F |
234 PRE_PHI_HUT | PRE_LIN_PSI |
235 FFTW_INIT | FFT_OUT_OF_PLACE,
236 FFTW_MEASURE | FFTW_DESTROY_INPUT);
237
238 /** realloc psi */
239 NFFT(free)(p.psi);
240 p.K = K;
241 p.psi = (R*) NFFT(malloc)((size_t)((p.K + 1) * p.d) * sizeof(R));
242
243 /** precomputation can be done before the nodes are initialised */
244 NFFT(precompute_one_psi)(&p);
245
246 /** init pseudo random nodes */
247 NFFT(vrand_shifted_unit_double)(p.x, p.d * p.M_total);
248
249 /** init pseudo random Fourier coefficients */
250 NFFT(vrand_unit_complex)(p.f_hat, p.N_total);
251
252 /** compute exact result */
253 CSWAP(p.f, swapndft);
254 NFFT(trafo_direct)(&p);
255 CSWAP(p.f, swapndft);
256
257 /** NFFT */
258 NFFT(trafo)(&p);
259 e = NFFT(error_l_2_complex)(swapndft, p.f, p.M_total);
260
261 // printf("%d\t%d\t%d\t%d\t%.2e\n",d,N,m,K,e);
262 printf("$%.1" __FES__ "$&\t", e);
263
264 fflush(stdout);
265
266 /** finalise */
267 NFFT(finalize)(&p);
268 NFFT(free)(swapndft);
269 }
270
main(int argc,char ** argv)271 int main(int argc, char **argv)
272 {
273 int l, trial;
274
275 if (argc <= 2)
276 {
277 fprintf(stderr, "flags type first last trials d m\n");
278 return EXIT_FAILURE;
279 }
280
281 if ((test == 0) && (atoi(argv[1]) < 2))
282 {
283 fprintf(stderr, "MEASURE_TIME in infft.h not set\n");
284 return EXIT_FAILURE;
285 }
286
287 fprintf(stderr, "Testing different precomputation schemes for the nfft.\n");
288 fprintf(stderr, "Columns: d, N=M, t_ndft, e_nfft, t_D, t_pre_phi_hut, ");
289 fprintf(stderr, "t_fftw, t_B, t_fg_psi, t_pre_lin_psi, t_pre_fg_psi, ");
290 fprintf(stderr, "t_pre_psi, t_pre_full_psi\n\n");
291
292 int arg2 = atoi(argv[2]);
293 int arg3 = atoi(argv[3]);
294 int arg4 = atoi(argv[4]);
295
296 /* time vs. N=M */
297 if (atoi(argv[1]) == 0)
298 {
299 int d = atoi(argv[5]);
300 int m = atoi(argv[6]);
301
302 for (l = arg2; l <= arg3; l++)
303 {
304 int N = (int)(1U << l);
305 int M = (int)(1U << (d * l));
306 for (trial = 0; trial < arg4; trial++)
307 {
308 time_accuracy(d, N, M, 2 * N, m, 0, 0);
309 }
310 }
311 }
312 else if (atoi(argv[1]) == 1) /* accuracy vs. time */
313 {
314 int d = atoi(argv[5]);
315 int N = atoi(argv[6]);
316 int m;
317
318 for (m = arg2; m <= arg3; m++)
319 {
320 for (trial = 0; trial < arg4; trial++)
321 {
322 time_accuracy(d, N, (int)(LRINT(POW((R)(N), (R)(d)))), 2 * N, m, 1, 1);
323 }
324 }
325 }
326 else if (atoi(argv[1]) == 2) /* accuracy vs. K for linear interpolation, assumes (m+1)|K */
327 {
328 int d = atoi(argv[5]);
329 int N = atoi(argv[6]);
330 int m = atoi(argv[7]);
331
332 printf("$\\log_2(K/(m+1))$&\t");
333
334 for (l = arg2; l < arg3; l++)
335 printf("$%d$&\t", l);
336
337 printf("$%d$\\\\\n", arg3);
338
339 printf("$\\tilde E_2$&\t");
340 for (l = arg2; l <= arg3; l++)
341 {
342 int x = (m + 1) * (int)(1U << l);
343 accuracy_pre_lin_psi(d, N, (int)(LRINT(POW((R)(N), (R)(d)))), 2 * N, m, x);
344 }
345
346 printf("\n");
347 }
348
349
350 return EXIT_SUCCESS;
351 }
352