Aufruf einer Fortran-Funktion von Julia, Rückgabe eines Arrays: unbekannte Funktion, segfault?

Question

Ich möchte Funktionen in meiner Fortran-Bibliothek von Julia aufrufen. In diesem Fall habe ich eine Funktion eye, die eine Ganzzahl annimmt und ein zweidimensionales Array von ganzen Zahlen zurückgibt.

Das Fortran-Modul wird in eine gemeinsame Bibliothek

$ gfortran -shared -fPIC -o matrix_routines.so matrix_routines.f90 

mit kompiliert und danach bin ich versucht, es von der interaktiven Julia Interpreter zu nennen wie das (Name von nm erhalten):

julia> n=5 
5 

julia> ccall((:__matrix_routines_MOD_eye, "/path/to/library/matrix_routines.so"), Array{Int64,2} , (Ptr{Int64},), &n) 

Dies führt jedoch sofort dazu, dass Julia bei mir einen Fehler verursacht:

signal (11): Segmentation fault 
__matrix_routines_MOD_eye at /path/to/library/matrix_routines.so (unknown line) 
anonymous at no file:0 
unknown function (ip: -1137818532) 
jl_f_top_eval at /usr/bin/../lib/julia/libjulia.so (unknown line) 
eval_user_input at REPL.jl:53 
jlcall_eval_user_input_19998 at (unknown line) 
jl_apply_generic at /usr/bin/../lib/julia/libjulia.so (unknown line) 
anonymous at task.jl:95 
jl_handle_stack_switch at /usr/bin/../lib/julia/libjulia.so (unknown line) 
julia_trampoline at /usr/bin/../lib/julia/libjulia.so (unknown line) 
unknown function (ip: 4199613) 
__libc_start_main at /usr/bin/../lib/libc.so.6 (unknown line) 
unknown function (ip: 4199667) 
unknown function (ip: 0) 
zsh: segmentation fault (core dumped) julia 

Rufe ich die Funktion falsch an? Wie lautet der korrekte Name der Funktion? (Es scheint nicht nur eye zu sein, da das auch nicht funktioniert.)

Als eine zusätzliche Frage: tut Julia irgendetwas mit der Speicherorientierung der resultierenden Arrays? Fortran und Julia sind beide Spalte-Major, aber ich frage mich, ob wegen ccall() Julia denken könnte, dass es sie transponieren sollte?

module matrix_routines 
    implicit none 

    private 

    public :: eye 

    contains 

     pure function eye(n,offset) result(um) !{{{ 
      integer, intent(in) :: n 
      integer, intent(in), optional :: offset 

      integer, dimension(n,n) :: um 

      integer :: i, l, u, os 

      um = 0 

      l = 1 
      u = n 
      os = 0 

      if (present(offset)) then 
       os = offset 
      end if 

      if (abs(os) < n) then 
       if (os > 0) then 
        u = n - os 
       else if (os < 0) then 
        l = 1 - os 
       end if 

       do i=l, u 
        um(i, i+os) = 1 
       end do 
      end if 

     end function eye !}}} 
end module matrix_routines 

Answer 1

Es gibt ein paar Probleme mit Ihrem Ansatz. Die direkte Rückgabe eines Arrays an Julia ist problematisch, da Fortran-Arrays nur dann mit C interoperabel sind, wenn bestimmte Bedingungen erfüllt sind. Wenn Sie das Array interoperabel zu machen (bind(C) zu Ihrer Prozedur hinzufügen und dem Array eines C-Typen geben) die Compiler (gfortran) wird ich beschweren:

Error: Return type of BIND(C) function 'um' at (1) cannot be an array 

Um dies zu umgehen wir das Array durch ein Scheinargument zurückkehren können. Sie sollten dies zu einem intent(inout)-Argument machen und das Array in julia konstruieren, um Speicher-/Oszilloskopprobleme beim Erstellen des Arrays in Fortran zu vermeiden.

Zweitens ist das optionale Argument problematisch und skimming die Julia Dokumente Ich bin mir nicht sicher, es wird sogar unterstützt. Beachten Sie, dass nicht einmal Fortran Fortran mit optionalen Argumenten ohne eine explizite Schnittstelle aufrufen kann und da Julia nicht mit der Datei .mod interagiert und die C-Methode zu erwarten scheint, wird es wahrscheinlich nicht funktionieren (und Fortran 2008 15.3.7 p2.6 scheint zu sagen, dass es nicht unterstützt wird). Es gibt jedoch Problemumgehungen - Sie können mehrere Fortran-Prozeduren mit einer unterschiedlichen Anzahl von Argumenten erstellen und dann die Prozedur mit optionalen Argumenten von ihnen aufrufen.

Zuerst diese Fortran-Modul betrachten, die mit Ihrem Beispiel begonnen, ist aber nur abgespeckte, was notwendig ist, die Interop zu demonstrieren:

module matrix_routines 
    implicit none 

    private 
    public :: eye 

contains 

    pure subroutine eye(n,um) bind(C,name="eye") !{{{ 
    use, intrinsic :: iso_c_binding, only: c_int 
    implicit none 
    integer(c_int), intent(in) :: n 
    integer(c_int), intent(inout), dimension(n,n) :: um 

    integer :: i, j 

    do j=1,n 
     do i=1,n 
      um(i,j) = i+j 
     end do 
    end do 

    end subroutine eye !}}} 
end module matrix_routines 

Bitte beachte, dass ich um zu sein ein inout Scheinargument verschoben habe Da wir keinen Wert zurückgeben, wurde die Prozedur in eine Subroutine geändert. Ich habe auch das optionale Argument entfernt. Ich habe auch C-Interop-Typen verwendet und einen C-Namen an die Prozedur gebunden. Sie können dies wie in Ihrer Frage zusammenstellen.

In Julia können Sie jetzt wie folgt vorgehen:

julia> n = 2 
2 

julia> um = zeros(Int32, n, n) 
2x2 Array{Int32,2}: 
0 0 
0 0 

julia> ccall((:eye, "matrix_routines.so"), Void, (Ptr{Int32}, Ptr{Array{Int32,2}}), &n, um) 

julia> um 
2x2 Array{Int32,2}: 
2 3 
3 4 

julia> n = 4 
4 

julia> um = zeros(Int32, n, n) 
4x4 Array{Int32,2}: 
0 0 0 0 
0 0 0 0 
0 0 0 0 
0 0 0 0 

julia> ccall((:eye, "matrix_routines.so"), Void, (Ptr{Int32}, Ptr{Array{Int32,2}}), &n, um) 

julia> um 
4x4 Array{Int32,2}: 
2 3 4 5 
3 4 5 6 
4 5 6 7 
5 6 7 8 

Beachten Sie, dass wir die Funktion aufrufen kann als nur :eye da wir die bind(C,name="eye") C Interop in unserem Fortran verwendet.

Und schließlich, wenn wir die do-Schleife in meinem Fortran Beispiel ändern um(i,j) = i*10+j zu sein, können wir sehen, dass keine Umsetzung in der Anordnung passiert:

julia> um 
3x3 Array{Int32,2}: 
11 12 13 
21 22 23 
31 32 33 

---

Der besondere Grund für Ihre segfault haben könnte war eine Reihe von Dingen - nicht übereinstimmende Datentypen, Probleme mit dem Rückgabetyp, Probleme mit dem optionalen Argument oder Nichtübereinstimmung im tatsächlichen Aufruf und Variablen bestanden.