函数c ++中的gsl集成

Question

我只是想在 C++ 中集成一个函数。我一直在尝试使用 gsl，因为我在网上看到了这个推荐。我遵循了 gsl 示例，但收效甚微。

这是我的 C++ 代码：

double inverseE(double z){
   double inverseE = 1.0/(std::sqrt(Om0*std::pow(1.0+z,3.0)+1.0-Om0));
   return inverseE;
}

double comoving_distance(double z){
   gsl_integration_workspace * w
     = gsl_integration_workspace_alloc (1000);

   double result, error;

   gsl_function iE;
   iE.function = &inverseE;

   gsl_integration_qags (&iE, 0, z, 0, 1e-7, 1000,
                           w, &result, &error);

   gsl_integration_workspace_free (w);

   cout << result << endl;

   return 0;
}

为了澄清，Python 中的相同代码（有效）如下所示：

def iE(z):
     return 1/(np.sqrt(Om0*np.power(1+z,3)+1-Om0))

def comoving_distance(z):
     return (c/H0)*quad(iE,0,z)[0]

quad 执行集成的地方（它是一个 scipy 模块）。

我收到两条错误消息：

ISO C++ 禁止使用未限定或带括号的非静态成员函数的地址来构成指向成员函数的指针。说‘&cosmo::inverseE’[-fpermissive]

无法将'double (cosmo::)(double)'转换为'double ()(double, void)'任务* cosmo 是包含这两个函数的 class 的名称。

看来这应该不是什么难事。关于我哪里出错的建议将不胜感激！

编辑：class

#include <iostream>    // using IO functions
#include <string>      // using string
#include <gsl/gsl_integration.h>
#include <cmath>
using namespace std;



class cosmo {
private:
   double H0;
   double Om0;
   double Ob0;
   double c;
   double pi;
   double G;

public:
   // Constructor with default values for data members
   cosmo(double Hubble0 = 70, double OmegaM0 = 0.3,
           double OmegaB0 = 0.05) {
         H0 = Hubble0;
         Om0 = OmegaM0;
         Ob0 = OmegaB0;
         c = 3e8;
         pi = 3.141592653589793;
         G = 6.67408e-11;

}

double inverseE(double z){
   double inverseE = 1.0/(std::sqrt(Om0*std::pow(1.0+z,3.0)+1.0-Om0));
   return inverseE;
}

double comoving_distance(double z){
   gsl_integration_workspace * w
     = gsl_integration_workspace_alloc (1000);

   double result, error;

   gsl_function iE;
   iE.function = &inverseE;

   gsl_integration_qags (&iE, 0, z, 0, 1e-7, 1000,
                           w, &result, &error);

   gsl_integration_workspace_free (w);

   cout << result << endl;

   return 0;
}


};

Answer 1

我看到两个问题：

1/ GSL 期望您的 inverseE 函数具有以下原型

double inverseE(double z,void *use_data);

在您的代码中您已声明：

double inverseE(double z);

2/ 就像您的代码在 C++ 中并且 GSL 是一个 C 库，您必须使您的 C++ 函数可从 C 调用。

解决方案是声明您的 inverseE 函数如下：

extern "C" {
  double inverseE(double z,void *) {
      double inverseE = 1.0/(std::sqrt(Om0*std::pow(1.0+z,3.0)+1.0-Om0));
      return inverseE;
   }
}

这 extern "C" 使您的 C++ 函数与 C 调用约定二进制兼容。

经过这两个修改，我认为你的代码应该没问题。

更新：2/

在我的回答中，我在 2/ 中认为 inverseE 是一个函数。这里我考虑的是你的 class.

的 method 的情况

这是 void *user_data 派上用场的示例：

声明这个包装器：

  extern "C" {
  double YourClass_f_wrap(double z, void *user_data)
  {
    YourClass *this_ptr = (YourClass *)user_data;
    return this_ptr->f(z);
  }
  }

那么YourClass定义如下：

class YourClass
{
 public:
  struct IntegrationResult
  {
    double result, error;
    size_t n_intervals;
  };

 public:
  double f(double x) const; // defines f to integrate
  IntegrationResult integrate_f() const; // integrates f using GSL lib

  ...
};

正如您在评论中提到的，必须谨慎对待前向声明。为了清楚起见，请在下面找到一个完整的运行示例，该示例重现了 the GSL official doc 的结果，但使用 C++ class 和方法 f()

完成运行代码：

可以编译为：

g++ gslIntegrationExample.cpp -lgsl -lcblas -o gslIntegrationExample

代码：

#include <cassert>
#include <cstdio>
#include <gsl/gsl_integration.h>

namespace details
{
  extern "C" {
  double YourClass_f_wrap(double z, void *user_data);
  }
}

class YourClass
{
 public:
  struct IntegrationResult
  {
    double result, error;
    size_t n_intervals;
  };

 public:
  double f(double x) const
  {
    return std::log(alpha_ * x) / std::sqrt(x);
  }

  IntegrationResult integrate_f() const
  {
    gsl_integration_workspace *w =
        gsl_integration_workspace_alloc(1000);

    assert(w != nullptr);

    IntegrationResult toReturn;

    gsl_function F;
    F.function = &details::YourClass_f_wrap;
    F.params = (void *)this;

    gsl_integration_qags(&F, 0, 1, 0, 1e-7, 1000, w, &toReturn.result,
                         &toReturn.error);

    toReturn.n_intervals = w->size;

    gsl_integration_workspace_free(w);

    return toReturn;
  }

 protected:
  double alpha_ = 1;
};

namespace details
{
  extern "C" {
  double YourClass_f_wrap(double z, void *user_data)
  {
    YourClass *this_ptr = (YourClass *)user_data;
    return this_ptr->f(z);
  }
  }
}

int main()
{
  YourClass yourClass;

  auto integrationResult = yourClass.integrate_f();

  double expected = -4.0;
  std::printf("result          = % .18f\n", integrationResult.result);
  std::printf("exact result    = % .18f\n", expected);
  std::printf("estimated error = % .18f\n", integrationResult.error);
  std::printf("actual error    = % .18f\n",
              integrationResult.result - expected);
  std::printf("intervals       = %zu\n",
              integrationResult.n_intervals);

  return EXIT_SUCCESS;
}

在我的电脑上我得到：

result          = -4.000000000000085265
exact result    = -4.000000000000000000
estimated error =  0.000000000000135447
actual error    = -0.000000000000085265
intervals       = 8

Answer 2

阅读一些评论，我终于 post 另一个包含较少样板代码的版本。

它的工作原理如下：

有一个 extern "C" gsl_function_cpp_wrapper_helper 函数来确保 C ABI（应用程序二进制接口）。这是 gsl 的子例程将使用的函数。
由于 std::function，有一个 C++ class gsl_function_cpp_wrapper 包装了 C++ function/class 方法。此包装器可转换为使用 gsl_function_cpp_wrapper_helper 作为 C 函数的 gsl_function_struct 指针，由于 gsl 的 void *params 参数
[=40，std::function 对象被传输=]
演示代码显示 两个示例，一个带有 自由函数，另一个带有 class方法。

完整的工作代码，用

编译

g++ -Wpedantic -std=c++11 testGsl.cpp -lgsl -lcblas

如下：

#include <gsl/gsl_integration.h>
#include <cassert>
#include <cstdio>
#include <functional>

namespace details
{
  extern "C" {
  double gsl_function_cpp_wrapper_helper(double z, void *cpp_f)
  {
    const auto p_cpp_f =
        reinterpret_cast<std::function<double(double)> *>(cpp_f);
    return (*p_cpp_f)(z);
  }
  }
}

class gsl_function_cpp_wrapper
{
  using Cpp_F = std::function<double(double)>;

 public:
  operator const gsl_function_struct *() const { return &gsl_f_; };

  gsl_function_cpp_wrapper(Cpp_F &&cpp_f) : cpp_f_(std::move(cpp_f))
  {
    gsl_f_.function = details::gsl_function_cpp_wrapper_helper;
    gsl_f_.params = &cpp_f_;
  }

  gsl_function_cpp_wrapper(double(f)(double))
      : gsl_function_cpp_wrapper(Cpp_F(f))
  {
  }

  template <typename OBJ>
  gsl_function_cpp_wrapper(double (OBJ::*method)(double) const, const OBJ *obj)
      : gsl_function_cpp_wrapper(
            Cpp_F(std::bind(method, obj, std::placeholders::_1)))
  {
  }
  template <typename OBJ>
  gsl_function_cpp_wrapper(double (OBJ::*method)(double), OBJ *obj)
      : gsl_function_cpp_wrapper(
            Cpp_F(std::bind(method, obj, std::placeholders::_1)))
  {
  }

 protected:
  Cpp_F cpp_f_;
  gsl_function_struct gsl_f_;
};

//----------------

class Class_Example
{
 public:
  double f(double x) const { return std::log(alpha_ * x) / std::sqrt(x); }

 protected:
  double alpha_ = 1;
};

double free_function_example(double x)
{
  const double alpha = 1;
  return std::log(alpha * x) / std::sqrt(x);
}

//----------------

int main()
{
  double result, error;
  gsl_integration_workspace *const w = gsl_integration_workspace_alloc(1000);
  assert(w != nullptr);

  //----------------

  Class_Example class_example;

  gsl_function_cpp_wrapper wrapper_1(&Class_Example::f, &class_example);

  gsl_integration_qags(wrapper_1, 0, 1, 0, 1e-7, 1000, w, &result, &error);

  std::printf("result          = % .18f\n", result);
  std::printf("estimated error = % .18f\n", error);

  //----------------

  gsl_function_cpp_wrapper wrapper_2(free_function_example);

  gsl_integration_qags(wrapper_2, 0, 1, 0, 1e-7, 1000, w, &result, &error);

  std::printf("result          = % .18f\n", result);
  std::printf("estimated error = % .18f\n", error);

  //----------------

  gsl_integration_workspace_free(w);

  return EXIT_SUCCESS;
}

函数c ++中的gsl集成

gsl integration within a function c++

c++

gsl