Fortranにおける「転置」の深淵：組込み関数とキャッシュ効率の冷徹な関係
数値計算の現場で、我々はしばしば「なぜこの程度の行列演算で、理論性能の10%も出ないのか？」という壁に突き当たる。そのボトルネックの正体は、多くの場合、CPUではなく「メモリ階層の無慈悲な構造」にある。
今回は、Fortranの組込み関数 `TRANSPOSE` が裏で何をしているのか、そして大規模シミュレーションにおいて我々が書くべき「魂のコード」とは何か、その極限の知見を共有しよう。
—
1. `TRANSPOSE` 関数の甘い罠
Fortranは列優先（Column-major）ストレージだ。メモリ上では `A(1,1), A(2,1), A(3,1)…` の順でデータが並んでいる。この並び順を維持したままループを回せば、CPUキャッシュは最大効率で機能する。
しかし、`TRANSPOSE(A)` を呼ぶと何が起きるか。コンパイラは内部的に一時的な作業領域（テンポラリ配列）を動的に確保し、そこへ要素をコピーして転置行列を構築する。
- メリット: コードが極めて簡潔になり、可読性が高い。
- デメリット: メモリ確保のオーバーヘッドと、コピー操作によるメモリ帯域の浪費。
数千×数千の行列でこれをループ内で叩けば、キャッシュミスが多発し、メモリコントローラが悲鳴を上げる。我々のような大規模シミュレーション屋が、クリティカルパスで安易に組込み関数を多用してはならない理由はここにある。
—
2. キャッシュを味方につける：ブロッキング（タイル化）の実装
メモリの物理配置を意識するなら、ループを「タイル」に分割し、キャッシュライン（通常64バイト）に収まる範囲でデータを処理する「ブロッキング」が必須だ。
以下に、コンパイラの自動ベクトル化を最大限に引き出す、堅牢な手動転置の実装例を示す。
subroutine optimized_transpose(n, A, B)
implicit none
integer, intent(in) :: n
real(8), intent(in) :: A(n, n)
real(8), intent(out) :: B(n, n)
! 現場の知見: ブロックサイズはCPUのL1/L2キャッシュサイズに合わせる
! 8×8や16×16のタイル化が、最近のIntel/AMDプロセッサでは安定して速い
integer, parameter :: block_size = 32
integer :: i, j, ii, jj
! 1. ループの入れ替えとタイリングにより、ストライドアクセスを抑制する
do jj = 1, n, block_size
do ii = 1, n, block_size
do j = jj, min(jj + block_size – 1, n)
! SIMDベクトル化を促すため、内側のループは連続メモリにアクセスさせる
do i = ii, min(ii + block_size – 1, n)
B(j, i) = A(i, j)
end do
end do
end do
end do
end subroutine optimized_transpose
このコードの肝は、`B(j, i) = A(i, j)` というアクセスにおいて、内側のループがメモリに対して可能な限り連続的なアクセスを維持するように組んでいる点だ。
—
3. コンパイラを「その気」にさせるビルド設定
せっかくコードを最適化しても、コンパイラが余計な安全チェックを入れたら台無しだ。特に大規模計算では、以下のフラグを使い分けるのが「プロ」の流儀である。
- Intel Fortran (ifort/ifx):
`-O3 -xHost -qopt-report=5 -qno-opt-dynamic-align`
- `-xHost`: そのCPUが持つAVX-512等の命令セットをフル活用する。
- `-qopt-report`: コンパイラがループをどう展開したか、ベクトル化に失敗した理由は何かが逐一出力される。これを見ないエンジニアはプロとは呼べない。
- GNU Fortran (gfortran):
`-O3 -march=native -ffast-math -funroll-loops`
- `-ffast-math`: 浮動小数点の厳密な結合法則を無視して良いなら必須。これだけで速度が数倍になることもある。
—
4. シニアアーキテクトからの助言：設計ルール
1. 「転置してから計算」は悪手: 転置した行列が必要なのではなく、転置したアクセス順序が必要なだけなら、アルゴリズム側を書き換えて `A(i, j)` を `A(j, i)` にアクセスできるようにループ順序を入れ替えろ。行列を物理的に転置する必要は殆どない。
2. 配列の形状（Shape）を意識せよ: 構造体や派生型の中に巨大な行列を詰め込むな。メモリレイアウトが崩れ、コンパイラのベクトル化が阻害される。
3. 計測なき最適化は単なる勘: 必ず `gprof` や Intel VTune Profiler を使い、実際にどこでキャッシュミスが起きているかを特定してから手を動かせ。
Fortranは、メモリを物理的に支配できる最後の高級言語だ。`TRANSPOSE` という便利な黒魔術に頼る前に、一度そのメモリマップを想像してほしい。貴殿の書くコードが、計算機の性能を極限まで引き出すことを期待している。
コメント