31 static void simple_test_nfsoft(
int bw,
int M)
41 unsigned int flags = NFSOFT_MALLOC_X | NFSOFT_MALLOC_F | NFSOFT_MALLOC_F_HAT;
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);
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);
62 for (j = 0; j < plan_nfsoft.
M_total; j++)
64 d1 = ((double) rand()) / RAND_MAX - 0.5;
65 d2 = 0.5 * ((double) rand()) / RAND_MAX;
66 d3 = ((double) rand()) / RAND_MAX - 0.5;
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;
78 for (j = 0; j < (bw + 1) * (4* (bw +1)*(bw+1)-1)/3;j++)
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;
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");
91 printf(
"\n---------------------------------------------\n");
94 nfsoft_precompute(&plan_nfsoft);
95 nfsoft_precompute(&plan_ndsoft);
99 t0 = nfft_clock_gettime_seconds();
100 nfsoft_trafo(&plan_nfsoft);
101 t1 = nfft_clock_gettime_seconds();
107 printf(
" computed in %11le seconds\n",time);
110 t0 = nfft_clock_gettime_seconds();
111 nfsoft_trafo(&plan_ndsoft);
112 t1 = nfft_clock_gettime_seconds();
118 printf(
" computed in %11le seconds\n",time);
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);
124 printf(
"\n---------------------------------------------\n");
126 plan_nfsoft.
f[0]=1.0;
127 plan_ndsoft.
f[0]=1.0;
131 t0 = nfft_clock_gettime_seconds();
132 nfsoft_adjoint(&plan_nfsoft);
133 t1 = nfft_clock_gettime_seconds();
135 if ((bw+1)*(4*(bw+1)*(bw+1)-1)/3<=20)
139 printf(
" computed in %11le seconds\n",time);
142 t0 = nfft_clock_gettime_seconds();
143 nfsoft_adjoint(&plan_ndsoft);
144 t1 = nfft_clock_gettime_seconds();
146 if ((bw+1)*(4*(bw+1)*(bw+1)-1)/3<=20)
150 printf(
" computed in %11le seconds\n",time);
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);
157 printf(
"\n---------------------------------------------\n");
160 nfsoft_finalize(&plan_ndsoft);
161 nfsoft_finalize(&plan_nfsoft);
180 int main(
int argc,
char **argv)
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");
199 "computing an NDSOFT, an NFSOFT, an adjoint NDSOFT, and an adjoint NFSOFT\n\n");
201 simple_test_nfsoft(N, M);
fftw_complex * f_hat
Fourier coefficients.
NFFT_INT M_total
Total number of samples.
void nfft_vpr_complex(fftw_complex *x, const NFFT_INT n, const char *text)
Print complex vector to standard output.