Julia から Fortran の subroutine を呼び出す

はじめに

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

と下のリンクを読む.

使用例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{T}のように Ref で括ること.
この Ref{} がないと Julia が強制終了してしまうので注意.
Fortran の 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

Front page   Edit Diff Attach Copy Rename Reload   New List of pages Search Recent changes   Help   RSS of recent changes
Last-modified: 2020-01-10 (Fri) 06:41:30 (47d)