演算器の暴力とメモリの沈黙：MATMUL/DOT_PRODUCTの「その先」にある真実
スパコンの演算性能がエクサスケールに到達しようとする今、我々が直面しているのは「演算器の飽和」ではなく「メモリ帯域の枯渇」である。
君たちが書いている `DO` ループは、コンパイラが自動ベクトル化してくれると信じているかもしれない。だが、その期待は多くの場合、裏切られる。メモリのレイアウト、ストライドの不一致、そして何よりコンパイラが「安全側に倒した」エイリアス解析の結果によって、君たちのコードは本来の性能の10%も出せていないはずだ。
今日は、Fortranの組込み関数 `MATMUL` と `DOT_PRODUCT` を巡る、現場の泥臭い最適化について語ろう。
—
1. なぜ「手書き」より組込み関数なのか
「手書きの三重ループの方が速い」と信じている諸君へ。それは、君がBLAS（Basic Linear Algebra Subprograms）の内部実装を完全に再現できるほど暇な場合か、あるいはコンパイラがインライン展開を完璧に行える単純なケースに限られる。
モダンFortranの `MATMUL` や `DOT_PRODUCT` は、単なる関数ではない。多くの商用コンパイラ（Intel oneAPI, Cray, Fujitsu）において、これらはリンク時にBLASライブラリへの呼び出しへ動的に置換されるか、あるいはコンパイラ最適化パスの深い領域で、SIMD命令（AVX-512やSVE）を駆使したカーネルへ書き換えられる。
重要なのは「何を使うか」ではなく「どうメモリに配置するか」だ。
Fortranの列優先順位（Column-Major）を理解していないコードは、現代のキャッシュハイアラキーにおいて「死」を意味する。
—
2. キャッシュミスを制する：メモリアクセスの鉄則
行列演算において最も避けるべきは、L1キャッシュへのロードを阻害するストライドアクセスだ。
! 非推奨：列優先のFortranで、行方向にアクセスするのはキャッシュ汚染の元
do j = 1, n
do i = 1, m
a(i, j) = b(i, j) c(i, j)
end do
end do
! ↑ これは正解。インデックスの順序とメモリ配置が一致している。
! 最悪のケース
do i = 1, m
do j = 1, n
a(i, j) = … ! ここでストライドnのアクセスが発生し、キャッシュラインを無駄に消費する
end do
end do
`MATMUL` を使う際も、引数として渡す行列のメモリレイアウトが `contiguous`（連続）であることを保証せねばならない。`POINTER` 配列や `ALLOCATABLE` の不適切なスライス渡しは、コンパイラに「一時的なコピー（Temporary Array）」を作成させ、計算以前にメモリ帯域を食い尽くす。
—
3. 最適化の極致：コンパイラフラグと並列実行の調律
数万コア規模のHPC環境で、`gfortran -O3` で満足している者はいないだろう。実戦では以下のフラグが必須となる。
Intel oneAPI (ifort/ifx) の例
-xHost: 実行環境のCPUアーキテクチャで可能な限りSIMDを叩く
-qopt-report=5: コンパイラがどこでベクトル化を諦めたかを確認する
-mkl: 組込み関数をMKL(BLAS/LAPACK)へ直結させる
ifx -O3 -xHost -qopt-report=5 -mkl=parallel -qopenmp source.f90 -o simulation.x
ここで重要なのは、`qopt-report` を見て「Vectorization diagnostics」を読み解くことだ。もしループがベクトル化されていないなら、そこには「データ依存性」か「不規則なメモリ参照」が存在する。
—
4. ハイブリッド並列の罠：MPI + OpenMP
数万コアを動かす場合、純粋なMPIよりもハイブリッド並列が好まれる。しかし、`MATMUL` をOpenMPの並列領域内で呼び出すときは注意が必要だ。
! 並列領域内でのMATMUL呼び出し
!$OMP PARALLEL DO
do i = 1, n_blocks
call compute_block(my_matrix(:,:,i))
! 内部でMATMULが呼ばれる際、MKLがスレッドをさらに生成するとオーバーサブスクリプションが発生する
end do
!$OMP END PARALLEL DO
現場の知見:
MKLなどのライブラリを使う場合、`MKL_NUM_THREADS=1` に設定し、OpenMPの並列数はMPIプロセス側のノード内スレッド数に委ねるべきだ。ライブラリとユーザーコードでスレッドを競合させるのは、パフォーマンス低下の最大の要因となる。
—
5. 結論：計測なき最適化は勘違いの温床
プロファイラ（Intel VTuneやScalasca）を使わずに「速くなったはずだ」と断言するのは、目隠しをしてF1カーを運転するようなものだ。
1. ベースラインの測定: まずは素の `MATMUL` でボトルネックを特定する。
2. メモリの連続性確保: `CONTIGUOUS` 属性を付与し、配列の形状を最適化する。
3. カーネルの強制: 必要であれば `dgemm` (BLAS) を直接叩き、行列積の計算順序（Tiling）を調整する。
Fortranは死んでいない。物理現象を記述する言語として、これほど計算機資源の深淵にアクセスできる言語は他にない。コンパイラの挙動を理解し、ハードウェアの悲鳴を聞きながらコードを磨き上げる。それが、我々数値計算エンジニアの矜持だ。
次回の記事では、`ISO_C_BINDING` を用いて、FortranとC/C++のメモリ共有を行う際の「型安全性」と「ABIの衝突」について、さらに深掘りしていく。期待していてくれ。
コメント