【Fortran】数値計算ライブラリを使用する【LAPACK/BLAS】

はじめに

Fortranには長い歴史があり、先人たちが作成した効率の良いアルゴリズムで実装されたプログラムがライブラリとして提供されている場合があります。これらのライブラリを自分たちのプログラムで使わない手はありません。

今回は、LAPACK / BLASという線形代数の計算ライブラリを利用する方法と、これらのライブラリを使うといかに実行速度が速くなるかを示したいと思います。

動作検証済み環境

macOS Sonoma(14.7), gfortran (gcc version 14.1.0),
Processor: 2.4 GHz 8-Core Intel Core i9,
Memory: 32 GB 2400 MHz DDR4

Pythonで学ぶ偏微分方程式の数値シミュレーションの技術書を販売中

計算方法、グラフやアニメーションの作成方法、計算ログの残し方を惜しみなく解説しています！
理論と実装の溝を埋めてくれる良書です！

技術書ページへ

LAPACK / BLAS とは

LAPACK / BLAS (linear algebra package / basic linear algebra subprogram) とは、線形代数ライブラリです。線形代数の基本的な計算（連立一方程式の計算、固有値問題など）のルーチンをサポートしています。

オープンソースライブラリで、すべてのソースコードが開示されています。オリジナルのバージョンはFortran90で記述されています。

LAPACK / BLAS を使えるようにするには

LAPACとBLASは、mac の場合homebrewからインストールすることができます (BLAS はいくつかの実装がありますが、今回はOpenBLASを使います)。

$ brew install openblas
$ brew install lapack

ライブラリは、/usr/local/opt/openblas/lib および /usr/local/opt/lapack/lib 以下に置かれます。

BLASを使って行列の積を計算する

今回は、BLASを使って行列の積を計算してみます。

比較に使用する行列の積を計算する3つの方法

$l\times m$ 行列 $\mathsf{A}$ と $m\times n$ 行列 $\mathsf{B}$ の積を計算して、$l\times n$ 行列 $\mathsf{C} = \mathsf{A}\times B$ を得たいとします。このような最も簡単な行列にかんする計算を通して、BLASを使うことの有用性を示したいと思います。以下に、3つの方法で、倍精度実数の2次元配列 a(1:l, 1:m) と b(1:m, 1:n) の行列の積 c(1:l, 1:n) を計算する方法を提示します。

単純にループで実装

あまり何も考えずに、doループで実装すると

  do j = 1, n
    do i = 1, l
    c(i, j) = 0.0d0
      do k = 1, m
        c(i, j) = c(i, j) + a(i, k)*b(k, j)
      end do
    end do
  end do

という3重ループになります。

組み込み関数で実装

Fortranには組み込み関数 matmul が用意されていて

c = matmul(a, b)

とすることもできます。

BLASのサブルーチンで実装

BLASでは、dgemmというサブルーチンを使って $\mathsf{C} := \alpha \mathrm{op}(\mathsf{A})\mathsf{B} + \beta \mathsf{C}$ という計算が可能です。$\mathrm{op}(\mathsf{A})$ はサブルーチンの引数に指定がなされたときに $\mathsf{A}$ の転置を取る関数です。行列の積を計算したいときは、転置をとらず、$\alpha = 1, \beta = 0$ というふうに指定すればよいです。具体的には以下のように記述します：

  call dgemm('n', 'n', l, n, m, 1.0d0, a, l, &
  &          b, m, 0.0d0, c, l)

サブルーチンに与える引数は、こちらを見てください。

計算速度の比較

以上の3つの手法で計算速度を比較します。ランダムに行列の値を設定し、時間は組み込み関数の system_clock を使って計測します。

単純にループで実装

program mat_prod_simple
  implicit none

  integer(4) :: count1, count2, count_rate, count_max
  integer(4), parameter :: l=1000, m=2000, n=3000
  real(8) :: a(1:l, 1:m), b(1:m, 1:n), c(1:l, 1:n)
  integer(4) :: i, j, k

  call random_number(a)
  call random_number(b)

  call system_clock(count1)
  do j = 1, n
    do i = 1, l
    c(i, j) = 0.0d0
      do k = 1, m
        c(i, j) = c(i, j) + a(i, k)*b(k, j)
      end do
    end do
  end do
  call system_clock(count2, count_rate, count_max)

  write(*, *) ' cpu time = ', (count2-count1)/dble(count_rate), ' s'

end program mat_prod_simple

組み込み関数で実装

program mat_prod_matmul
  implicit none

  integer(4) :: count1, count2, count_rate, count_max
  integer(4), parameter :: l=1000, m=2000, n=3000
  real(8) :: a(1:l, 1:m), b(1:m, 1:n), c(1:l, 1:n)

  call random_number(a)
  call random_number(b)

  call system_clock(count1)
  c = matmul(a, b)
  call system_clock(count2, count_rate, count_max)

  write(*, *) ' cpu time = ', (count2-count1)/dble(count_rate), ' s'

end program mat_prod_matmul

BLASのサブルーチンで実装

program mat_prod_blas
  implicit none

  integer(4) :: count1, count2, count_rate, count_max
  integer(4), parameter :: l=1000, m=2000, n=3000
  real(8) :: a(1:l, 1:m), b(1:m, 1:n), c(1:l, 1:n)

  call random_number(a)
  call random_number(b)

  call system_clock(count1)
  call dgemm('n', 'n', l, n, m, 1.0d0, a, l, &
  &          b, m, 0.0d0, c, l)
  call system_clock(count2, count_rate, count_max)

  write(*, *) ' cpu time = ', (count2-count1)/dble(count_rate), ' s'
  
end program mat_prod_blas

それぞれのプログラムを、適当なファイル名で適当なディレクトリに保存します。例えば次にようにします。

単純ループ：~Desktop/LabCode/fortran/06_mat-prod/simple/main.f90
組み込み関数 (matmul)：~Desktop/LabCode/fortran/06_mat-prod/matmul/main.f90
BLAS (dgemm)：~Desktop/LabCode/fortran/06_mat-prod/blas/main.f90

コンパイルと実行

以上のプログラムをそれぞれ別々のフォルダに保存し、それぞれのフォルダでコンパイル、実行してみます。

単純ループと組み込み関数の場合は、それぞれプログラムを保存したディレクトリに移動し

$ gfortran -o main main.f90

としてコンパイルします。

BLASを使用するものは、リンクが必要です。

$ gfortran -L/usr/local/opt/openblas/lib -lopenblas -o main main.f90

として、コンパイルを実行してください。

ここで、-L/opt/local/lib -lopenblas の意味ですが /opt/local/lib 以下のディレクトリにあるlibopenblas.a というライブラリアーカイブ（サブルーチン群をおさめたオブジェクトファイルの集合）をリンクするという意味です。

指定の仕方にややクセがありますが、-L<ディレクトリのパス>、l<libを除いたライブラリのファイル名> でディレクトリパスとライブラリを指定します。

実行結果

筆者の実行環境で実行した場合の実行時間（5回実行した際の平均値）は

単純ループ：42.3 秒
組み込み関数 (matmul)：0.3696秒
BLAS (dgemm)：0.0994秒

となりました。BLASを利用するとかなり高速に計算ができることがわかります！

最後に

今回は、ライブラリの利用方法について紹介しました。数値計算では最終的には線形代数の基本的な計算に還元されることがほとんどです。このような計算が利用できる場合には、偉大な先人たちが作ったライブラリの力を借りて無駄な時間を費やさないようにしましょう。

また、今回使用したライブラリの実装は実行環境に特別にチューニングされたものではなく、あくまで手に入りやすいものを使用しただけです。より速さを追求するにはハードウェアに対してチューニングされたものを利用できる場合もありますので、ぜひ調べてみてください。

参考文献

LAPACK documentation：https://www.netlib.org/lapack/explore-html/index.html
幸野智紀、LAPACK/BLAS 入門、森北出版
OpenBLAS：http://www.openmathlib.org/OpenBLAS/