I have been writing a Fortran and C interface, but there always appears errors. I'm not sure why there is the problem. The original C code is following:
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>
void gsl_vector_setvalue(gsl_vector* v, size_t index, double value)
{
gsl_vector_set(v, index-1, value);
}
double gsl_vector_getvalue(gsl_vector* v, size_t index)
{
return gsl_vector_get(v, index-1);
}
int rosenbrock_f (const gsl_vector * x, void * params, gsl_vector * f)
{
double x0 = gsl_vector_get (x, 0);
double x1 = gsl_vector_get (x, 1);
double y0 = 1. * (1 - x0);
double y1 = 10. * (x1 - x0 * x0);
gsl_vector_set (f, 0, y0);
gsl_vector_set (f, 1, y1);
return GSL_SUCCESS;
}
int print_state (size_t iter, gsl_multiroot_fsolver * s)
{
printf ("iter = %3zu x = % .3f % .3f " "f(x) = % .3e % .3e\n", iter, gsl_vector_get (s->x, 0), gsl_vector_get (s->x, 1), gsl_vector_get (s->f, 0), gsl_vector_get (s->f, 1));
}
void gsl_multiroots(int rosenbrock())
{
const gsl_multiroot_fsolver_type *T;
gsl_multiroot_fsolver *s;
int status;
size_t i, iter = 0;
const size_t n = 2;
//struct rparams p = {1.0, 10.0};
struct rparams *p=NULL;
gsl_multiroot_function f = {rosenbrock, n, p};
double x_init[2] = {-10.0, -5.0};
gsl_vector *x = gsl_vector_alloc (n);
gsl_vector_set (x, 0, x_init[0]);
gsl_vector_set (x, 1, x_init[1]);
T = gsl_multiroot_fsolver_hybrids;
s = gsl_multiroot_fsolver_alloc (T, 2);
gsl_multiroot_fsolver_set (s, &f, x);
print_state (iter, s);
do
{
iter++;
status = gsl_multiroot_fsolver_iterate (s);
print_state (iter, s);
if (status) break;
/* check if solver is stuck */
status = gsl_multiroot_test_residual (s->f, 1e-7);
}
while (status == GSL_CONTINUE && iter < 1000);
printf ("status = %s\n", gsl_strerror (status));
gsl_multiroot_fsolver_free (s);
gsl_vector_free (x);
}
int main()
{
gsl_multiroots(rosenbrock_f);
return 0;
}
And when Ｉ compile and link it, it can run and have the correct result. But when I write rosenbrock_f in a Fortran code and use the intrinsic mode iso_c_binding to make them interoperable, it can't have the result, only say
> Program received signal SIGSEGV: Segmentation fault - invalid memory
> reference.
>
> Backtrace for this error:
> #0 0x7FCE7F618777
> #1 0x7FCE7F618D7E
> #2 0x7FCE7F270D3F
> #3 0x7FCE7FCF68AD
> #4 0x400BAE in gsl_vector_getvalue
> #5 0x400E45 in rosenbrock.1876 at multiroot.f90:?
> #6 0x7FCE7FC4A923
> #7 0x400D6F in gsl_multiroots
> #8 0x400F0D in MAIN__ at multiroot.f90:? Segmentation fault (core dumped)
here is the C and Fortran code:
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>
void gsl_vector_setvalue(gsl_vector* v, size_t index, double value)
{
gsl_vector_set(v, index-1, value);
}
double gsl_vector_getvalue(gsl_vector* v, size_t index)
{
return gsl_vector_get(v, index-1);
}
int print_state (size_t iter, gsl_multiroot_fsolver * s)
{
printf ("iter = %3zu x = % .3f % .3f " "f(x) = % .3e % .3e\n", iter, gsl_vector_get (s->x, 0), gsl_vector_get (s->x, 1), gsl_vector_get (s->f, 0), gsl_vector_get (s->f, 1));
}
void gsl_multiroots(int rosenbrock())
{
const gsl_multiroot_fsolver_type *T;
gsl_multiroot_fsolver *s;
int status;
size_t i, iter = 0;
const size_t n = 2;
struct rparams *p=NULL;
gsl_multiroot_function f = {rosenbrock, n, p};
double x_init[2] = {-10.0, -5.0};
gsl_vector *x = gsl_vector_alloc (n);
gsl_vector_set (x, 0, x_init[0]);
gsl_vector_set (x, 1, x_init[1]);
T = gsl_multiroot_fsolver_hybrids;
s = gsl_multiroot_fsolver_alloc (T, 2);
gsl_multiroot_fsolver_set (s, &f, x);
print_state (iter, s);
do
{
iter++;
status = gsl_multiroot_fsolver_iterate (s);
print_state (iter, s);
if (status) break;
/* check if solver is stuck */
status = gsl_multiroot_test_residual (s->f, 1e-7);
}
while (status == GSL_CONTINUE && iter < 1000);
printf ("status = %s\n", gsl_strerror (status));
gsl_multiroot_fsolver_free (s);
gsl_vector_free (x);
}
fortran module:
MODULE gsl_interfaces
USE, INTRINSIC :: ISO_C_BINDING
IMPLICIT NONE
TYPE, BIND(C) :: gsl_block
INTEGER(C_SIZE_T) :: size
TYPE(C_PTR) :: data
END TYPE gsl_block
TYPE, BIND(C) :: gsl_vector
INTEGER(C_SIZE_T) :: size
INTEGER(C_SIZE_T) :: stride
REAL(C_DOUBLE) :: data
TYPE(C_PTR) :: block
INTEGER(C_INT) :: owner
END TYPE gsl_vector
TYPE, BIND(C) :: gsl_matrix
INTEGER(C_SIZE_T) :: size1
INTEGER(C_SIZE_T) :: size2
INTEGER(C_SIZE_T) :: tda
REAL(C_DOUBLE) :: data
TYPE(C_PTR) :: block
INTEGER(C_INT) :: owner
END TYPE gsl_matrix
TYPE, BIND(C) :: rparams
REAL(C_DOUBLE) :: A
real(c_double) :: b
END TYPE
INTERFACE
SUBROUTINE GSL_VECTOR_SETVALUE(V, ROW_INDEX, V_VALUE) BIND(C, NAME='gsl_vector_setvalue')
USE ISO_C_BINDING
import gsl_vector
TYPE(c_ptr), VALUE :: V
INTEGER :: ROW_INDEX
REAL(C_DOUBLE), VALUE :: V_VALUE
END SUBROUTINE GSL_VECTOR_SETVALUE
REAL(C_DOUBLE) FUNCTION GSL_VECTOR_GETVALUE(V, ROW_INDEX) BIND(C, NAME='gsl_vector_getvalue')
USE ISO_C_BINDING
import gsl_vector
TYPE(c_ptr), VALUE :: V
INTEGER :: ROW_INDEX
END FUNCTION GSL_VECTOR_GETVALUE
SUBROUTINE GSL_MULTIROOTS(ROSENBROCK) BIND(C, NAME='gsl_multiroots')
USE ISO_C_BINDING
type(c_funptr), value :: ROSENBROCK
! TYPE(c_ptr), value :: x
END SUBROUTINE
! INTEGER(C_INT) FUNCTION ROSENBROCK(x, params, f) BIND(C, NAME='rosenbrock_f')
! USE ISO_C_BINDING
! import gsl_vector
! IMPORT rparams
! TYPE(gsl_vector), target :: x, f
! TYPE(C_PTR), target :: params
! END FUNCTION ROSENBROCK
END INTERFACE
END MODULE gsl_interfaces
fortran program:
PROGRAM MULTIROOT
USE, INTRINSIC :: ISO_C_BINDING
USE GSL_INTERFACES
IMPLICIT NONE
integer(c_int), parameter :: GSL_SUCCESS=0
type(gsl_vector), target :: x, f
CALL GSL_MULTIROOTS(c_funloc(ROSENBROCK))
contains
integer(c_int) function rosenbrock(x, params, f) bind(c)
use iso_c_binding
use gsl_interfaces
type(c_ptr) :: x, f
type(c_ptr) :: params
real(c_double) :: x0, x1, y0, y1
!a=params%a
!b=params%b
x0 = gsl_vector_getvalue (x, 1)
x1 = gsl_vector_getvalue (x, 2)
y0 = 1. * (1 - x0)
y1 = 10. * (x1 - x0 * x0)
call gsl_vector_setvalue (f, 1, y0)
call gsl_vector_setvalue (f, 2, y1)
rosenbrock=GSL_SUCCESS
end function rosenbrock
END PROGRAM
Can someone help me out? I don't know why there is memory invalid reference.
