00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00030 #include <stdio.h>
00031 #include <math.h>
00032 #include <string.h>
00033 #include <stdlib.h>
00034 #include <complex.h>
00035
00036 #include "nfft3util.h"
00037 #include "nfft3.h"
00038
00039 typedef struct
00040 {
00041 nfft_plan p;
00043 int *idx0;
00045 double *deltax0;
00046 } taylor_plan;
00047
00059 void taylor_init(taylor_plan *ths, int N, int M, int n, int m)
00060 {
00061
00062 nfft_init_guru((nfft_plan*)ths, 1, &N, M, &n, m,
00063 MALLOC_X| MALLOC_F_HAT| MALLOC_F|
00064 FFTW_INIT| FFT_OUT_OF_PLACE,
00065 FFTW_ESTIMATE| FFTW_PRESERVE_INPUT);
00066
00067 ths->idx0=(int*)nfft_malloc(M*sizeof(int));
00068 ths->deltax0=(double*)nfft_malloc(M*sizeof(double));
00069 }
00070
00078 void taylor_precompute(taylor_plan *ths)
00079 {
00080 int j;
00081
00082 nfft_plan* cths=(nfft_plan*)ths;
00083
00084 for(j=0;j<cths->M_total;j++)
00085 {
00086 ths->idx0[j] = ((int)round((cths->x[j]+0.5)*cths->n[0]) +
00087 cths->n[0]/2)%cths->n[0];
00088 ths->deltax0[j] = cths->x[j] - (round((cths->x[j]+0.5)*cths->n[0]) /
00089 cths->n[0] - 0.5);
00090 }
00091 }
00092
00100 void taylor_finalize(taylor_plan *ths)
00101 {
00102 nfft_free(ths->deltax0);
00103 nfft_free(ths->idx0);
00104
00105 nfft_finalize((nfft_plan*)ths);
00106 }
00107
00118 void taylor_trafo(taylor_plan *ths)
00119 {
00120 int j,k,l,ll;
00121 double _Complex *f, *f_hat, *g1;
00122 double *deltax;
00123 int *idx;
00124
00125 nfft_plan *cths=(nfft_plan*)ths;
00126
00127 for(j=0, f=cths->f; j<cths->M_total; j++)
00128 *f++ = 0;
00129
00130 for(k=0; k<cths->n_total; k++)
00131 cths->g1[k]=0;
00132
00133 for(k=-cths->N_total/2, g1=cths->g1+cths->n_total-cths->N_total/2,
00134 f_hat=cths->f_hat; k<0; k++)
00135 (*g1++)=cpow( - 2*PI*_Complex_I*k,cths->m)* (*f_hat++);
00136
00137 cths->g1[0]=cths->f_hat[cths->N_total/2];
00138
00139 for(k=1, g1=cths->g1+1, f_hat=cths->f_hat+cths->N_total/2+1;
00140 k<cths->N_total/2; k++)
00141 (*g1++)=cpow( - 2*PI*_Complex_I*k,cths->m)* (*f_hat++);
00142
00143 for(l=cths->m-1; l>=0; l--)
00144 {
00145 for(k=-cths->N_total/2, g1=cths->g1+cths->n_total-cths->N_total/2;
00146 k<0; k++)
00147 (*g1++) /= (-2*PI*_Complex_I*k);
00148
00149 for(k=1, g1=cths->g1+1; k<cths->N_total/2; k++)
00150 (*g1++) /= (-2*PI*_Complex_I*k);
00151
00152 fftw_execute(cths->my_fftw_plan1);
00153
00154 ll=(l==0?1:l);
00155 for(j=0, f=cths->f, deltax=ths->deltax0, idx=ths->idx0; j<cths->M_total;
00156 j++, f++)
00157 (*f) = ((*f) * (*deltax++) + cths->g2[*idx++]) /ll;
00158 }
00159 }
00160
00174 void taylor_time_accuracy(int N, int M, int n, int m, int n_taylor,
00175 int m_taylor, unsigned test_accuracy)
00176 {
00177 int r;
00178 double t_ndft, t_nfft, t_taylor, t;
00179 double _Complex *swapndft;
00180
00181 taylor_plan tp;
00182 nfft_plan np;
00183
00184 printf("%d\t%d\t",N, M);
00185
00186 taylor_init(&tp,N,M,n_taylor,m_taylor);
00187
00188 nfft_init_guru(&np, 1, &N, M, &n, m,
00189 PRE_PHI_HUT| PRE_FG_PSI|
00190 FFTW_INIT| FFT_OUT_OF_PLACE,
00191 FFTW_ESTIMATE| FFTW_DESTROY_INPUT);
00192
00194 np.x=tp.p.x;
00195 np.f_hat=tp.p.f_hat;
00196 np.f=tp.p.f;
00197
00199 if(test_accuracy)
00200 swapndft=(double _Complex*)nfft_malloc(M*sizeof(double _Complex));
00201
00203 nfft_vrand_shifted_unit_double(np.x, np.M_total);
00204
00206 taylor_precompute(&tp);
00207
00209 if(np.nfft_flags & PRE_ONE_PSI)
00210 nfft_precompute_one_psi(&np);
00211
00213 nfft_vrand_unit_complex(np.f_hat, np.N_total);
00214
00216 if(test_accuracy)
00217 {
00218 NFFT_SWAP_complex(np.f,swapndft);
00219
00220 t_ndft=0;
00221 r=0;
00222 while(t_ndft<0.01)
00223 {
00224 r++;
00225 t=nfft_second();
00226 ndft_trafo(&np);
00227 t=nfft_second()-t;
00228 t_ndft+=t;
00229 }
00230 t_ndft/=r;
00231
00232 NFFT_SWAP_complex(np.f,swapndft);
00233 printf("%.2e\t",t_ndft);
00234 }
00235 else
00236 printf("nan\t\t");
00237
00239 t_nfft=0;
00240 r=0;
00241 while(t_nfft<0.01)
00242 {
00243 r++;
00244 t=nfft_second();
00245 nfft_trafo(&np);
00246 t=nfft_second()-t;
00247 t_nfft+=t;
00248 }
00249 t_nfft/=r;
00250
00251 printf("%.2f\t%d\t%.2e\t",((double)n)/N, m, t_nfft);
00252
00253 if(test_accuracy)
00254 printf("%.2e\t",nfft_error_l_infty_complex(swapndft, np.f, np.M_total));
00255 else
00256 printf("nan\t\t");
00257
00259 t_taylor=0;
00260 r=0;
00261 while(t_taylor<0.01)
00262 {
00263 r++;
00264 t=nfft_second();
00265 taylor_trafo(&tp);
00266 t=nfft_second()-t;
00267 t_taylor+=t;
00268 }
00269 t_taylor/=r;
00270
00271
00272 printf("%.2f\t%d\t%.2e\t",((double)n_taylor)/N,m_taylor,t_taylor);
00273
00274 if(test_accuracy)
00275 printf("%.2e\n",nfft_error_l_infty_complex(swapndft, np.f, np.M_total));
00276 else
00277 printf("nan\t\n");
00278
00279 fflush(stdout);
00280
00282 if(test_accuracy)
00283 nfft_free(swapndft);
00284
00285 nfft_finalize(&np);
00286 taylor_finalize(&tp);
00287 }
00288
00289 int main(int argc,char **argv)
00290 {
00291 int l,m,trial,N;
00292
00293 if(argc<=2)
00294 {
00295 fprintf(stderr,"taylor_nfft type first last trials sigma_nfft sigma_taylor.\n");
00296 return -1;
00297 }
00298
00299 fprintf(stderr,"Testing the Nfft & a Taylor expansion based version.\n\n");
00300 fprintf(stderr,"Columns: N, M, t_ndft, sigma_nfft, m_nfft, t_nfft, e_nfft");
00301 fprintf(stderr,", sigma_taylor, m_taylor, t_taylor, e_taylor\n");
00302
00303
00304 if(atoi(argv[1])==0)
00305 {
00306 fprintf(stderr,"Fixed target accuracy, timings.\n\n");
00307 for(l=atoi(argv[2]); l<=atoi(argv[3]); l++)
00308 for(trial=0; trial<atoi(argv[4]); trial++)
00309 if(l<=10)
00310 taylor_time_accuracy((1U<< l), (1U<< l), (int)(atof(argv[5])*
00311 (1U<< l)), 6, (int)(atof(argv[6])*(1U<< l)),
00312 6, 1);
00313 else
00314 taylor_time_accuracy((1U<< l), (1U<< l), (int)(atof(argv[5])*
00315 (1U<< l)), 6, (int)(atof(argv[6])*(1U<< l)),
00316 6, 0);
00317 }
00318
00319
00320 if(atoi(argv[1])==1)
00321 {
00322 N=atoi(argv[7]);
00323 fprintf(stderr,"Fixed N=M=%d, error vs. m.\n\n",N);
00324 for(m=atoi(argv[2]); m<=atoi(argv[3]); m++)
00325 for(trial=0; trial<atoi(argv[4]); trial++)
00326 taylor_time_accuracy(N,N, (int)(atof(argv[5])*N), m,
00327 (int)(atof(argv[6])*N), m, 1);
00328 }
00329
00330 return 1;
00331 }