NFFT  3.3.0
ndft_fast.c
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2002, 2015 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 /* $Id$ */
20 
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 static void ndft_horner_trafo(NFFT(plan) *ths)
41 {
42  INT j, k;
43  C *f_hat_k, *f_j;
44  C exp_omega_0;
45 
46  for (j = 0, f_j = ths->f; j < ths->M_total; j++, f_j++)
47  (*f_j) = K(0.0);
48 
49  for (j = 0, f_j = ths->f; j < ths->M_total; j++, f_j++)
50  {
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++)
53  {
54  (*f_j) += (*f_hat_k);
55  (*f_j) *= exp_omega_0;
56  }
57  (*f_j) *= CEXP(-KPI * II * (R)(ths->N[0]) * ths->x[j]);
58  }
59 } /* ndft_horner_trafo */
60 
61 static void ndft_pre_full_trafo(NFFT(plan) *ths, C *A)
62 {
63  INT j, k;
64  C *f_hat_k, *f_j;
65  C *A_local;
66 
67  for (j = 0, f_j = ths->f; j < ths->M_total; j++, f_j++)
68  (*f_j) = K(0.0);
69 
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);
73 } /* ndft_pre_full_trafo */
74 
75 static void ndft_pre_full_init(NFFT(plan) *ths, C *A)
76 {
77  INT j, k;
78  C *f_hat_k, *f_j, *A_local;
79 
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++)
82  (*A_local) = CEXP(
83  -K2PI * II * ((R) (k) - (R) (ths->N[0]) / K(2.0)) * ths->x[j]);
84 
85 } /* ndft_pre_full_init */
86 
87 static void ndft_time(int N, int M, unsigned test_ndft, unsigned test_pre_full)
88 {
89  int r;
90  R t, t_ndft, t_horner, t_pre_full, t_nfft;
91  C *A = NULL;
92  ticks t0, t1;
93 
94  NFFT(plan) np;
95 
96  printf("%d\t%d\t", N, M);
97 
98  NFFT(init_1d)(&np, N, M);
99 
100  /* Initialize pseudo random nodes. */
101  NFFT(vrand_shifted_unit_double)(np.x, np.M_total);
102 
103  if (test_pre_full)
104  {
105  A = (C*) NFFT(malloc)((size_t)(N * M) * sizeof(R));
106  ndft_pre_full_init(&np, A);
107  }
108 
109  /* Initialize pseudo random Fourier coefficients. */
110  NFFT(vrand_unit_complex)(np.f_hat, np.N_total);
111 
112  /* NDFT */
113  if (test_ndft)
114  {
115  t_ndft = K(0.0);
116  r = 0;
117  while (t_ndft < K(0.1))
118  {
119  r++;
120  t0 = getticks();
121  NFFT(trafo_direct)(&np);
122  t1 = getticks();
123  t = NFFT(elapsed_seconds)(t1, t0);
124  t_ndft += t;
125  }
126  t_ndft /= (R) (r);
127 
128  printf("%.2" __FES__ "\t", t_ndft);
129  }
130  else
131  printf("N/A\t\t");
132 
133  /* Horner NDFT */
134  t_horner = K(0.0);
135  r = 0;
136  while (t_horner < K(0.1))
137  {
138  r++;
139  t0 = getticks();
140  ndft_horner_trafo(&np);
141  t1 = getticks();
142  t = NFFT(elapsed_seconds)(t1, t0);
143  t_horner += t;
144  }
145  t_horner /= (R)(r);
146 
147  printf("%.2" __FES__ "\t", t_horner);
148 
149  /* Fully precomputed NDFT */
150  if (test_pre_full)
151  {
152  t_pre_full = K(0.0);
153  r = 0;
154  while (t_pre_full < K(0.1))
155  {
156  r++;
157  t0 = getticks();
158  ndft_pre_full_trafo(&np, A);
159  t1 = getticks();
160  t = NFFT(elapsed_seconds)(t1, t0);
161  t_pre_full += t;
162  }
163  t_pre_full /= (R)(r);
164 
165  printf("%.2" __FES__ "\t", t_pre_full);
166  }
167  else
168  printf("N/A\t\t");
169 
170  t_nfft = K(0.0);
171  r = 0;
172  while (t_nfft < K(0.1))
173  {
174  r++;
175  t0 = getticks();
176  NFFT(trafo)(&np);
177  t1 = getticks();
178  t = NFFT(elapsed_seconds)(t1, t0);
179  t_nfft += t;
180  }
181  t_nfft /= (R)(r);
182 
183  printf("%.2" __FES__ "\n", t_nfft);
184 
185  fflush(stdout);
186 
187  if (test_pre_full)
188  NFFT(free)(A);
189 
190  NFFT(finalize)(&np);
191 }
192 
193 int main(int argc, char **argv)
194 {
195  int l, trial;
196 
197  if (argc < 4)
198  {
199  fprintf(stderr, "ndft_fast type first last trials\n");
200  return EXIT_FAILURE;
201  }
202  else
203  {
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");
209 
210  /* time vs. N=M */
211  if (atoi(argv[1]) == 0)
212  {
213  for (l = arg2; l <= arg3; l++)
214  {
215  for (trial = 0; trial < arg4; trial++)
216  {
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);
220  }
221  }
222  }
223  }
224 
225  return EXIT_SUCCESS;
226 }