00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
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 }
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 }
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 }
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
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 }