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 #ifdef _OPENMP
00038 #include <omp.h>
00039 #endif
00040
00041 #include "fastsum.h"
00042 #include "kernels.h"
00043 #include "infft.h"
00044
00051 int main(int argc, char **argv)
00052 {
00053 int j,k,t;
00054 int d;
00055 int N;
00056 int M;
00057 int n;
00058 int m;
00059 int p;
00060 char *s;
00061 double _Complex (*kernel)(double , int , const double *);
00062 double c;
00063 fastsum_plan my_fastsum_plan;
00064 double _Complex *direct;
00065 ticks t0, t1;
00066 double time;
00067 double error=0.0;
00068 double eps_I;
00069 double eps_B;
00071 if (argc!=11)
00072 {
00073 printf("\nfastsum_test d N M n m p kernel c eps_I eps_B\n\n");
00074 printf(" d dimension \n");
00075 printf(" N number of source nodes \n");
00076 printf(" M number of target nodes \n");
00077 printf(" n expansion degree \n");
00078 printf(" m cut-off parameter \n");
00079 printf(" p degree of smoothness \n");
00080 printf(" kernel kernel function (e.g., gaussian)\n");
00081 printf(" c kernel parameter \n");
00082 printf(" eps_I inner boundary \n");
00083 printf(" eps_B outer boundary \n\n");
00084 exit(-1);
00085 }
00086 else
00087 {
00088 d=atoi(argv[1]);
00089 N=atoi(argv[2]); c=1.0/pow((double)N,1.0/(double)d);
00090 M=atoi(argv[3]);
00091 n=atoi(argv[4]);
00092 m=atoi(argv[5]);
00093 p=atoi(argv[6]);
00094 s=argv[7];
00095 c=atof(argv[8]);
00096 eps_I=atof(argv[9]);
00097 eps_B=atof(argv[10]);
00098 if (strcmp(s,"gaussian")==0)
00099 kernel = gaussian;
00100 else if (strcmp(s,"multiquadric")==0)
00101 kernel = multiquadric;
00102 else if (strcmp(s,"inverse_multiquadric")==0)
00103 kernel = inverse_multiquadric;
00104 else if (strcmp(s,"logarithm")==0)
00105 kernel = logarithm;
00106 else if (strcmp(s,"thinplate_spline")==0)
00107 kernel = thinplate_spline;
00108 else if (strcmp(s,"one_over_square")==0)
00109 kernel = one_over_square;
00110 else if (strcmp(s,"one_over_modulus")==0)
00111 kernel = one_over_modulus;
00112 else if (strcmp(s,"one_over_x")==0)
00113 kernel = one_over_x;
00114 else if (strcmp(s,"inverse_multiquadric3")==0)
00115 kernel = inverse_multiquadric3;
00116 else if (strcmp(s,"sinc_kernel")==0)
00117 kernel = sinc_kernel;
00118 else if (strcmp(s,"cosc")==0)
00119 kernel = cosc;
00120 else if (strcmp(s,"cot")==0)
00121 kernel = kcot;
00122 else
00123 {
00124 s="multiquadric";
00125 kernel = multiquadric;
00126 }
00127 }
00128 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);
00129 #ifdef NF_KUB
00130 printf("nearfield correction using piecewise cubic Lagrange interpolation\n");
00131 #elif defined(NF_QUADR)
00132 printf("nearfield correction using piecewise quadratic Lagrange interpolation\n");
00133 #elif defined(NF_LIN)
00134 printf("nearfield correction using piecewise linear Lagrange interpolation\n");
00135 #endif
00136
00137 #ifdef _OPENMP
00138 #pragma omp parallel
00139 {
00140 #pragma omp single
00141 {
00142 printf("nthreads=%d\n", omp_get_max_threads());
00143 }
00144 }
00145
00146 fftw_init_threads();
00147 #endif
00148
00150 fastsum_init_guru(&my_fastsum_plan, d, N, M, kernel, &c, 0, n, m, p, eps_I, eps_B);
00151
00152
00153 if (my_fastsum_plan.flags & NEARFIELD_BOXES)
00154 printf("determination of nearfield candidates based on partitioning into boxes\n");
00155 else
00156 printf("determination of nearfield candidates based on search tree\n");
00157
00159 k = 0;
00160 while (k < N)
00161 {
00162 double r_max = 0.25 - my_fastsum_plan.eps_B/2.0;
00163 double r2 = 0.0;
00164
00165 for (j=0; j<d; j++)
00166 my_fastsum_plan.x[k*d+j] = 2.0 * r_max * (double)rand()/(double)RAND_MAX - r_max;
00167
00168 for (j=0; j<d; j++)
00169 r2 += my_fastsum_plan.x[k*d+j] * my_fastsum_plan.x[k*d+j];
00170
00171 if (r2 >= r_max * r_max)
00172 continue;
00173
00174 k++;
00175 }
00176
00177 for (k=0; k<N; k++)
00178 {
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192 my_fastsum_plan.alpha[k] = (double)rand()/(double)RAND_MAX + _Complex_I*(double)rand()/(double)RAND_MAX;
00193 }
00194
00196 k = 0;
00197 while (k < M)
00198 {
00199 double r_max = 0.25 - my_fastsum_plan.eps_B/2.0;
00200 double r2 = 0.0;
00201
00202 for (j=0; j<d; j++)
00203 my_fastsum_plan.y[k*d+j] = 2.0 * r_max * (double)rand()/(double)RAND_MAX - r_max;
00204
00205 for (j=0; j<d; j++)
00206 r2 += my_fastsum_plan.y[k*d+j] * my_fastsum_plan.y[k*d+j];
00207
00208 if (r2 >= r_max * r_max)
00209 continue;
00210
00211 k++;
00212 }
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00230 printf("direct computation: "); fflush(NULL);
00231 t0 = getticks();
00232 fastsum_exact(&my_fastsum_plan);
00233 t1 = getticks();
00234 time=nfft_elapsed_seconds(t1,t0);
00235 printf("%fsec\n",time);
00236
00238 direct = (double _Complex *)nfft_malloc(my_fastsum_plan.M_total*(sizeof(double _Complex)));
00239 for (j=0; j<my_fastsum_plan.M_total; j++)
00240 direct[j]=my_fastsum_plan.f[j];
00241
00243 printf("pre-computation: "); fflush(NULL);
00244 t0 = getticks();
00245 fastsum_precompute(&my_fastsum_plan);
00246 t1 = getticks();
00247 time=nfft_elapsed_seconds(t1,t0);
00248 printf("%fsec\n",time);
00249
00251 printf("fast computation: "); fflush(NULL);
00252 t0 = getticks();
00253 fastsum_trafo(&my_fastsum_plan);
00254 t1 = getticks();
00255 time=nfft_elapsed_seconds(t1,t0);
00256 printf("%fsec\n",time);
00257
00259 error=0.0;
00260 for (j=0; j<my_fastsum_plan.M_total; j++)
00261 {
00262 if (cabs(direct[j]-my_fastsum_plan.f[j])/cabs(direct[j])>error)
00263 error=cabs(direct[j]-my_fastsum_plan.f[j])/cabs(direct[j]);
00264 }
00265 printf("max relative error: %e\n",error);
00266
00268 fastsum_finalize(&my_fastsum_plan);
00269
00270 return 0;
00271 }
00272