NFFT Logo 3.2.2
fastsum_matlab.c
Go to the documentation of this file.
00001 /*
00002  * Copyright (c) 2002, 2012 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 
00019 /* $Id: fastsum_matlab.c 3775 2012-06-02 16:39:48Z keiner $ */
00020 
00027 #include "config.h"
00028 
00029 #include <stdlib.h>
00030 #include <stdio.h>
00031 #include <string.h>
00032 #include <math.h>
00033 #ifdef HAVE_COMPLEX_H
00034   #include <complex.h>
00035 #endif
00036 
00037 #include "fastsum.h"
00038 #include "kernels.h"
00039 #include "infft.h"
00040 
00047 int main(int argc, char **argv)
00048 {
00049   int j,k,t;                                         
00050   int d;                                             
00051   int N;                                             
00052   int M;                                             
00053   int n;                                             
00054   int m;                                             
00055   int p;                                             
00056   char *s;                                           
00057   double _Complex (*kernel)(double , int , const double *);  
00058   double c;                                          
00059   fastsum_plan my_fastsum_plan;                      
00060   double _Complex *direct;                           
00061   ticks t0, t1;                                      
00062   double time;                                       
00063   double error=0.0;                                  
00064   double eps_I;                                      
00065   double eps_B;                                      
00066   FILE *fid1, *fid2;
00067   double temp;
00068 
00069   if (argc!=11)
00070   {
00071     printf("\nfastsum_test d N M n m p kernel c\n\n");
00072     printf("  d       dimension                 \n");
00073     printf("  N       number of source nodes    \n");
00074     printf("  M       number of target nodes    \n");
00075     printf("  n       expansion degree          \n");
00076     printf("  m       cut-off parameter         \n");
00077     printf("  p       degree of smoothness      \n");
00078     printf("  kernel  kernel function  (e.g., gaussian)\n");
00079     printf("  c       kernel parameter          \n");
00080     printf("  eps_I   inner boundary            \n");
00081     printf("  eps_B   outer boundary            \n\n");
00082     exit(-1);
00083   }
00084   else
00085   {
00086     d=atoi(argv[1]);
00087     N=atoi(argv[2]); c=1.0/pow((double)N,1.0/(double)d);
00088     M=atoi(argv[3]);
00089     n=atoi(argv[4]);
00090     m=atoi(argv[5]);
00091     p=atoi(argv[6]);
00092     s=argv[7];
00093     c=atof(argv[8]);
00094     eps_I=atof(argv[9]);
00095     eps_B=atof(argv[10]);
00096     if (strcmp(s,"gaussian")==0)
00097       kernel = gaussian;
00098     else if (strcmp(s,"multiquadric")==0)
00099       kernel = multiquadric;
00100     else if (strcmp(s,"inverse_multiquadric")==0)
00101       kernel = inverse_multiquadric;
00102     else if (strcmp(s,"logarithm")==0)
00103       kernel = logarithm;
00104     else if (strcmp(s,"thinplate_spline")==0)
00105       kernel = thinplate_spline;
00106     else if (strcmp(s,"one_over_square")==0)
00107       kernel = one_over_square;
00108     else if (strcmp(s,"one_over_modulus")==0)
00109       kernel = one_over_modulus;
00110     else if (strcmp(s,"one_over_x")==0)
00111       kernel = one_over_x;
00112     else if (strcmp(s,"inverse_multiquadric3")==0)
00113       kernel = inverse_multiquadric3;
00114     else if (strcmp(s,"sinc_kernel")==0)
00115       kernel = sinc_kernel;
00116     else if (strcmp(s,"cosc")==0)
00117       kernel = cosc;
00118     else if (strcmp(s,"cot")==0)
00119       kernel = kcot;
00120     else
00121     {
00122       s="multiquadric";
00123       kernel = multiquadric;
00124     }
00125   }
00126   printf("d=%d, N=%d, M=%d, n=%d, m=%d, p=%d, kernel=%s, c=%g, eps_I=%g, eps_B=%g \n",d,N,M,n,m,p,s,c,eps_I,eps_B);
00127 
00129   fastsum_init_guru(&my_fastsum_plan, d, N, M, kernel, &c, 0, n, m, p, eps_I, eps_B);
00130   /*fastsum_init_guru(&my_fastsum_plan, d, N, M, kernel, &c, EXACT_NEARFIELD, n, m, p);*/
00131 
00133   fid1=fopen("x.dat","r");
00134   fid2=fopen("alpha.dat","r");
00135   for (k=0; k<N; k++)
00136   {
00137     for (t=0; t<d; t++)
00138     {
00139       fscanf(fid1,"%le",&my_fastsum_plan.x[k*d+t]);
00140     }
00141     fscanf(fid2,"%le",&temp); my_fastsum_plan.alpha[k] = temp;
00142     fscanf(fid2,"%le",&temp); my_fastsum_plan.alpha[k] += temp*_Complex_I;
00143   }
00144   fclose(fid1);
00145   fclose(fid2);
00146 
00148   fid1=fopen("y.dat","r");
00149   for (j=0; j<M; j++)
00150   {
00151     for (t=0; t<d; t++)
00152     {
00153       fscanf(fid1,"%le",&my_fastsum_plan.y[j*d+t]);
00154     }
00155   }
00156   fclose(fid1);
00157 
00159   printf("direct computation: "); fflush(NULL);
00160   t0 = getticks();
00161   fastsum_exact(&my_fastsum_plan);
00162   t1 = getticks();
00163   time=nfft_elapsed_seconds(t1,t0);
00164   printf("%fsec\n",time);
00165 
00167   direct = (double _Complex *)nfft_malloc(my_fastsum_plan.M_total*(sizeof(double _Complex)));
00168   for (j=0; j<my_fastsum_plan.M_total; j++)
00169     direct[j]=my_fastsum_plan.f[j];
00170 
00172   printf("pre-computation:    "); fflush(NULL);
00173   t0 = getticks();
00174   fastsum_precompute(&my_fastsum_plan);
00175   t1 = getticks();
00176   time=nfft_elapsed_seconds(t1,t0);
00177   printf("%fsec\n",time);
00178 
00180   printf("fast computation:   "); fflush(NULL);
00181   t0 = getticks();
00182   fastsum_trafo(&my_fastsum_plan);
00183   t1 = getticks();
00184   time=nfft_elapsed_seconds(t1,t0);
00185   printf("%fsec\n",time);
00186 
00188   error=0.0;
00189   for (j=0; j<my_fastsum_plan.M_total; j++)
00190   {
00191    if (cabs(direct[j]-my_fastsum_plan.f[j])/cabs(direct[j])>error)
00192       error=cabs(direct[j]-my_fastsum_plan.f[j])/cabs(direct[j]);
00193   }
00194   printf("max relative error: %e\n",error);
00195 
00197   fid1=fopen("f.dat","w+");
00198   fid2=fopen("f_direct.dat","w+");
00199   if (fid1==NULL)
00200   {
00201     printf("Fehler!\n");
00202     exit(-1);
00203   }
00204   for (j=0; j<M; j++)
00205   {
00206     temp=creal(my_fastsum_plan.f[j]);
00207     fprintf(fid1,"  % .16e",temp);
00208     temp=cimag(my_fastsum_plan.f[j]);
00209     fprintf(fid1,"  % .16e\n",temp);
00210 
00211     temp=creal(direct[j]);
00212     fprintf(fid2,"  % .16e",temp);
00213     temp=cimag(direct[j]);
00214     fprintf(fid2,"  % .16e\n",temp);
00215   }
00216   fclose(fid1);
00217   fclose(fid2);
00218 
00220   fastsum_finalize(&my_fastsum_plan);
00221 
00222   return 0;
00223 }
00224 /* \} */

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409