如何将 std::valarray<double> 与 gsl 集成？答案

【问题标题】：How to integrate std::valarray<double> with gsl?如何将 std::valarray<double> 与 gsl 集成？
【发布时间】：2017-12-11 01:29:01
【问题描述】：

我对 C++ 比较陌生，但我有一些（很少）编码和数字经验。

我知道这个问题时不时会发布，你如何集成一个数组。在 MATLAB 中，你可以强制你的数组成为一个函数（我忘了怎么做，但我知道我以前做过）并将它发送给内置的积分器，所以我的问题是你如何在 C++ 中做到这一点。

我有这个积分：

I = integral(A(z)*sin(qz)*dz)

q 只是 double const，z 是积分变量，但 A(z) 是一个数组（我从现在开始将其称为实际函数），它与我的代码中的 z 轴具有相同的点数。积分边界是 z[0] 和 z[nz-1]。

我使用梯形规则计算了这个积分，对于 5000 个点的 z 轴，这需要 0.06 秒。我的问题是这个计算大约发生了 300 * 30 * 20 次（我有 3 个 for 循环），这个 0.06 秒的时间很快增长到模拟的 3 小时。而我的代码的整个瓶颈就是这种集成（我显然可以通过减少 z 来加快速度，但这不是重点。）

我知道库函数通常比用户编写的函数要好得多。我也知道我不能像辛普森规则那样使用更简单的东西，因为被积函数是高度振荡的，我想避免自己实现一些复杂的数值算法。

GSL 需要一个函数形式：

F = f(double x, void *params)

我或许可以使用来自 gsl 的 QAWO 自适应集成，但是如何使我的函数以将我的数组转换为函数的形式？

我的想法是：

F(double z, void *params)
{
  std::valarray<double> actualfunction = *(std::valarray<double> *) params;
  double dz = *(double *) params; // Pretty sure this is wrong
  unsigned int actual_index = z / dz; // crazy assumption (my z[0] was 0)
  return actualfunction[actual_index];
 }

这样的事情可能吗？我怀疑数值算法会使用与实际函数相同的空间差异，然后我应该以某种方式对实际函数进行插值吗？

还有比gsl更好的吗？

【问题讨论】：

标签： c++ integration gsl valarray

【解决方案1】：

template<class F>
struct func_ptr_helper {
  F f;
  void* pvoid(){ return std::addressof(f); }
  template<class R, class...Args>
  using before_ptr=R(*)(void*,Args...);
  template<class R, class...Args>
  using after_ptr=R(*)(Args...,void*);
  template<class R, class...Args>
  static before_ptr<R,Args...> before_func() {
    return [](void* p, Args...args)->R{
      return (*static_cast<F*>(p))(std::forward<Args>(args)...);
    };
  }
  template<class R, class...Args>
  static after_ptr<R,Args...> after_func() {
    return [](Args...args, void* p)->R{
      return (*static_cast<F*>(p))(std::forward<Args>(args)...);
    };
  }
};
template<class F>
func_ptr_helper<F> lambda_to_pfunc( F f ){ return {std::move(f)}; }

使用：

auto f = lambda_to_pfunc([&actualfunction, &dz](double z){
  unsigned int actual_index = z / dz; // crazy assumption (my z[0] was 0)
  return actualfunction[actual_index];
});

然后

void* pvoid - f.pvoid();
void(*pfun)(double, void*) = f.after_func();

你可以通过pfun和pvoid。

对任何错别字深表歉意。

我们的想法是我们编写一个 lambda 来做我们想做的事。然后 lambda_to_pfunc 将其包装起来，以便我们可以将其作为 void* 和指向 C 样式 API 的函数指针传递。

当然，你必须妥善管理所有事物的生命周期。

【讨论】：