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