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