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