NFFT Logo 3.2.2
ndft_fast.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: ndft_fast.c 3896 2012-10-10 12:19:26Z tovo $ */
00020 
00027 #include "config.h"
00028 
00029 #include <stdio.h>
00030 #include <math.h>
00031 #include <string.h>
00032 #include <stdlib.h>
00033 #ifdef HAVE_COMPLEX_H
00034 #include <complex.h>
00035 #endif
00036 
00037 #include "nfft3util.h"
00038 #include "nfft3.h"
00039 #include "infft.h"
00040 
00041 static void ndft_horner_trafo(nfft_plan *ths)
00042 {
00043   int j,k;
00044   double _Complex *f_hat_k, *f_j;
00045   double _Complex exp_omega_0;
00046 
00047   for(j=0, f_j=ths->f; j<ths->M_total; j++, f_j++)
00048     (*f_j) =0;
00049 
00050   for(j=0, f_j=ths->f; j<ths->M_total; j++, f_j++)
00051     {
00052       exp_omega_0 = cexp(+2*PI*_Complex_I*ths->x[j]);
00053       for(k=0, f_hat_k= ths->f_hat; k<ths->N[0]; k++, f_hat_k++)
00054         {
00055           (*f_j)+=(*f_hat_k);
00056           (*f_j)*=exp_omega_0;
00057   }
00058       (*f_j)*=cexp(-PI*_Complex_I*ths->N[0]*ths->x[j]);
00059     }
00060 } /* ndft_horner_trafo */
00061 
00062 static void ndft_pre_full_trafo(nfft_plan *ths, double _Complex *A)
00063 {
00064   int j,k;
00065   double _Complex *f_hat_k, *f_j;
00066   double _Complex *A_local;
00067 
00068   for(j=0, f_j=ths->f; j<ths->M_total; j++, f_j++)
00069     (*f_j) =0;
00070 
00071   for(j=0, f_j=ths->f, A_local=A; j<ths->M_total; j++, f_j++)
00072     for(k=0, f_hat_k= ths->f_hat; k<ths->N[0]; k++, f_hat_k++, A_local++)
00073       (*f_j) += (*f_hat_k)*(*A_local);
00074 } /* ndft_pre_full_trafo */
00075 
00076 static void ndft_pre_full_init(nfft_plan *ths, double _Complex *A)
00077 {
00078   int j,k;
00079   double _Complex *f_hat_k, *f_j, *A_local;
00080 
00081   for(j=0, f_j=ths->f, A_local=A; j<ths->M_total; j++, f_j++)
00082     for(k=0, f_hat_k= ths->f_hat; k<ths->N[0]; k++, f_hat_k++, A_local++)
00083       (*A_local) = cexp(-2*PI*_Complex_I*(k-ths->N[0]/2)*ths->x[j]);
00084 
00085 } /* ndft_pre_full_init */
00086 
00087 static void ndft_time(int N, int M, unsigned test_ndft, unsigned test_pre_full)
00088 {
00089   int r;
00090   double t, t_ndft, t_horner, t_pre_full, t_nfft;
00091   double _Complex *A = NULL;
00092   ticks t0, t1;
00093 
00094   nfft_plan np;
00095 
00096   printf("%d\t%d\t",N, M);
00097 
00098   nfft_init_1d(&np, N, M);
00099 
00101   nfft_vrand_shifted_unit_double(np.x, np.M_total);
00102 
00103   if(test_pre_full)
00104    {
00105      A=(double _Complex*)nfft_malloc(N*M*sizeof(double _Complex));
00106      ndft_pre_full_init(&np, A);
00107    }
00108 
00110   nfft_vrand_unit_complex(np.f_hat, np.N_total);
00111 
00113   if(test_ndft)
00114     {
00115       t_ndft=0;
00116       r=0;
00117       while(t_ndft<0.1)
00118         {
00119           r++;
00120           t0 = getticks();
00121           nfft_trafo_direct(&np);
00122           t1 = getticks();
00123           t = nfft_elapsed_seconds(t1,t0);
00124           t_ndft+=t;
00125         }
00126       t_ndft/=r;
00127 
00128       printf("%.2e\t",t_ndft);
00129     }
00130   else
00131     printf("nan\t\t");
00132 
00134   t_horner=0;
00135   r=0;
00136   while(t_horner<0.1)
00137     {
00138       r++;
00139       t0 = getticks();
00140       ndft_horner_trafo(&np);
00141       t1 = getticks();
00142       t = nfft_elapsed_seconds(t1,t0);
00143       t_horner+=t;
00144     }
00145   t_horner/=r;
00146 
00147   printf("%.2e\t", t_horner);
00148 
00150   if(test_pre_full)
00151     {
00152       t_pre_full=0;
00153       r=0;
00154       while(t_pre_full<0.1)
00155         {
00156           r++;
00157           t0 = getticks();
00158           ndft_pre_full_trafo(&np,A);
00159           t1 = getticks();
00160           t = nfft_elapsed_seconds(t1,t0);
00161           t_pre_full+=t;
00162         }
00163       t_pre_full/=r;
00164 
00165       printf("%.2e\t", t_pre_full);
00166     }
00167   else
00168     printf("nan\t\t");
00169 
00170   t_nfft=0;
00171   r=0;
00172   while(t_nfft<0.1)
00173     {
00174       r++;
00175       t0 = getticks();
00176       nfft_trafo(&np);
00177       t1 = getticks();
00178       t = nfft_elapsed_seconds(t1,t0);
00179       t_nfft+=t;
00180     }
00181   t_nfft/=r;
00182 
00183   printf("%.2e\n", t_nfft);
00184 
00185   fflush(stdout);
00186 
00187   if(test_pre_full)
00188     nfft_free(A);
00189 
00190   nfft_finalize(&np);
00191 }
00192 
00193 int main(int argc,char **argv)
00194 {
00195   int l,trial;
00196 
00197   if(argc<=2)
00198     {
00199       fprintf(stderr,"ndft_fast type first last trials\n");
00200       return -1;
00201     }
00202 
00203   fprintf(stderr,"Testing ndft, Horner-like ndft, fully precomputed ndft.\n");
00204   fprintf(stderr,"Columns: N, M, t_ndft, t_horner, t_pre_full, t_nfft\n\n");
00205 
00206   /* time vs. N=M */
00207   if(atoi(argv[1])==0)
00208     {
00209       for(l=atoi(argv[2]); l<=atoi(argv[3]); l++)
00210         for(trial=0; trial<atoi(argv[4]); trial++)
00211           if(l<=13)
00212             ndft_time((1U<< l), (1U<< l), 1, 1);
00213           else
00214             if(l<=15)
00215               ndft_time((1U<< l), (1U<< l), 1, 0);
00216             else
00217               ndft_time((1U<< l), (1U<< l), 0, 0);
00218     }
00219 
00220   return 1;
00221 }

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409