51 int main(
int argc,
char **argv)
61 C (*kernel)(R, int,
const R *);
73 printf(
"\nfastsum_test d N M n m p kernel c eps_I eps_B\n\n");
74 printf(
" d dimension \n");
75 printf(
" N number of source nodes \n");
76 printf(
" M number of target nodes \n");
77 printf(
" n expansion degree \n");
78 printf(
" m cut-off parameter \n");
79 printf(
" p degree of smoothness \n");
80 printf(
" kernel kernel function (e.g., gaussian)\n");
81 printf(
" c kernel parameter \n");
82 printf(
" eps_I inner boundary \n");
83 printf(
" eps_B outer boundary \n\n");
90 c = K(1.0) / POW((R)(N), K(1.0) / ((R)(d)));
96 c = (R)(atof(argv[8]));
97 eps_I = (R)(atof(argv[9]));
98 eps_B = (R)(atof(argv[10]));
99 if (strcmp(s,
"gaussian") == 0)
101 else if (strcmp(s,
"multiquadric") == 0)
102 kernel = multiquadric;
103 else if (strcmp(s,
"inverse_multiquadric") == 0)
104 kernel = inverse_multiquadric;
105 else if (strcmp(s,
"logarithm") == 0)
107 else if (strcmp(s,
"thinplate_spline") == 0)
108 kernel = thinplate_spline;
109 else if (strcmp(s,
"one_over_square") == 0)
110 kernel = one_over_square;
111 else if (strcmp(s,
"one_over_modulus") == 0)
112 kernel = one_over_modulus;
113 else if (strcmp(s,
"one_over_x") == 0)
115 else if (strcmp(s,
"inverse_multiquadric3") == 0)
116 kernel = inverse_multiquadric3;
117 else if (strcmp(s,
"sinc_kernel") == 0)
118 kernel = sinc_kernel;
119 else if (strcmp(s,
"cosc") == 0)
121 else if (strcmp(s,
"cot") == 0)
126 kernel = multiquadric;
130 "d=%d, N=%d, M=%d, n=%d, m=%d, p=%d, kernel=%s, c=%" __FGS__
", eps_I=%" __FGS__
", eps_B=%" __FGS__
" \n",
131 d, N, M, n, m, p, s, c, eps_I, eps_B);
133 printf(
"nearfield correction using piecewise cubic Lagrange interpolation\n");
134 #elif defined(NF_QUADR) 135 printf(
"nearfield correction using piecewise quadratic Lagrange interpolation\n");
136 #elif defined(NF_LIN) 137 printf(
"nearfield correction using piecewise linear Lagrange interpolation\n");
145 printf(
"nthreads=%d\n", omp_get_max_threads());
149 FFTW(init_threads)();
153 fastsum_init_guru(&my_fastsum_plan, d, N, M, kernel, &c, 0, n, m, p, eps_I,
157 if (my_fastsum_plan.
flags & NEARFIELD_BOXES)
159 "determination of nearfield candidates based on partitioning into boxes\n");
161 printf(
"determination of nearfield candidates based on search tree\n");
167 R r_max = K(0.25) - my_fastsum_plan.
eps_B / K(2.0);
170 for (j = 0; j < d; j++)
171 my_fastsum_plan.
x[k * d + j] = K(2.0) * r_max * NFFT(drand48)() - r_max;
173 for (j = 0; j < d; j++)
174 r2 += my_fastsum_plan.
x[k * d + j] * my_fastsum_plan.
x[k * d + j];
176 if (r2 >= r_max * r_max)
182 for (k = 0; k < N; k++)
197 my_fastsum_plan.
alpha[k] = NFFT(drand48)() + II * NFFT(drand48)();
204 R r_max = K(0.25) - my_fastsum_plan.
eps_B / K(2.0);
207 for (j = 0; j < d; j++)
208 my_fastsum_plan.
y[k * d + j] = K(2.0) * r_max * NFFT(drand48)() - r_max;
210 for (j = 0; j < d; j++)
211 r2 += my_fastsum_plan.
y[k * d + j] * my_fastsum_plan.
y[k * d + j];
213 if (r2 >= r_max * r_max)
235 printf(
"direct computation: ");
240 time = NFFT(elapsed_seconds)(t1, t0);
241 printf(__FI__
"sec\n", time);
244 direct = (C *) NFFT(malloc)((size_t)(my_fastsum_plan.
M_total) * (
sizeof(C)));
245 for (j = 0; j < my_fastsum_plan.
M_total; j++)
246 direct[j] = my_fastsum_plan.
f[j];
249 printf(
"pre-computation: ");
254 time = NFFT(elapsed_seconds)(t1, t0);
255 printf(__FI__
"sec\n", time);
258 printf(
"fast computation: ");
263 time = NFFT(elapsed_seconds)(t1, t0);
264 printf(__FI__
"sec\n", time);
268 for (j = 0; j < my_fastsum_plan.
M_total; j++)
270 if (CABS(direct[j] - my_fastsum_plan.
f[j]) / CABS(direct[j]) > error)
271 error = CABS(direct[j] - my_fastsum_plan.
f[j]) / CABS(direct[j]);
273 printf(
"max relative error: %" __FES__
"\n", error);
int M_total
number of target knots
Header file with predefined kernels for the fast summation algorithm.
plan for fast summation algorithm
Header file for the fast NFFT-based summation algorithm.
C * alpha
source coefficients
unsigned flags
flags precomp.
void fastsum_trafo(fastsum_plan *ths)
fast NFFT-based summation
void fastsum_precompute(fastsum_plan *ths)
precomputation for fastsum
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
R * y
target knots in d-ball with radius 1/4-eps_b/2
void fastsum_finalize(fastsum_plan *ths)
finalization of fastsum plan
void fastsum_exact(fastsum_plan *ths)
direct computation of sums
R * x
source knots in d-ball with radius 1/4-eps_b/2