NFFT  3.3.2alpha
fastsum_matlab.c
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2002, 2016 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 
45 int 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
120  {
121  s = "multiquadric";
122  kernel = multiquadric;
123  }
124  }
125  printf(
126  "d=%d, N=%d, M=%d, n=%d, m=%d, p=%d, kernel=%s, c=%" __FGS__ ", eps_I=%" __FGS__ ", eps_B=%" __FGS__ " \n",
127  d, N, M, n, m, p, s, c, eps_I, eps_B);
128 
130  fastsum_init_guru(&my_fastsum_plan, d, N, M, kernel, &c, 0, n, m, p, eps_I,
131  eps_B);
132  /*fastsum_init_guru(&my_fastsum_plan, d, N, M, kernel, &c, EXACT_NEARFIELD, n, m, p);*/
133 
135  fid1 = fopen("x.dat", "r");
136  fid2 = fopen("alpha.dat", "r");
137  for (k = 0; k < N; k++)
138  {
139  for (t = 0; t < d; t++)
140  {
141  fscanf(fid1, __FR__, &my_fastsum_plan.x[k * d + t]);
142  }
143  fscanf(fid2, __FR__, &temp);
144  my_fastsum_plan.alpha[k] = temp;
145  fscanf(fid2, __FR__, &temp);
146  my_fastsum_plan.alpha[k] += temp * II;
147  }
148  fclose(fid1);
149  fclose(fid2);
150 
152  fid1 = fopen("y.dat", "r");
153  for (j = 0; j < M; j++)
154  {
155  for (t = 0; t < d; t++)
156  {
157  fscanf(fid1, __FR__, &my_fastsum_plan.y[j * d + t]);
158  }
159  }
160  fclose(fid1);
161 
163  printf("direct computation: ");
164  fflush(NULL);
165  t0 = getticks();
166  fastsum_exact(&my_fastsum_plan);
167  t1 = getticks();
168  time = NFFT(elapsed_seconds)(t1, t0);
169  printf(__FI__ "sec\n", time);
170 
172  direct = (C *) NFFT(malloc)((size_t)(my_fastsum_plan.M_total) * (sizeof(C)));
173  for (j = 0; j < my_fastsum_plan.M_total; j++)
174  direct[j] = my_fastsum_plan.f[j];
175 
177  printf("pre-computation: ");
178  fflush(NULL);
179  t0 = getticks();
180  fastsum_precompute(&my_fastsum_plan);
181  t1 = getticks();
182  time = NFFT(elapsed_seconds)(t1, t0);
183  printf(__FI__ "sec\n", time);
184 
186  printf("fast computation: ");
187  fflush(NULL);
188  t0 = getticks();
189  fastsum_trafo(&my_fastsum_plan);
190  t1 = getticks();
191  time = NFFT(elapsed_seconds)(t1, t0);
192  printf(__FI__ "sec\n", time);
193 
195  error = K(0.0);
196  for (j = 0; j < my_fastsum_plan.M_total; j++)
197  {
198  if (CABS(direct[j] - my_fastsum_plan.f[j]) / CABS(direct[j]) > error)
199  error = CABS(direct[j] - my_fastsum_plan.f[j]) / CABS(direct[j]);
200  }
201  printf("max relative error: " __FE__ "\n", error);
202 
204  fid1 = fopen("f.dat", "w+");
205  fid2 = fopen("f_direct.dat", "w+");
206  if (fid1 == NULL)
207  {
208  printf("Fehler!\n");
209  exit(EXIT_FAILURE);
210  }
211  for (j = 0; j < M; j++)
212  {
213  temp = CREAL(my_fastsum_plan.f[j]);
214  fprintf(fid1, " % .16" __FES__ "", temp);
215  temp = CIMAG(my_fastsum_plan.f[j]);
216  fprintf(fid1, " % .16" __FES__ "\n", temp);
217 
218  temp = CREAL(direct[j]);
219  fprintf(fid2, " % .16" __FES__ "", temp);
220  temp = CIMAG(direct[j]);
221  fprintf(fid2, " % .16" __FES__ "\n", temp);
222  }
223  fclose(fid1);
224  fclose(fid2);
225 
227  fastsum_finalize(&my_fastsum_plan);
228 
229  return EXIT_SUCCESS;
230 }
231 /* \} */
int M_total
number of target knots
Definition: fastsum.h:79
Header file with predefined kernels for the fast summation algorithm.
plan for fast summation algorithm
Definition: fastsum.h:72
Header file for the fast NFFT-based summation algorithm.
C * alpha
source coefficients
Definition: fastsum.h:81
void fastsum_trafo(fastsum_plan *ths)
fast NFFT-based summation
Definition: fastsum.c:1055
void fastsum_precompute(fastsum_plan *ths)
precomputation for fastsum
Definition: fastsum.c:907
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:691
C * f
target evaluations
Definition: fastsum.h:82
R * y
target knots in d-ball with radius 1/4-eps_b/2
Definition: fastsum.h:85
void fastsum_finalize(fastsum_plan *ths)
finalization of fastsum plan
Definition: fastsum.c:844
void fastsum_exact(fastsum_plan *ths)
direct computation of sums
Definition: fastsum.c:877
R * x
source knots in d-ball with radius 1/4-eps_b/2
Definition: fastsum.h:84