26 #define CSWAP(x,y) {double _Complex * NFFT_SWAP_temp__; \
27 NFFT_SWAP_temp__=(x); (x)=(y); (y)=NFFT_SWAP_temp__;}
29 void simple_test_nnfft_1d(
void)
38 nnfft_init(&my_plan, 1, 3, 19, N);
41 for(j=0;j<my_plan.M_total;j++)
43 my_plan.x[j]=((double)rand())/((
double)RAND_MAX)-0.5;
46 for(j=0;j<my_plan.N_total;j++)
48 my_plan.v[j]=((double)rand())/((
double)RAND_MAX)-0.5;
64 nnfft_precompute_one_psi(&my_plan);
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);
71 nfft_vpr_complex(my_plan.f_hat,my_plan.N_total,
"given Fourier coefficients, vector f_hat");
74 nnfft_trafo_direct(&my_plan);
75 nfft_vpr_complex(my_plan.f,my_plan.M_total,
"nndft, vector f");
78 nnfft_trafo(&my_plan);
79 nfft_vpr_complex(my_plan.f,my_plan.M_total,
"nnfft, vector f");
82 nnfft_finalize(&my_plan);
85 static void simple_test_adjoint_nnfft_1d(
void)
94 nnfft_init(&my_plan, 1, 20, 33, N);
97 for(j=0;j<my_plan.M_total;j++)
99 my_plan.x[j]=((double)rand())/((
double)RAND_MAX)-0.5;
102 for(j=0;j<my_plan.N_total;j++)
104 my_plan.v[j]=((double)rand())/((
double)RAND_MAX)-0.5;
108 if(my_plan.nnfft_flags & PRE_PSI)
109 nnfft_precompute_psi(&my_plan);
111 if(my_plan.nnfft_flags & PRE_FULL_PSI)
112 nnfft_precompute_full_psi(&my_plan);
114 if(my_plan.nnfft_flags & PRE_LIN_PSI)
115 nnfft_precompute_lin_psi(&my_plan);
118 if(my_plan.nnfft_flags & PRE_PHI_HUT)
119 nnfft_precompute_phi_hut(&my_plan);
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);
125 nfft_vpr_complex(my_plan.f,my_plan.M_total,
"given Samples, vector f");
128 nnfft_adjoint_direct(&my_plan);
129 nfft_vpr_complex(my_plan.f_hat,my_plan.N_total,
"adjoint nndft, vector f_hat");
132 nnfft_adjoint(&my_plan);
133 nfft_vpr_complex(my_plan.f_hat,my_plan.N_total,
"adjoint nnfft, vector f_hat");
136 nnfft_finalize(&my_plan);
139 static void simple_test_nnfft_2d(
void)
149 nnfft_init(&my_plan, 2,12*14,19, N);
152 for(j=0;j<my_plan.M_total;j++)
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;
159 for(j=0;j<my_plan.N_total;j++)
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;
166 if(my_plan.nnfft_flags & PRE_PSI)
167 nnfft_precompute_psi(&my_plan);
169 if(my_plan.nnfft_flags & PRE_FULL_PSI)
170 nnfft_precompute_full_psi(&my_plan);
172 if(my_plan.nnfft_flags & PRE_LIN_PSI)
173 nnfft_precompute_lin_psi(&my_plan);
176 if(my_plan.nnfft_flags & PRE_PHI_HUT)
177 nnfft_precompute_phi_hut(&my_plan);
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);
183 nfft_vpr_complex(my_plan.f_hat,12,
184 "given Fourier coefficients, vector f_hat (first 12 entries)");
187 nnfft_trafo_direct(&my_plan);
188 nfft_vpr_complex(my_plan.f,my_plan.M_total,
"ndft, vector f");
191 nnfft_trafo(&my_plan);
192 nfft_vpr_complex(my_plan.f,my_plan.M_total,
"nfft, vector f");
195 nnfft_finalize(&my_plan);
198 static void simple_test_innfft_1d(
void)
202 solver_plan_complex my_iplan;
205 nnfft_init(&my_plan,1 ,8 ,8 ,&N);
208 solver_init_advanced_complex(&my_iplan,(nfft_mv_plan_complex*)(&my_plan),CGNR);
211 for(j=0;j<my_plan.M_total;j++)
212 my_plan.x[j]=((
double)rand())/((double)RAND_MAX)-0.5;
215 for(k=0;k<my_plan.N_total;k++)
216 my_plan.v[k]=((
double)rand())/((
double)RAND_MAX)-0.5;
219 if(my_plan.nnfft_flags & PRE_PSI)
220 nnfft_precompute_psi(&my_plan);
222 if(my_plan.nnfft_flags & PRE_FULL_PSI)
223 nnfft_precompute_full_psi(&my_plan);
226 if(my_plan.nnfft_flags & PRE_PHI_HUT)
227 nnfft_precompute_phi_hut(&my_plan);
230 for(j=0;j<my_plan.M_total;j++)
231 my_iplan.y[j] = ((
double)rand())/((double)RAND_MAX);
233 nfft_vpr_complex(my_iplan.y,my_plan.M_total,
"given data, vector given_f");
236 for(k=0;k<my_plan.N_total;k++)
237 my_iplan.f_hat_iter[k] = 0.0;
239 nfft_vpr_complex(my_iplan.f_hat_iter,my_plan.N_total,
240 "approximate solution, vector f_hat_iter");
243 solver_before_loop_complex(&my_iplan);
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");
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);
258 solver_finalize_complex(&my_iplan);
259 nnfft_finalize(&my_plan);
262 static void measure_time_nnfft_1d(
void)
270 for(my_N=16; my_N<=16384; my_N*=2)
272 nnfft_init(&my_plan,1,my_N,my_N,&N);
274 for(j=0;j<my_plan.M_total;j++)
275 my_plan.x[j]=((
double)rand())/((double)RAND_MAX)-0.5;
277 for(j=0;j<my_plan.N_total;j++)
278 my_plan.v[j]=((
double)rand())/((
double)RAND_MAX)-0.5;
280 if(my_plan.nnfft_flags & PRE_PSI)
281 nnfft_precompute_psi(&my_plan);
283 if(my_plan.nnfft_flags & PRE_FULL_PSI)
284 nnfft_precompute_full_psi(&my_plan);
286 if(my_plan.nnfft_flags & PRE_PHI_HUT)
287 nnfft_precompute_phi_hut(&my_plan);
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);
292 t0 = nfft_clock_gettime_seconds();
293 nnfft_trafo_direct(&my_plan);
294 t1 = nfft_clock_gettime_seconds();
296 printf(
"t_nndft=%e,\t",t);
298 t0 = nfft_clock_gettime_seconds();
299 nnfft_trafo(&my_plan);
300 t1 = nfft_clock_gettime_seconds();
302 printf(
"t_nnfft=%e\t",t);
304 printf(
"(N=M=%d)\n",my_N);
306 nnfft_finalize(&my_plan);
313 printf(
"1) computing a one dimensional nndft, nnfft\n\n");
314 simple_test_nnfft_1d();