线程安全函数指针gsl monte carlo集成

我想使用gsl monte carlo集成库，我目前正在运行线程安全问题，目前我还无法解决这个问题。线程安全函数指针gsl monte carlo集成

我有一个静态函数gsl_func_wrapper，它包含具有gsl（非成员函数指针）所需的正确类型。这是传递给gsl蒙特卡洛积分器内部的积分成员函数。

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需要这个非成员函数指针，所以我现在并没有真正看到如何让这个线程安全......有人能指出我正确的方向吗？ Thx

编辑：我添加了将由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; 
}

EDIT2： 我加入了呼吁的整体功能由Mike的要求。所以下面是进行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); 
}

来源

2016-02-19 steve

显示'Model2D :: eval'的代码。 –

你可以使用类作为函数指针来避免这个问题，但不是很好。 –

@Pete：Model2D :: eval是一个纯虚拟函数，但对于我来说它是一个非常复杂的函数，所以我不确定是否应该在这里发布所有这些...是否有关联？ – steve

假设Integral被称为在同一个线程，执行gsl_func_wrapper，您可以使用thread_local代替static作为current_model存储类。

这会给你一个单独的current_model变量为每个线程。

来源

2016-02-19 16:31:15 MikeMB

我测试了它，它似乎工作，我会明天测试一下更彻底，并给予更新。非常感谢迈克。 – steve

好的，所以它只是偶尔比平时跑得长一些。它仍然给我nan的积分一次偶尔的结果... – steve

@steve：那么问题很可能不是（只）在你告诉我们的代码中。请提供[mcve]（/ help/mcve）。你正在讨论多线程的问题，但是现在它完全不清楚，你在哪里产生这些线程以及每个线程访问什么数据 – MikeMB

好吧，我设法弄清楚了这个问题。这是这两个变量

const gsl_rng_type *T; gsl_rng *r; gsl_monte_vegas_state *s;

都是实例变量，因此以同样的方式使用的所有线程。如果我在每个线程的Integral（）函数内分别创建这些函数，那么每个线程都会运行良好...

来源

2016-02-21 12:03:10 steve

线程安全函数指针gsl monte carlo集成

回答

相关问题