「なぜFortranの行列計算は速いのか？」MATMULとBLASの境界線で見極める極限の最適化
こんにちは。長年、スパコンの冷却ファンの音を聞きながら、数千コアをぶん回して流体解析をしてきた「数値計算の現場」から来ました。
PythonやC言語からFortranの世界へ足を踏み入れた皆さん、ようこそ。Fortranは単なる古い言語ではありません。「コンパイラが計算機の限界性能を引き出すために、言語側が余計な足枷を外している」という、極めて戦闘力の高いツールです。
今日は、行列演算の二大巨頭である「組込み関数 `MATMUL`」と「外部ライブラリ `BLAS`」の使い分けについて、教科書には載っていない「現場の勘所」をお話しします。
—
1. そもそも、なぜ行列演算が「Fortranの主戦場」なのか？
Fortranには「配列のメモリ配置が列優先（Column-major）」であるという、C言語ユーザーには少し奇妙なルールがあります。しかし、これが線形代数計算において、CPUのキャッシュラインを極限まで使い切るための「最強の武器」になるんです。
まずは、最もシンプルな `MATMUL` を見てみましょう。
! シンプルな行列積の例
real(8), dimension(1000, 1000) :: A, B, C
! C = A B を計算
C = matmul(A, B)
これだけで動きます。非常に簡単ですよね。しかし、プロの現場では「これだけで終わらせない」理由があるのです。
—
2. MATMULの正体：コンパイラの健気な努力
`MATMUL` は組込み関数であり、コンパイラ（Intel Fortranやgfortran）はこれを「よし、可能な限り最適化してやろう！」と頑張ってくれます。ループのアンロールやベクトル化を自動的に適用し、近年のコンパイラは非常に優秀です。
しかし、`MATMUL` には限界があります。それは「汎用性」です。
`MATMUL` はあらゆるサイズ、あらゆる形状の行列に対して、コードが自動生成されます。一方で、プロの現場で使う BLAS (Basic Linear Algebra Subprograms) は、CPUメーカー（IntelならMKL、AMDならAOCLなど）が、そのプロセッサのキャッシュサイズ、パイプラインの深さまで知り尽くした上で、アセンブラレベルでチューニングされた「魔境の関数」です。
—
3. 実践：MATMUL vs BLAS (DGEMM)
行列サイズが小さいうちは `MATMUL` でも十分ですが、行列が $1000 \times 1000$ を超えたあたりから、性能差は絶望的に開きます。BLASの行列積関数 `DGEMM` を呼ぶ準備をしましょう。
! BLASのDGEMMを呼び出すためのインターフェース定義
! 実際にはモジュール化して使うのが定石です
interface
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
character, intent(in) :: transa, transb
integer, intent(in) :: m, n, k, lda, ldb, ldc
real(8), intent(in) :: alpha, beta
real(8), dimension(lda,), intent(in) :: a
real(8), dimension(ldb,), intent(in) :: b
real(8), dimension(ldc,), intent(inout) :: c
end subroutine
end interface
! 呼び出し側
! C = 1.0 A B + 0.0 C
call dgemm(‘N’, ‘N’, 1000, 1000, 1000, 1.0d0, A, 1000, B, 1000, 0.0d0, C, 1000)
なぜ `DGEMM` が速いのか？
1. ブロッキング（タイリング）: 行列を小さなブロックに分割し、CPUのL1/L2キャッシュに収まるように計算します。`MATMUL` も頑張りますが、BLASの手法はハードウェアの特性を完全に掌握しています。
2. マルチスレッド: MKLなどのライブラリは、内部でOpenMPを使い、行列の分割をコア数に合わせて最適化します。
3. SIMD命令の極限活用: AVX-512など、CPUの拡張命令を叩くために最適化されたコードが走ります。
—
4. 現場の結論：どっちを使うべき？
若手の皆さんがコードを書く際の「判断基準」を授けます。
- MATMULを使うべきケース:
- 行列サイズが小さい（数百要素程度）。
- 外部ライブラリ（MKL等）の依存関係を作れない、環境構築を極力シンプルにしたいとき。
- プロトタイプ作成段階。
- BLAS (DGEMM) を使うべきケース:
- 行列サイズが1000を超える計算が頻出する。
- シミュレーションの実行時間が「数時間〜数日」のレベルである。
- スパコンや高性能ワークステーションで動かすのが前提。
—
最後に：まずは「動かして測る」ことから
「最適化」とは、勘でコードを書き換えることではありません。まずは `MATMUL` でコードを書き上げ、`time` コマンドやプロファイラ（Intel VTuneなど）でボトルネックを特定してください。「ここが遅い」と分かってから、初めて `DGEMM` に差し替える。このステップを怠らないことが、プロの数値計算エンジニアへの第一歩です。
Fortranは、書いたコードがそのまま計算機の性能に直結する、非常に正直で面白い言語です。まずは小さな行列で、コンパイラの最適化オプション（`-O3 -march=native` など）をいじりながら、実行時間の変化を観察することから始めてみてください。
皆さんの計算機が、今日も最高のパフォーマンスを叩き出しますように！応援しています。
コメント