NFFT  3.3.0
nnfft/simple_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 <stdlib.h>
22 #include <math.h>
23 #include <complex.h>
24 
25 #include "nfft3.h"
26 
28 #define CSWAP(x,y) {double _Complex * NFFT_SWAP_temp__; \
29  NFFT_SWAP_temp__=(x); (x)=(y); (y)=NFFT_SWAP_temp__;}
30 
31 void simple_test_nnfft_1d(void)
32 {
33  int j,k;
34  nnfft_plan my_plan;
36  int N[1];
37  N[0]=10;
38 
40  nnfft_init(&my_plan, 1, 3, 19, N);
41 
43  for(j=0;j<my_plan.M_total;j++)
44  {
45  my_plan.x[j]=((double)rand())/((double)RAND_MAX)-0.5;
46  }
48  for(j=0;j<my_plan.N_total;j++)
49  {
50  my_plan.v[j]=((double)rand())/((double)RAND_MAX)-0.5;
51  }
52 
54 /* if(my_plan.nnfft_flags & PRE_PSI)
55  nnfft_precompute_psi(&my_plan);
56 
57  if(my_plan.nnfft_flags & PRE_FULL_PSI)
58  nnfft_precompute_full_psi(&my_plan);
59  if(my_plan.nnfft_flags & PRE_LIN_PSI)
60  nnfft_precompute_lin_psi(&my_plan);
62 /* if(my_plan.nnfft_flags & PRE_PHI_HUT)
63  nnfft_precompute_phi_hut(&my_plan);
64 */
65 
66 nnfft_precompute_one_psi(&my_plan);
67 
68 
70  for(k=0;k<my_plan.N_total;k++)
71  my_plan.f_hat[k] = ((double)rand())/((double)RAND_MAX) + _Complex_I*((double)rand())/((double)RAND_MAX);
72 
73  nfft_vpr_complex(my_plan.f_hat,my_plan.N_total,"given Fourier coefficients, vector f_hat");
74 
76  nnfft_trafo_direct(&my_plan);
77  nfft_vpr_complex(my_plan.f,my_plan.M_total,"nndft, vector f");
78 
80  nnfft_trafo(&my_plan);
81  nfft_vpr_complex(my_plan.f,my_plan.M_total,"nnfft, vector f");
82 
84  nnfft_finalize(&my_plan);
85 }
86 
87 static void simple_test_adjoint_nnfft_1d(void)
88 {
89  int j;
90  nnfft_plan my_plan;
92  int N[1];
93  N[0]=12;
94 
96  nnfft_init(&my_plan, 1, 20, 33, N);
97 
99  for(j=0;j<my_plan.M_total;j++)
100  {
101  my_plan.x[j]=((double)rand())/((double)RAND_MAX)-0.5;
102  }
104  for(j=0;j<my_plan.N_total;j++)
105  {
106  my_plan.v[j]=((double)rand())/((double)RAND_MAX)-0.5;
107  }
108 
110  if(my_plan.nnfft_flags & PRE_PSI)
111  nnfft_precompute_psi(&my_plan);
112 
113  if(my_plan.nnfft_flags & PRE_FULL_PSI)
114  nnfft_precompute_full_psi(&my_plan);
115 
116  if(my_plan.nnfft_flags & PRE_LIN_PSI)
117  nnfft_precompute_lin_psi(&my_plan);
118 
120  if(my_plan.nnfft_flags & PRE_PHI_HUT)
121  nnfft_precompute_phi_hut(&my_plan);
122 
124  for(j=0;j<my_plan.M_total;j++)
125  my_plan.f[j] = ((double)rand())/((double)RAND_MAX) + _Complex_I*((double)rand())/((double)RAND_MAX);
126 
127  nfft_vpr_complex(my_plan.f,my_plan.M_total,"given Samples, vector f");
128 
130  nnfft_adjoint_direct(&my_plan);
131  nfft_vpr_complex(my_plan.f_hat,my_plan.N_total,"adjoint nndft, vector f_hat");
132 
134  nnfft_adjoint(&my_plan);
135  nfft_vpr_complex(my_plan.f_hat,my_plan.N_total,"adjoint nnfft, vector f_hat");
136 
138  nnfft_finalize(&my_plan);
139 }
140 
141 static void simple_test_nnfft_2d(void)
142 {
143  int j,k;
144  nnfft_plan my_plan;
146  int N[2];
147  N[0]=12;
148  N[1]=14;
149 
151  nnfft_init(&my_plan, 2,12*14,19, N);
152 
154  for(j=0;j<my_plan.M_total;j++)
155  {
156  my_plan.x[2*j]=((double)rand())/((double)RAND_MAX)-0.5;
157  my_plan.x[2*j+1]=((double)rand())/((double)RAND_MAX)-0.5;
158  }
159 
161  for(j=0;j<my_plan.N_total;j++)
162  {
163  my_plan.v[2*j]=((double)rand())/((double)RAND_MAX)-0.5;
164  my_plan.v[2*j+1]=((double)rand())/((double)RAND_MAX)-0.5;
165  }
166 
168  if(my_plan.nnfft_flags & PRE_PSI)
169  nnfft_precompute_psi(&my_plan);
170 
171  if(my_plan.nnfft_flags & PRE_FULL_PSI)
172  nnfft_precompute_full_psi(&my_plan);
173 
174  if(my_plan.nnfft_flags & PRE_LIN_PSI)
175  nnfft_precompute_lin_psi(&my_plan);
176 
178  if(my_plan.nnfft_flags & PRE_PHI_HUT)
179  nnfft_precompute_phi_hut(&my_plan);
180 
182  for(k=0;k<my_plan.N_total;k++)
183  my_plan.f_hat[k] = ((double)rand())/((double)RAND_MAX) + _Complex_I*((double)rand())/((double)RAND_MAX);
184 
185  nfft_vpr_complex(my_plan.f_hat,12,
186  "given Fourier coefficients, vector f_hat (first 12 entries)");
187 
189  nnfft_trafo_direct(&my_plan);
190  nfft_vpr_complex(my_plan.f,my_plan.M_total,"ndft, vector f");
191 
193  nnfft_trafo(&my_plan);
194  nfft_vpr_complex(my_plan.f,my_plan.M_total,"nfft, vector f");
195 
197  nnfft_finalize(&my_plan);
198 }
199 
200 static void simple_test_innfft_1d(void)
201 {
202  int j,k,l,N=8;
203  nnfft_plan my_plan;
204  solver_plan_complex my_iplan;
207  nnfft_init(&my_plan,1 ,8 ,8 ,&N);
208 
210  solver_init_advanced_complex(&my_iplan,(nfft_mv_plan_complex*)(&my_plan),CGNR);
211 
213  for(j=0;j<my_plan.M_total;j++)
214  my_plan.x[j]=((double)rand())/((double)RAND_MAX)-0.5;
215 
217  for(k=0;k<my_plan.N_total;k++)
218  my_plan.v[k]=((double)rand())/((double)RAND_MAX)-0.5;
219 
221  if(my_plan.nnfft_flags & PRE_PSI)
222  nnfft_precompute_psi(&my_plan);
223 
224  if(my_plan.nnfft_flags & PRE_FULL_PSI)
225  nnfft_precompute_full_psi(&my_plan);
226 
228  if(my_plan.nnfft_flags & PRE_PHI_HUT)
229  nnfft_precompute_phi_hut(&my_plan);
230 
232  for(j=0;j<my_plan.M_total;j++)
233  my_iplan.y[j] = ((double)rand())/((double)RAND_MAX);
234 
235  nfft_vpr_complex(my_iplan.y,my_plan.M_total,"given data, vector given_f");
236 
238  for(k=0;k<my_plan.N_total;k++)
239  my_iplan.f_hat_iter[k] = 0.0;
240 
241  nfft_vpr_complex(my_iplan.f_hat_iter,my_plan.N_total,
242  "approximate solution, vector f_hat_iter");
243 
245  solver_before_loop_complex(&my_iplan);
246 
247  for(l=0;l<8;l++)
248  {
249  printf("iteration l=%d\n",l);
250  solver_loop_one_step_complex(&my_iplan);
251  nfft_vpr_complex(my_iplan.f_hat_iter,my_plan.N_total,
252  "approximate solution, vector f_hat_iter");
253 
254  CSWAP(my_iplan.f_hat_iter,my_plan.f_hat);
255  nnfft_trafo(&my_plan);
256  nfft_vpr_complex(my_plan.f,my_plan.M_total,"fitting the data, vector f");
257  CSWAP(my_iplan.f_hat_iter,my_plan.f_hat);
258  }
259 
260  solver_finalize_complex(&my_iplan);
261  nnfft_finalize(&my_plan);
262 }
263 
264 static void measure_time_nnfft_1d(void)
265 {
266  int j,k;
267  nnfft_plan my_plan;
268  int my_N,N=4;
269  double t;
270  double t0, t1;
271 
272  for(my_N=16; my_N<=16384; my_N*=2)
273  {
274  nnfft_init(&my_plan,1,my_N,my_N,&N);
275 
276  for(j=0;j<my_plan.M_total;j++)
277  my_plan.x[j]=((double)rand())/((double)RAND_MAX)-0.5;
278 
279  for(j=0;j<my_plan.N_total;j++)
280  my_plan.v[j]=((double)rand())/((double)RAND_MAX)-0.5;
281 
282  if(my_plan.nnfft_flags & PRE_PSI)
283  nnfft_precompute_psi(&my_plan);
284 
285  if(my_plan.nnfft_flags & PRE_FULL_PSI)
286  nnfft_precompute_full_psi(&my_plan);
287 
288  if(my_plan.nnfft_flags & PRE_PHI_HUT)
289  nnfft_precompute_phi_hut(&my_plan);
290 
291  for(k=0;k<my_plan.N_total;k++)
292  my_plan.f_hat[k] = ((double)rand())/((double)RAND_MAX) + _Complex_I*((double)rand())/((double)RAND_MAX);
293 
294  t0 = nfft_clock_gettime_seconds();
295  nnfft_trafo_direct(&my_plan);
296  t1 = nfft_clock_gettime_seconds();
297  t = t1 - t0;
298  printf("t_nndft=%e,\t",t);
299 
300  t0 = nfft_clock_gettime_seconds();
301  nnfft_trafo(&my_plan);
302  t1 = nfft_clock_gettime_seconds();
303  t = t1 - t0;
304  printf("t_nnfft=%e\t",t);
305 
306  printf("(N=M=%d)\n",my_N);
307 
308  nnfft_finalize(&my_plan);
309  }
310 }
311 
312 int main(void)
313 {
314  system("clear");
315  printf("1) computing a one dimensional nndft, nnfft SUSE\n\n");
316  simple_test_nnfft_1d();
317 
318  /*getc(stdin);
319 
320  system("clear");
321  printf("1a) computing a one dimensional adjoint nndft, nnfft\n\n");
322  simple_test_adjoint_nnfft_1d();
323 
324  getc(stdin);
325 
326  system("clear");
327  printf("2) computing a two dimensional nndft, nnfft\n\n");
328  simple_test_nnfft_2d();
329 
330  getc(stdin);
331 
332  system("clear");
333  printf("3) computing a one dimensional innfft\n\n");
334  simple_test_innfft_1d();
335 
336  getc(stdin);
337 
338  system("clear");
339  printf("4) computing times for one dimensional nnfft\n\n");
340  measure_time_nnfft_1d();
341 */
342  return 1;
343 }
fftw_complex * f_hat
Fourier coefficients.
Definition: nfft3.h:426
void nnfft_precompute_full_psi(nnfft_plan *ths_plan)
computes all entries of B explicitly
Definition: nnfft.c:426
unsigned nnfft_flags
flags for precomputation, malloc
Definition: nfft3.h:426
fftw_complex * f
Samples.
Definition: nfft3.h:426
double * v
nodes (in fourier domain)
Definition: nfft3.h:426
void nnfft_trafo(nnfft_plan *ths_plan)
user routines
Definition: nnfft.c:293
void nnfft_precompute_phi_hut(nnfft_plan *ths_plan)
initialisation of direct transform
Definition: nnfft.c:349
void nfft_vpr_complex(fftw_complex *x, const NFFT_INT n, const char *text)
Print complex vector to standard output.
NFFT_INT M_total
Total number of samples.
Definition: nfft3.h:426
data structure for an NNFFT (nonequispaced in time and frequency fast Fourier transform) plan with do...
Definition: nfft3.h:426
double * x
nodes (in time/spatial domain)
Definition: nfft3.h:426
void nnfft_precompute_lin_psi(nnfft_plan *ths_plan)
create a lookup table
Definition: nnfft.c:369
fftw_complex * y
right hand side, samples
Definition: nfft3.h:786
NFFT_INT N_total
Total number of Fourier coefficients.
Definition: nfft3.h:426
data structure for an inverse NFFT plan with double precision
Definition: nfft3.h:786
#define CSWAP(x, y)
Swap two vectors.
Definition: infft.h:152
fftw_complex * f_hat_iter
iterative solution
Definition: nfft3.h:786