00001
00008 #include <stdlib.h>
00009 #include <stdio.h>
00010 #include <string.h>
00011 #include <complex.h>
00012 #include <math.h>
00013
00014 #include "fastsum.h"
00015 #include "kernels.h"
00016
00023 int main(int argc, char **argv)
00024 {
00025 int j,k,t;
00026 int d;
00027 int N;
00028 int M;
00029 int n;
00030 int m;
00031 int p;
00032 char *s;
00033 complex (*kernel)(double , int , const double *);
00034 double c;
00035 fastsum_plan my_fastsum_plan;
00036 complex *direct;
00037 double time;
00038 double error=0.0;
00039 double eps_I;
00040 double eps_B;
00042 if (argc!=11)
00043 {
00044 printf("\nfastsum_test d N M n m p kernel c eps_I eps_B\n\n");
00045 printf(" d dimension \n");
00046 printf(" N number of source nodes \n");
00047 printf(" M number of target nodes \n");
00048 printf(" n expansion degree \n");
00049 printf(" m cut-off parameter \n");
00050 printf(" p degree of smoothness \n");
00051 printf(" kernel kernel function (e.g., gaussian)\n");
00052 printf(" c kernel parameter \n");
00053 printf(" eps_I inner boundary \n");
00054 printf(" eps_B outer boundary \n\n");
00055 exit(-1);
00056 }
00057 else
00058 {
00059 d=atoi(argv[1]);
00060 N=atoi(argv[2]); c=1.0/pow((double)N,1.0/(double)d);
00061 M=atoi(argv[3]);
00062 n=atoi(argv[4]);
00063 m=atoi(argv[5]);
00064 p=atoi(argv[6]);
00065 s=argv[7];
00066 c=atof(argv[8]);
00067 eps_I=atof(argv[9]);
00068 eps_B=atof(argv[10]);
00069 if (strcmp(s,"gaussian")==0)
00070 kernel = gaussian;
00071 else if (strcmp(s,"multiquadric")==0)
00072 kernel = multiquadric;
00073 else if (strcmp(s,"inverse_multiquadric")==0)
00074 kernel = inverse_multiquadric;
00075 else if (strcmp(s,"logarithm")==0)
00076 kernel = logarithm;
00077 else if (strcmp(s,"thinplate_spline")==0)
00078 kernel = thinplate_spline;
00079 else if (strcmp(s,"one_over_square")==0)
00080 kernel = one_over_square;
00081 else if (strcmp(s,"one_over_modulus")==0)
00082 kernel = one_over_modulus;
00083 else if (strcmp(s,"one_over_x")==0)
00084 kernel = one_over_x;
00085 else if (strcmp(s,"inverse_multiquadric3")==0)
00086 kernel = inverse_multiquadric3;
00087 else if (strcmp(s,"sinc_kernel")==0)
00088 kernel = sinc_kernel;
00089 else if (strcmp(s,"cosc")==0)
00090 kernel = cosc;
00091 else if (strcmp(s,"cot")==0)
00092 kernel = cot;
00093 else
00094 {
00095 s="multiquadric";
00096 kernel = multiquadric;
00097 }
00098 }
00099 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);
00100
00102 fastsum_init_guru(&my_fastsum_plan, d, N, M, kernel, &c, 0, n, m, p, eps_I, eps_B);
00103
00104
00106 for (k=0; k<N; k++)
00107 {
00108 double r=(0.25-my_fastsum_plan.eps_B/2.0)*pow((double)rand()/(double)RAND_MAX,1.0/d);
00109 my_fastsum_plan.x[k*d+0] = r;
00110 for (j=1; j<d; j++)
00111 {
00112 double phi=2.0*PI*(double)rand()/(double)RAND_MAX;
00113 my_fastsum_plan.x[k*d+j] = r;
00114 for (t=0; t<j; t++)
00115 {
00116 my_fastsum_plan.x[k*d+t] *= cos(phi);
00117 }
00118 my_fastsum_plan.x[k*d+j] *= sin(phi);
00119 }
00120
00121 my_fastsum_plan.alpha[k] = (double)rand()/(double)RAND_MAX + I*(double)rand()/(double)RAND_MAX;
00122 }
00123
00125 for (k=0; k<M; k++)
00126 {
00127 double r=(0.25-my_fastsum_plan.eps_B/2.0)*pow((double)rand()/(double)RAND_MAX,1.0/d);
00128 my_fastsum_plan.y[k*d+0] = r;
00129 for (j=1; j<d; j++)
00130 {
00131 double phi=2.0*PI*(double)rand()/(double)RAND_MAX;
00132 my_fastsum_plan.y[k*d+j] = r;
00133 for (t=0; t<j; t++)
00134 {
00135 my_fastsum_plan.y[k*d+t] *= cos(phi);
00136 }
00137 my_fastsum_plan.y[k*d+j] *= sin(phi);
00138 }
00139 }
00140
00142 printf("direct computation: "); fflush(NULL);
00143 time=nfft_second();
00144 fastsum_exact(&my_fastsum_plan);
00145 time=nfft_second()-time;
00146 printf("%fsec\n",time);
00147
00149 direct = (complex *)malloc(my_fastsum_plan.M_total*(sizeof(complex)));
00150 for (j=0; j<my_fastsum_plan.M_total; j++)
00151 direct[j]=my_fastsum_plan.f[j];
00152
00154 printf("pre-computation: "); fflush(NULL);
00155 time=nfft_second();
00156 fastsum_precompute(&my_fastsum_plan);
00157 time=nfft_second()-time;
00158 printf("%fsec\n",time);
00159
00161 printf("fast computation: "); fflush(NULL);
00162 time=nfft_second();
00163 fastsum_trafo(&my_fastsum_plan);
00164 time=nfft_second()-time;
00165 printf("%fsec\n",time);
00166
00168 error=0.0;
00169 for (j=0; j<my_fastsum_plan.M_total; j++)
00170 {
00171 if (cabs(direct[j]-my_fastsum_plan.f[j])/cabs(direct[j])>error)
00172 error=cabs(direct[j]-my_fastsum_plan.f[j])/cabs(direct[j]);
00173 }
00174 printf("max relative error: %e\n",error);
00175
00177 fastsum_finalize(&my_fastsum_plan);
00178
00179 return 0;
00180 }
00181
