00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00028 #include <stdlib.h>
00029 #include <stdio.h>
00030 #include <string.h>
00031 #include <math.h>
00032 #include <complex.h>
00033
00034 #include "fastsum.h"
00035 #include "kernels.h"
00036
00043 int main(int argc, char **argv)
00044 {
00045 int j,k,t;
00046 int d;
00047 int N;
00048 int M;
00049 int n;
00050 int m;
00051 int p;
00052 char *s;
00053 double _Complex (*kernel)(double , int , const double *);
00054 double c;
00055 fastsum_plan my_fastsum_plan;
00056 double _Complex *direct;
00057 double time;
00058 double error=0.0;
00059 double eps_I;
00060 double eps_B;
00061 FILE *fid1, *fid2;
00062 double temp;
00063
00064 if (argc!=11)
00065 {
00066 printf("\nfastsum_test d N M n m p kernel c\n\n");
00067 printf(" d dimension \n");
00068 printf(" N number of source nodes \n");
00069 printf(" M number of target nodes \n");
00070 printf(" n expansion degree \n");
00071 printf(" m cut-off parameter \n");
00072 printf(" p degree of smoothness \n");
00073 printf(" kernel kernel function (e.g., gaussian)\n");
00074 printf(" c kernel parameter \n");
00075 printf(" eps_I inner boundary \n");
00076 printf(" eps_B outer boundary \n\n");
00077 exit(-1);
00078 }
00079 else
00080 {
00081 d=atoi(argv[1]);
00082 N=atoi(argv[2]); c=1.0/pow((double)N,1.0/(double)d);
00083 M=atoi(argv[3]);
00084 n=atoi(argv[4]);
00085 m=atoi(argv[5]);
00086 p=atoi(argv[6]);
00087 s=argv[7];
00088 c=atof(argv[8]);
00089 eps_I=atof(argv[9]);
00090 eps_B=atof(argv[10]);
00091 if (strcmp(s,"gaussian")==0)
00092 kernel = gaussian;
00093 else if (strcmp(s,"multiquadric")==0)
00094 kernel = multiquadric;
00095 else if (strcmp(s,"inverse_multiquadric")==0)
00096 kernel = inverse_multiquadric;
00097 else if (strcmp(s,"logarithm")==0)
00098 kernel = logarithm;
00099 else if (strcmp(s,"thinplate_spline")==0)
00100 kernel = thinplate_spline;
00101 else if (strcmp(s,"one_over_square")==0)
00102 kernel = one_over_square;
00103 else if (strcmp(s,"one_over_modulus")==0)
00104 kernel = one_over_modulus;
00105 else if (strcmp(s,"one_over_x")==0)
00106 kernel = one_over_x;
00107 else if (strcmp(s,"inverse_multiquadric3")==0)
00108 kernel = inverse_multiquadric3;
00109 else if (strcmp(s,"sinc_kernel")==0)
00110 kernel = sinc_kernel;
00111 else if (strcmp(s,"cosc")==0)
00112 kernel = cosc;
00113 else if (strcmp(s,"cot")==0)
00114 kernel = kcot;
00115 else
00116 {
00117 s="multiquadric";
00118 kernel = multiquadric;
00119 }
00120 }
00121 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);
00122
00124 fastsum_init_guru(&my_fastsum_plan, d, N, M, kernel, &c, 0, n, m, p, eps_I, eps_B);
00125
00126
00128 fid1=fopen("x.dat","r");
00129 fid2=fopen("alpha.dat","r");
00130 for (k=0; k<N; k++)
00131 {
00132 for (t=0; t<d; t++)
00133 {
00134 fscanf(fid1,"%le",&my_fastsum_plan.x[k*d+t]);
00135 }
00136 fscanf(fid2,"%le",&temp); my_fastsum_plan.alpha[k] = temp;
00137 fscanf(fid2,"%le",&temp); my_fastsum_plan.alpha[k] += temp*_Complex_I;
00138 }
00139 fclose(fid1);
00140 fclose(fid2);
00141
00143 fid1=fopen("y.dat","r");
00144 for (j=0; j<M; j++)
00145 {
00146 for (t=0; t<d; t++)
00147 {
00148 fscanf(fid1,"%le",&my_fastsum_plan.y[j*d+t]);
00149 }
00150 }
00151 fclose(fid1);
00152
00154 printf("direct computation: "); fflush(NULL);
00155 time=nfft_second();
00156 fastsum_exact(&my_fastsum_plan);
00157 time=nfft_second()-time;
00158 printf("%fsec\n",time);
00159
00161 direct = (double _Complex *)nfft_malloc(my_fastsum_plan.M_total*(sizeof(double _Complex)));
00162 for (j=0; j<my_fastsum_plan.M_total; j++)
00163 direct[j]=my_fastsum_plan.f[j];
00164
00166 printf("pre-computation: "); fflush(NULL);
00167 time=nfft_second();
00168 fastsum_precompute(&my_fastsum_plan);
00169 time=nfft_second()-time;
00170 printf("%fsec\n",time);
00171
00173 printf("fast computation: "); fflush(NULL);
00174 time=nfft_second();
00175 fastsum_trafo(&my_fastsum_plan);
00176 time=nfft_second()-time;
00177 printf("%fsec\n",time);
00178
00180 error=0.0;
00181 for (j=0; j<my_fastsum_plan.M_total; j++)
00182 {
00183 if (cabs(direct[j]-my_fastsum_plan.f[j])/cabs(direct[j])>error)
00184 error=cabs(direct[j]-my_fastsum_plan.f[j])/cabs(direct[j]);
00185 }
00186 printf("max relative error: %e\n",error);
00187
00189 fid1=fopen("f.dat","w+");
00190 fid2=fopen("f_direct.dat","w+");
00191 if (fid1==NULL)
00192 {
00193 printf("Fehler!\n");
00194 exit(-1);
00195 }
00196 for (j=0; j<M; j++)
00197 {
00198 temp=creal(my_fastsum_plan.f[j]);
00199 fprintf(fid1," % .16e",temp);
00200 temp=cimag(my_fastsum_plan.f[j]);
00201 fprintf(fid1," % .16e\n",temp);
00202
00203 temp=creal(direct[j]);
00204 fprintf(fid2," % .16e",temp);
00205 temp=cimag(direct[j]);
00206 fprintf(fid2," % .16e\n",temp);
00207 }
00208 fclose(fid1);
00209 fclose(fid2);
00210
00212 fastsum_finalize(&my_fastsum_plan);
00213
00214 return 0;
00215 }
00216