#author("2020-01-14T09:44:42+09:00","default:Miyashita","Miyashita") #author("2022-04-13T00:25:17+09:00","default:Miyashita","Miyashita") *Julia から Fortran の subroutine を呼び出す [#x947bacd] #contents **はじめに [#r30b4559] ccall という関数で Fortran の subroutine と function を使うことができる．~ 使用法について，まずは公式ドキュメント -[[Calling C and Fortran Code>https://docs.julialang.org/en/v1/manual/calling-c-and-fortran-code/]]~ と下のリンクを読む． -[[JuliaからFortranのsubroutineを呼び出す - Qiita>https://qiita.com/cometscome_phys/items/d098b00ab97f44916783]] -[[JuliaからFortranの構造体を引数にもつsubroutineを呼び出す - Qiita>https://qiita.com/cometscome_phys/items/fdf5720764b579ae8640]] **使用例1 テスト [#qff22c28] 実行バージョンは Julia v1.3.1~ ***Fortran 側 [#gc22704d] 例として下にあるような簡単な Fortran の subroutine を使う． #codeprettify(lang-fortran){{ 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 とする)を，コンパイルして呼び出せるようにする．~ #codeprettify{{ ifort ccalltest.f90 -o ccalltest.so -fPIC -shared # gfortran でも同じオプション }} ***Julia 側 [#j0e7f576] この生成された ccalltest.so から test1 をJulia 上で呼ぶことを考える．~ Fortran の integer, real*8 型は それぞれ Julia の Int32, Float64 に該当するので，以下のように Fortran の integer, real*8 型は それぞれ Julia の Int32, Float64 に該当するので，以下のように書く． #codeprettify(lang-julia){{ 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 の引数は #codeprettify(lang-julia){{ ccall((function_name, library), returntype, (argtype1, ...), argvalue1, ...) }} であり，注意すべきは Fortran の subroutine では returntype を Nothing に，argtype については Ref{T}のように Ref で括ること．~ であり，注意すべきは Fortran の subroutine では returntype を Nothing に，argtype については Ref{Float64}のように Ref で変数型名を括る必要がある点．~ この Ref{} がないと Julia が強制終了してしまうので注意．~ Fortran の function を ccall する場合，returntype に Float64 等の型名を入力し，argtype については Ref{} なしで型名をそのまま入力．~ subroutine ではなく function を ccall で呼び出す場合，returntype に Float64 等の型名を入力し，argtype については Ref{} なしで型名をそのまま入力．~ ~ 上の Julia のコードを実行すると b = 32.0000000000000 c = 10.0000000000000 [32.0] [10.0] となる．~ ~ もう1つ気をつける点は，値が返される変数は配列でなくても Vector{Float64} のように1要素の配列を作る必要がある点．~ b, c を配列でなく #codeprettify(lang-julia){{ b = Float64(0.0) c = Float64(0.0) }} のようにしてしまうと，実行結果は b = 32.0000000000000 c = 10.0000000000000 0.0 0.0 となり，処理は問題なく実行できるにも関わらず，値が Julia の変数に返ってこない．理由は（ちゃんと調べてないから）わからない．~ ~ **使用例2 実践 [#n11c0adc] F77 形式で書かれている防災科研の DC3D モデル (Okada, 1992) を Julia で使うことを試みる．~ ソースコードは [[断層モデルによる地殻変動計算プログラム>http://www.bosai.go.jp/study/application/dc3d/DC3Dhtml_J.html]] に公開されている．~ 呼び出すためにまずはコンパイル．~ #codeprettify{{ ifort DC3Dfortran.f -fPIC -shared -o DC3Dfortran.f.so }} 次に，subroutine DC3D を呼び出す時の名前を，nm コマンドで確認する．~ #codeprettify{{ nm DC3Dfortran.f.so }} 実行すると 0000000000004640 T dc3d_ のように サブルーチン名にアンダースコア("_")がついた文字列が出てきて，これを Julia 側から呼び出す時に使う．~ FORTRAN のソースが単精度浮動小数点数なことに注意して，呼び出す function を作る．~ 引数が多いとだいぶ面倒くさい． #codeprettify{{ 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 を使ってみる． #codeprettify(lang-julia){{ 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)) }} #htmlinsert(Julia/DC3D.html)