從Fortran調用GSL例程CQUAD

Question

我正在嘗試從Fortran中調用GSL rountine CQUAD。我的想法是編寫一個調用gsl rountine的.c子例程，並依賴於函數和邊界。兩個問題：我對c和fortrans的iso_c_binding知之甚少。我嘗試如下：

一個簡單的調用程序（在No output from a Fortran program using the Gnu Scientific Library via a c wrapper類似MSB的職位）：

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到綜合功能：

測試：

> 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 

P.S.我還在結束之前刪除了您的stop和return聲明。不知何故，它總是讓我發瘋，但那可能只是我的強迫症。我很習慣從遠古時代的舊節目中看到它。

P.P.S：您可能希望看到由Vincius Miranda http://www.lrz.de/services/software/mathematik/gsl/fortran/鏈接的FGSL接口包。我知道這一點，但我主要嘗試指出錯誤，以便您可以自己製作類似的接口，而不需要現成的包裝。