NFFT  3.3.0
nfsoft/simple_test.c
1 /*
2  * Copyright (c) 2002, 2015 Jens Keiner, Stefan Kunis, Daniel Potts
3  *
4  * This program is free software; you can redistribute it and/or modify it under
5  * the terms of the GNU General Public License as published by the Free Software
6  * Foundation; either version 2 of the License, or (at your option) any later
7  * version.
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12  * details.
13  *
14  * You should have received a copy of the GNU General Public License along with
15  * this program; if not, write to the Free Software Foundation, Inc., 51
16  * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  */
18 
19 /* $Id$ */
20 
21 /* Include standard C headers. */
22 #include <stdio.h>
23 #include <math.h>
24 #include <string.h>
25 #include <stdlib.h>
26 #include <complex.h>
27 
28 /* Include NFFT3 library header. */
29 #include "nfft3.h"
30 
31 static void simple_test_nfsoft(int bw, int M)
32 {
33  nfsoft_plan plan_nfsoft;
34  nfsoft_plan plan_ndsoft;
36  double t0, t1;
37  int j;
38  int k, m;
39  double d1, d2, d3;
40  double time, error;
41  unsigned int flags = NFSOFT_MALLOC_X | NFSOFT_MALLOC_F | NFSOFT_MALLOC_F_HAT;
44  k = 1000;
45  m = 5;
54  nfsoft_init_guru(&plan_ndsoft, bw, M, flags | NFSOFT_USE_NDFT
55  | NFSOFT_USE_DPT, PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT
56  | MALLOC_F | FFTW_INIT | FFT_OUT_OF_PLACE, m, k);
57 
58  nfsoft_init_guru(&plan_nfsoft, bw, M, flags, PRE_PHI_HUT | PRE_PSI | MALLOC_X
59  | MALLOC_F_HAT | MALLOC_F | FFTW_INIT | FFT_OUT_OF_PLACE, m, k);
60 
62  for (j = 0; j < plan_nfsoft.M_total; j++)
63  {
64  d1 = ((double) rand()) / RAND_MAX - 0.5;
65  d2 = 0.5 * ((double) rand()) / RAND_MAX;
66  d3 = ((double) rand()) / RAND_MAX - 0.5;
67 
68  plan_nfsoft.x[3* j ] = d1;
69  plan_nfsoft.x[3* j + 1] = d2;
70  plan_nfsoft.x[3* j + 2] = d3;
72  plan_ndsoft.x[3* j ] = d1;
73  plan_ndsoft.x[3* j + 1] = d2;
74  plan_ndsoft.x[3* j + 2] = d3;
75  }
76 
78  for (j = 0; j < (bw + 1) * (4* (bw +1)*(bw+1)-1)/3;j++)
79  {
80  d1=((double)rand())/RAND_MAX - 0.5;
81  d2=((double)rand())/RAND_MAX - 0.5;
82  plan_nfsoft.f_hat[j]=d1 + I*d2;
83  plan_ndsoft.f_hat[j]=d1 + I*d2;
84  }
85 
86  if ((bw+1)*(4*(bw+1)*(bw+1)-1)/3<=20)
87  nfft_vpr_complex(plan_nfsoft.f_hat,(bw+1)*(4*(bw+1)*(bw+1)-1)/3,"randomly generated SO(3) Fourier coefficients");
88  else
89  nfft_vpr_complex(plan_ndsoft.f_hat,20,"1st-20th randomly generated SO(3) Fourier coefficient");
90 
91  printf("\n---------------------------------------------\n");
92 
94  nfsoft_precompute(&plan_nfsoft);
95  nfsoft_precompute(&plan_ndsoft);
96 
97 
99  t0 = nfft_clock_gettime_seconds();
100  nfsoft_trafo(&plan_nfsoft);
101  t1 = nfft_clock_gettime_seconds();
102  time = t1-t0;
103  if (plan_nfsoft.M_total<=20)
104  nfft_vpr_complex(plan_nfsoft.f,plan_nfsoft.M_total,"NFSOFT, function samples");
105  else
106  nfft_vpr_complex(plan_nfsoft.f,20,"NFSOFT, 1st-20th function sample");
107  printf(" computed in %11le seconds\n",time);
108 
110  t0 = nfft_clock_gettime_seconds();
111  nfsoft_trafo(&plan_ndsoft);
112  t1 = nfft_clock_gettime_seconds();
113  time = t1-t0;
114  if (plan_ndsoft.M_total<=20)
115  nfft_vpr_complex(plan_ndsoft.f,plan_ndsoft.M_total,"NDSOFT, function samples");
116  else
117  nfft_vpr_complex(plan_ndsoft.f,20,"NDSOFT, 1st-20th function sample");
118  printf(" computed in %11le seconds\n",time);
119 
121  error= nfft_error_l_infty_complex(plan_ndsoft.f,plan_nfsoft.f, plan_nfsoft.M_total);
122  printf("\n The NFSOFT of bandwidth=%d for %d rotations has infty-error %11le \n",bw, M,error);
123 
124  printf("\n---------------------------------------------\n");
125 
126  plan_nfsoft.f[0]=1.0;
127  plan_ndsoft.f[0]=1.0;
128  nfft_vpr_complex(plan_ndsoft.f,plan_ndsoft.M_total, "function samples to start adjoint trafo");
129 
131  t0 = nfft_clock_gettime_seconds();
132  nfsoft_adjoint(&plan_nfsoft);
133  t1 = nfft_clock_gettime_seconds();
134  time = t1-t0;
135  if ((bw+1)*(4*(bw+1)*(bw+1)-1)/3<=20)
136  nfft_vpr_complex(plan_nfsoft.f_hat,(bw+1)*(4*(bw+1)*(bw+1)-1)/3,"SO(3) Fourier coefficients");
137  else
138  nfft_vpr_complex(plan_nfsoft.f_hat,20,"adjoint NFSOFT, 1st-20th Fourier coefficient");
139  printf(" computed in %11le seconds\n",time);
140 
142  t0 = nfft_clock_gettime_seconds();
143  nfsoft_adjoint(&plan_ndsoft);
144  t1 = nfft_clock_gettime_seconds();
145  time = t1-t0;
146  if ((bw+1)*(4*(bw+1)*(bw+1)-1)/3<=20)
147  nfft_vpr_complex(plan_ndsoft.f_hat,(bw+1)*(4*(bw+1)*(bw+1)-1)/3,"SO(3) Fourier coefficients");
148  else
149  nfft_vpr_complex(plan_ndsoft.f_hat,20,"adjoint NDSOFT, 1st-20th Fourier coefficient");
150  printf(" computed in %11le seconds\n",time);
151 
152 
154  error=nfft_error_l_infty_complex(plan_ndsoft.f_hat,plan_nfsoft.f_hat, (bw+1)*(4*(bw+1)*(bw+1)-1)/3);
155  printf("\n The adjoint NFSOFT of bandwidth=%d for %d rotations has infty-error %11le \n",bw, M,error);
156 
157  printf("\n---------------------------------------------\n");
158 
160  nfsoft_finalize(&plan_ndsoft);
161  nfsoft_finalize(&plan_nfsoft);
162 }
163 
180 int main(int argc, char **argv)
181 {
182  int N;
183  int M;
185  if (argc < 2)
186  {
187  printf(
188  "This test programm computes the NFSOFT with maximum polynomial degree N at M input rotations\n");
189  printf("Usage: simple_test N M \n");
190  printf("e.g.: simple_test 8 64\n");
191  exit(0);
192  }
193 
195  N = atoi(argv[1]);
196  M = atoi(argv[2]);
197 
198  printf(
199  "computing an NDSOFT, an NFSOFT, an adjoint NDSOFT, and an adjoint NFSOFT\n\n");
200 
201  simple_test_nfsoft(N, M);
202 
203 
204 
205  /* Exit the program. */
206  return EXIT_SUCCESS;
207 
208 }
fftw_complex * f_hat
Fourier coefficients.
Definition: nfft3.h:686
fftw_complex * f
Samples.
Definition: nfft3.h:686
NFFT_INT M_total
Total number of samples.
Definition: nfft3.h:686
double * x
input nodes
Definition: nfft3.h:686
void nfft_vpr_complex(fftw_complex *x, const NFFT_INT n, const char *text)
Print complex vector to standard output.