从 Fortran 调用 GSL 例程 CQUAD

Question

我正在尝试从 Fortran 调用 GSL 常规 CQUAD。我的想法是编写一个 .c 子例程，它调用 gsl 例程并依赖于函数和边界。两个问题：我对 c 和 fortrans iso_c_binding 知之甚少。我的尝试如下：

一个简单的调用程序（类似于MSB在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

.c子程序（基本上取自CQUAD作者在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;
}

单独的 .c 子例程似乎工作正常。这可以通过以下方式进行测试：

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;
}

但是与 .f90 调用程序一起，在它编译的那一刻，但在运行时我得到了一个我不太明白的分段错误。

此外，拥有某种根据 fortran 类型函数创建 c 类型函数的包装器当然会很好。我正在考虑类似的事情：

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

但这并不涵盖虚拟变量。

提前感谢，非常抱歉代码量。

Answer 1

注意，根据 Fortran 标准，内部函数不应具有 bind(C) 属性。我将函数移到了一个模块中。

a 和 b 必须按值传递给 my_cquad 和 x 到集成函数：

module functions_to_integrate

   use iso_c_binding

contains
   function g(x) bind(C)
      real(kind=c_double) :: g
      real(kind=c_double), value :: x
      g = sin(x)/x
   end function g
end module

program test

   use iso_c_binding

   use functions_to_integrate

   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
               real(kind=c_double), value :: x 
            end function
         end interface

         real (kind=c_double), value :: 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

end program test

测试：

> gfortran my_cquad.c test_cquad.f90 -lgsl -lopenblas
> ./a.out 
  0.94608307036718275  

这是一个普通 Fortran 函数的包装示例（请不要使用kind=8，原因在其他问题中有解释）：

module functions_to_integrate
   use iso_fortran_env
   use iso_c_binding


   integer, parameter :: wp = real64
contains

    pure function g(x)
      real(kind=wp) :: g
      real(kind=wp), intent(in) :: x
      g = sin(x)/x
    end function g

    function g_to_c(x) bind(C)
      real(kind=c_double) :: g_to_c
      real(kind=c_double), value :: x
      g_to_c = real(g(x),kind=c_double)
    end function

end module

program test

    use iso_c_binding

    use functions_to_integrate

    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
                real(kind=c_double), value :: x 
            end function
          end interface

          real (kind=c_double), value :: a,b
      end function my_cquad
    end interface

    real (kind=c_double) :: y,a,b

    a = 0 ; b = 1
    y = my_cquad(g_to_c,a,b)

    print *,y

end program test

附：我还在结尾之前删除了您的stop 和return 声明。不知何故，它总是让我发疯，但这可能只是我的强迫症。我太习惯在古代的旧程序中看到它了。

P.P.S：您可能希望看到由 Vincius Miranda http://www.lrz.de/services/software/mathematik/gsl/fortran/ 链接的 FGSL 接口包。我知道那个，但我主要是想指出错误，以便您可以自己制作类似的界面，没有现成的包可用。