I'm trying to call the GSL rountine CQUAD from Fortran. My idea was to write a .c subroutine that calls the gsl rountine and depends on a function and bounds. Two problems: I have only very little idea about c and fortrans iso_c_binding. My attempt is as follows:
A simple calling program (similar to M.S.B's post in No output from a Fortran program using the Gnu Scientific Library via a c wrapper ):
program test
use iso_c_binding
interface
function my_cquad (f,a,b) bind(c)
import
real (kind=c_double) :: my_cquad
interface
function f(x) bind(c)
import
real(kind=c_double) :: f,x
end function
end interface
real (kind=c_double) :: a,b
end function my_cquad
end interface
real (kind=c_double) :: y,a,b
a=0. ; b=1.
y=my_cquad(g,a,b)
print*,y
stop
contains
function g(x) bind(C)
real(kind=c_double) :: g,x
g=sin(x)/x
return
end function g
end program test
The .c subroutine (basically taken from the example given by the author of CQUAD in https://scicomp.stackexchange.com/questions/4730/numerical-integration-handling-nans-c-fortran):
#include <stdio.h>
#include <gsl/gsl_integration.h>
double my_cquad ( double my_f() , double a , double b )
{
gsl_function f;
gsl_integration_cquad_workspace *ws = NULL;
double res, abserr;
size_t neval;
/* Prepare the function. */
f.function = my_f;
f.params = NULL;
/* Initialize the workspace. */
if ( ( ws = gsl_integration_cquad_workspace_alloc( 200 ) ) == NULL ) {
printf( "call to gsl_integration_cquad_workspace_alloc failed.\n" );
abort();
}
/* Call the integrator. */
if ( gsl_integration_cquad( &f, a , b , 1.0e-10 , 1.0e-10 , ws , &res , &abserr , &neval ) != 0 ) {
printf( "call to gsl_integration_cquad failed.\n" );
abort();
}
/* Free the workspace. */
gsl_integration_cquad_workspace_free( ws );
/* Bye. */
return res;
}
The .c subroutine alone seems to work fine. This can be tested with:
double g (double x)
{
return sin(x)/x;
}
int main () {
double y;
y=my_cquad(g,0.,1.);
printf("y: %2.18f\n", y);
return 0;
}
But together with the .f90 calling program, at the moment it compiles but at runtime I get a segmentation fault that I don't quite get.
Additionally, it would of course be good to have some kind of wrapper that creates a c-type function depending on a fortran type function. I'm thinking about something like:
function f_to_c(f) bind(c)
real(kind=c_double) :: f_to_c
real(kind=8) :: f
f_to_c=real(f,kind=c_double)
end function
But this desn't cover the dummy variables.
Thanks in advance and very sorry for the amount of code.
