NFFT  3.3.0
nsfft_test.c
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 
21 #include <stdio.h>
22 #include <math.h>
23 #include <string.h>
24 #include <stdlib.h>
25 #include <complex.h>
26 
27 #include "nfft3.h"
28 
30 #define CSWAP(x,y) {double _Complex * NFFT_SWAP_temp__; \
31  NFFT_SWAP_temp__=(x); (x)=(y); (y)=NFFT_SWAP_temp__;}
32 
33 static void accuracy_nsfft(int d, int J, int M, int m)
34 {
35  nsfft_plan p;
36  double _Complex *swap_sndft_trafo, *swap_sndft_adjoint;
37 
38  nsfft_init(&p, d, J, M, m, NSDFT);
39 
40  swap_sndft_trafo=(double _Complex*) nfft_malloc(p.M_total*
41  sizeof(double _Complex));
42  swap_sndft_adjoint=(double _Complex*) nfft_malloc(p.N_total*
43  sizeof(double _Complex));
44 
45  nsfft_init_random_nodes_coeffs(&p);
46 
48  nsfft_trafo_direct(&p);
49 
50  CSWAP(swap_sndft_trafo,p.f);
51 
53  nsfft_trafo(&p);
54 
55  printf("%5d\t %+.5E\t",J,
56  nfft_error_l_infty_1_complex(swap_sndft_trafo, p.f, p.M_total,
57  p.f_hat, p.N_total));
58  fflush(stdout);
59 
61 
63  nsfft_adjoint_direct(&p);
64 
65  CSWAP(swap_sndft_adjoint,p.f_hat);
66 
68  nsfft_adjoint(&p);
69 
70  printf("%+.5E\n",
71  nfft_error_l_infty_1_complex(swap_sndft_adjoint, p.f_hat,
72  p.N_total,
73  p.f, p.M_total));
74  fflush(stdout);
75 
76  nfft_free(swap_sndft_adjoint);
77  nfft_free(swap_sndft_trafo);
78 
80  nsfft_finalize(&p);
81 }
82 
83 static void time_nsfft(int d, int J, int M, unsigned test_nsdft, unsigned test_nfft)
84 {
85  int r, N[d], n[d];
86  double t, t_nsdft, t_nfft, t_nsfft;
87  double t0, t1;
88 
89  nsfft_plan p;
90  nfft_plan np;
91 
92  for(r=0;r<d;r++)
93  {
94  N[r]= nfft_exp2i(J+2);
95  n[r]=(3*N[r])/2;
96  /*n[r]=2*N[r];*/
97  }
98 
100  nsfft_init(&p, d, J, M, 4, NSDFT);
101  nsfft_init_random_nodes_coeffs(&p);
102 
103  /* transforms */
104  if(test_nsdft)
105  {
106  t_nsdft=0;
107  r=0;
108  while(t_nsdft<0.1)
109  {
110  r++;
111  t0 = nfft_clock_gettime_seconds();
112  nsfft_trafo_direct(&p);
113  t1 = nfft_clock_gettime_seconds();
114  t = (t1 - t0);
115  t_nsdft+=t;
116  }
117  t_nsdft/=r;
118  }
119  else
120  t_nsdft=nan("");
121 
122  if(test_nfft)
123  {
124  nfft_init_guru(&np,d,N,M,n,6, FG_PSI| MALLOC_F_HAT| MALLOC_F| FFTW_INIT,
125  FFTW_MEASURE);
126  np.x=p.act_nfft_plan->x;
127  if(np.flags & PRE_ONE_PSI)
128  nfft_precompute_one_psi(&np);
129  nsfft_cp(&p, &np);
130 
131  t_nfft=0;
132  r=0;
133  while(t_nfft<0.1)
134  {
135  r++;
136  t0 = nfft_clock_gettime_seconds();
137  nfft_trafo(&np);
138  t1 = nfft_clock_gettime_seconds();
139  t = (t1 - t0);
140  t_nfft+=t;
141  }
142  t_nfft/=r;
143 
144  nfft_finalize(&np);
145  }
146  else
147  {
148  t_nfft=nan("");
149  }
150 
151  t_nsfft=0;
152  r=0;
153  while(t_nsfft<0.1)
154  {
155  r++;
156  t0 = nfft_clock_gettime_seconds();
157  nsfft_trafo(&p);
158  t1 = nfft_clock_gettime_seconds();
159  t = (t1 - t0);
160  t_nsfft+=t;
161  }
162  t_nsfft/=r;
163 
164  printf("%d\t%.2e\t%.2e\t%.2e\n", J, t_nsdft, t_nfft, t_nsfft);
165 
166  fflush(stdout);
167 
169  nsfft_finalize(&p);
170 }
171 
172 
173 int main(int argc,char **argv)
174 {
175  int d, J, M;
176 
177  if(argc<=2)
178  {
179  fprintf(stderr,"nsfft_test type d [first last trials]\n");
180  return -1;
181  }
182 
183  d=atoi(argv[2]);
184  fprintf(stderr,"Testing the nfft on the hyperbolic cross (nsfft).\n");
185 
186  if(atoi(argv[1])==1)
187  {
188  fprintf(stderr,"Testing the accuracy of the nsfft vs. nsdft\n");
189  fprintf(stderr,"Columns: d, E_{1,\\infty}(trafo) E_{1,\\infty}(adjoint)\n\n");
190  for(J=1; J<10; J++)
191  accuracy_nsfft(d, J, 1000, 6);
192  }
193 
194  if(atoi(argv[1])==2)
195  {
196  fprintf(stderr,"Testing the computation time of the nsdft, nfft, and nsfft\n");
197  fprintf(stderr,"Columns: d, J, M, t_nsdft, t_nfft, t_nsfft\n\n");
198  for(J=atoi(argv[3]); J<=atoi(argv[4]); J++)
199  {
200  if(d==2)
201  M=(J+4)*nfft_exp2i(J+1);
202  else
203  M=6*nfft_exp2i(J)*(nfft_exp2i((J+1)/2+1)-1)+nfft_exp2i(3*(J/2+1));
204 
205  if(d*(J+2)<=24)
206  time_nsfft(d, J, M, 1, 1);
207  else
208  if(d*(J+2)<=24)
209  time_nsfft(d, J, M, 0, 1);
210  else
211  time_nsfft(d, J, M, 0, 0);
212  }
213  }
214 
215  return 1;
216 }
NFFT_INT N_total
Total number of Fourier coefficients.
Definition: nfft3.h:477
NFFT_INT M_total
Total number of samples.
Definition: nfft3.h:477
void nfft_free(void *p)
fftw_complex * f
Samples.
Definition: nfft3.h:477
data structure for an NFFT (nonequispaced fast Fourier transform) plan with double precision ...
Definition: nfft3.h:194
void nfft_vrand_unit_complex(fftw_complex *x, const NFFT_INT n)
Inits a vector of random complex numbers in .
void * nfft_malloc(size_t n)
double * x
Nodes in time/spatial domain, size is doubles.
Definition: nfft3.h:194
unsigned flags
Flags for precomputation, (de)allocation, and FFTW usage, default setting is PRE_PHI_HUT | PRE_PSI | ...
Definition: nfft3.h:194
data structure for an NSFFT (nonequispaced sparse fast Fourier transform) plan with double precision ...
Definition: nfft3.h:477
nfft_plan * act_nfft_plan
current nfft block
Definition: nfft3.h:477
#define CSWAP(x, y)
Swap two vectors.
Definition: infft.h:152
fftw_complex * f_hat
Fourier coefficients.
Definition: nfft3.h:477