How to integrate std::valarray with gsl?

Question

I am relatively new to C++, but I have some (scarce) coding and numerical experience.

I know that this question gets posted every now and then, how do you integrate an array. In MATLAB you can force your array to be a function (I forgot how, but I know I did it before) and send it to inbuilt integrators, so my question is how do you do it in C++.

I have this integral:

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

q is just double const, z is the integrating variable, but A(z) is an array (I'll call it actualfunction from now on) that has same number of points as z axis in my code. Integrating boundaries are z[0] and z[nz-1].

I calculated this integral by using trapezium rule, and for z-axis of 5000 points this takes 0.06 sec. My problem is that this calculation occurs roughly 300 * 30 * 20 times (I have 3 for loops), and this 0.06 sec grows very quickly to 3 hours of the simulation. And the entire bottleneck of my code is this integration (I can obviously speed up by reducing z, but that's not the point.)

I know that library functions are usually much better then user-written ones. I also know that I can't use something simpler as Simpson's rule, because the integrand is highly oscillatory, and I want to avoid my own implementation of some complicated numerical alghoritm.

GSL needs a function in form:

F = f(double x, void *params)

and I can probably use QAWO adaptive integration from gsl, but how do I make my function in form that turns my array into function?

I am thinking something as:

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

Is something like this possible? I doubt that numerical algorithm will use same spatial difference as actualfunction had, should I then somehow do interpolation of actualfunction or something?

Is there something better then gsl?

Yakk - Adam Nevraumont · Accepted Answer · 2021-08-30T23:42:22.877

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

use:

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

then

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

and you can pass pfun and pvoid through.

Apologies for any typos.

The idea is we write a lambda that does what we want. Then lambda_to_pfunc wraps it up so we can pass it as a void* and function pointer to C style APIs.

You'll have to properly manage the lifetime of everything of course.

How to integrate std::valarray with gsl?

1 Answers1