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;
00062 if (argc!=11)
00063 {
00064 printf("\nfastsum_test d N M n m p kernel c eps_I eps_B\n\n");
00065 printf(" d dimension \n");
00066 printf(" N number of source nodes \n");
00067 printf(" M number of target nodes \n");
00068 printf(" n expansion degree \n");
00069 printf(" m cut-off parameter \n");
00070 printf(" p degree of smoothness \n");
00071 printf(" kernel kernel function (e.g., gaussian)\n");
00072 printf(" c kernel parameter \n");
00073 printf(" eps_I inner boundary \n");
00074 printf(" eps_B outer boundary \n\n");
00075 exit(-1);
00076 }
00077 else
00078 {
00079 d=atoi(argv[1]);
00080 N=atoi(argv[2]); c=1.0/pow((double)N,1.0/(double)d);
00081 M=atoi(argv[3]);
00082 n=atoi(argv[4]);
00083 m=atoi(argv[5]);
00084 p=atoi(argv[6]);
00085 s=argv[7];
00086 c=atof(argv[8]);
00087 eps_I=atof(argv[9]);
00088 eps_B=atof(argv[10]);
00089 if (strcmp(s,"gaussian")==0)
00090 kernel = gaussian;
00091 else if (strcmp(s,"multiquadric")==0)
00092 kernel = multiquadric;
00093 else if (strcmp(s,"inverse_multiquadric")==0)
00094 kernel = inverse_multiquadric;
00095 else if (strcmp(s,"logarithm")==0)
00096 kernel = logarithm;
00097 else if (strcmp(s,"thinplate_spline")==0)
00098 kernel = thinplate_spline;
00099 else if (strcmp(s,"one_over_square")==0)
00100 kernel = one_over_square;
00101 else if (strcmp(s,"one_over_modulus")==0)
00102 kernel = one_over_modulus;
00103 else if (strcmp(s,"one_over_x")==0)
00104 kernel = one_over_x;
00105 else if (strcmp(s,"inverse_multiquadric3")==0)
00106 kernel = inverse_multiquadric3;
00107 else if (strcmp(s,"sinc_kernel")==0)
00108 kernel = sinc_kernel;
00109 else if (strcmp(s,"cosc")==0)
00110 kernel = cosc;
00111 else if (strcmp(s,"cot")==0)
00112 kernel = kcot;
00113 else
00114 {
00115 s="multiquadric";
00116 kernel = multiquadric;
00117 }
00118 }
00119 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);
00120
00122 fastsum_init_guru(&my_fastsum_plan, d, N, M, kernel, &c, 0, n, m, p, eps_I, eps_B);
00123
00124
00126 for (k=0; k<N; k++)
00127 {
00128 double r=(0.25-my_fastsum_plan.eps_B/2.0)*pow((double)rand()/(double)RAND_MAX,1.0/d);
00129 my_fastsum_plan.x[k*d+0] = r;
00130 for (j=1; j<d; j++)
00131 {
00132 double phi=2.0*PI*(double)rand()/(double)RAND_MAX;
00133 my_fastsum_plan.x[k*d+j] = r;
00134 for (t=0; t<j; t++)
00135 {
00136 my_fastsum_plan.x[k*d+t] *= cos(phi);
00137 }
00138 my_fastsum_plan.x[k*d+j] *= sin(phi);
00139 }
00140
00141 my_fastsum_plan.alpha[k] = (double)rand()/(double)RAND_MAX + _Complex_I*(double)rand()/(double)RAND_MAX;
00142 }
00143
00145 for (k=0; k<M; k++)
00146 {
00147 double r=(0.25-my_fastsum_plan.eps_B/2.0)*pow((double)rand()/(double)RAND_MAX,1.0/d);
00148 my_fastsum_plan.y[k*d+0] = r;
00149 for (j=1; j<d; j++)
00150 {
00151 double phi=2.0*PI*(double)rand()/(double)RAND_MAX;
00152 my_fastsum_plan.y[k*d+j] = r;
00153 for (t=0; t<j; t++)
00154 {
00155 my_fastsum_plan.y[k*d+t] *= cos(phi);
00156 }
00157 my_fastsum_plan.y[k*d+j] *= sin(phi);
00158 }
00159 }
00160
00162 printf("direct computation: "); fflush(NULL);
00163 time=nfft_second();
00164 fastsum_exact(&my_fastsum_plan);
00165 time=nfft_second()-time;
00166 printf("%fsec\n",time);
00167
00169 direct = (double _Complex *)nfft_malloc(my_fastsum_plan.M_total*(sizeof(double _Complex)));
00170 for (j=0; j<my_fastsum_plan.M_total; j++)
00171 direct[j]=my_fastsum_plan.f[j];
00172
00174 printf("pre-computation: "); fflush(NULL);
00175 time=nfft_second();
00176 fastsum_precompute(&my_fastsum_plan);
00177 time=nfft_second()-time;
00178 printf("%fsec\n",time);
00179
00181 printf("fast computation: "); fflush(NULL);
00182 time=nfft_second();
00183 fastsum_trafo(&my_fastsum_plan);
00184 time=nfft_second()-time;
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 fastsum_finalize(&my_fastsum_plan);
00198
00199 return 0;
00200 }
00201