Calling Fortran subroutine from Julia passing array arguments

General Usage
1

Fortran subroutine looks like this:

 subroutine MatMult(a, b, m, n, k, answer)
        !DEC$ ATTRIBUTES DLLEXPORT :: MatMult
        !DEC$ ATTRIBUTES ALIAS: 'MatMult' :: MatMult
        !DEC$ ATTRIBUTES REFERENCE :: a
        !DEC$ ATTRIBUTES REFERENCE :: b
        !DEC$ ATTRIBUTES REFERENCE :: answer
        !DEC$ ATTRIBUTES VALUE :: m, n, k
        integer, intent(in) :: m, n, k
        real(8), intent(in) :: a(m, n)
        real(8), intent(in) :: b(n, k)
        real(8), intent(out) :: answer(m, k)
        !this is just a test
        answer = transpose(matmul(transpose(a), transpose(b)));
    end subroutine

Fortran’s code is compiled from f90 into a dll.
In Julia I call MatMult passing my arrays

function mult()
    matrixSize = 2
    mA = [[2 3] ; [4 5]] 
    mB =  [[2 1] ; [4 1]] 
    mC = zeros(Float64, matrixSize, matrixSize)
    ccall((:MatMult, "FortranWorker.dll"), Cvoid, (Ptr{Matrix{Float64}}, Ptr{Matrix{Float64}}, Int64, Int64, Int64,Ptr{Matrix{Float64}}), mA, mB, matrixSize, matrixSize, matrixSize, mC)
end

The array returned contains all zeroes. Unfortunately I can’t step into fortran code when debugging in VS Code. I’ve tried compiling fortran code as .so library - not the dll - but then Julia pretends not see my .so library. Could be because fortran is in f90 format and not in f95?
Any suggestions are appreciated!!!

2

This works:

% gfortran teste.f90 -shared -o test.so

julia> a = rand(3,3); b = rand(3,3); 

julia> answer = zeros(3,3);

julia> ccall((:matmult_, "./test.so"), Nothing, (Ref{Float64}, Ref{Float64}, Ref{Int}, Ref{Int}, Ref{Int}, Ref{Float64}), a, b, 3, 3, 3, answer)

julia> answer
3×3 Matrix{Float64}:
 0.181291  0.829994  1.25334
 0.311001  0.804451  0.874293
 0.386725  0.91879   0.93275

I’m not sure if this is a toy example, but of course you can do in Julia:

julia> transpose(a) * transpose(b)
3×3 Matrix{Float64}:
 0.181291  0.311001  0.386725
 0.829994  0.804451  0.91879
 1.25334   0.874293  0.93275

#or

julia> using LinearAlgebra

julia> mul!(answer, transpose(a), transpose(b))
3×3 Matrix{Float64}:
 0.181291  0.311001  0.386725
 0.829994  0.804451  0.91879
 1.25334   0.874293  0.93275

ps: I didn’t transpose the answer in my examples, and I noticed not that you do it in the fortran code.

3

Where is your test.so located? I put mine in the same folder as Julia source code - but it tells me it cannot find that file

4

Here in the same directory where I started Julia. Note the ./ in ./test.so, without it it does not find the library.

5

julia> ccall((:MatMult, “./test.so”), Nothing, (Ref{Float64}, Ref{Float64}, Ref{Int}, Ref{Int}, Ref{Int}, Ref{Float64}), a, b, 3, 3, 3, answer)
ERROR: could not load library “./test.so”
The specified module could not be found.

So sorry for taking your time on this - not sure why it’s not finding that library. Could be smth wrong with the lib itself? I did get compilation warning when creating the .so from fortran source: Warning: Legacy Extension: REAL array index at (1)

6

I think you have not started Julia in the same directory where the library is.

Also, you need to use :matmult_ (like I posted), otherwise it does not work either (don’t ask me why one needs that underline at the end).