Chiamare una funzione Fortran da Julia, restituendo un array: funzione sconosciuta, segfault?

Question

Voglio chiamare le funzioni nella mia libreria Fortran da Julia. In questo caso, ho una funzione eye che accetta un intero e restituisce un array bidimensionale di numeri interi.

Il modulo Fortran viene compilato in una libreria condivisa utilizzando

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

E da allora in poi sto cercando di chiamare da l'interprete interattivo Julia genere (nome ottenuto dal nm):

julia> n=5 
5 

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

Ciò, tuttavia, si traduce immediatamente in Julia lanciando un segfault a me:

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 

Sto chiamando la funzione nel modo sbagliato? Qual è il nome corretto della funzione? (Non sembra essere solo eye, perché non funziona neanche.)

Come domanda aggiuntiva: Julia fa qualcosa con l'orientamento della memoria degli array risultanti? Fortran e Julia sono entrambi importanti per la colonna, ma mi chiedo se a causa di ccall() Julia potrebbe pensare che dovrebbe metterli in pratica?

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

Ci sono un paio di problemi con il vostro approccio. Restituire un array direttamente su julia è problematico perché gli array Fortran non sono interoperabili con C a meno che non siano soddisfatte condizioni specifiche. Quando si effettua la matrice di interoperabilità (aggiungi bind(C) per la vostra procedura e dare la matrice di tipo C) il compilatore (gfortran) si lamenta:

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

Per ovviare a questo siamo in grado di restituire la matrice attraverso un argomento fittizio. Dovrai renderlo un argomento intent(inout) e costruire l'array in julia per evitare problemi di memoria/ambito con la creazione dell'array in Fortran.

In secondo luogo, l'argomento facoltativo è problematico e il passaggio dei documenti di Julia non sono sicuro che sia supportato. Nota che nemmeno Fortran può chiamare Fortran con argomenti opzionali senza un'interfaccia esplicita e poiché Julia non interagisce con il file .mod e sembra aspettarsi il modo C di fare le cose, probabilmente non funzionerà (e Fortran 2008 15.3.7 p2.6 sembra dire che non è supportato). Esistono tuttavia soluzioni alternative: è possibile creare più procedure Fortran con un numero variabile di argomenti e quindi richiamare la procedura con argomenti facoltativi.

primo luogo, considerare questo modulo Fortran, che è iniziato con il vostro esempio ma viene ridotto al solo ciò che è necessario per dimostrare l'interoperabilità:

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 

Si noti che ho spostato um ad essere un argomento fittizio inout e dal momento che non restituiamo un valore, la procedura è cambiata in una subroutine. Ho anche rimosso l'argomento opzionale. Ho anche usato i tipi di interoperabilità C e ho associato un nome C alla procedura. Puoi compilarlo come nella tua domanda.

In Julia ora è possibile effettuare le seguenti operazioni:

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 

Nota che possiamo chiamare la funzione come solo :eye dato che abbiamo usato l'interoperabilità bind(C,name="eye") C nel nostro Fortran.

E, infine, se cambiamo il ciclo do nel mio esempio Fortran essere um(i,j) = i*10+j, possiamo vedere che nessuna trasposizione avviene nella matrice:

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

---

La ragione particolare per il vostro segfault potrebbe avere sono state un numero di cose: tipi di dati non corrispondenti, problemi con il tipo restituito, problemi con l'argomento opzionale o mancata corrispondenza nella chiamata effettiva e variabili passate.