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