00001
00009 #include <stdio.h>
00010 #include <math.h>
00011 #include <string.h>
00012 #include <stdlib.h>
00013 #include "util.h"
00014 #include "nfft3.h"
00015 #include "options.h"
00016
00017 #ifdef GAUSSIAN
00018 unsigned test_fg=1;
00019 #else
00020 unsigned test_fg=0;
00021 #endif
00022
00023 #ifdef MEASURE_TIME_FFTW
00024 unsigned test_fftw=1;
00025 #else
00026 unsigned test_fftw=0;
00027 #endif
00028
00029 #ifdef MEASURE_TIME
00030 unsigned test=1;
00031 #else
00032 unsigned test=0;
00033 #endif
00034
00035 void flags_cp(nfft_plan *dst, nfft_plan *src)
00036 {
00037 dst->x=src->x;
00038 dst->f_hat=src->f_hat;
00039 dst->f=src->f;
00040 dst->g1=src->g1;
00041 dst->g2=src->g2;
00042 dst->my_fftw_plan1=src->my_fftw_plan1;
00043 dst->my_fftw_plan2=src->my_fftw_plan2;
00044 }
00045
00046 void time_accuracy(int d, int N, int M, int n, int m, unsigned test_ndft,
00047 unsigned test_pre_full_psi)
00048 {
00049 int r, NN[d], nn[d];
00050 double t_ndft, t, e;
00051 double complex *swapndft;
00052
00053 nfft_plan p;
00054 nfft_plan p_pre_phi_hut;
00055 nfft_plan p_fg_psi;
00056 nfft_plan p_pre_lin_psi;
00057 nfft_plan p_pre_fg_psi;
00058 nfft_plan p_pre_psi;
00059 nfft_plan p_pre_full_psi;
00060
00061 printf("%d\t%d\t", d, N);
00062
00063 for(r=0; r<d; r++)
00064 {
00065 NN[r]=N;
00066 nn[r]=n;
00067 }
00068
00070 if(test_ndft)
00071 swapndft=(double complex*)fftw_malloc(M*sizeof(double complex));
00072
00073 nfft_init_guru(&p, d, NN, M, nn, m,
00074 MALLOC_X| MALLOC_F_HAT| MALLOC_F|
00075 FFTW_INIT| FFT_OUT_OF_PLACE,
00076 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00077
00079 nfft_vrand_shifted_unit_double(p.x, p.d*p.M_total);
00080
00081 nfft_init_guru(&p_pre_phi_hut, d, NN, M, nn, m, PRE_PHI_HUT,0);
00082 flags_cp(&p_pre_phi_hut, &p);
00083 nfft_precompute_one_psi(&p_pre_phi_hut);
00084
00085 if(test_fg)
00086 {
00087 nfft_init_guru(&p_fg_psi, d, NN, M, nn, m, FG_PSI,0);
00088 flags_cp(&p_fg_psi, &p);
00089 nfft_precompute_one_psi(&p_fg_psi);
00090 }
00091
00092 nfft_init_guru(&p_pre_lin_psi, d, NN, M, nn, m, PRE_LIN_PSI,0);
00093 flags_cp(&p_pre_lin_psi, &p);
00094 nfft_precompute_one_psi(&p_pre_lin_psi);
00095
00096 if(test_fg)
00097 {
00098 nfft_init_guru(&p_pre_fg_psi, d, NN, M, nn, m, PRE_FG_PSI,0);
00099 flags_cp(&p_pre_fg_psi, &p);
00100 nfft_precompute_one_psi(&p_pre_fg_psi);
00101 }
00102
00103 nfft_init_guru(&p_pre_psi, d, NN, M, nn, m, PRE_PSI,0);
00104 flags_cp(&p_pre_psi, &p);
00105 nfft_precompute_one_psi(&p_pre_psi);
00106
00107 if(test_pre_full_psi)
00108 {
00109 nfft_init_guru(&p_pre_full_psi, d, NN, M, nn, m, PRE_FULL_PSI,0);
00110 flags_cp(&p_pre_full_psi, &p);
00111 nfft_precompute_one_psi(&p_pre_full_psi);
00112 }
00113
00115 nfft_vrand_unit_complex(p.f_hat, p.N_total);
00116
00118 if(test_ndft)
00119 {
00120 NFFT_SWAP_complex(p.f,swapndft);
00121
00122 t_ndft=0;
00123 r=0;
00124 while(t_ndft<0.01)
00125 {
00126 r++;
00127 t=nfft_second();
00128 ndft_trafo(&p);
00129 t=nfft_second()-t;
00130 t_ndft+=t;
00131 }
00132 t_ndft/=r;
00133
00134 NFFT_SWAP_complex(p.f,swapndft);
00135 }
00136 else
00137 t_ndft=nan("");
00138
00140 nfft_trafo(&p);
00141 nfft_trafo(&p_pre_phi_hut);
00142 if(test_fg)
00143 nfft_trafo(&p_fg_psi);
00144 else
00145 p_fg_psi.MEASURE_TIME_t[2]=nan("");
00146 nfft_trafo(&p_pre_lin_psi);
00147 if(test_fg)
00148 nfft_trafo(&p_pre_fg_psi);
00149 else
00150 p_pre_fg_psi.MEASURE_TIME_t[2]=nan("");
00151 nfft_trafo(&p_pre_psi);
00152 if(test_pre_full_psi)
00153 nfft_trafo(&p_pre_full_psi);
00154 else
00155 p_pre_full_psi.MEASURE_TIME_t[2]=nan("");
00156
00157 if(test_fftw==0)
00158 p.MEASURE_TIME_t[1]=nan("");
00159
00160 if(test_ndft)
00161 e=nfft_error_l_2_complex(swapndft, p.f, p.M_total);
00162 else
00163 e=nan("");
00164
00165 printf("%.2e\t%d\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
00166 t_ndft,
00167 m,
00168 e,
00169 p.MEASURE_TIME_t[0],
00170 p_pre_phi_hut.MEASURE_TIME_t[0],
00171
00172 p.MEASURE_TIME_t[1],
00173
00174 p.MEASURE_TIME_t[2],
00175 p_fg_psi.MEASURE_TIME_t[2],
00176 p_pre_lin_psi.MEASURE_TIME_t[2],
00177 p_pre_fg_psi.MEASURE_TIME_t[2],
00178 p_pre_psi.MEASURE_TIME_t[2],
00179 p_pre_full_psi.MEASURE_TIME_t[2]);
00180
00181 fflush(stdout);
00182
00184 if(test_pre_full_psi)
00185 nfft_finalize(&p_pre_full_psi);
00186 nfft_finalize(&p_pre_psi);
00187 if(test_fg)
00188 nfft_finalize(&p_pre_fg_psi);
00189 nfft_finalize(&p_pre_lin_psi);
00190 if(test_fg)
00191 nfft_finalize(&p_fg_psi);
00192 nfft_finalize(&p_pre_phi_hut);
00193 nfft_finalize(&p);
00194
00195 if(test_ndft)
00196 fftw_free(swapndft);
00197 }
00198
00199 void accuracy_pre_lin_psi(int d, int N, int M, int n, int m, int K)
00200 {
00201 int r, NN[d], nn[d];
00202 double e;
00203 double complex *swapndft;
00204
00205 nfft_plan p;
00206
00207 for(r=0; r<d; r++)
00208 {
00209 NN[r]=N;
00210 nn[r]=n;
00211 }
00212
00214 swapndft=(double complex*)fftw_malloc(M*sizeof(double complex));
00215
00216 nfft_init_guru(&p, d, NN, M, nn, m,
00217 MALLOC_X| MALLOC_F_HAT| MALLOC_F|
00218 PRE_PHI_HUT| PRE_LIN_PSI|
00219 FFTW_INIT| FFT_OUT_OF_PLACE,
00220 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00221
00223 fftw_free(p.psi);
00224 p.K=K;
00225 p.psi=(double*) fftw_malloc((p.K+1)*p.d*sizeof(double));
00226
00228 nfft_precompute_one_psi(&p);
00229
00231 nfft_vrand_shifted_unit_double(p.x, p.d*p.M_total);
00232
00234 nfft_vrand_unit_complex(p.f_hat, p.N_total);
00235
00237 NFFT_SWAP_complex(p.f,swapndft);
00238 ndft_trafo(&p);
00239 NFFT_SWAP_complex(p.f,swapndft);
00240
00242 nfft_trafo(&p);
00243 e=nfft_error_l_2_complex(swapndft, p.f, p.M_total);
00244
00245
00246 printf("$%.1e$&\t",e);
00247
00248 fflush(stdout);
00249
00251 nfft_finalize(&p);
00252 fftw_free(swapndft);
00253 }
00254
00255
00256 int main(int argc,char **argv)
00257 {
00258 int l,m,d,trial,N;
00259
00260 if(argc<=2)
00261 {
00262 fprintf(stderr,"flags type first last trials d m\n");
00263 return -1;
00264 }
00265
00266 if((test==0)&&(atoi(argv[1])<2))
00267 {
00268 fprintf(stderr,"MEASURE_TIME in options.h not set\n");
00269 return -1;
00270 }
00271
00272 fprintf(stderr,"Testing different precomputation schemes for the nfft.\n");
00273 fprintf(stderr,"Columns: d, N=M, t_ndft, e_nfft, t_D, t_pre_phi_hut, ");
00274 fprintf(stderr,"t_fftw, t_B, t_fg_psi, t_pre_lin_psi, t_pre_fg_psi, ");
00275 fprintf(stderr,"t_pre_psi, t_pre_full_psi\n\n");
00276
00277 d=atoi(argv[5]);
00278 m=atoi(argv[6]);
00279
00280
00281 if(atoi(argv[1])==0)
00282 for(l=atoi(argv[2]); l<=atoi(argv[3]); l++)
00283 for(trial=0; trial<atoi(argv[4]); trial++)
00284 {
00285 if(l<=20)
00286 time_accuracy(d, (1U<< l), (1U<< (d*l)), (1U<< (l+1)), m, 0, 0);
00287 else
00288 time_accuracy(d, (1U<< l), (1U<< (d*l)), (1U<< (l+1)), m, 0, 0);
00289 }
00290
00291 d=atoi(argv[5]);
00292 N=atoi(argv[6]);
00293
00294
00295 if(atoi(argv[1])==1)
00296 for(m=atoi(argv[2]); m<=atoi(argv[3]); m++)
00297 for(trial=0; trial<atoi(argv[4]); trial++)
00298 time_accuracy(d, N, (int)pow(N,d), 2*N, m, 1, 1);
00299
00300 d=atoi(argv[5]);
00301 N=atoi(argv[6]);
00302 m=atoi(argv[7]);
00303
00304
00305 if(atoi(argv[1])==2)
00306 {
00307 printf("$\\log_2(K/(m+1))$&\t");
00308 for(l=atoi(argv[2]); l<atoi(argv[3]); l++)
00309 printf("$%d$&\t",l);
00310
00311 printf("$%d$\\\\\n",atoi(argv[3]));
00312
00313 printf("$\\tilde E_2$&\t");
00314 for(l=atoi(argv[2]); l<=atoi(argv[3]); l++)
00315 accuracy_pre_lin_psi(d, N, (int)pow(N,d), 2*N, m, (m+1)*(1U<< l));
00316
00317 printf("\n");
00318 }
00319
00320 return 1;
00321 }