45 int main(
int argc,
char **argv)
55 C (*kernel)(R, int,
const R *);
69 printf(
"\nfastsum_test d N M n m p kernel c\n\n");
70 printf(
" d dimension \n");
71 printf(
" N number of source nodes \n");
72 printf(
" M number of target nodes \n");
73 printf(
" n expansion degree \n");
74 printf(
" m cut-off parameter \n");
75 printf(
" p degree of smoothness \n");
76 printf(
" kernel kernel function (e.g., gaussian)\n");
77 printf(
" c kernel parameter \n");
78 printf(
" eps_I inner boundary \n");
79 printf(
" eps_B outer boundary \n\n");
86 c = K(1.0) / POW((R)(N), K(1.0) / ((R)(d)));
92 c = (R)(atof(argv[8]));
93 eps_I = (R)(atof(argv[9]));
94 eps_B = (R)(atof(argv[10]));
95 if (strcmp(s,
"gaussian") == 0)
97 else if (strcmp(s,
"multiquadric") == 0)
99 else if (strcmp(s,
"inverse_multiquadric") == 0)
101 else if (strcmp(s,
"logarithm") == 0)
103 else if (strcmp(s,
"thinplate_spline") == 0)
105 else if (strcmp(s,
"one_over_square") == 0)
107 else if (strcmp(s,
"one_over_modulus") == 0)
109 else if (strcmp(s,
"one_over_x") == 0)
111 else if (strcmp(s,
"inverse_multiquadric3") == 0)
113 else if (strcmp(s,
"sinc_kernel") == 0)
115 else if (strcmp(s,
"cosc") == 0)
117 else if (strcmp(s,
"cot") == 0)
119 else if (strcmp(s,
"one_over_cube") == 0)
121 else if (strcmp(s,
"log_sin") == 0)
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);
134 fastsum_init_guru(&my_fastsum_plan, d, N, M, kernel, &c, 0, n, m, p, eps_I,
139 fid1 = fopen(
"x.dat",
"r");
140 fid2 = fopen(
"alpha.dat",
"r");
141 for (k = 0; k < N; k++)
143 for (t = 0; t < d; t++)
145 fscanf(fid1, __FR__, &my_fastsum_plan.
x[k * d + t]);
147 fscanf(fid2, __FR__, &temp);
148 my_fastsum_plan.
alpha[k] = temp;
149 fscanf(fid2, __FR__, &temp);
150 my_fastsum_plan.
alpha[k] += temp * II;
156 fid1 = fopen(
"y.dat",
"r");
157 for (j = 0; j < M; j++)
159 for (t = 0; t < d; t++)
161 fscanf(fid1, __FR__, &my_fastsum_plan.
y[j * d + t]);
167 printf(
"direct computation: ");
172 time = NFFT(elapsed_seconds)(t1, t0);
173 printf(__FI__
"sec\n", time);
176 direct = (C *) NFFT(malloc)((size_t)(my_fastsum_plan.
M_total) * (
sizeof(C)));
177 for (j = 0; j < my_fastsum_plan.
M_total; j++)
178 direct[j] = my_fastsum_plan.
f[j];
181 printf(
"pre-computation: ");
186 time = NFFT(elapsed_seconds)(t1, t0);
187 printf(__FI__
"sec\n", time);
190 printf(
"fast computation: ");
195 time = NFFT(elapsed_seconds)(t1, t0);
196 printf(__FI__
"sec\n", time);
200 for (j = 0; j < my_fastsum_plan.
M_total; j++)
202 if (CABS(direct[j] - my_fastsum_plan.
f[j]) / CABS(direct[j]) > error)
203 error = CABS(direct[j] - my_fastsum_plan.
f[j]) / CABS(direct[j]);
205 printf(
"max relative error: " __FE__
"\n", error);
208 fid1 = fopen(
"f.dat",
"w+");
209 fid2 = fopen(
"f_direct.dat",
"w+");
215 for (j = 0; j < M; j++)
217 temp = CREAL(my_fastsum_plan.
f[j]);
218 fprintf(fid1,
" % .16" __FES__
"", temp);
219 temp = CIMAG(my_fastsum_plan.
f[j]);
220 fprintf(fid1,
" % .16" __FES__
"\n", temp);
222 temp = CREAL(direct[j]);
223 fprintf(fid2,
" % .16" __FES__
"", temp);
224 temp = CIMAG(direct[j]);
225 fprintf(fid2,
" % .16" __FES__
"\n", temp);
int M_total
number of target knots
C inverse_multiquadric3(R x, int der, const R *param)
K(x) = 1/sqrt(x^2+c^2)^3.
C one_over_cube(R x, int der, const R *param)
K(x) = 1/x^3.
C gaussian(R x, int der, const R *param)
K(x)=exp(-x^2/c^2)
Header file with predefined kernels for the fast summation algorithm.
C cosc(R x, int der, const R *param)
K(x) = cos(cx)/x.
C one_over_x(R x, int der, const R *param)
K(x) = 1/x.
C thinplate_spline(R x, int der, const R *param)
K(x) = x^2 log |x|.
C sinc_kernel(R x, int der, const R *param)
K(x) = sin(cx)/x.
plan for fast summation algorithm
Header file for the fast NFFT-based summation algorithm.
C * alpha
source coefficients
void fastsum_trafo(fastsum_plan *ths)
fast NFFT-based summation
void fastsum_precompute(fastsum_plan *ths)
precomputation for fastsum
C one_over_square(R x, int der, const R *param)
K(x) = 1/x^2.
C one_over_modulus(R x, int der, const R *param)
K(x) = 1/|x|.
C logarithm(R x, int der, const R *param)
K(x)=log |x|.
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
C inverse_multiquadric(R x, int der, const R *param)
K(x)=1/sqrt(x^2+c^2)
R * y
target knots in d-ball with radius 1/4-eps_b/2
C kcot(R x, int der, const R *param)
K(x) = cot(cx)
C multiquadric(R x, int der, const R *param)
K(x)=sqrt(x^2+c^2)
void fastsum_finalize(fastsum_plan *ths)
finalization of fastsum plan
void fastsum_exact(fastsum_plan *ths)
direct computation of sums
C log_sin(R x, int der, const R *param)
K(x) = log(|sin(cx)|)
R * x
source knots in d-ball with radius 1/4-eps_b/2