Possible to do C++ interoperation with std::function?

My main program is in fortran. It calls a C++ library fn (that I wrote) and passes to it a fortran function as one of the args. This works well when the function is declared as double(*f1)((void)). I want to declare this function as std::function<double((void))>f1 (for perhaps questionable reasons) and found that this causes a segfault at the point where this function is called. Is there something that can be done to make interop work with std::function type functions?

Additional background: The reason I thought of using std::function is to be able to call this library fn from other C++ codes and have the option of giving it a C++ lambda. The callback function is a supposed to be a sampling function that can can have varying number/types parameters and “source distributions” (uniform, log-normal etc. etc.). I use a module and make these parameters available to the function through the module and, where the samples are needed this function is called with no arguments. This worked well with the callback declared as double(*f1)((void)) in the C++ side.

Now, I wanted to have the flexibility of calling the library from a C++ driver. To have the same argument-less calling behavior I defined the parameters in the body of the caller and used the C++ anonymous function feature to capture the relevant parameters and use this lambda in the library call . This did not work with the function pointer declaration. But, it works with the std::function declaration. Later on I realized I could use global variables in the C++ driver to get the desired outcome i.e. argument-less callback that works with both fortran and C++ drivers using the double(*f1)((void))-type declaration. This is what I ended up doing. Finally, I note that there are additional nice features one could use when dealing with an std::function-type function.

Minimal reproducers along with build command are given below:

gfortran -ggdb -O0  -fcheck=all -fbacktrace -c main.f90
g++ -ggdb -O0  -c test.cxx
gfortran -ggdb -O0  -fcheck=all -fbacktrace -o test main.o test.o -lstdc++
main.f90
module vars
   use iso_c_binding
   implicit none
   integer, parameter::WP=c_double

   type, bind(C):: part_info
      real(wp) :: d_mean
      real(wp) :: d_std
      real(wp) :: d_min
      real(wp) :: d_max
   endtype part_info
 
   type(part_info) :: pinfo

   !interface to c++ function
   interface
   subroutine test(npart, part, f1,f2) bind(C)
       use iso_c_binding
       integer(C_INT), intent(out) :: npart
       type(c_ptr), intent(out) :: part
       type(C_FUNPTR), value:: f1,f2
   endsubroutine test
   endinterface
  
   contains
   function dsamp() bind(C)
      implicit none
       real(wp) :: dsamp
      dsamp = 0.0_wp
      !get a random number between pinfo%d_min and pinfo%d_max
      call random_number( dsamp )
      dsamp = dsamp * ( pinfo%d_max - pinfo%d_min ) + pinfo%d_min
   endfunction dsamp

   function dist(x) bind(C)
      implicit none
      real(wp), intent(in) :: x(3)
      real(wp) :: dist
      dist = 1.0 - sum( x**2 )
   endfunction dist
 
endmodule vars

program main
  use vars
  implicit none

  procedure(dsamp), pointer :: dsampler_fptr=>null()
  procedure(dist), pointer :: dist_fptr=>null()
  type(C_FUNPTR) :: dsampler_cptr, dist_cptr
  type(c_ptr):: part_cptr
  integer(C_INT) :: npart
  real(wp), pointer :: part(:,:)=>null()
  integer i

  !initialize some parameters
  pinfo%d_mean = 0.0_wp
  pinfo%d_std  = 0.0_wp
  pinfo%d_min  = 0.01_wp
  pinfo%d_max  = 0.05_wp

  dsampler_fptr => dsamp
  dist_fptr => dist

  dsampler_cptr = c_funloc( dsampler_fptr )
  dist_cptr = c_funloc( dist_fptr )

  call test( npart, part_cptr, dsampler_cptr, dist_cptr )
  call c_f_pointer( part_cptr, part, [4,npart] )
  !print first 10
  do i=1,min(10,npart)
     write(*,"(4(g16.8,x))") part(:,i)
  enddo
endprogram main

test.cxx
#include <functional>
// #include <iostream>
extern "C"
{
// void test(int &npart, double *&part, std::function<double((void))>f1,std::function<double(double*)>f2)
void test(int &npart, double *&part, double(*f1)(void), double (*f2)(double*))
{
  npart = 12;
  part = new double[4 * npart];

  for (int i = 0; i < npart; i++)
  {
    part[4*i]   = 0.0;
    part[4*i+1] = 0.0;
    part[4*i+2] = 0.0;
    part[4*i+3] = f1();
    double *x = &part[4*i];
    // std::cout << "f2(x) = " << f2(x) << std::endl;
  }
}
}
main.cxx (standalone ver.)
#include <iostream>
#include <random>
#include <cmath>

using namespace std;
extern "C"
{
  typedef struct
  {
    double d_mean;
    double d_std;
    double d_min;
    double d_max;
  } part_info;
}

part_info pinfo;

double dia_sampler()
{
  std::random_device rd;
  std::mt19937 gen(rd());
  std::lognormal_distribution<> d(pinfo.d_mean, pinfo.d_std);
  // only return values between d_min and d_max
  double d_val = d(gen);
  while (d_val < pinfo.d_min || d_val > pinfo.d_max)
  {
    d_val = d(gen);
  }
  return d_val;
};

double dist_fn(double *x)
{
  double x_c[3] = {0.5, 0.5, 0.5}; // center of cylinder
  double r = 0.5;                  // radius of cylinder
  double h = 1.0;                  // height of cylinder
  double d = sqrt((x[0] - x_c[0]) * (x[0] - x_c[0]) + (x[1] - x_c[1]) * (x[1] - x_c[1])) - r;
  double d1 = x[2] - x_c[2] - h / 2.0;
  double d2 = x_c[2] - x[2] - h / 2.0;
  double d3 = max(d1, d2);
  return -max(d, d3);
};

#include <functional>
// #include <iostream>
extern "C"
{
  // void test(int &npart, double *&part, std::function<double((void))>f1,std::function<double(double*)>f2)
  void test(int &npart, double *&part, double (*f1)(void), double (*f2)(double *))
  {
    npart = 12;
    part = new double[4 * npart];

    for (int i = 0; i < npart; i++)
    {
      part[4 * i] = 0.0;
      part[4 * i + 1] = 0.0;
      part[4 * i + 2] = 0.0;
      part[4 * i + 3] = f1();
      double *x = &part[4 * i];
      // std::cout << "f2(x) = " << f2(x) << std::endl;
    }
  }
}

int main()
{
  int npart;
  double *part;
  pinfo.d_mean = 0.0;
  pinfo.d_std = 4.0;
  pinfo.d_min = 0.02;
  pinfo.d_max = 0.05;

  double (*f1)() = dia_sampler;
  double (*f2)(double *) = dist_fn;

  test(npart, part, f1, f2);

  for (int i = 0; i < min(10,npart); i++)
  {
    std::cout << "part[" << i << "] = " << part[4 * i] << ", " << part[4 * i + 1] << ", " << part[4 * i + 2] << ", " << part[4 * i + 3] << std::endl;
  }
}

No, not if you are thinking of standard Fortran with its interoperability facility which is actually with a C companion processor only.

You will note std::function in C++ is a special thing, technically it is a class template for a wrapper type with polymorphism for general-purpose handling of various forms of functions including, as you noted, lamdas. Note standard Fortran does not know of basic C++ itself, let alone any interoperability with such modern features in C++.

If you remain interested in standard Fortran, your best bet is to keep it simple using C at the Fortran interface. And if you see value, consider suitable wrapping with modern C++ of this C bridge toward consumption in C++.

1 Like

I have a few comments related to your C++ code:

  • The array allocation in test feels like the wrong approach; forgetting to free part will lead to a memory leak. Because of such perils, there is a lot of momentum at the moment to move away from C and C++ (see Is it time to retire C and C++ for Rust in new programs? • The Register) to something else (Rust/Go/Swift/Carbon/…).
  • If npart is a constant equal to 12, it can be made part of the “contract” instead of a hidden parameter.
  • When the calling code is C++ it may make sense to just return a std::vector which is memory-safe (it has ownership of the data). But if you’d like to call the routine from C and Fortran, it could be easier to just decouple memory ownership from the behavior/algorithm.
  • The way your dia_sampler is written currently, the PRNG devices are initialized repeatedly at each call. You can add the static specifer to initialize them only the first time. One challenge of this approach is modifying/setting the seed value. A function-object would be better IMO, and would work with your proposed std::function-based interface.
  • Drop the d_ prefix in part_info (I’m guessing this stands for double). Such prefixing is also known as Hungarian notation in the sarcastic essay How to Write Unmaintainable Code.

Here are some excerpts of what the modifications could look like (C++20 required):

C++ Code
struct part_info {
  double mean;
  double std;
  double min;
  double max;
};

template<typename T = std::mt19937>
class dia_sampler {
public:

  // We should be able to seed the sampler explicitly for reproducibility    
  explicit dia_sampler(part_info info, int seed) : 
    prng_(seed),
    d_(info.mean,info.std),
    min_(info.min), 
    max_(info.max) {}

  // For convenience we can provide a default seed
  explicit dia_sampler(part_info pinfo) :
    dia_sampler(pinfo,std::random_device{}()) {}

  double operator()() {
    double r = d_(prng_);
    while (r < min_ || r > max_) {
        r = d_(prng_);
    }
    return r;
  }

private:
  T prng_;
  std::lognormal_distribution<double> d_;
  double min_, max_;
};

Some minor nit-picking in the distance function (use of const and std::hypot):

double dist_fn(const double x[3])
{
  const double x_c[3] = {0.5, 0.5, 0.5}; // center of cylinder
  const double r = 0.5;                  // radius of cylinder
  const double h = 1.0;                  // height of cylinder

  const double d = std::hypot(x[0] - x_c[0], x[1] - x_c[1]) - r;
  const double d1 = x[2] - x_c[2] - h / 2.0;
  const double d2 = x_c[2] - x[2] - h / 2.0;
  const double d3 = std::max(d1, d2);

  return -std::max(d, d3);
};

The core function test can be rewritten to return a std::vector:

using samplefunc = std::function<double(void)>;
using distfunc   = std::function<double(const double x[3])>;

std::vector<double> test(samplefunc sample, distfunc dist)
{
  const int npart = 12;
  std::vector<double> part(4*npart);
  for (int i = 0; i < npart; i++) {
    part[4 * i] = 0.0;
    part[4 * i + 1] = 0.0;
    part[4 * i + 2] = 0.0;
    part[4 * i + 3] = sample();
  }
  
  // ... dist() ...

  return part;
}

Alternatively (my preference), it could become a function-object too:

class TestFunctor {
public:
  
  using SampleFunction   = std::function<double(void)>;
  using DistanceFunction = std::function<double(const double [3])>;

    // The argument to DistanceFunction could also be std::span<double,3> or
    // std::array<double,3>. For easier interop with C, we will keep it
    // as a fixed-size buffer (the length is part of the contract).

  explicit TestFunctor(SampleFunction smpl, DistanceFunction dist) :
    smpl_(std::move(smpl)), dist_(std::move(dist)) {}

  void operator()(std::span<double> part) const {
    assert(part.size() % 4 == 0);        // use -DNDEBUG to turn off
    const int npart = part.size()/4;
    //assert(npart == 12);                   // ?
    for (int i = 0; i < npart; i++) {
      part[4 * i] = 0.0;
      part[4 * i + 1] = 0.0;
      part[4 * i + 2] = 0.0;
      part[4 * i + 3] = smpl_();
      // ... dist_() ... ;
    }
  }

private:
  SampleFunction smpl_;
  DistanceFunction dist_;
};

Here is the main program, demonstrating both approaches:

int main(int argc, char const *argv[])
{
  part_info info = {
    .mean = 0.0, .std = 4.0,
    .min = 0.02, .max = 0.05
  };                          // NB: C++20 designated initializers (DI)

  // Test the sampling function-object works
  dia_sampler sampler(info);
  std::cout << "sample: " << sampler() << '\n'; 

  //
  // Approach A, result returned as std::vector (automatic clean-up)
  //
  auto part = test(sampler,dist_fn);
  int npart = part.size() / 4;       // number of parts must be inferred

  //
  // Approach B, behavior and memory separated
  //
  auto tester = TestFunctor(
    dia_sampler(info),dist_fn
  );

  // Helper function to print some items
  auto display = [](std::span<double> part) {
    for (int i = 0; i < 4; ++i) {
      std::cout << "part[" << i << "] = " << part[i] << '\n';
    }
    std::cout << "----" << '\n';
  };

  // By using a std::span we are free to pick the container
  // of our liking. We do need to adjust the size ourselves. 

  const int num_parts = 12;   // state the number of items
  const int N = 4*num_parts;  // array length is 4Ă— of that

  // Vector container
  std::vector<double> vpart(N);
  tester(vpart);
  display(vpart);

  // Array container
  std::array<double,N> apart;
  tester(apart);
  display(apart);

  // Static array
  double spart[N];
  tester(spart);
  display(spart);

  // Dynamically-allocated array
  double *cpart = new double[N];
  tester({cpart,N});  // We need to state the size manually in this case
  display({cpart,N});
  delete[] cpart;     // (!) don't forget the [] when deallocating an array

  return 0;

}

To call the routine from Fortran you will need to write an adapter routine with a C-conforming function prototype. Using the function-object, which doesn’t require handling of the std::vector, this could be something simple like:

#include "test.h"    // Public C++ API
#include "c_test.h"  // C17-compatible public header

extern "C" {

void c_test(int npart, double* part, 
            const double(*smpl)(void), 
            const double(*dist)(const double [3]))
{
    // Forward the C-compatible callbacks as lambdas
    auto tester = TestFunctor(
        [=](void) { return smpl(); },
        [=](const double x[3]) { return dist(x); }
    );

    // Passing them directly also appears to work in practice.
    // With the lambdas you can pass/capture extra data variables if needed.
    // There could be some subtle issues with extern "C" and the binding
    // expected for the callbacks, but I haven't investigated it.
    //
    // auto tester = TestFunctor(smpl,dist);

    const size_t N = 4 * npart;
    std::span s{part,N};

    tester(s);
}

}

The wrapper code in Fortran is straightforward. You could use internal procedures to adapt the callback interface to your liking (classic Adapter pattern).

module test_m

  use, intrinsic :: iso_c_binding
  implicit none
  private

  public :: test
  public :: smplfunc, distfunc
  public :: wp

  integer, parameter :: wp = c_double

  abstract interface
    function smplfunc() bind(c)
      import c_double
      real(c_double) :: smplfunc
    end function
    function distfunc(x) bind(c)
      import c_double
      real(c_double), intent(in) :: x(3)
      real(c_double) :: distfunc
    end function
  end interface

contains

  subroutine test(npart,part,smpl,dist)
    integer, intent(in) :: npart
    real(c_double), intent(out), target :: part(4,npart)
    procedure(smplfunc) :: smpl
    procedure(distfunc) :: dist

    interface
      subroutine c_test(npart,part,smpl,dist) bind(c,name="c_test")
        import c_int, c_double, c_funptr
        integer(c_int), intent(in), value :: npart
        real(c_double), intent(out) :: part(4*npart)
        type(c_funptr), value :: smpl
        type(c_funptr), value :: dist
      end subroutine
    end interface
    real(c_double), pointer :: ptr_part(:) => null()

    ptr_part(1:4*npart) => part

    call c_test(npart,ptr_part,c_funloc(smpl),c_funloc(dist))

  end subroutine

end module

Finally, in the main program you can just give the locations of your callbacks directly:

program main
  use test_m
  implicit none
  integer, parameter :: npart = 12
  real(wp), allocatable :: part(:,:)
  type :: part_info
    real(wp) :: mean, std, min, max
  end type
  type(part_info) :: pinfo
  
  pinfo = part_info( &
    mean = 0.0_wp, &
    std = 0.0_wp, &
    min = 0.01_wp, &
    max = 0.05_wp)
  
  allocate(part(4,npart))
  call test(npart,part,dsamp,dist)
  write(*,"(4(g16.8,2X))") part(:,1)

contains
  !> pinfo taken from host scope
  function dsamp() bind(C)
    real(wp) :: dsamp
    dsamp = 0.0_wp
    call random_number( dsamp )
    dsamp = dsamp * ( pinfo%max - pinfo%min ) + pinfo%min
  end function dsamp
  function dist(x) bind(C)
    implicit none
    real(wp), intent(in) :: x(3)
    real(wp) :: dist
    dist = 1.0_wp - sum( x**2 )
  end function
end program

If you want to test the modified version, I’ve placed all the needed files in cppexample.tar.gz.txt (4.4 KB). After downloaing the text file, all you need to do is:

$ mv cppexample.tar.gz.txt cppexample.tar.gz
$ tar -xzvf cppexample.tar.gz
$ cd cppexample/
$ make CXX=g++ FC=gfortran
$ make run
1 Like

Hi Ivan. Thank you for the detailed explanation and suggestions! I am half-wishing I had posted the full C++ code :slight_smile:. I appreciate the detailed review.

The C++ library code generates random sphere packings. So, npart (and the size of the part array) is not known beforehand. (the d_ prefix was referring to diameter. Not as bad as the Hungarian notation but could have been more descriptive)
I will take a closer look at your changes. I hope to be able to incorporate your modifications.
Thanks again!