Jupyter 上の IJulia 経由の Fortran 呼び出しで図を描く

能書き

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 関数の引数リストはタプルなので、変数１個の時にはコンマで終わる「( 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言語に否定的な方々の意見　１年以上前のものですが
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 代替のニッチ向きとのこと。

私は配列１始まりで、ブロックを end で閉じる方が落ち着くんですけど。