我一直在写Fortran 和C 接口(interface)，但总是出现错误。我不确定为什么会出现问题。原始C代码如下:
#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;
}
并且在Ｉ编译链接时，可以运行并得到正确的结果。 但是当我用 Fortran 代码写 rosenbrock_f 并使用内在模式 iso_c_binding 让它们互操作时，它不能有结果，只能说
> 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)
这是 C 和 Fortran 代码:
#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 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 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
有人能帮帮我吗？不知道为什么会出现内存无效引用。
最佳答案
注意:M. S. B.指出存在一个Fortran wrapper已经用于 GNU 科学图书馆。
以下是适用于我的 gcc 5.2 的内容。
我修改了您的 C 代码以与 gsl documentation 中给出的示例保持一致，因为我不确定，您的功能描述的顺序是否与界面匹配。此外，我更改了函数指针声明以匹配 gsl 的期望。
#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)(gsl_vector *x, void *p, gsl_vector *f))
{
const gsl_multiroot_fsolver_type *T;
gsl_multiroot_fsolver *s;
int status;
size_t i, iter = 0;
const size_t n = 2;
void *p = NULL;
gsl_multiroot_function f;
f.f = rosenbrock;
f.n = n;
f.params = &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);
}
在接口(interface)中，我为整数添加了 value 属性，这些整数应该作为按值参数调用:
MODULE gsl_interfaces
USE, INTRINSIC :: ISO_C_BINDING
IMPLICIT NONE
INTERFACE
SUBROUTINE GSL_VECTOR_SETVALUE(V, ROW_INDEX, V_VALUE) BIND(C, NAME='gsl_vector_setvalue')
USE ISO_C_BINDING
TYPE(c_ptr), VALUE :: V
INTEGER,value :: 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
TYPE(c_ptr), VALUE :: V
INTEGER,value :: 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
END SUBROUTINE
END INTERFACE
END MODULE gsl_interfaces
在 rosenbrock 函数中，我为 c_ptr 参数添加了缺失值属性:
PROGRAM MULTIROOT
USE, INTRINSIC :: ISO_C_BINDING
USE GSL_INTERFACES
IMPLICIT NONE
integer(c_int), parameter :: GSL_SUCCESS=0
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),value :: x, f
type(c_ptr),value :: 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
运行它时，我得到:
iter = 0 x = -10.000 -5.000 f(x) = 1.100e+01 -1.050e+03
iter = 1 x = -10.000 -5.000 f(x) = 1.100e+01 -1.050e+03
iter = 2 x = -3.976 24.827 f(x) = 4.976e+00 9.020e+01
iter = 3 x = -3.976 24.827 f(x) = 4.976e+00 9.020e+01
iter = 4 x = -3.976 24.827 f(x) = 4.976e+00 9.020e+01
iter = 5 x = -1.274 -5.680 f(x) = 2.274e+00 -7.302e+01
iter = 6 x = -1.274 -5.680 f(x) = 2.274e+00 -7.302e+01
iter = 7 x = 0.249 0.298 f(x) = 7.511e-01 2.359e+00
iter = 8 x = 0.249 0.298 f(x) = 7.511e-01 2.359e+00
iter = 9 x = 1.000 0.878 f(x) = -2.653e-10 -1.218e+00
iter = 10 x = 1.000 0.989 f(x) = -2.353e-11 -1.080e-01
iter = 11 x = 1.000 1.000 f(x) = 0.000e+00 0.000e+00
status = success
我希望这对您有所帮助，可以帮助您实现您真正想要实现的目标。如果可能的话，我会远离为 C 结构声明 Fortran 类型。如果您可以只使用 c_ptr 作为调用方的句柄，那应该会少很多痛苦。
关于c - Fortran 和 c 的互操作性，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/32374612/