NFFT  3.3.0
fastsum_benchomp_detail.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: simple_test.c 3372 2009-10-21 06:04:05Z skunis $ */
20 
21 #include <stdio.h>
22 #include <math.h>
23 #include <string.h>
24 #include <stdlib.h>
25 #include <complex.h>
26 
27 #include "nfft3.h"
28 #include "infft.h"
29 #include "fastsum.h"
30 #include "kernels.h"
31 #ifdef _OPENMP
32 #include <omp.h>
33 #endif
34 
35 int bench_openmp(FILE *infile, int n, int m, int p,
36  C (*kernel)(R, int, const R *), R c, R eps_I, R eps_B)
37 {
38  fastsum_plan my_fastsum_plan;
39  int d, L, M;
40  int t, j;
41  R re, im;
42  R r_max = K(0.25) - my_fastsum_plan.eps_B / K(2.0);
43  ticks t0, t1;
44  R tt_total;
45 
46  fscanf(infile, "%d %d %d", &d, &L, &M);
47 
48 #ifdef _OPENMP
49  FFTW(import_wisdom_from_filename)("fastsum_benchomp_detail_threads.plan");
50 #else
51  FFTW(import_wisdom_from_filename)("fastsum_benchomp_detail_single.plan");
52 #endif
53 
54  fastsum_init_guru(&my_fastsum_plan, d, L, M, kernel, &c, NEARFIELD_BOXES, n,
55  m, p, eps_I, eps_B);
56 
57 #ifdef _OPENMP
58  FFTW(export_wisdom_to_filename)("fastsum_benchomp_detail_threads.plan");
59 #else
60  FFTW(export_wisdom_to_filename)("fastsum_benchomp_detail_single.plan");
61 #endif
62 
63  for (j = 0; j < L; j++)
64  {
65  for (t = 0; t < d; t++)
66  {
67  R v;
68  fscanf(infile, __FR__, &v);
69  my_fastsum_plan.x[d * j + t] = v * r_max;
70  }
71  }
72 
73  for (j = 0; j < L; j++)
74  {
75  fscanf(infile, __FR__ " " __FR__, &re, &im);
76  my_fastsum_plan.alpha[j] = re + II * im;
77  }
78 
79  for (j = 0; j < M; j++)
80  {
81  for (t = 0; t < d; t++)
82  {
83  R v;
84  fscanf(infile, __FR__, &v);
85  my_fastsum_plan.y[d * j + t] = v * r_max;
86  }
87  }
88 
90  t0 = getticks();
91  fastsum_precompute(&my_fastsum_plan);
92 
94  fastsum_trafo(&my_fastsum_plan);
95  t1 = getticks();
96  tt_total = NFFT(elapsed_seconds)(t1, t0);
97 
98 #ifndef MEASURE_TIME
99  my_fastsum_plan.MEASURE_TIME_t[0] = K(0.0);
100  my_fastsum_plan.MEASURE_TIME_t[1] = K(0.0);
101  my_fastsum_plan.MEASURE_TIME_t[2] = K(0.0);
102  my_fastsum_plan.MEASURE_TIME_t[3] = K(0.0);
103  my_fastsum_plan.MEASURE_TIME_t[4] = K(0.0);
104  my_fastsum_plan.MEASURE_TIME_t[5] = K(0.0);
105  my_fastsum_plan.MEASURE_TIME_t[6] = K(0.0);
106  my_fastsum_plan.MEASURE_TIME_t[7] = K(0.0);
107  my_fastsum_plan.mv1.MEASURE_TIME_t[0] = K(0.0);
108  my_fastsum_plan.mv1.MEASURE_TIME_t[2] = K(0.0);
109  my_fastsum_plan.mv2.MEASURE_TIME_t[0] = K(0.0);
110  my_fastsum_plan.mv2.MEASURE_TIME_t[2] = K(0.0);
111 #endif
112 #ifndef MEASURE_TIME_FFTW
113  my_fastsum_plan.mv1.MEASURE_TIME_t[1] = K(0.0);
114  my_fastsum_plan.mv2.MEASURE_TIME_t[1] = K(0.0);
115 #endif
116 
117  printf(
118  "%.6" __FES__ " %.6" __FES__ " %.6" __FES__ " %6" __FES__ " %.6" __FES__ " %.6" __FES__ " %.6" __FES__ " %.6" __FES__ " %.6" __FES__ " %6" __FES__ " %.6" __FES__ " %.6" __FES__ " %6" __FES__ " %.6" __FES__ " %.6" __FES__ " %6" __FES__ "\n",
119  my_fastsum_plan.MEASURE_TIME_t[0], my_fastsum_plan.MEASURE_TIME_t[1],
120  my_fastsum_plan.MEASURE_TIME_t[2], my_fastsum_plan.MEASURE_TIME_t[3],
121  my_fastsum_plan.MEASURE_TIME_t[4], my_fastsum_plan.MEASURE_TIME_t[5],
122  my_fastsum_plan.MEASURE_TIME_t[6], my_fastsum_plan.MEASURE_TIME_t[7],
123  tt_total - my_fastsum_plan.MEASURE_TIME_t[0]
124  - my_fastsum_plan.MEASURE_TIME_t[1]
125  - my_fastsum_plan.MEASURE_TIME_t[2]
126  - my_fastsum_plan.MEASURE_TIME_t[3]
127  - my_fastsum_plan.MEASURE_TIME_t[4]
128  - my_fastsum_plan.MEASURE_TIME_t[5]
129  - my_fastsum_plan.MEASURE_TIME_t[6]
130  - my_fastsum_plan.MEASURE_TIME_t[7], tt_total,
131  my_fastsum_plan.mv1.MEASURE_TIME_t[0],
132  my_fastsum_plan.mv1.MEASURE_TIME_t[1],
133  my_fastsum_plan.mv1.MEASURE_TIME_t[2],
134  my_fastsum_plan.mv2.MEASURE_TIME_t[0],
135  my_fastsum_plan.mv2.MEASURE_TIME_t[1],
136  my_fastsum_plan.mv2.MEASURE_TIME_t[2]);
137 
138  fastsum_finalize(&my_fastsum_plan);
139 
140  return 0;
141 }
142 
143 int main(int argc, char **argv)
144 {
145  int n;
146  int m;
147  int p;
148  char *s;
149  C (*kernel)(R, int, const R *);
150  R c;
151  R eps_I;
152  R eps_B;
154 #ifdef _OPENMP
155  int nthreads;
156 
157  if (argc != 9)
158  return EXIT_FAILURE;
159 
160  nthreads = atoi(argv[8]);
161  FFTW(init_threads)();
162  omp_set_num_threads(nthreads);
163 #else
164  if (argc != 8)
165  return EXIT_FAILURE;
166 #endif
167 
168  n = atoi(argv[1]);
169  m = atoi(argv[2]);
170  p = atoi(argv[3]);
171  s = argv[4];
172  c = atof(argv[5]);
173  eps_I = (R)(atof(argv[6]));
174  eps_B = (R)(atof(argv[7]));
175  if (strcmp(s, "gaussian") == 0)
176  kernel = gaussian;
177  else if (strcmp(s, "multiquadric") == 0)
178  kernel = multiquadric;
179  else if (strcmp(s, "inverse_multiquadric") == 0)
180  kernel = inverse_multiquadric;
181  else if (strcmp(s, "logarithm") == 0)
182  kernel = logarithm;
183  else if (strcmp(s, "thinplate_spline") == 0)
184  kernel = thinplate_spline;
185  else if (strcmp(s, "one_over_square") == 0)
186  kernel = one_over_square;
187  else if (strcmp(s, "one_over_modulus") == 0)
188  kernel = one_over_modulus;
189  else if (strcmp(s, "one_over_x") == 0)
190  kernel = one_over_x;
191  else if (strcmp(s, "inverse_multiquadric3") == 0)
192  kernel = inverse_multiquadric3;
193  else if (strcmp(s, "sinc_kernel") == 0)
194  kernel = sinc_kernel;
195  else if (strcmp(s, "cosc") == 0)
196  kernel = cosc;
197  else if (strcmp(s, "cot") == 0)
198  kernel = kcot;
199  else
200  {
201  s = "multiquadric";
202  kernel = multiquadric;
203  }
204 
205  bench_openmp(stdin, n, m, p, kernel, c, eps_I, eps_B);
206 
207  return EXIT_SUCCESS;
208 }
209 
R eps_B
outer boundary
Definition: fastsum.h:105
Header file with predefined kernels for the fast summation algorithm.
plan for fast summation algorithm
Definition: fastsum.h:74
Header file for the fast NFFT-based summation algorithm.
C * alpha
source coefficients
Definition: fastsum.h:83
void fastsum_trafo(fastsum_plan *ths)
fast NFFT-based summation
Definition: fastsum.c:1057
void fastsum_precompute(fastsum_plan *ths)
precomputation for fastsum
Definition: fastsum.c:909
R MEASURE_TIME_t[8]
Measured time for each step if MEASURE_TIME is set.
Definition: fastsum.h:123
void fastsum_init_guru(fastsum_plan *ths, int d, int N_total, int M_total, kernel k, R *param, unsigned flags, int nn, int m, int p, R eps_I, R eps_B)
initialization of fastsum plan
Definition: fastsum.c:693
R * y
target knots in d-ball with radius 1/4-eps_b/2
Definition: fastsum.h:87
void fastsum_finalize(fastsum_plan *ths)
finalization of fastsum plan
Definition: fastsum.c:846
R * x
source knots in d-ball with radius 1/4-eps_b/2
Definition: fastsum.h:86