NFFT Logo 3.2.2
taylor_nfft.c
Go to the documentation of this file.
00001 /*
00002  * Copyright (c) 2002, 2012 Jens Keiner, Stefan Kunis, Daniel Potts
00003  *
00004  * This program is free software; you can redistribute it and/or modify it under
00005  * the terms of the GNU General Public License as published by the Free Software
00006  * Foundation; either version 2 of the License, or (at your option) any later
00007  * version.
00008  *
00009  * This program is distributed in the hope that it will be useful, but WITHOUT
00010  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00011  * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
00012  * details.
00013  *
00014  * You should have received a copy of the GNU General Public License along with
00015  * this program; if not, write to the Free Software Foundation, Inc., 51
00016  * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
00017  */
00018 
00019 /* $Id: taylor_nfft.c 3896 2012-10-10 12:19:26Z tovo $ */
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   /* Note: no nfft precomputation! */
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   /* time vs. N=M */
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   /* error vs. m */
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 }

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409