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