NFFT Logo 3.2.2
nfsoft/simple_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: simple_test.c 3896 2012-10-10 12:19:26Z tovo $ */
00020 
00021 #include "config.h"
00022 
00023 /* Include standard C headers. */
00024 #include <stdio.h>
00025 #include <math.h>
00026 #include <string.h>
00027 #include <stdlib.h>
00028 #ifdef HAVE_COMPLEX_H
00029 #include <complex.h>
00030 #endif
00031 
00032 /* Include NFFT 3 utilities headers. */
00033 #include "nfft3util.h"
00034 /* Include NFFT3 library header. */
00035 #include "nfft3.h"
00036 #include "infft.h"
00037 
00038 static void simple_test_nfsoft(int bw, int M)
00039 {
00040   nfsoft_plan plan_nfsoft; 
00041   nfsoft_plan plan_ndsoft; 
00043   ticks t0, t1;
00044   int j; 
00045   int k, m; 
00046   double d1, d2, d3; 
00047   double time, error; 
00048   unsigned int flags = NFSOFT_MALLOC_X | NFSOFT_MALLOC_F | NFSOFT_MALLOC_F_HAT; 
00051   k = 1000; 
00052   m = 5; 
00061   nfsoft_init_guru(&plan_ndsoft, bw, M, flags | NFSOFT_USE_NDFT
00062       | NFSOFT_USE_DPT, PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT
00063       | MALLOC_F | FFTW_INIT | FFT_OUT_OF_PLACE, m, k);
00064 
00065   nfsoft_init_guru(&plan_nfsoft, bw, M, flags, PRE_PHI_HUT | PRE_PSI | MALLOC_X
00066       | MALLOC_F_HAT | MALLOC_F | FFTW_INIT | FFT_OUT_OF_PLACE, m, k);
00067 
00069   for (j = 0; j < plan_nfsoft.M_total; j++)
00070   {
00071     d1 = ((R) rand()) / RAND_MAX - 0.5;
00072     d2 = 0.5 * ((R) rand()) / RAND_MAX;
00073     d3 = ((R) rand()) / RAND_MAX - 0.5;
00074 
00075     plan_nfsoft.x[3* j ] = d1; 
00076     plan_nfsoft.x[3* j + 1] = d2; 
00077     plan_nfsoft.x[3* j + 2] = d3; 
00079     plan_ndsoft.x[3* j ] = d1; 
00080     plan_ndsoft.x[3* j + 1] = d2; 
00081     plan_ndsoft.x[3* j + 2] = d3; 
00082   }
00083 
00085   for (j = 0; j < (bw + 1) * (4* (bw +1)*(bw+1)-1)/3;j++)
00086   {
00087     d1=((R)rand())/RAND_MAX - 0.5;
00088     d2=((R)rand())/RAND_MAX - 0.5;
00089     plan_nfsoft.f_hat[j]=d1 + I*d2;
00090     plan_ndsoft.f_hat[j]=d1 + I*d2;
00091   }
00092 
00093   if ((bw+1)*(4*(bw+1)*(bw+1)-1)/3<=20)
00094   nfft_vpr_complex(plan_nfsoft.f_hat,(bw+1)*(4*(bw+1)*(bw+1)-1)/3,"randomly generated SO(3) Fourier coefficients");
00095   else
00096   nfft_vpr_complex(plan_ndsoft.f_hat,20,"1st-20th randomly generated SO(3) Fourier coefficient");
00097 
00098   printf("\n---------------------------------------------\n");
00099 
00101   nfsoft_precompute(&plan_nfsoft);
00102   nfsoft_precompute(&plan_ndsoft);
00103 
00104 
00106   t0 = getticks();
00107   nfsoft_trafo(&plan_nfsoft);
00108   t1 = getticks();
00109   time = nfft_elapsed_seconds(t1,t0);
00110   if (plan_nfsoft.M_total<=20)
00111     nfft_vpr_complex(plan_nfsoft.f,plan_nfsoft.M_total,"NFSOFT, function samples");
00112   else
00113     nfft_vpr_complex(plan_nfsoft.f,20,"NFSOFT, 1st-20th function sample");
00114   printf(" computed in %11le seconds\n",time);
00115 
00117   t0 = getticks();
00118   nfsoft_trafo(&plan_ndsoft);
00119   t1 = getticks();
00120 time = nfft_elapsed_seconds(t1,t0);
00121   if (plan_ndsoft.M_total<=20)
00122     nfft_vpr_complex(plan_ndsoft.f,plan_ndsoft.M_total,"NDSOFT, function samples");
00123   else
00124     nfft_vpr_complex(plan_ndsoft.f,20,"NDSOFT, 1st-20th function sample");
00125   printf(" computed in %11le seconds\n",time);
00126 
00128   error= X(error_l_infty_complex)(plan_ndsoft.f,plan_nfsoft.f, plan_nfsoft.M_total);
00129   printf("\n The NFSOFT of bandwidth=%d for %d rotations has infty-error %11le \n",bw, M,error);
00130 
00131   printf("\n---------------------------------------------\n");
00132 
00133   plan_nfsoft.f[0]=1.0;
00134   plan_ndsoft.f[0]=1.0;
00135   nfft_vpr_complex(plan_ndsoft.f,plan_ndsoft.M_total, "function samples to start adjoint trafo");
00136 
00138   t0 = getticks();
00139   nfsoft_adjoint(&plan_nfsoft);
00140   t1 = getticks();
00141 time = nfft_elapsed_seconds(t1,t0);
00142   if ((bw+1)*(4*(bw+1)*(bw+1)-1)/3<=20)
00143      nfft_vpr_complex(plan_nfsoft.f_hat,(bw+1)*(4*(bw+1)*(bw+1)-1)/3,"SO(3) Fourier coefficients");
00144   else
00145      nfft_vpr_complex(plan_nfsoft.f_hat,20,"adjoint NFSOFT, 1st-20th Fourier coefficient");
00146   printf(" computed in %11le seconds\n",time);
00147 
00149   t0 = getticks();
00150   nfsoft_adjoint(&plan_ndsoft);
00151   t1 = getticks();
00152 time = nfft_elapsed_seconds(t1,t0);
00153   if ((bw+1)*(4*(bw+1)*(bw+1)-1)/3<=20)
00154   nfft_vpr_complex(plan_ndsoft.f_hat,(bw+1)*(4*(bw+1)*(bw+1)-1)/3,"SO(3) Fourier coefficients");
00155   else
00156     nfft_vpr_complex(plan_ndsoft.f_hat,20,"adjoint NDSOFT, 1st-20th Fourier coefficient");
00157   printf(" computed in %11le seconds\n",time);
00158 
00159 
00161   error=X(error_l_infty_complex)(plan_ndsoft.f_hat,plan_nfsoft.f_hat, (bw+1)*(4*(bw+1)*(bw+1)-1)/3);
00162   printf("\n The adjoint NFSOFT of bandwidth=%d for %d rotations has infty-error %11le \n",bw, M,error);
00163 
00164   printf("\n---------------------------------------------\n");
00165 
00167   nfsoft_finalize(&plan_ndsoft);
00168   nfsoft_finalize(&plan_nfsoft);
00169 }
00170 
00187 int main(int argc, char **argv)
00188 {
00189   int N; 
00190   int M; 
00192   if (argc < 2)
00193   {
00194     printf(
00195         "This test programm computes the NFSOFT with maximum polynomial degree N at M input rotations\n");
00196     printf("Usage: simple_test N M \n");
00197     printf("e.g.: simple_test 8 64\n");
00198     exit(0);
00199   }
00200 
00202   N = atoi(argv[1]);
00203   M = atoi(argv[2]);
00204 
00205   printf(
00206       "computing an NDSOFT, an NFSOFT, an adjoint NDSOFT, and an adjoint NFSOFT\n\n");
00207 
00208   simple_test_nfsoft(N, M);
00209 
00210 
00211 
00212   /* Exit the program. */
00213   return EXIT_SUCCESS;
00214 
00215 }

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409