組み込み `MATMUL` か、それともBLASか？ ― 現場のシニアが教える計算機資源の「極限利用」術
数値計算の最前線でコードを叩いている諸君。君たちが書いた `C = MATMUL(A, B)` という一行。コンパイラが自動で高速化してくれると思ったら大間違いだ。
多くのジュニアエンジニアが陥る罠は、MATMULとBLASの境界を「単なる好み」で選んでしまうことにある。しかし、我々のような大規模シミュレーションを回す人間にとって、この選択は計算時間数週間分の差を分かつ致命的な分岐点となる。今回は、この「見えないコスト」の正体と、現場で生き残るための最適化戦略を伝授しよう。
—
1. `MATMUL` の正体と限界：コンパイラの限界値
`MATMUL` はFortranの組み込み関数であり、極めてポータブルだ。しかし、注意が必要だ。
- 自動ベクトル化の皮肉: コンパイラは `MATMUL` を見ると、ループ展開やベクトル化を試みる。しかし、これは「汎用的なループ」に対する最適化に過ぎない。
- キャッシュ効率の欠如: `MATMUL` は、行列のサイズが小さい（例えば100×100以下）場合には、関数呼び出しのオーバーヘッドが少なく非常に速い。しかし、キャッシュラインを意識した「ブロッキング（タイリング）」処理においては、手書きのループやBLASには到底及ばない。
つまり、`MATMUL` は「小規模な行列を頻繁に計算する」場合には最適だが、大規模な科学技術計算には不向きなのだ。
2. BLAS（MKL/OpenBLAS）が圧倒的に速い理由
なぜベンダー提供のBLAS（`dgemm`など）は速いのか？ それは、彼らが「CPUのアーキテクチャそのもの」を熟知しているからだ。
- キャッシュ・ブロッキング: 行列をCPUキャッシュサイズに合わせて細かく切り分け（タイリング）、L1/L2/L3キャッシュのヒット率を極限まで高めている。
- NUMA対応: マルチソケット環境において、メモリがどのノードにあるかを考慮した並列化を行っている。
- SIMD命令の極限活用: AVX-512等の最新命令セットを、アセンブリレベルでチューニングして叩いている。
これらを自前で実装するのは愚策だ。我々は「車輪の再発明」をせず、これらのライブラリを「正しく呼ぶ」ことに専念すべきである。
—
3. 実践：セキュアかつ高速な実装パターン
現場で私が愛用している、BLASを叩く際のラッパーパターンを紹介しよう。Fortranの `ISO_C_BINDING` を活用し、インターフェースを明確に定義することで、型不一致によるセグメンテーションフォールトを未然に防ぐ。
module mat_math_module
use, intrinsic :: iso_fortran_env, only: real64
implicit none
private
public :: fast_matmul
! BLASのdgemmをインターフェース宣言（型安全を担保）
interface
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc) &
bind(c, name=”dgemm_”)
import :: real64
character(kind=c_char), intent(in) :: transa, transb
integer(c_int), intent(in) :: m, n, k, lda, ldb, ldc
real(real64), intent(in) :: alpha, beta
real(real64), intent(in) :: a(lda,), b(ldb,)
real(real64), intent(out) :: c(ldc,)
end subroutine dgemm
end interface
contains
! 安全なラッパー関数
subroutine fast_matmul(A, B, C)
real(real64), intent(in) :: A(:,:), B(:,:)
real(real64), intent(out) :: C(:,:)
! BLAS呼び出しのための引数
integer :: m, n, k
m = size(A, 1)
k = size(A, 2)
n = size(B, 2)
! DGEMM: C = alphaAB + betaC
! ここでメモリレイアウト（列優先）を意識した引数を渡す
call dgemm(‘N’, ‘N’, m, n, k, 1.0_real64, A, m, B, k, 0.0_real64, C, m)
end subroutine fast_matmul
end module mat_math_module
なぜこの実装が良いのか？
1. インターフェース定義: `bind(c)` を使うことで、リンクエラーや引数の型不一致をコンパイル時に検知できる。
2. LDA/LDBの明示: 行列の先頭アドレスと実際のメモリ上の幅（Leading Dimension）を正しく渡すことで、サブ行列の計算にも対応できる設計にしている。
3. 可読性: 呼び出し側は `fast_matmul` を呼ぶだけでよく、複雑なBLASの引数に悩まされることがない。
—
4. シニアからの鉄則：ビルド設定と最適化フラグ
コードがどれほど洗練されていても、コンパイラの設定が甘ければ全て台無しだ。以下はIntel Fortran (ifort/ifx) を想定した現場の標準的なコンパイル設定だ。
最適化の基本
FC=ifx
FFLAGS=”-O3 -march=native -qopenmp -ipo”
LIBS=”-lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -liomp5 -lpthread”
-O3: 最強の最適化
-march=native: 実行環境のCPU命令をフル活用
-ipo: プロシージャ間最適化（モジュールをまたいだインライン化）
結び：エンジニアとして持つべき視点
`MATMUL` は「プロトタイピングと小規模行列」のために使い、`BLAS` は「プロダクション環境の数値シミュレーション」のために使え。
計算機科学の真髄は、「ハードウェアがどうデータを読み込み、どう演算器に流し込んでいるか」を想像することにある。メモリの列優先順位（Column-major order）を意識し、キャッシュラインを汚さないアクセスを心がければ、君たちのシミュレーションは確実に加速する。
次は、このBLAS呼び出しを「OpenMP」でどうラップし、複数ノード間（MPI）へとスケーリングさせていくかについて話をしよう。コードの先にある物理現象が見えてきたとき、それが本物のエンジニアリングだ。
コメント