00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include "config.h"
00021
00022 #include <stdlib.h>
00023 #include <math.h>
00024 #ifdef HAVE_COMPLEX_H
00025 #include <complex.h>
00026 #endif
00027
00028 #include "nfft3util.h"
00029 #include "nfft3.h"
00030 #include "infft.h"
00031
00032 void simple_test_nnfft_1d(void)
00033 {
00034 int j,k;
00035 nnfft_plan my_plan;
00037 int N[1];
00038 N[0]=12;
00039
00041 nnfft_init(&my_plan, 1, 3, 19, N);
00042
00044 for(j=0;j<my_plan.M_total;j++)
00045 {
00046 my_plan.x[j]=((double)rand())/((double)RAND_MAX)-0.5;
00047 }
00049 for(j=0;j<my_plan.N_total;j++)
00050 {
00051 my_plan.v[j]=((double)rand())/((double)RAND_MAX)-0.5;
00052 }
00053
00055 if(my_plan.nnfft_flags & PRE_PSI)
00056 nnfft_precompute_psi(&my_plan);
00057
00058 if(my_plan.nnfft_flags & PRE_FULL_PSI)
00059 nnfft_precompute_full_psi(&my_plan);
00060
00061 if(my_plan.nnfft_flags & PRE_LIN_PSI)
00062 nnfft_precompute_lin_psi(&my_plan);
00063
00065 if(my_plan.nnfft_flags & PRE_PHI_HUT)
00066 nnfft_precompute_phi_hut(&my_plan);
00067
00069 for(k=0;k<my_plan.N_total;k++)
00070 my_plan.f_hat[k] = ((double)rand())/((double)RAND_MAX) + _Complex_I*((double)rand())/((double)RAND_MAX);
00071
00072 nfft_vpr_complex(my_plan.f_hat,my_plan.N_total,"given Fourier coefficients, vector f_hat");
00073
00075 nnfft_trafo_direct(&my_plan);
00076 nfft_vpr_complex(my_plan.f,my_plan.M_total,"nndft, vector f");
00077
00079 nnfft_trafo(&my_plan);
00080 nfft_vpr_complex(my_plan.f,my_plan.M_total,"nnfft, vector f");
00081
00083 nnfft_finalize(&my_plan);
00084 }
00085
00086 static void simple_test_adjoint_nnfft_1d(void)
00087 {
00088 int j;
00089 nnfft_plan my_plan;
00091 int N[1];
00092 N[0]=12;
00093
00095 nnfft_init(&my_plan, 1, 20, 33, N);
00096
00098 for(j=0;j<my_plan.M_total;j++)
00099 {
00100 my_plan.x[j]=((double)rand())/((double)RAND_MAX)-0.5;
00101 }
00103 for(j=0;j<my_plan.N_total;j++)
00104 {
00105 my_plan.v[j]=((double)rand())/((double)RAND_MAX)-0.5;
00106 }
00107
00109 if(my_plan.nnfft_flags & PRE_PSI)
00110 nnfft_precompute_psi(&my_plan);
00111
00112 if(my_plan.nnfft_flags & PRE_FULL_PSI)
00113 nnfft_precompute_full_psi(&my_plan);
00114
00115 if(my_plan.nnfft_flags & PRE_LIN_PSI)
00116 nnfft_precompute_lin_psi(&my_plan);
00117
00119 if(my_plan.nnfft_flags & PRE_PHI_HUT)
00120 nnfft_precompute_phi_hut(&my_plan);
00121
00123 for(j=0;j<my_plan.M_total;j++)
00124 my_plan.f[j] = ((double)rand())/((double)RAND_MAX) + _Complex_I*((double)rand())/((double)RAND_MAX);
00125
00126 nfft_vpr_complex(my_plan.f,my_plan.M_total,"given Samples, vector f");
00127
00129 nnfft_adjoint_direct(&my_plan);
00130 nfft_vpr_complex(my_plan.f_hat,my_plan.N_total,"adjoint nndft, vector f_hat");
00131
00133 nnfft_adjoint(&my_plan);
00134 nfft_vpr_complex(my_plan.f_hat,my_plan.N_total,"adjoint nnfft, vector f_hat");
00135
00137 nnfft_finalize(&my_plan);
00138 }
00139
00140 static void simple_test_nnfft_2d(void)
00141 {
00142 int j,k;
00143 nnfft_plan my_plan;
00145 int N[2];
00146 N[0]=12;
00147 N[1]=14;
00148
00150 nnfft_init(&my_plan, 2,12*14,19, N);
00151
00153 for(j=0;j<my_plan.M_total;j++)
00154 {
00155 my_plan.x[2*j]=((double)rand())/((double)RAND_MAX)-0.5;
00156 my_plan.x[2*j+1]=((double)rand())/((double)RAND_MAX)-0.5;
00157 }
00158
00160 for(j=0;j<my_plan.N_total;j++)
00161 {
00162 my_plan.v[2*j]=((double)rand())/((double)RAND_MAX)-0.5;
00163 my_plan.v[2*j+1]=((double)rand())/((double)RAND_MAX)-0.5;
00164 }
00165
00167 if(my_plan.nnfft_flags & PRE_PSI)
00168 nnfft_precompute_psi(&my_plan);
00169
00170 if(my_plan.nnfft_flags & PRE_FULL_PSI)
00171 nnfft_precompute_full_psi(&my_plan);
00172
00173 if(my_plan.nnfft_flags & PRE_LIN_PSI)
00174 nnfft_precompute_lin_psi(&my_plan);
00175
00177 if(my_plan.nnfft_flags & PRE_PHI_HUT)
00178 nnfft_precompute_phi_hut(&my_plan);
00179
00181 for(k=0;k<my_plan.N_total;k++)
00182 my_plan.f_hat[k] = ((double)rand())/((double)RAND_MAX) + _Complex_I*((double)rand())/((double)RAND_MAX);
00183
00184 nfft_vpr_complex(my_plan.f_hat,12,
00185 "given Fourier coefficients, vector f_hat (first 12 entries)");
00186
00188 nnfft_trafo_direct(&my_plan);
00189 nfft_vpr_complex(my_plan.f,my_plan.M_total,"ndft, vector f");
00190
00192 nnfft_trafo(&my_plan);
00193 nfft_vpr_complex(my_plan.f,my_plan.M_total,"nfft, vector f");
00194
00196 nnfft_finalize(&my_plan);
00197 }
00198
00199 static void simple_test_innfft_1d(void)
00200 {
00201 int j,k,l,N=8;
00202 nnfft_plan my_plan;
00203 solver_plan_complex my_iplan;
00206 nnfft_init(&my_plan,1 ,8 ,8 ,&N);
00207
00209 solver_init_advanced_complex(&my_iplan,(nfft_mv_plan_complex*)(&my_plan),CGNR);
00210
00212 for(j=0;j<my_plan.M_total;j++)
00213 my_plan.x[j]=((double)rand())/((double)RAND_MAX)-0.5;
00214
00216 for(k=0;k<my_plan.N_total;k++)
00217 my_plan.v[k]=((double)rand())/((double)RAND_MAX)-0.5;
00218
00220 if(my_plan.nnfft_flags & PRE_PSI)
00221 nnfft_precompute_psi(&my_plan);
00222
00223 if(my_plan.nnfft_flags & PRE_FULL_PSI)
00224 nnfft_precompute_full_psi(&my_plan);
00225
00227 if(my_plan.nnfft_flags & PRE_PHI_HUT)
00228 nnfft_precompute_phi_hut(&my_plan);
00229
00231 for(j=0;j<my_plan.M_total;j++)
00232 my_iplan.y[j] = ((double)rand())/((double)RAND_MAX);
00233
00234 nfft_vpr_complex(my_iplan.y,my_plan.M_total,"given data, vector given_f");
00235
00237 for(k=0;k<my_plan.N_total;k++)
00238 my_iplan.f_hat_iter[k] = 0.0;
00239
00240 nfft_vpr_complex(my_iplan.f_hat_iter,my_plan.N_total,
00241 "approximate solution, vector f_hat_iter");
00242
00244 solver_before_loop_complex(&my_iplan);
00245
00246 for(l=0;l<8;l++)
00247 {
00248 printf("iteration l=%d\n",l);
00249 solver_loop_one_step_complex(&my_iplan);
00250 nfft_vpr_complex(my_iplan.f_hat_iter,my_plan.N_total,
00251 "approximate solution, vector f_hat_iter");
00252
00253 NFFT_SWAP_complex(my_iplan.f_hat_iter,my_plan.f_hat);
00254 nnfft_trafo(&my_plan);
00255 nfft_vpr_complex(my_plan.f,my_plan.M_total,"fitting the data, vector f");
00256 NFFT_SWAP_complex(my_iplan.f_hat_iter,my_plan.f_hat);
00257 }
00258
00259 solver_finalize_complex(&my_iplan);
00260 nnfft_finalize(&my_plan);
00261 }
00262
00263 static void measure_time_nnfft_1d(void)
00264 {
00265 int j,k;
00266 nnfft_plan my_plan;
00267 int my_N,N=4;
00268 double t;
00269 ticks t0, t1;
00270
00271 for(my_N=16; my_N<=16384; my_N*=2)
00272 {
00273 nnfft_init(&my_plan,1,my_N,my_N,&N);
00274
00275 for(j=0;j<my_plan.M_total;j++)
00276 my_plan.x[j]=((double)rand())/((double)RAND_MAX)-0.5;
00277
00278 for(j=0;j<my_plan.N_total;j++)
00279 my_plan.v[j]=((double)rand())/((double)RAND_MAX)-0.5;
00280
00281 if(my_plan.nnfft_flags & PRE_PSI)
00282 nnfft_precompute_psi(&my_plan);
00283
00284 if(my_plan.nnfft_flags & PRE_FULL_PSI)
00285 nnfft_precompute_full_psi(&my_plan);
00286
00287 if(my_plan.nnfft_flags & PRE_PHI_HUT)
00288 nnfft_precompute_phi_hut(&my_plan);
00289
00290 for(k=0;k<my_plan.N_total;k++)
00291 my_plan.f_hat[k] = ((double)rand())/((double)RAND_MAX) + _Complex_I*((double)rand())/((double)RAND_MAX);
00292
00293 t0 = getticks();
00294 nnfft_trafo_direct(&my_plan);
00295 t1 = getticks();
00296 t = nfft_elapsed_seconds(t1,t0);
00297 printf("t_nndft=%e,\t",t);
00298
00299 t0 = getticks();
00300 nnfft_trafo(&my_plan);
00301 t1 = getticks();
00302 t = nfft_elapsed_seconds(t1,t0);
00303 printf("t_nnfft=%e\t",t);
00304
00305 printf("(N=M=%d)\n",my_N);
00306
00307 nnfft_finalize(&my_plan);
00308 }
00309 }
00310
00311 int main(void)
00312 {
00313 system("clear");
00314 printf("1) computing a one dimensional nndft, nnfft\n\n");
00315 simple_test_nnfft_1d();
00316
00317 getc(stdin);
00318
00319 system("clear");
00320 printf("1a) computing a one dimensional adjoint nndft, nnfft\n\n");
00321 simple_test_adjoint_nnfft_1d();
00322
00323 getc(stdin);
00324
00325 system("clear");
00326 printf("2) computing a two dimensional nndft, nnfft\n\n");
00327 simple_test_nnfft_2d();
00328
00329 getc(stdin);
00330
00331 system("clear");
00332 printf("3) computing a one dimensional innfft\n\n");
00333 simple_test_innfft_1d();
00334
00335 getc(stdin);
00336
00337 system("clear");
00338 printf("4) computing times for one dimensional nnfft\n\n");
00339 measure_time_nnfft_1d();
00340
00341 return 1;
00342 }