NFFT 3.5.3alpha
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
49int 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/* \} */
Header file for the fast NFFT-based summation algorithm.
void fastsum_precompute(fastsum_plan *ths)
precomputation for fastsum
Definition fastsum.c:1173
C inverse_multiquadric(R x, int der, const R *param)
K(x)=1/sqrt(x^2+c^2)
Definition kernels.c:90
int M_total
number of target knots
Definition fastsum.h:89
C logarithm(R x, int der, const R *param)
K(x)=log |x|.
Definition kernels.c:116
C multiquadric(R x, int der, const R *param)
K(x)=sqrt(x^2+c^2)
Definition kernels.c:64
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
R * x
source knots in d-ball with radius 1/4-eps_b/2
Definition fastsum.h:94
C one_over_cube(R x, int der, const R *param)
K(x) = 1/x^3.
Definition kernels.c:374
C one_over_square(R x, int der, const R *param)
K(x) = 1/x^2.
Definition kernels.c:177
C sinc_kernel(R x, int der, const R *param)
K(x) = sin(cx)/x.
Definition kernels.c:287
C kcot(R x, int der, const R *param)
K(x) = cot(cx)
Definition kernels.c:346
C one_over_modulus(R x, int der, const R *param)
K(x) = 1/|x|.
Definition kernels.c:205
C * f
target evaluations
Definition fastsum.h:92
C * alpha
source coefficients
Definition fastsum.h:91
R eps_B
outer boundary
Definition fastsum.h:114
C thinplate_spline(R x, int der, const R *param)
K(x) = x^2 log |x|.
Definition kernels.c:149
void fastsum_trafo(fastsum_plan *ths)
fast NFFT-based summation
Definition fastsum.c:1180
void fastsum_exact(fastsum_plan *ths)
direct computation of sums
Definition fastsum.c:1056
unsigned flags
flags precomp.
Definition fastsum.h:100
void fastsum_finalize(fastsum_plan *ths)
finalization of fastsum plan
Definition fastsum.c:1048
C inverse_multiquadric3(R x, int der, const R *param)
K(x) = 1/sqrt(x^2+c^2)^3.
Definition kernels.c:261
C gaussian(R x, int der, const R *param)
K(x)=exp(-x^2/c^2)
Definition kernels.c:38
C one_over_x(R x, int der, const R *param)
K(x) = 1/x.
Definition kernels.c:233
C cosc(R x, int der, const R *param)
K(x) = cos(cx)/x.
Definition kernels.c:314
R * y
target knots in d-ball with radius 1/4-eps_b/2
Definition fastsum.h:95
C log_sin(R x, int der, const R *param)
K(x) = log(|sin(cx)|)
Definition kernels.c:402
C laplacian_rbf(R x, int der, const R *param)
K(x) = exp(-|x|/c)
Definition kernels.c:417
Internal header file for auxiliary definitions and functions.
Header file with predefined kernels for the fast summation algorithm.
plan for fast summation algorithm
Definition fastsum.h:83