|
|
|
00001 #include "util.h" 00002 #include "nfft3.h" 00003 #include "stdlib.h" 00004 00005 void simple_test_nnfft_1d() 00006 { 00007 int j,k; 00008 nnfft_plan my_plan; 00010 int N[1]; 00011 N[0]=12; 00012 00014 nnfft_init(&my_plan, 1, 12, 19, N); 00015 00017 for(j=0;j<my_plan.M_total;j++) 00018 { 00019 my_plan.x[j]=((double)rand())/RAND_MAX-0.5; 00020 } 00022 for(j=0;j<my_plan.N_total;j++) 00023 { 00024 my_plan.v[j]=((double)rand())/RAND_MAX-0.5; 00025 } 00026 00028 if(my_plan.nnfft_flags & PRE_PSI) 00029 nnfft_precompute_psi(&my_plan); 00030 00031 if(my_plan.nnfft_flags & PRE_FULL_PSI) 00032 nnfft_precompute_full_psi(&my_plan); 00033 00034 if(my_plan.nnfft_flags & PRE_LIN_PSI) 00035 nnfft_precompute_lin_psi(&my_plan); 00036 00038 if(my_plan.nnfft_flags & PRE_PHI_HUT) 00039 nnfft_precompute_phi_hut(&my_plan); 00040 00042 for(k=0;k<my_plan.N_total;k++) 00043 my_plan.f_hat[k] = ((double)rand())/RAND_MAX +I*((double)rand())/RAND_MAX; 00044 00045 nfft_vpr_complex(my_plan.f_hat,my_plan.N_total,"given Fourier coefficients, vector f_hat"); 00046 00048 nndft_trafo(&my_plan); 00049 nfft_vpr_complex(my_plan.f,my_plan.M_total,"nndft, vector f"); 00050 00052 nnfft_trafo(&my_plan); 00053 nfft_vpr_complex(my_plan.f,my_plan.M_total,"nnfft, vector f"); 00054 00056 nnfft_finalize(&my_plan); 00057 } 00058 00059 void simple_test_nnfft_2d() 00060 { 00061 int j,k; 00062 nnfft_plan my_plan; 00064 int N[2]; 00065 N[0]=12; 00066 N[1]=14; 00067 00069 nnfft_init(&my_plan, 2,12*14,19, N); 00070 00072 for(j=0;j<my_plan.M_total;j++) 00073 { 00074 my_plan.x[2*j]=((double)rand())/RAND_MAX-0.5; 00075 my_plan.x[2*j+1]=((double)rand())/RAND_MAX-0.5; 00076 } 00077 00079 for(j=0;j<my_plan.N_total;j++) 00080 { 00081 my_plan.v[2*j]=((double)rand())/RAND_MAX-0.5; 00082 my_plan.v[2*j+1]=((double)rand())/RAND_MAX-0.5; 00083 } 00084 00086 if(my_plan.nnfft_flags & PRE_PSI) 00087 nnfft_precompute_psi(&my_plan); 00088 00089 if(my_plan.nnfft_flags & PRE_FULL_PSI) 00090 nnfft_precompute_full_psi(&my_plan); 00091 00092 if(my_plan.nnfft_flags & PRE_LIN_PSI) 00093 nnfft_precompute_lin_psi(&my_plan); 00094 00096 if(my_plan.nnfft_flags & PRE_PHI_HUT) 00097 nnfft_precompute_phi_hut(&my_plan); 00098 00100 for(k=0;k<my_plan.N_total;k++) 00101 my_plan.f_hat[k] = ((double)rand())/RAND_MAX + I*((double)rand())/RAND_MAX; 00102 00103 nfft_vpr_complex(my_plan.f_hat,12, 00104 "given Fourier coefficients, vector f_hat (first 12 entries)"); 00105 00107 nndft_trafo(&my_plan); 00108 nfft_vpr_complex(my_plan.f,my_plan.M_total,"ndft, vector f"); 00109 00111 nnfft_trafo(&my_plan); 00112 nfft_vpr_complex(my_plan.f,my_plan.M_total,"nfft, vector f"); 00113 00115 nnfft_finalize(&my_plan); 00116 } 00117 00118 void simple_test_innfft_1d() 00119 { 00120 int j,k,l,N=8; 00121 nnfft_plan my_plan; 00122 innfft_plan my_iplan; 00125 nnfft_init(&my_plan,1 ,8 ,8 ,&N); 00126 00128 innfft_init_advanced(&my_iplan,&my_plan,CGNR); 00129 00131 for(j=0;j<my_plan.M_total;j++) 00132 my_plan.x[j]=((double)rand())/RAND_MAX-0.5; 00133 00135 for(j=0;j<my_plan.N_total;j++) 00136 my_plan.v[j]=((double)rand())/RAND_MAX-0.5; 00137 00139 if(my_plan.nnfft_flags & PRE_PSI) 00140 nnfft_precompute_psi(&my_plan); 00141 00142 if(my_plan.nnfft_flags & PRE_FULL_PSI) 00143 nnfft_precompute_full_psi(&my_plan); 00144 00146 if(my_plan.nnfft_flags & PRE_PHI_HUT) 00147 nnfft_precompute_phi_hut(&my_plan); 00148 00150 for(j=0;j<my_plan.M_total;j++) 00151 my_iplan.y[j] = ((double)rand())/RAND_MAX; 00152 00153 nfft_vpr_complex(my_iplan.y,my_plan.M_total,"given data, vector given_f"); 00154 00156 for(k=0;k<my_plan.N_total;k++) 00157 my_iplan.f_hat_iter[k] = 0.0; 00158 00159 nfft_vpr_complex(my_iplan.f_hat_iter,my_plan.N_total, 00160 "approximate solution, vector f_hat_iter"); 00161 00163 innfft_before_loop(&my_iplan); 00164 for(l=0;l<8;l++) 00165 { 00166 printf("iteration l=%d\n",l); 00167 innfft_loop_one_step(&my_iplan); 00168 nfft_vpr_complex(my_iplan.f_hat_iter,my_plan.N_total, 00169 "approximate solution, vector f_hat_iter"); 00170 00171 NFFT_SWAP_complex(my_iplan.f_hat_iter,my_plan.f_hat); 00172 nnfft_trafo(&my_plan); 00173 nfft_vpr_complex(my_plan.f,my_plan.M_total,"fitting the data, vector f"); 00174 NFFT_SWAP_complex(my_iplan.f_hat_iter,my_plan.f_hat); 00175 } 00176 00177 innfft_finalize(&my_iplan); 00178 nnfft_finalize(&my_plan); 00179 } 00180 00181 void measure_time_nnfft_1d() 00182 { 00183 int j,k; 00184 nnfft_plan my_plan; 00185 int my_N,N=4; 00186 double t; 00187 00188 for(my_N=16; my_N<=16384; my_N*=2) 00189 { 00190 nnfft_init(&my_plan,1,my_N,my_N,&N); 00191 00192 for(j=0;j<my_plan.M_total;j++) 00193 my_plan.x[j]=((double)rand())/RAND_MAX-0.5; 00194 00195 for(j=0;j<my_plan.N_total;j++) 00196 my_plan.v[j]=((double)rand())/RAND_MAX-0.5; 00197 00198 if(my_plan.nnfft_flags & PRE_PSI) 00199 nnfft_precompute_psi(&my_plan); 00200 00201 if(my_plan.nnfft_flags & PRE_FULL_PSI) 00202 nnfft_precompute_full_psi(&my_plan); 00203 00204 if(my_plan.nnfft_flags & PRE_PHI_HUT) 00205 nnfft_precompute_phi_hut(&my_plan); 00206 00207 for(k=0;k<my_plan.N_total;k++) 00208 my_plan.f_hat[k] = ((double)rand())/RAND_MAX + I*((double)rand())/RAND_MAX; 00209 00210 t=nfft_second(); 00211 nndft_trafo(&my_plan); 00212 t=nfft_second()-t; 00213 printf("t_nndft=%e,\t",t); 00214 00215 t=nfft_second(); 00216 nnfft_trafo(&my_plan); 00217 t=nfft_second()-t; 00218 printf("t_nnfft=%e\t",t); 00219 00220 printf("(N=M=%d)\n",my_N); 00221 00222 nnfft_finalize(&my_plan); 00223 } 00224 } 00225 00226 int main() 00227 { 00228 system("clear"); 00229 printf("1) computing an one dimensional nndft, nnfft\n\n"); 00230 simple_test_nnfft_1d(); 00231 00232 getc(stdin); 00233 00234 system("clear"); 00235 printf("2) computing a two dimensional nndft, nfft\n\n"); 00236 simple_test_nnfft_2d(); 00237 00238 getc(stdin); 00239 00240 system("clear"); 00241 printf("3) computing an one dimensional infft\n\n"); 00242 simple_test_innfft_1d(); 00243 00244 getc(stdin); 00245 00246 system("clear"); 00247 printf("4) computing times for one dimensional nnfft\n\n"); 00248 measure_time_nnfft_1d(); 00249 00250 return 1; 00251 }