JuliaからFortranのsubroutineを呼び出す

Julia 1.0を使っているときに、Fortranで作ったルーチンを呼び出したいこともあるかと思います。
そのような時にどうするか、ということをまとめました。例えば、MKLを使ったFortranの高速なコードをわざわざJuliaに移植せずにそのまま使いたい、などがあるかもしれません。

まず、以下のようなFortran90のmoduleがあるとしましょう。

module multi
  contains
    subroutine mat(val,col,row,val_l,col_l,row_l,x,n,b) !複素数行列の行列ベクトル積を計算
      implicit none
      integer(8),intent(in)::val_l,col_l,row_l
      complex(8),intent(in)::val(1:val_l)
      integer(8),intent(in)::col(1:col_l),row(1:row_l)
      integer(8),intent(in)::n
      complex(8),intent(in)::x(1:n)
      complex(8),intent(out)::b(1:n)
      integer::i,j
      b = 0d0

      call mkl_zcsrgemv("T", n, val, col, row, x, b) !MKLのスパースBLAS レベル2を使用

      return
    end subroutine mat
end module multi

このmoduleの中のsubroutineのmatでは、MKLのスパースBLAS レベル2を利用することによって、疎行列ベクトル積を計算しています。最近はMKLの利用は無料となっています。
https://software.intel.com/en-us/mkl
このコードをコンパイルするために、ifortを使っても良いのですが、ifortは有料ですので、gfortranを使います。

gfortran test.f90 -o test.so -shared -fPIC -lmkl_intel_ilp64  -lmkl_sequential -lmkl_core -lpthread -lm -ldl

ここで、このmoduleを共有ライブラリ化するために、-shared -fPIC というオプションをつけています。後ろのmkl関連は、OSなどによって異なると思います。今回はMacのgfortranでMKLを呼んでいます。なお、ilpにするとinteger(8)、lpにするとintegerをMKLのルーチンに入れることになりますので、MKLのルーチンに入れる整数の引数に合わせてオプションを指定しないと変な結果が返ってきますので注意してください。
ここでできたtest.soは

nm test.so 

とすることによって、中に入っているsubroutineの名前を見ることができます。gfortranかifortかで変わります。
このコードをJulia側から呼び出すには、

using SparseArrays
using Random

Random.seed!(1234) #乱数シードを固定
n = 10
A = spzeros(ComplexF64,n,n)
for i=1:n
    j = rand(1:n)
    A[i,j] = rand()+im*rand()
end
#10x10の適当な行列を用意した。

col = A.colptr #疎行列の列
val = A.nzval #疎行列の値
row = A.rowval #疎行列の行
#行列ベクトル積 b = A*xを実行したい
x = rand(Float64,n)+im*rand(Float64,n) 
b = zeros(ComplexF64,n) 
val_l = length(val)
row_l = length(row)
col_l = length(col)
println("Ax ", A*x)

ccall((:__multi_MOD_mat,"./test.so"),Nothing, #gfortranでコンパイルした場合
#ccall((:multi_mp_mat_,"./test.so"),Nothing, #ifortだと:multi_mp_mat_となる。
    (Ref{ComplexF64}, #val
    Ref{Int64}, #col
    Ref{Int64}, #row
    Ref{Int64}, #val_l
    Ref{Int64}, #col_l
    Ref{Int64}, #row_l
    Ref{ComplexF64}, #x
    Ref{Int64}, #n
    Ref{ComplexF64}), #b
    val,col,row,val_l,col_l,row_l,x,n,b)
println("Ax from Fortran ", b)

とすればよいです。
ここで、Fortranを呼び出す時には、Ref{Float64}やRef{ComplexF64}などとします。ここでのFloat64は引数が配列であれば配列の要素の型となります。Int64はFortranではinteger(8)となることに注意してください。どの型であれRef{}の形で定義しておけば、問題なく引数を渡すことができます。

Nothingは古いバージョンのJuliaではVoidでした。そのため、呼び出しの説明でVoidを使っている場合がありますが、Julia 1.0以降ではNothingを使います。これは、関数の返り値の型です。

疎行列格納形式の注意

今回は疎行列を演算するFortranコードを呼び出していますので、Julia側の疎行列の中身をval,row,colと取り出しています。ここで注意しなければならないのは、FortranでのMKLの疎行列サブルーチンは行格納型（CSR）ですが、Juliaでは列格納型(CSC)であることです。しかし、これは特に大きな問題になりません。というのは、FortranのMKLのルーチンmkl_zcsrgemvは行列の転置の演算に対応しているからです。よって、mkl_zcsrgemvの引数のrowとcolの場所を入れ替えて"T"を引数にすれば、Juliaでの疎行列をそのまま入れることができます。

追記

cure_honeyさんに教えていただきましたが、ccallで呼び出すときのFortranのコードはFortran2003以降に対応していれば（最近のgfortranで使える事を確認しました）、

subroutine mat(val,col,row,val_l,col_l,row_l,x,n,b) bind(c, name = 'mat')

とbindを使えば名前をつけられるようです。これを使えば、呼び出しの際に、:__multi_MOD_matの代わりにmatを使って呼び出すことができるようです。