c++ - 线程安全函数指针 gsl 蒙特卡洛积分

Question

我想使用 gsl monte carlo 集成库，但我目前遇到了迄今为​​止无法解决的线程安全问题。

我有一个静态函数 gsl_func_wrapper，它包装了 gsl 所需的正确类型(非成员函数指针)。它被传递给 Integral 成员函数内的 gsl monte carlo 积分器。

class IntegralStrategyGSL2D: public IntegralStrategy2D {
  static Model2D *current_model;

  const gsl_rng_type *T;
  gsl_rng *r;
  gsl_monte_vegas_state *s;
  size_t calls;
  size_t maxcalls;

  static double gsl_func_wrapper(double *x, size_t dim, void *params) {
    return current_model->eval(x);
  }
...
}

IntegralStrategyGSL2D::IntegralStrategyGSL2D() :
    T(gsl_rng_default) {
  calls = 500;    // keep calls low at first
  maxcalls = 100000;
  gsl_rng_env_setup();
  r = gsl_rng_alloc(T);
  s = gsl_monte_vegas_alloc(2);
}

double IntegralStrategyGSL2D::Integral(Model2D *model2d, double xlow,
    double xhigh, double ylow, double yhigh, double precision) {
  current_model = model2d;
  double result, error;

  double xl[2] = { xlow, ylow };
  double xu[2] = { xhigh, yhigh }

  gsl_monte_function G = { &gsl_func_wrapper, 2, 0 };

  gsl_monte_vegas_init(s);

  while (true) {
    gsl_monte_vegas_integrate(&G, xl, xu, 2, calls, r, s, &result, &error);
    ...
  }
  ...
}

如果我使用 1 个线程运行，所有这些都有效，如果我使用多个线程，我偶尔会得到 nan 积分。由于 gsl 需要这个非成员函数指针，我现在真的不知道如何使这个线程安全......有人能指出我正确的方向吗？谢谢

编辑:我添加了将由 gsl_func_wrapper 调用的 eval 函数。

double GaussianModel2D::eval(const double *x) const {
  // see wikipedia definition

  double normalization = 1.0
        / (2.0 * M_PI * gauss_sigma_var1->getValue()
                * gauss_sigma_var2->getValue()
                * sqrt(1 - pow(gauss_rho->getValue(), 2.0)));

  double exp_value = exp(
        -(pow(x[0] - gauss_mean_var1->getValue(), 2.0)
                / pow(gauss_sigma_var1->getValue(), 2.0)
                + pow(x[1] - gauss_mean_var2->getValue(), 2.0)
                        / pow(gauss_sigma_var2->getValue(), 2.0)
                + 2.0 * gauss_rho->getValue()
                        * (x[0] - gauss_mean_var1->getValue())
                        * (x[1] - gauss_mean_var2->getValue())
                        / (gauss_sigma_var1->getValue()
                                * gauss_sigma_var2->getValue()))
                / (2.0 * (1 - pow(gauss_rho->getValue(), 2.0))));

  return gauss_amplitude->getValue() * normalization * exp_value;
}

编辑2:我按照 Mike 的要求添加了对 Integral 函数的调用。下面是进行 Integral 调用的两个相关代码块。我真的不明白为什么 current_model 变量的 thread_local 关键字可以解决问题。虽然代码有点难看，但 current_model 变量不应该总是指向同一个 model2d 指针吗？因此，即使两个线程可能有写入冲突，写入的值始终相同，所以它应该无关紧要吗？

double smearing_probability = divergence_model->Integral(int_range, 1e-4);

double Model2D::Integral(const std::vector<DataStructs::DimensionRange &ranges,
    double precision) {
    return integral_strategy->Integral(this, ranges[0].range_low,
        ranges[0].range_high, ranges[1].range_low, ranges[1].range_high,
        precision);
}

Answer 1

假设 Integral 在执行 gsl_func_wrapper 的同一线程中调用，您可以使用 thread_local而不是 static 作为 current_model 的存储类。

这将为每个线程提供一个单独的 current_model 变量。