31 void bench_openmp_readfile(FILE *infile,
int *trafo_adjoint,
int *N,
int *M,
double **x, C **f_hat, C **f)
37 fscanf(infile,
"%d %d %d", trafo_adjoint, N, M);
38 *x = (
double *)nfft_malloc(2*(*M)*
sizeof(double));
39 *f_hat = (C*)nfft_malloc((2*(*N)+2) * (2*(*N)+2) *
sizeof(C));
40 *f = (C*)nfft_malloc((*M)*
sizeof(C));
42 memset(*f_hat,0U,(2*(*N)+2) * (2*(*N)+2) *
sizeof(C));
43 memset(*f,0U,(*M)*
sizeof(C));
46 fftw_import_wisdom_from_filename(
"nfsft_benchomp_detail_threads.plan");
48 fftw_import_wisdom_from_filename(
"nfsft_benchomp_detail_single.plan");
51 nfsft_init_guru(&plan, *N, *M, NFSFT_MALLOC_X | NFSFT_MALLOC_F |
52 NFSFT_MALLOC_F_HAT | NFSFT_NORMALIZED | NFSFT_PRESERVE_F_HAT,
53 PRE_PHI_HUT | FFTW_INIT | FFT_OUT_OF_PLACE, 6);
56 fftw_export_wisdom_to_filename(
"nfsft_benchomp_detail_threads.plan");
58 fftw_export_wisdom_to_filename(
"nfsft_benchomp_detail_single.plan");
61 for (j=0; j < *M; j++)
63 fscanf(infile,
"%lg", (*x)+2*j+t);
67 for (k = 0; k <= *N; k++)
68 for (n = -k; n <= k; n++)
70 fscanf(infile,
"%lg %lg", &re, &im);
71 (*f_hat)[NFSFT_INDEX(k,n,&plan)] = re + _Complex_I * im;
76 for (j=0; j < *M; j++)
78 fscanf(infile,
"%lg %lg", &re, &im);
79 (*f)[j] = re + _Complex_I * im;
83 nfsft_finalize(&plan);
86 void bench_openmp(
int trafo_adjoint,
int N,
int M,
double *x, C *f_hat, C *f,
int m,
int nfsft_flags,
int psi_flags)
93 double tt_total, tt_pre;
112 nfsft_init_guru(&plan, N, M, nfsft_flags | NFSFT_MALLOC_X | NFSFT_MALLOC_F |
113 NFSFT_MALLOC_F_HAT | NFSFT_NORMALIZED | NFSFT_PRESERVE_F_HAT,
114 PRE_PHI_HUT | psi_flags | FFTW_INIT | FFT_OUT_OF_PLACE, m);
122 for (j=0; j < plan.
M_total; j++)
124 for (t=0; t < 2; t++)
126 plan.
x[2*j+t] = x[2*j+t];
129 if (trafo_adjoint==0)
131 memset(plan.
f_hat,0U,plan.
N_total*
sizeof(
double _Complex));
132 for (k = 0; k <= plan.
N; k++)
133 for (n = -k; n <= k; n++)
137 plan.
f_hat[NFSFT_INDEX(k,n,&plan)] = f_hat[NFSFT_INDEX(k,n,&plan)];
142 for (j=0; j < plan.
M_total; j++)
149 memset(plan.
f_hat,0U,plan.
N_total*
sizeof(
double _Complex));
154 nfsft_precompute_x(&plan);
158 if (trafo_adjoint==0)
161 nfsft_adjoint(&plan);
170 #ifndef MEASURE_TIME_FFTW
177 nfsft_finalize(&plan);
180 int main(
int argc,
char **argv)
182 int m, nfsft_flags, psi_flags;
184 int trafo_adjoint, N, M, r;
193 nthreads = atoi(argv[5]);
195 omp_set_num_threads(nthreads);
202 nfsft_flags = atoi(argv[2]);
203 psi_flags = atoi(argv[3]);
204 nrepeat = atoi(argv[4]);
206 bench_openmp_readfile(stdin, &trafo_adjoint, &N, &M, &x, &f_hat, &f);
209 nfsft_precompute(N,1000.0,0U,0U);
211 for (r = 0; r < nrepeat; r++)
212 bench_openmp(trafo_adjoint, N, M, x, f_hat, f, m, nfsft_flags, psi_flags);
double * x
the nodes for ,
fftw_complex * f_hat
Fourier coefficients.
R nfft_elapsed_seconds(ticks t1, ticks t0)
Return number of elapsed seconds between two time points.
data structure for an NFSFT (nonequispaced fast spherical Fourier transform) plan with double precisi...
NFFT_INT M_total
Total number of samples.
double MEASURE_TIME_t[3]
Measured time for each step if MEASURE_TIME is set.
NFFT_INT N_total
Total number of Fourier coefficients.