能書き
IPython がいつのころか Jupyter に名前を変えて、Python に限定されないプラットフォームに進化しております。
図表と記号処理と数値計算が混じり合った処理は、旧来 Maple/Mathematica/Matlab 等で行なわれてきましたが、年々上がる年貢の賦課に耐えかねて、これらに変わって Jupyter 上の寄せ集めでやろうという動きもあるようです。
一方モダン Fortran を使うには、IPython は微妙に Python 経由での Fortran 呼び出しが面倒でしたが、Jupyter 上の IJulia 経由だと、Fortran 呼び出しが楽だという噂を小耳にはさんだので、挑戦してみます。
以下のサイトを参考にしました。
divrot.hatenablog.com
Calling C and Fortran Code — Julia Language 0.4.8-pre documentation
下準備
上記の参考サイトを見ればわかるように、根本的なからくりはダイナミック・ライブラリを作ってそれを呼び出す形になるため、Fortran コンパイラの吐き出すライブラリがちゃんと呼び出せる必要があります。
Windows上では、DLL の形式が MSC のリンカ―のバージョンなどに依存して、いらぬ気苦労をしそうなので、敬して遠ざけることとします。そこで ubuntu on bash on windows 上で Jupyter/IPython/IJulia 環境を作って、bash on windows 側の gfortran を使うことにしました。
ネット情報を頼りにやりましたが、何ヶ所か躓きました。ひどく時間はかかりましたが、エラーメッセージを我慢して読むことで、なんとか解決できました。
重症例
- ZMQ ライブラリ構築失敗 ⇒ g++ がなかったため。
- PyPlot が IPython では作画可能ながら IJulia で失敗 ⇒ OpenGL がらみのライブラリが足りなかった。エラーで無いとい言われたものを検索、Ubuntu apt-get でパッケージをインストール。
[追記:sudo apt install libzmq5
pip uninstall matplotlib
pip install matplotlib
matplotlib のバージョン違いで問題になったりする。
しかし環境依存。
]
Fortran
とりあえず、エラトステネスのふるいで素数を求めて個数をプロットしてみることにします。Fortran 側の関数で個数を求めて、IJulia 側では、関数呼び出しと作図のみを行うことにします。
Fortran 側でも、デバッグの為メインプログラムで個数を求めてみます。
F2003 以降の C言語互換用の BIND を使うことで名前問題は解決するので、関数は Module 内においても大丈夫。またダイナミック・ライブラリにするので、メイン・プログラム付きのままコンパイルしても大丈夫(のはず)です。関数の引数は値呼び出しにしました。
module m_prime implicit none contains integer function nprim(n) bind(c, name = 'nprim') integer, value :: n logical(1), allocatable :: tab(:) ! Sieve of Eratosthenes integer :: i allocate(tab(n)) tab = .true. tab(1) = .false. do i = 1, int(sqrt(real(n))) if (tab(i)) tab(i**2::i) = .false. end do nprim = count(tab) end function nprim end module m_prime program eratos3 use m_prime implicit none integer :: i do i = 1, 8 print *, nprim(10**i) end do end program eratos3
これは windows 側の visual studio で作りました。
これを、bash on windows 側に送って、gfortran でコンパイルします。
gfortran prime.f90 -o prime.so -shared -fPIC
これで、ダイナミック・リンク用のライブラリが出来ました。
IJulia 側での処理
とりえあえず、関数を呼び出してみます。Fortran の引数は基本的に参照呼び出しなのですが、ここでは値呼び出しにすることで厄介をさけています。
よく分からないでの、ライブラリは絶対 path を与えておきました。
呼び出し時の fortran 関数の引数リストはタプルなので、変数1個の時にはコンマで終わる「( Int32, )」ということにしばらく気づきませんでした。
n = 100
t = ccall( (:nprim, "/mnt/d/home/prime.so"), Int32, ( Int32, ), n )
これで百までの素数の個数である 25 が返されます。
次にこれをもとに Fortran の関数を呼び出す Julia の関数を定義します。
function np(i) n = Int32(10^i) return ccall( (:nprim, "/mnt/d/home/prime.so"), Int32, ( Int32, ), n ) end
10 から 10^9 までの 10のべき乗個以下の素数の個数を求め配列に入れます。
x = 1:9
y = np.(x)
print(y)
出てきた結果は以下のようになります。
Int32[4, 25, 168, 1229, 9592, 78498, 664579, 5761455, 50847534]
出た結果を対数グラフに出力します。
using PyPlot
grid("on")
xlabel("x")
ylabel("Number of Primes pi(10^x)")
ax = gca()
ax[:set_yscale]("log")
grid("on")
plot(x, y)
今回は引数の受け渡しなどが素朴だったので、比較的簡単に実行できました。
付録
Julia言語に否定的な方々の意見 1年以上前のものですが
Giving up on Julia
その中での引用先
http://www.davideaversa.it/2015/12/the-most-promising-languages-of-2016/
A review of the Julia programming language
遅く、バグが多くかつ放置気味で、言語仕様が MATLAB 的でダサいようです。
MATLAB 代替のニッチ向きとのこと。
私は配列1始まりで、ブロックを end で閉じる方が落ち着くんですけど。