NFFT  3.3.1
fastsum_matlab.c
Go to the documentation of this file.
00001 /*
00002  * Copyright (c) 2002, 2016 Jens Keiner, Stefan Kunis, Daniel Potts
00003  *
00004  * This program is free software; you can redistribute it and/or modify it under
00005  * the terms of the GNU General Public License as published by the Free Software
00006  * Foundation; either version 2 of the License, or (at your option) any later
00007  * version.
00008  *
00009  * This program is distributed in the hope that it will be useful, but WITHOUT
00010  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00011  * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
00012  * details.
00013  *
00014  * You should have received a copy of the GNU General Public License along with
00015  * this program; if not, write to the Free Software Foundation, Inc., 51
00016  * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
00017  */
00018 
00025 #include "config.h"
00026 
00027 #include <stdlib.h>
00028 #include <stdio.h>
00029 #include <string.h>
00030 #include <math.h>
00031 #ifdef HAVE_COMPLEX_H
00032   #include <complex.h>
00033 #endif
00034 
00035 #include "fastsum.h"
00036 #include "kernels.h"
00037 #include "infft.h"
00038 
00045 int main(int argc, char **argv)
00046 {
00047   int j, k, t; 
00048   int d; 
00049   int N; 
00050   int M; 
00051   int n; 
00052   int m; 
00053   int p; 
00054   const char *s; 
00055   C (*kernel)(R, int, const R *); 
00056   R c; 
00057   fastsum_plan my_fastsum_plan; 
00058   C *direct; 
00059   ticks t0, t1; 
00060   R time; 
00061   R error = K(0.0); 
00062   R eps_I; 
00063   R eps_B; 
00064   FILE *fid1, *fid2;
00065   R temp;
00066 
00067   if (argc != 11)
00068   {
00069     printf("\nfastsum_test d N M n m p kernel c\n\n");
00070     printf("  d       dimension                 \n");
00071     printf("  N       number of source nodes    \n");
00072     printf("  M       number of target nodes    \n");
00073     printf("  n       expansion degree          \n");
00074     printf("  m       cut-off parameter         \n");
00075     printf("  p       degree of smoothness      \n");
00076     printf("  kernel  kernel function  (e.g., gaussian)\n");
00077     printf("  c       kernel parameter          \n");
00078     printf("  eps_I   inner boundary            \n");
00079     printf("  eps_B   outer boundary            \n\n");
00080     exit(-1);
00081   }
00082   else
00083   {
00084     d = atoi(argv[1]);
00085     N = atoi(argv[2]);
00086     c = K(1.0) / POW((R)(N), K(1.0) / ((R)(d)));
00087     M = atoi(argv[3]);
00088     n = atoi(argv[4]);
00089     m = atoi(argv[5]);
00090     p = atoi(argv[6]);
00091     s = argv[7];
00092     c = (R)(atof(argv[8]));
00093     eps_I = (R)(atof(argv[9]));
00094     eps_B = (R)(atof(argv[10]));
00095     if (strcmp(s, "gaussian") == 0)
00096       kernel = gaussian;
00097     else if (strcmp(s, "multiquadric") == 0)
00098       kernel = multiquadric;
00099     else if (strcmp(s, "inverse_multiquadric") == 0)
00100       kernel = inverse_multiquadric;
00101     else if (strcmp(s, "logarithm") == 0)
00102       kernel = logarithm;
00103     else if (strcmp(s, "thinplate_spline") == 0)
00104       kernel = thinplate_spline;
00105     else if (strcmp(s, "one_over_square") == 0)
00106       kernel = one_over_square;
00107     else if (strcmp(s, "one_over_modulus") == 0)
00108       kernel = one_over_modulus;
00109     else if (strcmp(s, "one_over_x") == 0)
00110       kernel = one_over_x;
00111     else if (strcmp(s, "inverse_multiquadric3") == 0)
00112       kernel = inverse_multiquadric3;
00113     else if (strcmp(s, "sinc_kernel") == 0)
00114       kernel = sinc_kernel;
00115     else if (strcmp(s, "cosc") == 0)
00116       kernel = cosc;
00117     else if (strcmp(s, "cot") == 0)
00118       kernel = kcot;
00119     else
00120     {
00121       s = "multiquadric";
00122       kernel = multiquadric;
00123     }
00124   }
00125   printf(
00126       "d=%d, N=%d, M=%d, n=%d, m=%d, p=%d, kernel=%s, c=%" __FGS__ ", eps_I=%" __FGS__ ", eps_B=%" __FGS__ " \n",
00127       d, N, M, n, m, p, s, c, eps_I, eps_B);
00128 
00130   fastsum_init_guru(&my_fastsum_plan, d, N, M, kernel, &c, 0, n, m, p, eps_I,
00131       eps_B);
00132   /*fastsum_init_guru(&my_fastsum_plan, d, N, M, kernel, &c, EXACT_NEARFIELD, n, m, p);*/
00133 
00135   fid1 = fopen("x.dat", "r");
00136   fid2 = fopen("alpha.dat", "r");
00137   for (k = 0; k < N; k++)
00138   {
00139     for (t = 0; t < d; t++)
00140     {
00141       fscanf(fid1, __FR__, &my_fastsum_plan.x[k * d + t]);
00142     }
00143     fscanf(fid2, __FR__, &temp);
00144     my_fastsum_plan.alpha[k] = temp;
00145     fscanf(fid2, __FR__, &temp);
00146     my_fastsum_plan.alpha[k] += temp * II;
00147   }
00148   fclose(fid1);
00149   fclose(fid2);
00150 
00152   fid1 = fopen("y.dat", "r");
00153   for (j = 0; j < M; j++)
00154   {
00155     for (t = 0; t < d; t++)
00156     {
00157       fscanf(fid1, __FR__, &my_fastsum_plan.y[j * d + t]);
00158     }
00159   }
00160   fclose(fid1);
00161 
00163   printf("direct computation: ");
00164   fflush(NULL);
00165   t0 = getticks();
00166   fastsum_exact(&my_fastsum_plan);
00167   t1 = getticks();
00168   time = NFFT(elapsed_seconds)(t1, t0);
00169   printf(__FI__ "sec\n", time);
00170 
00172   direct = (C *) NFFT(malloc)((size_t)(my_fastsum_plan.M_total) * (sizeof(C)));
00173   for (j = 0; j < my_fastsum_plan.M_total; j++)
00174     direct[j] = my_fastsum_plan.f[j];
00175 
00177   printf("pre-computation:    ");
00178   fflush(NULL);
00179   t0 = getticks();
00180   fastsum_precompute(&my_fastsum_plan);
00181   t1 = getticks();
00182   time = NFFT(elapsed_seconds)(t1, t0);
00183   printf(__FI__ "sec\n", time);
00184 
00186   printf("fast computation:   ");
00187   fflush(NULL);
00188   t0 = getticks();
00189   fastsum_trafo(&my_fastsum_plan);
00190   t1 = getticks();
00191   time = NFFT(elapsed_seconds)(t1, t0);
00192   printf(__FI__ "sec\n", time);
00193 
00195   error = K(0.0);
00196   for (j = 0; j < my_fastsum_plan.M_total; j++)
00197   {
00198     if (CABS(direct[j] - my_fastsum_plan.f[j]) / CABS(direct[j]) > error)
00199       error = CABS(direct[j] - my_fastsum_plan.f[j]) / CABS(direct[j]);
00200   }
00201   printf("max relative error: " __FE__ "\n", error);
00202 
00204   fid1 = fopen("f.dat", "w+");
00205   fid2 = fopen("f_direct.dat", "w+");
00206   if (fid1 == NULL)
00207   {
00208     printf("Fehler!\n");
00209     exit(EXIT_FAILURE);
00210   }
00211   for (j = 0; j < M; j++)
00212   {
00213     temp = CREAL(my_fastsum_plan.f[j]);
00214     fprintf(fid1, "  % .16" __FES__ "", temp);
00215     temp = CIMAG(my_fastsum_plan.f[j]);
00216     fprintf(fid1, "  % .16" __FES__ "\n", temp);
00217 
00218     temp = CREAL(direct[j]);
00219     fprintf(fid2, "  % .16" __FES__ "", temp);
00220     temp = CIMAG(direct[j]);
00221     fprintf(fid2, "  % .16" __FES__ "\n", temp);
00222   }
00223   fclose(fid1);
00224   fclose(fid2);
00225 
00227   fastsum_finalize(&my_fastsum_plan);
00228 
00229   return EXIT_SUCCESS;
00230 }
00231 /* \} */