NFFT  3.5.0
fastsum_test.c
Go to the documentation of this file.
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 
25 #include "config.h"
26 
27 #include <stdlib.h>
28 #include <stdio.h>
29 #include <string.h>
30 #include <math.h>
31 #ifdef HAVE_COMPLEX_H
32  #include <complex.h>
33 #endif
34 
35 #ifdef _OPENMP
36  #include <omp.h>
37 #endif
38 
39 #include "fastsum.h"
40 #include "kernels.h"
41 #include "infft.h"
42 
49 int main(int argc, char **argv)
50 {
51  int j, k;
52  int d;
53  int N;
54  int M;
55  int n;
56  int m;
57  int p;
58  const char *s;
59  C (*kernel)(R, int, const R *);
60  R c;
61  fastsum_plan my_fastsum_plan;
62  C *direct;
63  ticks t0, t1;
64  R time;
65  R error = K(0.0);
66  R eps_I;
67  R eps_B;
69  if (argc != 11)
70  {
71  printf("\nfastsum_test d N M n m p kernel c eps_I eps_B\n\n");
72  printf(" d dimension \n");
73  printf(" N number of source nodes \n");
74  printf(" M number of target nodes \n");
75  printf(" n expansion degree \n");
76  printf(" m cut-off parameter \n");
77  printf(" p degree of smoothness \n");
78  printf(" kernel kernel function (e.g., gaussian)\n");
79  printf(" c kernel parameter \n");
80  printf(" eps_I inner boundary \n");
81  printf(" eps_B outer boundary \n\n");
82  exit(EXIT_FAILURE);
83  }
84  else
85  {
86  d = atoi(argv[1]);
87  N = atoi(argv[2]);
88  c = K(1.0) / POW((R)(N), K(1.0) / ((R)(d)));
89  M = atoi(argv[3]);
90  n = atoi(argv[4]);
91  m = atoi(argv[5]);
92  p = atoi(argv[6]);
93  s = argv[7];
94  c = (R)(atof(argv[8]));
95  eps_I = (R)(atof(argv[9]));
96  eps_B = (R)(atof(argv[10]));
97  if (strcmp(s, "gaussian") == 0)
98  kernel = gaussian;
99  else if (strcmp(s, "multiquadric") == 0)
100  kernel = multiquadric;
101  else if (strcmp(s, "inverse_multiquadric") == 0)
102  kernel = inverse_multiquadric;
103  else if (strcmp(s, "logarithm") == 0)
104  kernel = logarithm;
105  else if (strcmp(s, "thinplate_spline") == 0)
106  kernel = thinplate_spline;
107  else if (strcmp(s, "one_over_square") == 0)
108  kernel = one_over_square;
109  else if (strcmp(s, "one_over_modulus") == 0)
110  kernel = one_over_modulus;
111  else if (strcmp(s, "one_over_x") == 0)
112  kernel = one_over_x;
113  else if (strcmp(s, "inverse_multiquadric3") == 0)
114  kernel = inverse_multiquadric3;
115  else if (strcmp(s, "sinc_kernel") == 0)
116  kernel = sinc_kernel;
117  else if (strcmp(s, "cosc") == 0)
118  kernel = cosc;
119  else if (strcmp(s, "cot") == 0)
120  kernel = kcot;
121  else if (strcmp(s, "one_over_cube") == 0)
122  kernel = one_over_cube;
123  else if (strcmp(s, "log_sin") == 0)
124  kernel = log_sin;
125  else if (strcmp(s, "laplacian_rbf") == 0)
126  kernel = laplacian_rbf;
127  else
128  {
129  s = "multiquadric";
130  kernel = multiquadric;
131  }
132  }
133  printf(
134  "d=%d, N=%d, M=%d, n=%d, m=%d, p=%d, kernel=%s, c=%" __FGS__ ", eps_I=%" __FGS__ ", eps_B=%" __FGS__ " \n",
135  d, N, M, n, m, p, s, c, eps_I, eps_B);
136 #ifdef NF_KUB
137  printf("nearfield correction using piecewise cubic Lagrange interpolation\n");
138 #elif defined(NF_QUADR)
139  printf("nearfield correction using piecewise quadratic Lagrange interpolation\n");
140 #elif defined(NF_LIN)
141  printf("nearfield correction using piecewise linear Lagrange interpolation\n");
142 #endif
143 
144 #ifdef _OPENMP
145 #pragma omp parallel
146  {
147 #pragma omp single
148  {
149  printf("nthreads=%d\n", omp_get_max_threads());
150  }
151  }
152 
153  FFTW(init_threads)();
154 #endif
155 
157  fastsum_init_guru(&my_fastsum_plan, d, N, M, kernel, &c, 0, n, m, p, eps_I,
158  eps_B);
159  //fastsum_init_guru(&my_fastsum_plan, d, N, M, kernel, &c, NEARFIELD_BOXES, n, m, p, eps_I, eps_B);
160 
161  if (my_fastsum_plan.flags & NEARFIELD_BOXES)
162  printf(
163  "determination of nearfield candidates based on partitioning into boxes\n");
164  else
165  printf("determination of nearfield candidates based on search tree\n");
166 
168  k = 0;
169  while (k < N)
170  {
171  R r_max = K(0.25) - my_fastsum_plan.eps_B / K(2.0);
172  R r2 = K(0.0);
173 
174  for (j = 0; j < d; j++)
175  my_fastsum_plan.x[k * d + j] = K(2.0) * r_max * NFFT(drand48)() - r_max;
176 
177  for (j = 0; j < d; j++)
178  r2 += my_fastsum_plan.x[k * d + j] * my_fastsum_plan.x[k * d + j];
179 
180  if (r2 >= r_max * r_max)
181  continue;
182 
183  k++;
184  }
185 
186  for (k = 0; k < N; k++)
187  {
188  /* R r=(0.25-my_fastsum_plan.eps_B/2.0)*pow((R)rand()/(R)RAND_MAX,1.0/d);
189  my_fastsum_plan.x[k*d+0] = r;
190  for (j=1; j<d; j++)
191  {
192  R phi=2.0*KPI*(R)rand()/(R)RAND_MAX;
193  my_fastsum_plan.x[k*d+j] = r;
194  for (t=0; t<j; t++)
195  {
196  my_fastsum_plan.x[k*d+t] *= cos(phi);
197  }
198  my_fastsum_plan.x[k*d+j] *= sin(phi);
199  }
200  */
201  my_fastsum_plan.alpha[k] = NFFT(drand48)() + II * NFFT(drand48)();
202  }
203 
205  k = 0;
206  while (k < M)
207  {
208  R r_max = K(0.25) - my_fastsum_plan.eps_B / K(2.0);
209  R r2 = K(0.0);
210 
211  for (j = 0; j < d; j++)
212  my_fastsum_plan.y[k * d + j] = K(2.0) * r_max * NFFT(drand48)() - r_max;
213 
214  for (j = 0; j < d; j++)
215  r2 += my_fastsum_plan.y[k * d + j] * my_fastsum_plan.y[k * d + j];
216 
217  if (r2 >= r_max * r_max)
218  continue;
219 
220  k++;
221  }
222  /* for (k=0; k<M; k++)
223  {
224  R r=(0.25-my_fastsum_plan.eps_B/2.0)*pow((R)rand()/(R)RAND_MAX,1.0/d);
225  my_fastsum_plan.y[k*d+0] = r;
226  for (j=1; j<d; j++)
227  {
228  R phi=2.0*KPI*(R)rand()/(R)RAND_MAX;
229  my_fastsum_plan.y[k*d+j] = r;
230  for (t=0; t<j; t++)
231  {
232  my_fastsum_plan.y[k*d+t] *= cos(phi);
233  }
234  my_fastsum_plan.y[k*d+j] *= sin(phi);
235  }
236  } */
237 
239  printf("direct computation: ");
240  fflush(NULL);
241  t0 = getticks();
242  fastsum_exact(&my_fastsum_plan);
243  t1 = getticks();
244  time = NFFT(elapsed_seconds)(t1, t0);
245  printf(__FI__ "sec\n", time);
246 
248  direct = (C *) NFFT(malloc)((size_t)(my_fastsum_plan.M_total) * (sizeof(C)));
249  for (j = 0; j < my_fastsum_plan.M_total; j++)
250  direct[j] = my_fastsum_plan.f[j];
251 
253  printf("pre-computation: ");
254  fflush(NULL);
255  t0 = getticks();
256  fastsum_precompute(&my_fastsum_plan);
257  t1 = getticks();
258  time = NFFT(elapsed_seconds)(t1, t0);
259  printf(__FI__ "sec\n", time);
260 
262  printf("fast computation: ");
263  fflush(NULL);
264  t0 = getticks();
265  fastsum_trafo(&my_fastsum_plan);
266  t1 = getticks();
267  time = NFFT(elapsed_seconds)(t1, t0);
268  printf(__FI__ "sec\n", time);
269 
271  error = K(0.0);
272  for (j = 0; j < my_fastsum_plan.M_total; j++)
273  {
274  if (CABS(direct[j] - my_fastsum_plan.f[j]) / CABS(direct[j]) > error)
275  error = CABS(direct[j] - my_fastsum_plan.f[j]) / CABS(direct[j]);
276  }
277  printf("max relative error: %" __FES__ "\n", error);
278 
280  fastsum_finalize(&my_fastsum_plan);
281 
282  return EXIT_SUCCESS;
283 }
284 /* \} */
int M_total
number of target knots
Definition: fastsum.h:89
C inverse_multiquadric3(R x, int der, const R *param)
K(x) = 1/sqrt(x^2+c^2)^3.
Definition: kernels.c:261
C one_over_cube(R x, int der, const R *param)
K(x) = 1/x^3.
Definition: kernels.c:374
R eps_B
outer boundary
Definition: fastsum.h:114
C gaussian(R x, int der, const R *param)
K(x)=exp(-x^2/c^2)
Definition: kernels.c:38
Header file with predefined kernels for the fast summation algorithm.
C cosc(R x, int der, const R *param)
K(x) = cos(cx)/x.
Definition: kernels.c:314
C one_over_x(R x, int der, const R *param)
K(x) = 1/x.
Definition: kernels.c:233
C thinplate_spline(R x, int der, const R *param)
K(x) = x^2 log |x|.
Definition: kernels.c:149
C sinc_kernel(R x, int der, const R *param)
K(x) = sin(cx)/x.
Definition: kernels.c:287
plan for fast summation algorithm
Definition: fastsum.h:82
Header file for the fast NFFT-based summation algorithm.
C * alpha
source coefficients
Definition: fastsum.h:91
unsigned flags
flags precomp.
Definition: fastsum.h:100
void fastsum_trafo(fastsum_plan *ths)
fast NFFT-based summation
Definition: fastsum.c:1180
void fastsum_precompute(fastsum_plan *ths)
precomputation for fastsum
Definition: fastsum.c:1173
C one_over_square(R x, int der, const R *param)
K(x) = 1/x^2.
Definition: kernels.c:177
C one_over_modulus(R x, int der, const R *param)
K(x) = 1/|x|.
Definition: kernels.c:205
C logarithm(R x, int der, const R *param)
K(x)=log |x|.
Definition: kernels.c:116
C laplacian_rbf(R x, int der, const R *param)
K(x) = exp(-|x|/c)
Definition: kernels.c:417
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:987
C * f
target evaluations
Definition: fastsum.h:92
C inverse_multiquadric(R x, int der, const R *param)
K(x)=1/sqrt(x^2+c^2)
Definition: kernels.c:90
R * y
target knots in d-ball with radius 1/4-eps_b/2
Definition: fastsum.h:95
C kcot(R x, int der, const R *param)
K(x) = cot(cx)
Definition: kernels.c:346
C multiquadric(R x, int der, const R *param)
K(x)=sqrt(x^2+c^2)
Definition: kernels.c:64
void fastsum_finalize(fastsum_plan *ths)
finalization of fastsum plan
Definition: fastsum.c:1048
void fastsum_exact(fastsum_plan *ths)
direct computation of sums
Definition: fastsum.c:1056
C log_sin(R x, int der, const R *param)
K(x) = log(|sin(cx)|)
Definition: kernels.c:402
R * x
source knots in d-ball with radius 1/4-eps_b/2
Definition: fastsum.h:94