gsl 集成在函数 c++ 中

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 是包含这两个函数的类的名称。

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

编辑：类

#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一个函数。在这里我考虑它是你的类的方法的情况。

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

  ...
};

正如您在评论中提到的，必须注意前瞻性声明。为了清楚起见，请在下面找到一个完整的运行示例，该示例重现了GSL 官方文档的结果，但使用带有方法的 C++ 类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