Fortran, Python, Juliaの計算速度比較:
モンテカルロ法を例に

1. 問題設定: モンテカルロ法を用いた円周率計算

Fortran, Python, Juliaの計算速度を比較するために、(0,1)×(0,1)の領域にランダムに点を打っていき、原点(0,0)からの距離が1以下の点を数えることで円周率を評価する比較的単純なモンテカルロ計算を用います。打った点の総数を\(N_{total}\)、円の中の点の総数を\(N_{inside}\)とすると、円周率(π)は次のように近似されます。 \[ \pi \approx \frac{4N_{inside}}{N_{total}}. \tag {1} \] 非常にたくさんの点をランダムに打っていくことで、より精度よく円周率を求めることが出来ます。

今回の計算速度の比較では、異なる言語の速度を公平に比較するための、共通の乱数生成ルーチンを用います。具体的には、線形合同法を用いた乱数生成ルーチンを用意してモンテカルロ計算を行います。

2. Fortran, Python, Juliaのコード

Fortran, Python, Juliaの計算速度を比較するために、モンテカルロ計算により円周率の近似値を求めます。今回の計算では、それぞれのコードで共通の乱数生成ルーチンとして線形合同法を実装して、10⁹個の乱数を生成して計算を行っています。計算の詳細はそれぞれのコードを確認してください。

Fortranコード

(mc.f90)

program main
  implicit none
  integer :: seed, i, n
  integer :: num_inside, num_points
  real(8) :: x, y, r2

  seed = 20231226
  num_points = 1000000000
              
  num_inside = 0
  do i = 1, num_points
    call LCGs(seed, x)
    call LCGs(seed, y)
    r2 = x**2 + y**2
    if(r2<1d0)num_inside = num_inside + 1
  end do

  write(*,*)"pi=",4*dble(num_inside)/num_points


  contains
! Linear congruential generators (LCGs)
! Parameters are provided by Park and Miller
! See https://c-faq.com/lib/rand.html
    subroutine LCGs(seed, rand_num)
      implicit none
      integer,parameter :: a = 48271
      integer,parameter :: m = 2147483647
      integer,parameter :: q = m/a
      integer,parameter :: r = mod(m,a)
      integer,intent(inout) :: seed
      real(8),intent(out) :: rand_num
      integer :: hi, lo, test

      hi = seed/q
      lo = mod(seed, q)
      test = a * lo - r * hi

      if(test>0)then
        seed = test
      else
        seed = test + m
      end if

      rand_num = dble(seed)/m

    end subroutine LCGs
end program main

Pythonコード

(mc.py)

import numpy as np

# Linear congruential generators (LCGs)
# Parameters are provided by Park and Miller
# See https://c-faq.com/lib/rand.html
def LCGs(seed):
    a = 48271
    m = 2147483647
    q = m // a
    r = m % a

    hi = seed // q
    lo = seed % q
    test = a * lo - r * hi

    if test > 0:
        seed = test
    else:
        seed = test + m

    rand_num = float(seed) / m
    return seed, rand_num


seed = 20231226
num_points = 1000000000
num_inside = 0

for i in range(num_points):
    seed, x = LCGs(seed)
    seed, y = LCGs(seed)
    r2 = x**2 + y**2
    if r2 < 1.0:
        num_inside += 1

print("pi =", 4.0 * float(num_inside) / num_points)

Juliaコード

(mc.jl)

# Linear congruential generators (LCGs)
# Parameters are provided by Park and Miller
# See https://c-faq.com/lib/rand.html
function LCGs(seed::Integer)
    a = 48271
    m = 2147483647
    q = m ÷ a
    r = m % a

    hi = seed ﾃｷ q
    lo = seed % q
    test = a * lo - r * hi

    if test > 0
        seed = test
    else
        seed = test + m
    end

    return seed, Float64(seed) / m
end

function main()
    seed = 20231226
    num_points = 1000000000
    num_inside = 0

    for i in 1:num_points
        seed, x = LCGs(seed)
        seed, y = LCGs(seed)
        r2 = x^2 + y^2
        if r2 < 1.0
            num_inside += 1
        end
    end

    println("pi=", 4.0 * num_inside / num_points)
end

main()

3. コンパイル、及び計算の実行

今回の比較では、Fortranのコンパイラとしてgfortran (GNU Fortran)とifx (Intel Fortran; 旧ifort)を用いてます。GNU Fortranは、コンパイル時に計算速度を向上させるための最適化オプションを指定することが出来ますが、今回はO0, O1, O2, 及びO3という4つの最適化オプションについてそれぞれ計算速度を計測することで、コンパイルオプションによる計算速度の影響も調べることにします。また、Intel Fortranを用いる場合は、これらのコンパイルオプションに加えて-xHOSTという実行プロセッサに応じた最適化を行うオプションを追加して検証を行います。

また、Pythonコード(mc.py)はターミナル上でpython mc.pyとして実行し、Juliaコード(mc.jl)についても同様にターミナル上でjulia mc.jlとして実行します。JuliaとFortranの計算速度の比較の際に、Juliaのコンパイル時間について気になる方がいるかもしれませんが、この点については今回の比較結果には影響がないことを確認したので、その点については後述します。

今回用意したFortranコードmc.f90, Pytohnコードmc.py, 及びJuliaコードmc.jlを実行すると、10⁹個の乱数を生成してモンテカルロサンプリングを行うことで円周率の近似値を求め、画面に得られた値がそれぞれ出力されます。3つのコードの実行結果を比較することで、それぞれのコードが正しく動作していることを確認することが出来ます。

4. Fortran, Python, Juliaの計算速度の比較

Fortran, Python, Juliaの計算速度の比較を行うために、上記の3つの計算コードをそれぞれ5回実行し、各回毎に計算時間の計測を行いました。以下の表に、今回の計測によって得られた実行時間の平均値と標準誤差を示します。

上の表で、まずFortranのコンパイラ、及びコンパイルオプションごとの結果を比較してみましょう。今回、比較に用いるFortranコンパイラとして、GNU Fortranコンパイラ (gfortran) とIntel Fortran コンパイラ (ifx; 旧ifort)を採用しています。GNU Fortranの場合、-O3のオプションを指定した場合の結果が最速で、実行時間は7.91秒となりました。これに対して、Intel Fortranの場合は-O3のオプションを指定した場合の結果が最速で、実行時間は8.39秒でした。興味深いことに、今回の線形合同法を用いたモンテカルロ計算の場合、Intel FortranコンパイラよりもGNU Fortranコンパイラの方が計算速度が6%程速くなっています。これは、今回用いた線形合同法が整数演算を主とする計算内容となっており、GNU FortranがIntel Fortranに比べて、整数演算に対する最適化が進んでいることを示唆していると考えられます。

次に、Fortran (gfortran -O3)とPythonの実行時間を比較してみると、Fortranの方がPythonと比べて115倍ほど速いです。また、Fortran (gfortran -O3)とJuliaの実行時間を比較してみると、Fortranの方がJuliaと比べて1.4倍ほど速いです。さらに、JuliaとPythonの実行時間を比較すると、Juliaの方がPythonと比べて79.9倍ほど速いです。

5. 速度比較の感想とコメント

今回の比較では、GNU FortranがIntel Fortranよりも高速になる場合があるという興味深い結果を得ることが出来ました。それぞれのコンパイラがどのような演算に重きを置いているのかが異なることで、このような結果が得られたと考えられます。例えば、Intel Fortranの場合は、科学技術計算で良く用いられる浮動小数演算に重きを置いているために、整数演算の性能がGNU Fortranの性能に追いついていない等のシナリオが考えられます。様々な条件で比較を行うことで、それぞれのコンパイラの特徴が明らかになるかもしれません。

また、今回のようなモンテカルロ計算の場合、JuliaよりもFortranの方が1.4倍程高速に計算できるということが分かりました。この1.4倍という数字をどう評価するかは難しいですが、仮に、あなた自身の研究・業務で用いているコードが1.4倍の実行時間が掛かるようになった場合にどう感じるかを想像するとJuliaの善し悪しが見えてくるかもしれません。

6. FortranとJuliaのもう少し真面目な比較

上記の比較では、Juliaコードをjulia mc.jlとして実行して実行時間を計測しているので、コンパイル時間が実行時間に含まれています。したがって、「JuliaとFortranの正確が出来ていない」と感じる人もいるかもしれません。そこで、追加の解析としてモンテカルロ法のサンプル数を2倍にした計算で実行時間を計測することで、実質的な計算時間を2倍にし、相対的にJuliaのコンパイル時間の影響を抑えた調査を行いました。また、Fortranの方はgfortran -O3のコンパイルコマンドを用いて比較を行います。

追加の解析の結果を下の表に示しました。実質的な計算量を2倍にした計算の計測結果を用いてFortranとJuliaの計算速度を比較しても、Fortranの方がJuliaよりも1.4倍ほど速いということが確認できました。

7. 比較の実行環境について

今回の検証は、2023年12月28日に以下のような環境で実行しています。

• CPU: Intel(R) Core(TM) i5-1145G7 @ 2.60GHz
• OS: Ubuntu 20.04.6
• GNU Fortran: gfortran version 9.4.0
• Intel Fortran: ifx 2024.0.2 20231213
• Python: Python 3.8.10
• Julia: Julia version 1.7.3