NFFT Logo 3.2.2
fastsum_test.c
Go to the documentation of this file.
00001 /*
00002  * Copyright (c) 2002, 2012 Jens Keiner, Stefan Kunis, Daniel Potts
00003  *
00004  * This program is free software; you can redistribute it and/or modify it under
00005  * the terms of the GNU General Public License as published by the Free Software
00006  * Foundation; either version 2 of the License, or (at your option) any later
00007  * version.
00008  *
00009  * This program is distributed in the hope that it will be useful, but WITHOUT
00010  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00011  * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
00012  * details.
00013  *
00014  * You should have received a copy of the GNU General Public License along with
00015  * this program; if not, write to the Free Software Foundation, Inc., 51
00016  * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
00017  */
00018 
00019 /* $Id: fastsum_test.c 3896 2012-10-10 12:19:26Z tovo $ */
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   //fastsum_init_guru(&my_fastsum_plan, d, N, M, kernel, &c, NEARFIELD_BOXES, n, m, p, eps_I, eps_B);
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 /*    double r=(0.25-my_fastsum_plan.eps_B/2.0)*pow((double)rand()/(double)RAND_MAX,1.0/d);
00180     my_fastsum_plan.x[k*d+0] = r;
00181     for (j=1; j<d; j++)
00182     {
00183       double phi=2.0*PI*(double)rand()/(double)RAND_MAX;
00184       my_fastsum_plan.x[k*d+j] = r;
00185       for (t=0; t<j; t++)
00186       {
00187         my_fastsum_plan.x[k*d+t] *= cos(phi);
00188       }
00189       my_fastsum_plan.x[k*d+j] *= sin(phi);
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 /*  for (k=0; k<M; k++)
00214   {
00215     double r=(0.25-my_fastsum_plan.eps_B/2.0)*pow((double)rand()/(double)RAND_MAX,1.0/d);
00216     my_fastsum_plan.y[k*d+0] = r;
00217     for (j=1; j<d; j++)
00218     {
00219       double phi=2.0*PI*(double)rand()/(double)RAND_MAX;
00220       my_fastsum_plan.y[k*d+j] = r;
00221       for (t=0; t<j; t++)
00222       {
00223         my_fastsum_plan.y[k*d+t] *= cos(phi);
00224       }
00225       my_fastsum_plan.y[k*d+j] *= sin(phi);
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 /* \} */

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409