Julia/Fortranのsubroutine呼び出し

Julia から Fortran の subroutine を呼び出す †

Julia から Fortran の subroutine を呼び出す

↑

はじめに †

ccall という関数で Fortran の subroutine と function を使うことができる．
使用法について，まずは公式ドキュメント

Calling C and Fortran Code

と下のリンクを読む．

↑

使用例1 テスト †

実行バージョンは Julia v1.3.1

↑

Fortran 側 †

例として下にあるような簡単な Fortran の subroutine を使う．

subroutine test1(n,a,b,c) bind(c, name="test1")
  implicit none
  integer , intent(in) :: n
  real*8, intent(in) :: a
  real*8, intent(out) :: b, c

  b = a**n
  c = n*a

  ! 確認
  print *, "b = ",b
  print *, "c = ",c

  return
end subroutine test1

整数 n と倍精度実数 a を入力して，倍精度実数 b, c を返すという subroutine．
このコードが書いてあるファイル( ccalltest.f90 とする)を，コンパイルして呼び出せるようにする．

ifort ccalltest.f90 -o ccalltest.so -fPIC -shared    # gfortran でも同じオプション

↑

Julia 側 †

この生成された ccalltest.so から test1 をJulia 上で呼ぶことを考える．
Fortran の integer, real*8 型はそれぞれ Julia の Int32, Float64 に該当するので，以下のように書く．

n = Int32(5)
a = Float64(2.0)
b = [Float64(0.0)] # 配列にする b = zeros(Float64, 1) 等でも可
c = [Float64(0.0)] # 配列にする

ccall((:test1, "./ccalltest.so"), Nothing, (Ref{Int32}, Ref{Float64}, Ref{Float64}, Ref{Float64}), n, a, b, c)

ここで最後の行に着目すると，ccall の引数は

ccall((function_name, library), returntype, (argtype1, ...), argvalue1, ...)

であり，注意すべきは Fortran の subroutine では returntype を Nothing に，argtype については Ref{Float64}のように Ref で変数型名を括る必要がある点．
この Ref{} がないと Julia が強制終了してしまうので注意．
subroutine ではなく function を ccall で呼び出す場合，returntype に Float64 等の型名を入力し，argtype については Ref{} なしで型名をそのまま入力．

上の Julia のコードを実行すると

b =    32.0000000000000     
c =    10.0000000000000     
[32.0]
[10.0]

となる．

もう1つ気をつける点は，値が返される変数は配列でなくても Vector{Float64} のように1要素の配列を作る必要がある点．
b, c を配列でなく

b = Float64(0.0)
c = Float64(0.0)

のようにしてしまうと，実行結果は

b =    32.0000000000000     
c =    10.0000000000000     
0.0
0.0

となり，処理は問題なく実行できるにも関わらず，値が Julia の変数に返ってこない．理由は（ちゃんと調べてないから）わからない．

↑

使用例2 実践 †

F77 形式で書かれている防災科研の DC3D モデル (Okada, 1992) を Julia で使うことを試みる．
ソースコードは断層モデルによる地殻変動計算プログラムに公開されている．
呼び出すためにまずはコンパイル．

ifort DC3Dfortran.f -fPIC -shared -o DC3Dfortran.f.so

次に，subroutine DC3D を呼び出す時の名前を，nm コマンドで確認する．

nm DC3Dfortran.f.so

実行すると

0000000000004640 T dc3d_

のようにサブルーチン名にアンダースコア("_")がついた文字列が出てきて，これを Julia 側から呼び出す時に使う．
FORTRAN のソースが単精度浮動小数点数なことに注意して，呼び出す function を作る．
引数が多いとだいぶ面倒くさい．

function DC3Df(x::T, y::T, z::T, depth::T, dip::T,
               L1::T, L2::T, W1::T, W2::T,
               Δu₁::T, Δu₂::T, Δu₃::T; α::T=2/3) where T <: Float64

    # convert
    α, x, y, z, depth, dip, L1, L2, W1, W2, Δu₁, Δu₂, Δu₃ = map(A->convert(Float32,A),
    [α, x, y, z, depth, dip, L1, L2, W1, W2, Δu₁, Δu₂, Δu₃])

    # args
    ux, uy, uz, uxx, uyx, uzx, uxy, uyy, uzy, uxz, uyz, uzz = [zeros(Float32,1) for i=1:12]
    IRET = zeros(Int32,1)

    # call
    ccall((:dc3d_, "./DC3Dfortran.f.so"), Nothing,
          ( # argtype
           Ref{Float32}, Ref{Float32}, Ref{Float32},
           Ref{Float32}, Ref{Float32}, Ref{Float32},
           Ref{Float32}, Ref{Float32}, Ref{Float32},
           Ref{Float32}, Ref{Float32}, Ref{Float32}, Ref{Float32}, # input
           Ref{Float32}, Ref{Float32}, Ref{Float32},
           Ref{Float32}, Ref{Float32}, Ref{Float32},
           Ref{Float32}, Ref{Float32}, Ref{Float32},
           Ref{Float32}, Ref{Float32}, Ref{Float32}, Ref{Int32} # output
           ),
           α, x, y, z, depth, dip, L1, L2, W1, W2, Δu₁, Δu₂, Δu₃,
           ux, uy, uz, uxx, uyx, uzx, uxy, uyy, uzy, uxz, uyz, uzz, IRET
          )


    # convert
    ux, uy, uz, uxx, uyx, uzx, uxy, uyy, uzy, uxz, uyz, uzz = map(A->convert(Float64,A[1]),
    [ux, uy, uz, uxx, uyx, uzx, uxy, uyy, uzy, uxz, uyz, uzz])

    # return
    return [ux, uy, uz, uxx, uyx, uzx, uxy, uyy, uzy, uxz, uyz, uzz]
end

適当な値を入力してこの function を使ってみる．

X = -99.0:2.0:99.0
Y = X
nx = length(X)
ny = length(Y)

L1 = -50.0
L2 = 50.0
W1 = -50.0
W2 = 50.0

z = -30.0
depth = 50.0
dip =70.0
Δu₁ = 200.0
Δu₂ = -150.0
Δu₃ = 100.0

uz = [DC3Df(X[j], Y[i], z, depth, dip, L1, L2, W1, W2, Δu₁, Δu₂, Δu₃)[3] for i=1:ny, j=1:nx]

# plot
using Plots
plotlyjs()
plt = plot(X, Y, uz, linetype=:surface, c=:balance, clims=(-50.0,50.0))

Plots.jl

Takuya Miyashita