NFFT 3.5.3alpha
nsfft_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 <stdio.h>
20#include <math.h>
21#include <string.h>
22#include <stdlib.h>
23#include <complex.h>
24
25#include "nfft3.h"
26
28#define CSWAP(x,y) {double _Complex * NFFT_SWAP_temp__; \
29 NFFT_SWAP_temp__=(x); (x)=(y); (y)=NFFT_SWAP_temp__;}
30
31static void accuracy_nsfft(int d, int J, int M, int m)
32{
33 nsfft_plan p;
34 double _Complex *swap_sndft_trafo, *swap_sndft_adjoint;
35
36 nsfft_init(&p, d, J, M, m, NSDFT);
37
38 swap_sndft_trafo=(double _Complex*) nfft_malloc(p.M_total*
39 sizeof(double _Complex));
40 swap_sndft_adjoint=(double _Complex*) nfft_malloc(p.N_total*
41 sizeof(double _Complex));
42
43 nsfft_init_random_nodes_coeffs(&p);
44
46 nsfft_trafo_direct(&p);
47
48 CSWAP(swap_sndft_trafo,p.f);
49
51 nsfft_trafo(&p);
52
53 printf("%5d\t %+.5E\t",J,
54 nfft_error_l_infty_1_complex(swap_sndft_trafo, p.f, p.M_total,
55 p.f_hat, p.N_total));
56 fflush(stdout);
57
59
61 nsfft_adjoint_direct(&p);
62
63 CSWAP(swap_sndft_adjoint,p.f_hat);
64
66 nsfft_adjoint(&p);
67
68 printf("%+.5E\n",
69 nfft_error_l_infty_1_complex(swap_sndft_adjoint, p.f_hat,
70 p.N_total,
71 p.f, p.M_total));
72 fflush(stdout);
73
74 nfft_free(swap_sndft_adjoint);
75 nfft_free(swap_sndft_trafo);
76
78 nsfft_finalize(&p);
79}
80
81static void time_nsfft(int d, int J, int M, unsigned test_nsdft, unsigned test_nfft)
82{
83 int r, N[d], n[d];
84 double t, t_nsdft, t_nfft, t_nsfft;
85 double t0, t1;
86
87 nsfft_plan p;
88 nfft_plan np;
89
90 for(r=0;r<d;r++)
91 {
92 N[r]= nfft_exp2i(J+2);
93 n[r]=(3*N[r])/2;
94 /*n[r]=2*N[r];*/
95 }
96
98 nsfft_init(&p, d, J, M, 4, NSDFT);
99 nsfft_init_random_nodes_coeffs(&p);
100
101 /* transforms */
102 if(test_nsdft)
103 {
104 t_nsdft=0;
105 r=0;
106 while(t_nsdft<0.1)
107 {
108 r++;
109 t0 = nfft_clock_gettime_seconds();
110 nsfft_trafo_direct(&p);
111 t1 = nfft_clock_gettime_seconds();
112 t = (t1 - t0);
113 t_nsdft+=t;
114 }
115 t_nsdft/=r;
116 }
117 else
118 t_nsdft=nan("");
119
120 if(test_nfft)
121 {
122 nfft_init_guru(&np,d,N,M,n,6, FG_PSI| MALLOC_F_HAT| MALLOC_F| FFTW_INIT,
123 FFTW_MEASURE);
124 np.x=p.act_nfft_plan->x;
125 if(np.flags & PRE_ONE_PSI)
126 nfft_precompute_one_psi(&np);
127 nsfft_cp(&p, &np);
128
129 t_nfft=0;
130 r=0;
131 while(t_nfft<0.1)
132 {
133 r++;
134 t0 = nfft_clock_gettime_seconds();
135 nfft_trafo(&np);
136 t1 = nfft_clock_gettime_seconds();
137 t = (t1 - t0);
138 t_nfft+=t;
139 }
140 t_nfft/=r;
141
142 nfft_finalize(&np);
143 }
144 else
145 {
146 t_nfft=nan("");
147 }
148
149 t_nsfft=0;
150 r=0;
151 while(t_nsfft<0.1)
152 {
153 r++;
154 t0 = nfft_clock_gettime_seconds();
155 nsfft_trafo(&p);
156 t1 = nfft_clock_gettime_seconds();
157 t = (t1 - t0);
158 t_nsfft+=t;
159 }
160 t_nsfft/=r;
161
162 printf("%d\t%.2e\t%.2e\t%.2e\n", J, t_nsdft, t_nfft, t_nsfft);
163
164 fflush(stdout);
165
167 nsfft_finalize(&p);
168}
169
170
171int main(int argc,char **argv)
172{
173 int d, J, M;
174
175 if(argc<=2)
176 {
177 fprintf(stderr,"nsfft_test type d [first last trials]\n");
178 return -1;
179 }
180
181 d=atoi(argv[2]);
182 fprintf(stderr,"Testing the nfft on the hyperbolic cross (nsfft).\n");
183
184 if(atoi(argv[1])==1)
185 {
186 fprintf(stderr,"Testing the accuracy of the nsfft vs. nsdft\n");
187 fprintf(stderr,"Columns: d, E_{1,\\infty}(trafo) E_{1,\\infty}(adjoint)\n\n");
188 for(J=1; J<10; J++)
189 accuracy_nsfft(d, J, 1000, 6);
190 }
191
192 if(atoi(argv[1])==2)
193 {
194 fprintf(stderr,"Testing the computation time of the nsdft, nfft, and nsfft\n");
195 fprintf(stderr,"Columns: d, J, M, t_nsdft, t_nfft, t_nsfft\n\n");
196 for(J=atoi(argv[3]); J<=atoi(argv[4]); J++)
197 {
198 if(d==2)
199 M=(J+4)*nfft_exp2i(J+1);
200 else
201 M=6*nfft_exp2i(J)*(nfft_exp2i((J+1)/2+1)-1)+nfft_exp2i(3*(J/2+1));
202
203 if(d*(J+2)<=24)
204 time_nsfft(d, J, M, 1, 1);
205 else
206 if(d*(J+2)<=24)
207 time_nsfft(d, J, M, 0, 1);
208 else
209 time_nsfft(d, J, M, 0, 0);
210 }
211 }
212
213 return 1;
214}
#define FG_PSI
Definition nfft3.h:188
#define MALLOC_F_HAT
Definition nfft3.h:194
#define PRE_ONE_PSI
Definition nfft3.h:200
#define MALLOC_F
Definition nfft3.h:195
#define FFTW_INIT
Definition nfft3.h:197
#define CSWAP(x, y)
Swap two vectors.
Definition infft.h:139
#define NSDFT
Definition nfft3.h:482
Header file for the nfft3 library.
void nfft_vrand_unit_complex(fftw_complex *x, const NFFT_INT n)
Inits a vector of random complex numbers in .
void * nfft_malloc(size_t n)
void nfft_free(void *p)
data structure for an NSFFT (nonequispaced sparse fast Fourier transform) plan with double precision
Definition nfft3.h:478
NFFT_INT M_total
Total number of samples.
Definition nfft3.h:478
fftw_complex * f
Samples.
Definition nfft3.h:478
nfft_plan * act_nfft_plan
current nfft block
Definition nfft3.h:478
fftw_complex * f_hat
Fourier coefficients.
Definition nfft3.h:478
NFFT_INT N_total
Total number of Fourier coefficients.
Definition nfft3.h:478