Fortran 和 c 互操作性
Posted
技术标签:
【中文标题】Fortran 和 c 互操作性【英文标题】：Fortran and c interoperability 【发布时间】：2015-09-03 11:40:05 【问题描述】：我一直在写一个Fortran和C接口，但是总是出现错误。我不确定为什么会出现问题。原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
有人可以帮帮我吗？不知道为什么会出现内存无效引用。
【问题讨论】：
虽然不是很确定，但是写“void gsl_multiroots(int rosenbrock())”而不是“void gsl_multiroots(int (* rosenbrock)())”是不是没问题...？跨度> 我试过你的答案，但它不起作用。谢谢你回答我的问题。我想整数（c_int）函数rosenbrock中的参数有一些错误。 再试一次：如何将 VALUE 属性附加到 rosenbrock() 中的所有虚拟参数（x、f 和 params）......？ 接口块中的ROW_INDEX可能还需要VALUE...嗯，很复杂。 GSL 有一个公共的 Fortran 接口：lrz.de/services/software/mathematik/gsl/fortran。查看他们的实施可能会有所帮助。 【参考方案1】：注意：M. S. B. 指出，GNU 科学图书馆已经存在 Fortran wrapper。
以下是 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);
在接口中，我为整数添加了 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 作为句柄，那应该不会那么痛苦。
【讨论】：
这个结果正是我想要的。但我很困惑为什么 x、params 和 f 应该具有价值属性。它们是 c 代码中的指针。和 x & f 是 gsl_vector 类型，它们不是值。 @zmwang 在你的函数声明中它们是 type(c_ptr)，因此它们需要具有 value 属性，因为 c 接口也“只有”一个 c 指针而不是指向指针的指针.它与 GSL_VECTOR_SET|GETVALUE 中的 v 的逻辑相同，您正确地进行了此声明。以上是关于Fortran 和 c 互操作性的主要内容，如果未能解决你的问题，请参考以下文章
在默认类型和 C 互操作类型之间转换 Fortran 字符和逻辑数组