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