00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include "config.h"
00021
00022 #include <stdio.h>
00023 #include <math.h>
00024 #include <string.h>
00025 #include <stdlib.h>
00026 #ifdef HAVE_COMPLEX_H
00027 #include <complex.h>
00028 #endif
00029
00030 #include "nfft3util.h"
00031 #include "nfft3.h"
00032 #include "infft.h"
00033
00034 static void accuracy_nsfft(int d, int J, int M, int m)
00035 {
00036 nsfft_plan p;
00037 double _Complex *swap_sndft_trafo, *swap_sndft_adjoint;
00038
00039 nsfft_init(&p, d, J, M, m, NSDFT);
00040
00041 swap_sndft_trafo=(double _Complex*) nfft_malloc(p.M_total*
00042 sizeof(double _Complex));
00043 swap_sndft_adjoint=(double _Complex*) nfft_malloc(p.N_total*
00044 sizeof(double _Complex));
00045
00046 nsfft_init_random_nodes_coeffs(&p);
00047
00049 nsfft_trafo_direct(&p);
00050
00051 NFFT_SWAP_complex(swap_sndft_trafo,p.f);
00052
00054 nsfft_trafo(&p);
00055
00056 printf("%5d\t %+.5E\t",J,
00057 X(error_l_infty_1_complex)(swap_sndft_trafo, p.f, p.M_total,
00058 p.f_hat, p.N_total));
00059 fflush(stdout);
00060
00061 nfft_vrand_unit_complex(p.f, p.M_total);
00062
00064 nsfft_adjoint_direct(&p);
00065
00066 NFFT_SWAP_complex(swap_sndft_adjoint,p.f_hat);
00067
00069 nsfft_adjoint(&p);
00070
00071 printf("%+.5E\n",
00072 X(error_l_infty_1_complex)(swap_sndft_adjoint, p.f_hat,
00073 p.N_total,
00074 p.f, p.M_total));
00075 fflush(stdout);
00076
00077 nfft_free(swap_sndft_adjoint);
00078 nfft_free(swap_sndft_trafo);
00079
00081 nsfft_finalize(&p);
00082 }
00083
00084 static void time_nsfft(int d, int J, int M, unsigned test_nsdft, unsigned test_nfft)
00085 {
00086 int r, N[d], n[d];
00087 double t, t_nsdft, t_nfft, t_nsfft;
00088 ticks t0, t1;
00089
00090 nsfft_plan p;
00091 nfft_plan np;
00092
00093 for(r=0;r<d;r++)
00094 {
00095 N[r]= X(exp2i)(J+2);
00096 n[r]=(3*N[r])/2;
00097
00098 }
00099
00101 nsfft_init(&p, d, J, M, 4, NSDFT);
00102 nsfft_init_random_nodes_coeffs(&p);
00103
00104
00105 if(test_nsdft)
00106 {
00107 t_nsdft=0;
00108 r=0;
00109 while(t_nsdft<0.1)
00110 {
00111 r++;
00112 t0 = getticks();
00113 nsfft_trafo_direct(&p);
00114 t1 = getticks();
00115 t = nfft_elapsed_seconds(t1,t0);
00116 t_nsdft+=t;
00117 }
00118 t_nsdft/=r;
00119 }
00120 else
00121 t_nsdft=nan("");
00122
00123 if(test_nfft)
00124 {
00125 nfft_init_guru(&np,d,N,M,n,6, FG_PSI| MALLOC_F_HAT| MALLOC_F| FFTW_INIT,
00126 FFTW_MEASURE);
00127 np.x=p.act_nfft_plan->x;
00128 if(np.nfft_flags & PRE_ONE_PSI)
00129 nfft_precompute_one_psi(&np);
00130 nsfft_cp(&p, &np);
00131
00132 t_nfft=0;
00133 r=0;
00134 while(t_nfft<0.1)
00135 {
00136 r++;
00137 t0 = getticks();
00138 nfft_trafo(&np);
00139 t1 = getticks();
00140 t = nfft_elapsed_seconds(t1,t0);
00141 t_nfft+=t;
00142 }
00143 t_nfft/=r;
00144
00145 nfft_finalize(&np);
00146 }
00147 else
00148 {
00149 t_nfft=nan("");
00150 }
00151
00152 t_nsfft=0;
00153 r=0;
00154 while(t_nsfft<0.1)
00155 {
00156 r++;
00157 t0 = getticks();
00158 nsfft_trafo(&p);
00159 t1 = getticks();
00160 t = nfft_elapsed_seconds(t1,t0);
00161 t_nsfft+=t;
00162 }
00163 t_nsfft/=r;
00164
00165 printf("%d\t%.2e\t%.2e\t%.2e\n", J, t_nsdft, t_nfft, t_nsfft);
00166
00167 fflush(stdout);
00168
00170 nsfft_finalize(&p);
00171 }
00172
00173
00174 int main(int argc,char **argv)
00175 {
00176 int d, J, M;
00177
00178 if(argc<=2)
00179 {
00180 fprintf(stderr,"nsfft_test type d [first last trials]\n");
00181 return -1;
00182 }
00183
00184 d=atoi(argv[2]);
00185 fprintf(stderr,"Testing the nfft on the hyperbolic cross (nsfft).\n");
00186
00187 if(atoi(argv[1])==1)
00188 {
00189 fprintf(stderr,"Testing the accuracy of the nsfft vs. nsdft\n");
00190 fprintf(stderr,"Columns: d, E_{1,\\infty}(trafo) E_{1,\\infty}(adjoint)\n\n");
00191 for(J=1; J<10; J++)
00192 accuracy_nsfft(d, J, 1000, 6);
00193 }
00194
00195 if(atoi(argv[1])==2)
00196 {
00197 fprintf(stderr,"Testing the computation time of the nsdft, nfft, and nsfft\n");
00198 fprintf(stderr,"Columns: d, J, M, t_nsdft, t_nfft, t_nsfft\n\n");
00199 for(J=atoi(argv[3]); J<=atoi(argv[4]); J++)
00200 {
00201 if(d==2)
00202 M=(J+4)*X(exp2i)(J+1);
00203 else
00204 M=6*X(exp2i)(J)*(X(exp2i)((J+1)/2+1)-1)+X(exp2i)(3*(J/2+1));
00205
00206 if(d*(J+2)<=24)
00207 time_nsfft(d, J, M, 1, 1);
00208 else
00209 if(d*(J+2)<=24)
00210 time_nsfft(d, J, M, 0, 1);
00211 else
00212 time_nsfft(d, J, M, 0, 0);
00213 }
00214 }
00215
00216 return 1;
00217 }