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