NFFT 3.5.3alpha
fastsum_matlab.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#include "fastsum.h"
36#include "kernels.h"
37#include "infft.h"
38
45int main(int argc, char **argv)
46{
47 int j, k, t;
48 int d;
49 int N;
50 int M;
51 int n;
52 int m;
53 int p;
54 const char *s;
55 C (*kernel)(R, int, const R *);
56 R c;
57 fastsum_plan my_fastsum_plan;
58 C *direct;
59 ticks t0, t1;
60 R time;
61 R error = K(0.0);
62 R eps_I;
63 R eps_B;
64 FILE *fid1, *fid2;
65 R temp;
66
67 if (argc != 11)
68 {
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");
80 exit(-1);
81 }
82 else
83 {
84 d = atoi(argv[1]);
85 N = atoi(argv[2]);
86 c = K(1.0) / POW((R)(N), K(1.0) / ((R)(d)));
87 M = atoi(argv[3]);
88 n = atoi(argv[4]);
89 m = atoi(argv[5]);
90 p = atoi(argv[6]);
91 s = argv[7];
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)
96 kernel = gaussian;
97 else if (strcmp(s, "multiquadric") == 0)
98 kernel = multiquadric;
99 else if (strcmp(s, "inverse_multiquadric") == 0)
100 kernel = inverse_multiquadric;
101 else if (strcmp(s, "logarithm") == 0)
102 kernel = logarithm;
103 else if (strcmp(s, "thinplate_spline") == 0)
104 kernel = thinplate_spline;
105 else if (strcmp(s, "one_over_square") == 0)
106 kernel = one_over_square;
107 else if (strcmp(s, "one_over_modulus") == 0)
108 kernel = one_over_modulus;
109 else if (strcmp(s, "one_over_x") == 0)
110 kernel = one_over_x;
111 else if (strcmp(s, "inverse_multiquadric3") == 0)
112 kernel = inverse_multiquadric3;
113 else if (strcmp(s, "sinc_kernel") == 0)
114 kernel = sinc_kernel;
115 else if (strcmp(s, "cosc") == 0)
116 kernel = cosc;
117 else if (strcmp(s, "cot") == 0)
118 kernel = kcot;
119 else if (strcmp(s, "one_over_cube") == 0)
120 kernel = one_over_cube;
121 else if (strcmp(s, "log_sin") == 0)
122 kernel = log_sin;
123 else if (strcmp(s, "laplacian_rbf") == 0)
124 kernel = laplacian_rbf;
125 else
126 {
127 printf("Unrecognized kernel function!\n");
128 exit(EXIT_FAILURE);
129 }
130 }
131 printf(
132 "d=%d, N=%d, M=%d, n=%d, m=%d, p=%d, kernel=%s, c=%" __FGS__ ", eps_I=%" __FGS__ ", eps_B=%" __FGS__ " \n",
133 d, N, M, n, m, p, s, c, eps_I, eps_B);
134
136 fastsum_init_guru(&my_fastsum_plan, d, N, M, kernel, &c, 0, n, m, p, eps_I,
137 eps_B);
138 /*fastsum_init_guru(&my_fastsum_plan, d, N, M, kernel, &c, EXACT_NEARFIELD, n, m, p);*/
139
141 fid1 = fopen("x.dat", "r");
142 fid2 = fopen("alpha.dat", "r");
143 for (k = 0; k < N; k++)
144 {
145 for (t = 0; t < d; t++)
146 {
147 fscanf(fid1, __FR__, &my_fastsum_plan.x[k * d + t]);
148 }
149 fscanf(fid2, __FR__, &temp);
150 my_fastsum_plan.alpha[k] = temp;
151 fscanf(fid2, __FR__, &temp);
152 my_fastsum_plan.alpha[k] += temp * II;
153 }
154 fclose(fid1);
155 fclose(fid2);
156
158 fid1 = fopen("y.dat", "r");
159 for (j = 0; j < M; j++)
160 {
161 for (t = 0; t < d; t++)
162 {
163 fscanf(fid1, __FR__, &my_fastsum_plan.y[j * d + t]);
164 }
165 }
166 fclose(fid1);
167
169 printf("direct computation: ");
170 fflush(NULL);
171 t0 = getticks();
172 fastsum_exact(&my_fastsum_plan);
173 t1 = getticks();
174 time = NFFT(elapsed_seconds)(t1, t0);
175 printf(__FI__ "sec\n", time);
176
178 direct = (C *) NFFT(malloc)((size_t)(my_fastsum_plan.M_total) * (sizeof(C)));
179 for (j = 0; j < my_fastsum_plan.M_total; j++)
180 direct[j] = my_fastsum_plan.f[j];
181
183 printf("pre-computation: ");
184 fflush(NULL);
185 t0 = getticks();
186 fastsum_precompute(&my_fastsum_plan);
187 t1 = getticks();
188 time = NFFT(elapsed_seconds)(t1, t0);
189 printf(__FI__ "sec\n", time);
190
192 printf("fast computation: ");
193 fflush(NULL);
194 t0 = getticks();
195 fastsum_trafo(&my_fastsum_plan);
196 t1 = getticks();
197 time = NFFT(elapsed_seconds)(t1, t0);
198 printf(__FI__ "sec\n", time);
199
201 error = K(0.0);
202 for (j = 0; j < my_fastsum_plan.M_total; j++)
203 {
204 if (CABS(direct[j] - my_fastsum_plan.f[j]) / CABS(direct[j]) > error)
205 error = CABS(direct[j] - my_fastsum_plan.f[j]) / CABS(direct[j]);
206 }
207 printf("max relative error: " __FE__ "\n", error);
208
210 fid1 = fopen("f.dat", "w+");
211 fid2 = fopen("f_direct.dat", "w+");
212 if (fid1 == NULL)
213 {
214 printf("Error writing to file f.dat!\n");
215 exit(EXIT_FAILURE);
216 }
217 for (j = 0; j < M; j++)
218 {
219 temp = CREAL(my_fastsum_plan.f[j]);
220 fprintf(fid1, " % .16" __FES__ "", temp);
221 temp = CIMAG(my_fastsum_plan.f[j]);
222 fprintf(fid1, " % .16" __FES__ "\n", temp);
223
224 temp = CREAL(direct[j]);
225 fprintf(fid2, " % .16" __FES__ "", temp);
226 temp = CIMAG(direct[j]);
227 fprintf(fid2, " % .16" __FES__ "\n", temp);
228 }
229 fclose(fid1);
230 fclose(fid2);
231
233 fastsum_finalize(&my_fastsum_plan);
234
235 return EXIT_SUCCESS;
236}
237/* \} */
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
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
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