00001
00010 #include <stdio.h>
00011 #include <math.h>
00012 #include <string.h>
00013 #include <stdlib.h>
00014 #include "util.h"
00015 #include "nfft3.h"
00016
00017 typedef struct
00018 {
00019 nfft_plan p;
00021 int *idx0;
00023 double *deltax0;
00024 } taylor_plan;
00025
00037 void taylor_init(taylor_plan *ths, int N, int M, int n, int m)
00038 {
00039
00040 nfft_init_guru((nfft_plan*)ths, 1, &N, M, &n, m,
00041 MALLOC_X| MALLOC_F_HAT| MALLOC_F|
00042 FFTW_INIT| FFT_OUT_OF_PLACE,
00043 FFTW_ESTIMATE| FFTW_PRESERVE_INPUT);
00044
00045 ths->idx0=(int*)fftw_malloc(M*sizeof(int));
00046 ths->deltax0=(double*)fftw_malloc(M*sizeof(double));
00047 }
00048
00056 void taylor_precompute(taylor_plan *ths)
00057 {
00058 int j;
00059
00060 nfft_plan* cths=(nfft_plan*)ths;
00061
00062 for(j=0;j<cths->M_total;j++)
00063 {
00064 ths->idx0[j] = ((int)round((cths->x[j]+0.5)*cths->n[0]) +
00065 cths->n[0]/2)%cths->n[0];
00066 ths->deltax0[j] = cths->x[j] - (round((cths->x[j]+0.5)*cths->n[0]) /
00067 cths->n[0] - 0.5);
00068 }
00069 }
00070
00078 void taylor_finalize(taylor_plan *ths)
00079 {
00080 fftw_free(ths->deltax0);
00081 fftw_free(ths->idx0);
00082
00083 nfft_finalize((nfft_plan*)ths);
00084 }
00085
00096 void taylor_trafo(taylor_plan *ths)
00097 {
00098 int j,k,l,ll;
00099 double complex *f, *f_hat, *g1;
00100 double *deltax;
00101 int *idx;
00102
00103 nfft_plan *cths=(nfft_plan*)ths;
00104
00105 for(j=0, f=cths->f; j<cths->M_total; j++)
00106 *f++ = 0;
00107
00108 for(k=0; k<cths->n_total; k++)
00109 cths->g1[k]=0;
00110
00111 for(k=-cths->N_total/2, g1=cths->g1+cths->n_total-cths->N_total/2,
00112 f_hat=cths->f_hat; k<0; k++)
00113 (*g1++)=cpow( - 2*PI*I*k,cths->m)* (*f_hat++);
00114
00115 cths->g1[0]=cths->f_hat[cths->N_total/2];
00116
00117 for(k=1, g1=cths->g1+1, f_hat=cths->f_hat+cths->N_total/2+1;
00118 k<cths->N_total/2; k++)
00119 (*g1++)=cpow( - 2*PI*I*k,cths->m)* (*f_hat++);
00120
00121 for(l=cths->m-1; l>=0; l--)
00122 {
00123 for(k=-cths->N_total/2, g1=cths->g1+cths->n_total-cths->N_total/2;
00124 k<0; k++)
00125 (*g1++) /= (-2*PI*I*k);
00126
00127 for(k=1, g1=cths->g1+1; k<cths->N_total/2; k++)
00128 (*g1++) /= (-2*PI*I*k);
00129
00130 fftw_execute(cths->my_fftw_plan1);
00131
00132 ll=(l==0?1:l);
00133 for(j=0, f=cths->f, deltax=ths->deltax0, idx=ths->idx0; j<cths->M_total;
00134 j++, f++)
00135 (*f) = ((*f) * (*deltax++) + cths->g2[*idx++]) /ll;
00136 }
00137 }
00138
00152 void taylor_time_accuracy(int N, int M, int n, int m, int n_taylor,
00153 int m_taylor, unsigned test_accuracy)
00154 {
00155 int r;
00156 double t_ndft, t_nfft, t_taylor, t;
00157 double complex *swapndft;
00158
00159 taylor_plan tp;
00160 nfft_plan np;
00161
00162 printf("%d\t%d\t",N, M);
00163
00164 taylor_init(&tp,N,M,n_taylor,m_taylor);
00165
00166 nfft_init_guru(&np, 1, &N, M, &n, m,
00167 PRE_PHI_HUT| PRE_FG_PSI|
00168 FFTW_INIT| FFT_OUT_OF_PLACE,
00169 FFTW_ESTIMATE| FFTW_DESTROY_INPUT);
00170
00172 np.x=tp.p.x;
00173 np.f_hat=tp.p.f_hat;
00174 np.f=tp.p.f;
00175
00177 if(test_accuracy)
00178 swapndft=(double complex*)fftw_malloc(M*sizeof(double complex));
00179
00181 nfft_vrand_shifted_unit_double(np.x, np.M_total);
00182
00184 taylor_precompute(&tp);
00185
00187 if(np.nfft_flags & PRE_ONE_PSI)
00188 nfft_precompute_one_psi(&np);
00189
00191 nfft_vrand_unit_complex(np.f_hat, np.N_total);
00192
00194 if(test_accuracy)
00195 {
00196 NFFT_SWAP_complex(np.f,swapndft);
00197
00198 t_ndft=0;
00199 r=0;
00200 while(t_ndft<0.01)
00201 {
00202 r++;
00203 t=nfft_second();
00204 ndft_trafo(&np);
00205 t=nfft_second()-t;
00206 t_ndft+=t;
00207 }
00208 t_ndft/=r;
00209
00210 NFFT_SWAP_complex(np.f,swapndft);
00211 printf("%.2e\t",t_ndft);
00212 }
00213 else
00214 printf("nan\t\t");
00215
00217 t_nfft=0;
00218 r=0;
00219 while(t_nfft<0.01)
00220 {
00221 r++;
00222 t=nfft_second();
00223 nfft_trafo(&np);
00224 t=nfft_second()-t;
00225 t_nfft+=t;
00226 }
00227 t_nfft/=r;
00228
00229 printf("%.2f\t%d\t%.2e\t",((double)n)/N, m, t_nfft);
00230
00231 if(test_accuracy)
00232 printf("%.2e\t",nfft_error_l_infty_complex(swapndft, np.f, np.M_total));
00233 else
00234 printf("nan\t\t");
00235
00237 t_taylor=0;
00238 r=0;
00239 while(t_taylor<0.01)
00240 {
00241 r++;
00242 t=nfft_second();
00243 taylor_trafo(&tp);
00244 t=nfft_second()-t;
00245 t_taylor+=t;
00246 }
00247 t_taylor/=r;
00248
00249
00250 printf("%.2f\t%d\t%.2e\t",((double)n_taylor)/N,m_taylor,t_taylor);
00251
00252 if(test_accuracy)
00253 printf("%.2e\n",nfft_error_l_infty_complex(swapndft, np.f, np.M_total));
00254 else
00255 printf("nan\t\n");
00256
00257 fflush(stdout);
00258
00260 if(test_accuracy)
00261 fftw_free(swapndft);
00262
00263 nfft_finalize(&np);
00264 taylor_finalize(&tp);
00265 }
00266
00267 int main(int argc,char **argv)
00268 {
00269 int l,m,trial,N;
00270
00271 if(argc<=2)
00272 {
00273 fprintf(stderr,"taylor_nfft type first last trials sigma_nfft sigma_taylor.\n");
00274 return -1;
00275 }
00276
00277 fprintf(stderr,"Testing the Nfft & a Taylor expansion based version.\n\n");
00278 fprintf(stderr,"Columns: N, M, t_ndft, sigma_nfft, m_nfft, t_nfft, e_nfft");
00279 fprintf(stderr,", sigma_taylor, m_taylor, t_taylor, e_taylor\n");
00280
00281
00282 if(atoi(argv[1])==0)
00283 {
00284 fprintf(stderr,"Fixed target accuracy, timings.\n\n");
00285 for(l=atoi(argv[2]); l<=atoi(argv[3]); l++)
00286 for(trial=0; trial<atoi(argv[4]); trial++)
00287 if(l<=10)
00288 taylor_time_accuracy((1U<< l), (1U<< l), (int)(atof(argv[5])*
00289 (1U<< l)), 6, (int)(atof(argv[6])*(1U<< l)),
00290 6, 1);
00291 else
00292 taylor_time_accuracy((1U<< l), (1U<< l), (int)(atof(argv[5])*
00293 (1U<< l)), 6, (int)(atof(argv[6])*(1U<< l)),
00294 6, 0);
00295 }
00296
00297
00298 if(atoi(argv[1])==1)
00299 {
00300 N=atoi(argv[7]);
00301 fprintf(stderr,"Fixed N=M=%d, error vs. m.\n\n",N);
00302 for(m=atoi(argv[2]); m<=atoi(argv[3]); m++)
00303 for(trial=0; trial<atoi(argv[4]); trial++)
00304 taylor_time_accuracy(N,N, (int)(atof(argv[5])*N), m,
00305 (int)(atof(argv[6])*N), m, 1);
00306 }
00307
00308 return 1;
00309 }