00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
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
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