NFFT Logo 3.2.2
flags.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: flags.c 3896 2012-10-10 12:19:26Z tovo $ */
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   //  printf("%d\t%d\t%d\t%d\t%.2e\n",d,N,m,K,e);
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   /* time vs. N=M */
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   /* accuracy vs. time */
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   /* accuracy vs. K for linear interpolation, assumes (m+1)|K */
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 }

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409