NFFT 3.5.3alpha
nfsoft/simple_test.c
1/*
2 * Copyright (c) 2002, 2017 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/* Include standard C headers. */
20#include <stdio.h>
21#include <math.h>
22#include <string.h>
23#include <stdlib.h>
24#include <complex.h>
25
26/* Include NFFT3 library header. */
27#include "nfft3.h"
28
29static void simple_test_nfsoft(int bw, int M)
30{
31 nfsoft_plan plan_nfsoft;
32 nfsoft_plan plan_ndsoft;
34 double t0, t1;
35 int j;
36 int k, m;
37 double d1, d2, d3;
38 double time, error;
39 unsigned int flags = NFSOFT_MALLOC_X | NFSOFT_MALLOC_F | NFSOFT_MALLOC_F_HAT;
42 k = 1000;
43 m = 5;
52 nfsoft_init_guru(&plan_ndsoft, bw, M, flags | NFSOFT_USE_NDFT
55
56 nfsoft_init_guru(&plan_nfsoft, bw, M, flags, PRE_PHI_HUT | PRE_PSI | MALLOC_X
58
60 for (j = 0; j < plan_nfsoft.M_total; j++)
61 {
62 d1 = ((double) rand()) / RAND_MAX - 0.5;
63 d2 = 0.5 * ((double) rand()) / RAND_MAX;
64 d3 = ((double) rand()) / RAND_MAX - 0.5;
65
66 plan_nfsoft.x[3* j ] = d1;
67 plan_nfsoft.x[3* j + 1] = d2;
68 plan_nfsoft.x[3* j + 2] = d3;
70 plan_ndsoft.x[3* j ] = d1;
71 plan_ndsoft.x[3* j + 1] = d2;
72 plan_ndsoft.x[3* j + 2] = d3;
73 }
74
76 for (j = 0; j < (bw + 1) * (4* (bw +1)*(bw+1)-1)/3;j++)
77 {
78 d1=((double)rand())/RAND_MAX - 0.5;
79 d2=((double)rand())/RAND_MAX - 0.5;
80 plan_nfsoft.f_hat[j]=d1 + I*d2;
81 plan_ndsoft.f_hat[j]=d1 + I*d2;
82 }
83
84 if ((bw+1)*(4*(bw+1)*(bw+1)-1)/3<=20)
85 nfft_vpr_complex(plan_nfsoft.f_hat,(bw+1)*(4*(bw+1)*(bw+1)-1)/3,"randomly generated SO(3) Fourier coefficients");
86 else
87 nfft_vpr_complex(plan_ndsoft.f_hat,20,"1st-20th randomly generated SO(3) Fourier coefficient");
88
89 printf("\n---------------------------------------------\n");
90
92 nfsoft_precompute(&plan_nfsoft);
93 nfsoft_precompute(&plan_ndsoft);
94
95
97 t0 = nfft_clock_gettime_seconds();
98 nfsoft_trafo(&plan_nfsoft);
99 t1 = nfft_clock_gettime_seconds();
100 time = t1-t0;
101 if (plan_nfsoft.M_total<=20)
102 nfft_vpr_complex(plan_nfsoft.f,plan_nfsoft.M_total,"NFSOFT, function samples");
103 else
104 nfft_vpr_complex(plan_nfsoft.f,20,"NFSOFT, 1st-20th function sample");
105 printf(" computed in %11le seconds\n",time);
106
108 t0 = nfft_clock_gettime_seconds();
109 nfsoft_trafo(&plan_ndsoft);
110 t1 = nfft_clock_gettime_seconds();
111 time = t1-t0;
112 if (plan_ndsoft.M_total<=20)
113 nfft_vpr_complex(plan_ndsoft.f,plan_ndsoft.M_total,"NDSOFT, function samples");
114 else
115 nfft_vpr_complex(plan_ndsoft.f,20,"NDSOFT, 1st-20th function sample");
116 printf(" computed in %11le seconds\n",time);
117
119 error= nfft_error_l_infty_complex(plan_ndsoft.f,plan_nfsoft.f, plan_nfsoft.M_total);
120 printf("\n The NFSOFT of bandwidth=%d for %d rotations has infty-error %11le \n",bw, M,error);
121
122 printf("\n---------------------------------------------\n");
123
124 plan_nfsoft.f[0]=1.0;
125 plan_ndsoft.f[0]=1.0;
126 nfft_vpr_complex(plan_ndsoft.f,plan_ndsoft.M_total, "function samples to start adjoint trafo");
127
129 t0 = nfft_clock_gettime_seconds();
130 nfsoft_adjoint(&plan_nfsoft);
131 t1 = nfft_clock_gettime_seconds();
132 time = t1-t0;
133 if ((bw+1)*(4*(bw+1)*(bw+1)-1)/3<=20)
134 nfft_vpr_complex(plan_nfsoft.f_hat,(bw+1)*(4*(bw+1)*(bw+1)-1)/3,"SO(3) Fourier coefficients");
135 else
136 nfft_vpr_complex(plan_nfsoft.f_hat,20,"adjoint NFSOFT, 1st-20th Fourier coefficient");
137 printf(" computed in %11le seconds\n",time);
138
140 t0 = nfft_clock_gettime_seconds();
141 nfsoft_adjoint(&plan_ndsoft);
142 t1 = nfft_clock_gettime_seconds();
143 time = t1-t0;
144 if ((bw+1)*(4*(bw+1)*(bw+1)-1)/3<=20)
145 nfft_vpr_complex(plan_ndsoft.f_hat,(bw+1)*(4*(bw+1)*(bw+1)-1)/3,"SO(3) Fourier coefficients");
146 else
147 nfft_vpr_complex(plan_ndsoft.f_hat,20,"adjoint NDSOFT, 1st-20th Fourier coefficient");
148 printf(" computed in %11le seconds\n",time);
149
150
152 error=nfft_error_l_infty_complex(plan_ndsoft.f_hat,plan_nfsoft.f_hat, (bw+1)*(4*(bw+1)*(bw+1)-1)/3);
153 printf("\n The adjoint NFSOFT of bandwidth=%d for %d rotations has infty-error %11le \n",bw, M,error);
154
155 printf("\n---------------------------------------------\n");
156
158 nfsoft_finalize(&plan_ndsoft);
159 nfsoft_finalize(&plan_nfsoft);
160}
161
178int main(int argc, char **argv)
179{
180 int N;
181 int M;
183 if (argc < 2)
184 {
185 printf(
186 "This test programm computes the NFSOFT with maximum polynomial degree N at M input rotations\n");
187 printf("Usage: simple_test N M \n");
188 printf("e.g.: simple_test 8 64\n");
189 exit(0);
190 }
191
193 N = atoi(argv[1]);
194 M = atoi(argv[2]);
195
196 printf(
197 "computing an NDSOFT, an NFSOFT, an adjoint NDSOFT, and an adjoint NFSOFT\n\n");
198
199 simple_test_nfsoft(N, M);
200
201
202
203 /* Exit the program. */
204 return EXIT_SUCCESS;
205
206}
#define MALLOC_F_HAT
Definition nfft3.h:194
#define MALLOC_X
Definition nfft3.h:193
#define FFT_OUT_OF_PLACE
Definition nfft3.h:196
#define PRE_PSI
Definition nfft3.h:191
#define MALLOC_F
Definition nfft3.h:195
#define FFTW_INIT
Definition nfft3.h:197
#define PRE_PHI_HUT
Definition nfft3.h:187
#define NFSOFT_USE_NDFT
Definition nfft3.h:704
#define NFSOFT_USE_DPT
Definition nfft3.h:705
#define NFSOFT_MALLOC_F_HAT
Definition nfft3.h:708
#define NFSOFT_MALLOC_X
Definition nfft3.h:706
#define NFSOFT_MALLOC_F
Definition nfft3.h:709
Header file for the nfft3 library.
void nfft_vpr_complex(fftw_complex *x, const NFFT_INT n, const char *text)
Print complex vector to standard output.
fftw_complex * f
Samples.
Definition nfft3.h:699
fftw_complex * f_hat
Fourier coefficients.
Definition nfft3.h:699
NFFT_INT M_total
Total number of samples.
Definition nfft3.h:699
double * x
input nodes
Definition nfft3.h:699