NFFT  3.3.1
simple_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 standard C headers. */
00020 #include <stdio.h>
00021 #include <math.h>
00022 #include <string.h>
00023 #include <stdlib.h>
00024 #include <complex.h>
00025 
00026 /* Include NFFT3 library header. */
00027 #include "nfft3.h"
00028 
00029 static void simple_test_nfsoft(int bw, int M)
00030 {
00031   nfsoft_plan plan_nfsoft; 
00032   nfsoft_plan plan_ndsoft; 
00034   double t0, t1;
00035   int j; 
00036   int k, m; 
00037   double d1, d2, d3; 
00038   double time, error; 
00039   unsigned int flags = NFSOFT_MALLOC_X | NFSOFT_MALLOC_F | NFSOFT_MALLOC_F_HAT; 
00042   k = 1000; 
00043   m = 5; 
00052   nfsoft_init_guru(&plan_ndsoft, bw, M, flags | NFSOFT_USE_NDFT
00053       | NFSOFT_USE_DPT, PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT
00054       | MALLOC_F | FFTW_INIT | FFT_OUT_OF_PLACE, m, k);
00055 
00056   nfsoft_init_guru(&plan_nfsoft, bw, M, flags, PRE_PHI_HUT | PRE_PSI | MALLOC_X
00057       | MALLOC_F_HAT | MALLOC_F | FFTW_INIT | FFT_OUT_OF_PLACE, m, k);
00058 
00060   for (j = 0; j < plan_nfsoft.M_total; j++)
00061   {
00062     d1 = ((double) rand()) / RAND_MAX - 0.5;
00063     d2 = 0.5 * ((double) rand()) / RAND_MAX;
00064     d3 = ((double) rand()) / RAND_MAX - 0.5;
00065 
00066     plan_nfsoft.x[3* j ] = d1; 
00067     plan_nfsoft.x[3* j + 1] = d2; 
00068     plan_nfsoft.x[3* j + 2] = d3; 
00070     plan_ndsoft.x[3* j ] = d1; 
00071     plan_ndsoft.x[3* j + 1] = d2; 
00072     plan_ndsoft.x[3* j + 2] = d3; 
00073   }
00074 
00076   for (j = 0; j < (bw + 1) * (4* (bw +1)*(bw+1)-1)/3;j++)
00077   {
00078     d1=((double)rand())/RAND_MAX - 0.5;
00079     d2=((double)rand())/RAND_MAX - 0.5;
00080     plan_nfsoft.f_hat[j]=d1 + I*d2;
00081     plan_ndsoft.f_hat[j]=d1 + I*d2;
00082   }
00083 
00084   if ((bw+1)*(4*(bw+1)*(bw+1)-1)/3<=20)
00085   nfft_vpr_complex(plan_nfsoft.f_hat,(bw+1)*(4*(bw+1)*(bw+1)-1)/3,"randomly generated SO(3) Fourier coefficients");
00086   else
00087   nfft_vpr_complex(plan_ndsoft.f_hat,20,"1st-20th randomly generated SO(3) Fourier coefficient");
00088 
00089   printf("\n---------------------------------------------\n");
00090 
00092   nfsoft_precompute(&plan_nfsoft);
00093   nfsoft_precompute(&plan_ndsoft);
00094 
00095 
00097   t0 = nfft_clock_gettime_seconds();
00098   nfsoft_trafo(&plan_nfsoft);
00099   t1 = nfft_clock_gettime_seconds();
00100   time = t1-t0;
00101   if (plan_nfsoft.M_total<=20)
00102     nfft_vpr_complex(plan_nfsoft.f,plan_nfsoft.M_total,"NFSOFT, function samples");
00103   else
00104     nfft_vpr_complex(plan_nfsoft.f,20,"NFSOFT, 1st-20th function sample");
00105   printf(" computed in %11le seconds\n",time);
00106 
00108   t0 = nfft_clock_gettime_seconds();
00109   nfsoft_trafo(&plan_ndsoft);
00110   t1 = nfft_clock_gettime_seconds();
00111   time = t1-t0;
00112   if (plan_ndsoft.M_total<=20)
00113     nfft_vpr_complex(plan_ndsoft.f,plan_ndsoft.M_total,"NDSOFT, function samples");
00114   else
00115     nfft_vpr_complex(plan_ndsoft.f,20,"NDSOFT, 1st-20th function sample");
00116   printf(" computed in %11le seconds\n",time);
00117 
00119   error= nfft_error_l_infty_complex(plan_ndsoft.f,plan_nfsoft.f, plan_nfsoft.M_total);
00120   printf("\n The NFSOFT of bandwidth=%d for %d rotations has infty-error %11le \n",bw, M,error);
00121 
00122   printf("\n---------------------------------------------\n");
00123 
00124   plan_nfsoft.f[0]=1.0;
00125   plan_ndsoft.f[0]=1.0;
00126   nfft_vpr_complex(plan_ndsoft.f,plan_ndsoft.M_total, "function samples to start adjoint trafo");
00127 
00129   t0 = nfft_clock_gettime_seconds();
00130   nfsoft_adjoint(&plan_nfsoft);
00131   t1 = nfft_clock_gettime_seconds();
00132   time = t1-t0;
00133   if ((bw+1)*(4*(bw+1)*(bw+1)-1)/3<=20)
00134      nfft_vpr_complex(plan_nfsoft.f_hat,(bw+1)*(4*(bw+1)*(bw+1)-1)/3,"SO(3) Fourier coefficients");
00135   else
00136      nfft_vpr_complex(plan_nfsoft.f_hat,20,"adjoint NFSOFT, 1st-20th Fourier coefficient");
00137   printf(" computed in %11le seconds\n",time);
00138 
00140   t0 = nfft_clock_gettime_seconds();
00141   nfsoft_adjoint(&plan_ndsoft);
00142   t1 = nfft_clock_gettime_seconds();
00143   time = t1-t0;
00144   if ((bw+1)*(4*(bw+1)*(bw+1)-1)/3<=20)
00145   nfft_vpr_complex(plan_ndsoft.f_hat,(bw+1)*(4*(bw+1)*(bw+1)-1)/3,"SO(3) Fourier coefficients");
00146   else
00147     nfft_vpr_complex(plan_ndsoft.f_hat,20,"adjoint NDSOFT, 1st-20th Fourier coefficient");
00148   printf(" computed in %11le seconds\n",time);
00149 
00150 
00152   error=nfft_error_l_infty_complex(plan_ndsoft.f_hat,plan_nfsoft.f_hat, (bw+1)*(4*(bw+1)*(bw+1)-1)/3);
00153   printf("\n The adjoint NFSOFT of bandwidth=%d for %d rotations has infty-error %11le \n",bw, M,error);
00154 
00155   printf("\n---------------------------------------------\n");
00156 
00158   nfsoft_finalize(&plan_ndsoft);
00159   nfsoft_finalize(&plan_nfsoft);
00160 }
00161 
00178 int main(int argc, char **argv)
00179 {
00180   int N; 
00181   int M; 
00183   if (argc < 2)
00184   {
00185     printf(
00186         "This test programm computes the NFSOFT with maximum polynomial degree N at M input rotations\n");
00187     printf("Usage: simple_test N M \n");
00188     printf("e.g.: simple_test 8 64\n");
00189     exit(0);
00190   }
00191 
00193   N = atoi(argv[1]);
00194   M = atoi(argv[2]);
00195 
00196   printf(
00197       "computing an NDSOFT, an NFSOFT, an adjoint NDSOFT, and an adjoint NFSOFT\n\n");
00198 
00199   simple_test_nfsoft(N, M);
00200 
00201 
00202 
00203   /* Exit the program. */
00204   return EXIT_SUCCESS;
00205 
00206 }