1

我正在尝试编写一个包装器以将 gsl 库与 Fortran 一起使用。我设法让一个简单的包装器工作 - 来自http://www.helsinki.fi/~fyl_tlpk/luento/ohj-13-GSL-e.html的示例

Fortran 代码

program gsltest
    implicit none

    real(kind=selected_real_kind(12)) :: a = 0.11, res
    external :: glsgateway

    call gslgateway(a,res)
    write(*,*) 'x', a, 'atanh(x)', res

end program gsltest

函数

#include <gsl/gsl_math.h>

void gslgateway_(double *x, double *res){
   *res = gsl_atanh(*x);
}

这一切都很好。但是，我在使用更复杂的包装器时遇到了问题。我从http://apwillis.staff.shef.ac.uk/aco/freesoftware.html的示例修改了以下代码

c 包装器 (rng_initialise.c)

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

static gsl_rng* r;

void rng_initialise__(int* s) {
   r = gsl_rng_alloc(gsl_rng_taus); 
   gsl_rng_set(r, (unsigned long int)(*s)); 
}

Fortran 主程序 (main.f90)

PROGRAM main
    integer seed

    call system_clock(seed)
    WRITE (*,*) 'calling rng_initialise'
    call rng_initialise(seed)

END PROGRAM main

然后我编译并链接

gcc -c rng_initialise.c
g95 -c main.f90
g95 -o main main.o rng_initialise.o -L/usr/libs -lgsl

当我运行这个程序时，我没有得到任何输出。但是，如果我注释掉 rng_initialise 中的行

...
void rng_initialise__(int* s) {
   // r = gsl_rng_alloc(gsl_rng_taus);  
   // gsl_rng_set(r, (unsigned long int)(*s)); 
}

然后我从 Fortran 代码中获得输出（它将“calling_rng_initialise”写入 STDOUT）。

所以，问题似乎是对 gsl_rng_alloc 和 gsl_rng_set 的调用。但是我没有收到任何错误消息，我不知道他们为什么会阻止 Fortran 代码做任何事情。有任何想法吗？

4

3 回答 3

4

如前所述，最好的方法是使用 Fortran ISO C 绑定，因为它将指示 Fortran 使用 C 调用约定来匹配 GSL 库的 C 例程。ISO C 绑定是 Fortran 2003 标准的一部分，并且多年来一直在许多 Fortran 95 编译器中可用。作为标准的一部分，它使 Fortran 和 C 的接口在两个方向上都具有可移植性，并且比过去需要的操作系统和编译器相关的 hack 更容易。我建议忽略早于 ISO C 绑定的说明已过时。

通常，您不需要编写任何 C 包装代码来调用 GSL 库，只需使用 Fortran 规范语句来描述 Fortran 语法中的 C 例程接口即可。这是一个调用 GSL 例程 gsl_cdf_chisq_Q 的简单示例

program cum_chisq_prob

   use iso_c_binding

   interface GSL_CummulativeChiSq_Prob_Upper   

      function gsl_cdf_chisq_Q  ( x, nu )  bind ( C, name="gsl_cdf_chisq_Q" )

         import

         real (kind=c_double) :: gsl_cdf_chisq_Q

         real (kind=c_double), VALUE, intent (in) :: x
         real (kind=c_double), VALUE, intent (in) :: nu

      end function gsl_cdf_chisq_Q

   end interface GSL_CummulativeChiSq_Prob_Upper

   real (kind=c_double) :: chisq_value
   real (kind=c_double) :: DoF
   real (kind=c_double) :: Upper_Tail_Prob    

   write ( *, '( / "Calculates cumulative upper-tail probability for Chi-Square distribution." )' )
   write ( *, '( "Input Chisq Value, Degrees of Freedom: " )', advance='no' )
   read ( *, * ) chisq_value, DoF

   Upper_Tail_Prob = gsl_cdf_chisq_Q  ( chisq_value, DoF )

   write ( *, '( "Probability is:", 1PG17.10 )' )  Upper_Tail_Prob

   stop

end program cum_chisq_prob

更简单：您可以在http://www.lrz.de/services/software/mathematik/gsl/fortran/找到一个预先编写的库，允许您从 Fortran 调用 GSL

于 2011-09-08T17:11:29.547 回答
1

很可能您在某种程度上将这两个例程之间的联系弄错了。如果在您通过该界面时未正确处理堆栈，那么任何事情都可能发生。

我没有注意到 Fortran 端或 C 端的任何代码指定了对方的调用约定。我不是 Gnu Fortran 的专家，但我知道大多数编译器都需要某种注释，说明它们应该使用另一个编译器的调用约定，否则可能会发生坏事。

通过网络搜索，我看到G95 Fortran 手册 (PDF)有一个很长的部分，标题为“与 G95 程序的接口”，似乎对此进行了详细介绍。仅从略读来看，您似乎应该在BIND(C)该 C 例程的 Fortran 函数声明中使用该属性。

于 2011-09-08T12:45:04.527 回答
0

问题出在文件中static gsl_rng* r;定义的问题上。C但我不知道/理解为什么标准不允许这样做。fgsl在稍微研究了包的源文件后，我发现了一个可行的调整。fortran文件_random_f.f90

module fgsl
    use, intrinsic :: iso_c_binding
    implicit none

    type, bind(C) :: fgsl_rng_type
        type(c_ptr) :: gsl_rng_type_ptr = c_null_ptr
    end type fgsl_rng_type

    type, bind(C) :: fgsl_rng
        type(c_ptr) :: gsl_rng_ptr = c_null_ptr
    end type fgsl_rng
end module fgsl

PROGRAM call_fgsl_rndm
    use, intrinsic :: iso_c_binding
    use fgsl
    implicit none

    interface
        subroutine srndm(seed, t, r) bind(C)
            import
            integer(C_INT) :: seed
            type(fgsl_rng) :: r
            type(fgsl_rng_type) :: t
        end subroutine srndm

        function rndm(r) bind(C)
            import
            real(C_DOUBLE) :: rndm
            type(fgsl_rng) :: r
        end function rndm
    end interface

    type(fgsl_rng) :: r
    type(fgsl_rng_type) :: t

    integer(C_INT) :: seed
    real(C_DOUBLE) :: xi
    seed = 1
    call srndm(seed, t, r)
    xi = rndm(r)
    print *, xi
    xi = rndm(r)
    print *, xi
    xi = rndm(r)
    print *, xi
END PROGRAM

和C文件random_c.c

#include <gsl/gsl_rng.h>

typedef struct{
    gsl_rng *gsl_rng_ptr;
} fgsl_rng;

typedef struct{
    gsl_rng_type *gsl_rng_type_ptr;
} fgsl_rng_type;

void srndm(int *seed, fgsl_rng_type* t, fgsl_rng* r) {
    t->gsl_rng_type_ptr = (gsl_rng_type *) gsl_rng_mt19937;  // cast to remove the const qualifier
    r->gsl_rng_ptr = gsl_rng_alloc(gsl_rng_mt19937);
    gsl_rng_set(r->gsl_rng_ptr, *seed);
}

double rndm(fgsl_rng* r) {
    return gsl_rng_uniform(r->gsl_rng_ptr);
}

虽然只用到了结构体中的指针，但引入fgsl_rngandfgsl_rng_type是必要的。否则，程序将失败。不幸的是，我仍然不清楚为什么它必须以这种方式工作。

于 2018-10-18T09:19:16.410 回答