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