NFFT Logo 3.2.2
nsfft_test.c
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: nsfft_test.c 3896 2012-10-10 12:19:26Z tovo $ */
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     /*n[r]=2*N[r];*/
00098   }
00099 
00101   nsfft_init(&p, d, J, M, 4, NSDFT);
00102   nsfft_init_random_nodes_coeffs(&p);
00103 
00104   /* transforms */
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 }

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409