c - 蒙特卡罗积分优化

Question

我正在使用 gsl 库中实现的 montecarlo 方法。我需要计算这个积分的多次重复，改变被积函数中的参数。所以我需要让我的子程序更快。看来最耗时的部分是在随机点处评估被积函数。在我的具体情况下如何才能更快地进行评估？这是一个最小的例子:

#include <gsl/gsl_rng.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_monte_vegas.h>

double q=0.0;
double mu=0.001;
double eta=0.1;
double kF=1.0;
double Kcut=10;
long int Nmax=10000000;
int Nwu=1000000;
double w=1;


struct my_f_params { double y;};

double 
g (double *k, size_t dim, void *p)
{
  double A;
  struct my_f_params * fp = (struct my_f_params *)p;

  double PQ=q*q+k[1]*k[1]-2*q*k[1]*cos(k[3])+mu;
  double QK=k[0]*k[0]+k[1]*k[1]-2*k[0]*k[1]* (cos(k[2])*cos(k[3])+cos(k[4])*sin(k[2])*sin(k[3]))+mu;
  double KPQ=q*q+k[0]*k[0]+k[1]*k[1]+2*k[0]*cos(k[2])*(q-k[1]*cos(k[3]))+2*k[1]*       (q*cos(k[3])+k[0]*cos(k[4])*sin(k[2])*sin(k[3]));
  double denFreq=fp->y-0.5*(k[0]*k[0]+k[1]*k[1]+KPQ);
  double vol=k[0]*k[0]*k[1]*k[1]*sin(k[2])*sin(k[3]);
  if (sqrt(KPQ) < kF) {
    A = vol*denFreq*(1/QK-1/PQ)/(QK*(pow(denFreq,2)+eta*eta));
  }
  else {
    A = 0;
  }

  return A; 
}


int
main (void)
{
  double res, err;
  double xl[5] = {0, kF, 0, 0, 0};
  double xu[5] = {kF, Kcut, M_PI, M_PI, 2*M_PI};
  const gsl_rng_type *T;
  gsl_rng *r;
  gsl_monte_function G;
  size_t calls = Nmax;

  gsl_rng_env_setup ();
  struct my_f_params params;
  T = gsl_rng_default;
  r = gsl_rng_alloc (T);

  params.y=w;
  G.f=&g;
  G.dim=5;
  G.params=&params;

  {
    gsl_monte_vegas_state *s = gsl_monte_vegas_alloc (5);
    gsl_monte_vegas_integrate (&G, xl, xu, 5, Nwu, r, s,&res, &err);
    do
      {
        gsl_monte_vegas_integrate (&G, xl, xu, 5, calls/5, r, s,&res, &err);
      }
    while (fabs (gsl_monte_vegas_chisq (s) - 1.0) > 0.5);
    gsl_monte_vegas_free (s);
  }

  printf ("%.6f   %.6f   %.6f\n", w,res,err);
  gsl_rng_free (r);



  return 0;
}

Answer 1

扩展 Bob__ 的评论，您可以使用 sincos 来计算同一参数的 sin 和 cos (k[2 ] 和 k[3])，并定义一个要在 main 中初始化的 kF_sqr 并在 g 函数中使用它以避免 sqrt 调用。通过这些优化，在我的机器上进行的快速而肮脏的测试显示，与您的代码相比，速度提高了约 5%。