gsl integration within a function c++

Question

I am just trying integrate over a function in C++. I have been trying to use gsl as I have seen this recommended online. I followed the gsl example with little success.

This is my C++ code:

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;
}

For clarification the same code in Python (which works) looks like this:

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]

Where quad performs the integration (it's a scipy module).

I get two error messages:

ISO C++ forbids taking the address of an unqualified or parenthesized non-static member function to form a pointer to member function. Say ‘&cosmo::inverseE’ [-fpermissive]

cannot convert ‘double (cosmo::)(double)’ to ‘double ()(double, void)’ in assignment* cosmo is the name of the class which contains both of these functions.

It seems to be that this should not be a difficult thing to do. Advice as to where I am going wrong would be much appreciated!

EDIT: 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

I see two problems:

1/ GSL expects your inverseE function to have the following prototype

double inverseE(double z,void *use_data);

in your code you have declared:

double inverseE(double z);

2/ Like your code is in C++ and GSL is a C library, you have to make your C++ function callable from C.

The solution is to declare your inverseE function as follow:

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;
   }
}

This extern "C" makes your C++ function binary compatible with C call convention.

With these two modifications I think your code should be ok.

---

UPDATE: 2/

In my answer I considered in 2/ that inverseE was a function. Here I consider the case where it is a method of your class.

This is an example where void *user_data comes to the rescue:

Declare this wrapper:

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

Then YourClass is defined as follow:

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

  ...
};

As you mentioned in your comment some care must be taken concerning forward declaration. To be clear, please find below a complete runnning example that reproduce the result of the GSL official doc but using a C++ class with a method f()

Complete running code:

Can be compiled with:

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

Code:

#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;
}

On my computer I get:

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

Answer 2

Reading some comments I finally post another version which involves less boilerplate code.

It works as follow:

• there is an extern "C" gsl_function_cpp_wrapper_helper function to insure C ABI (application binary interface). This is the function the gsl's subroutines will use.

• there is a C++ class gsl_function_cpp_wrapper that wraps the C++ function/class method thanks to a std::function. This wrapper is convertible to a gsl_function_struct pointer that uses gsl_function_cpp_wrapper_helper as C function, the std::function object is transmitted thank to the gsl's void *params parameter

• the demonstration code shows two examples, one with a free function, the other with a class method.

The complete working code, to be compiled with

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

is as follows:

#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;
}