疎行列を「殺す」な：PACK/UNPACKによるメモリアクセスの極致
大規模シミュレーションの現場において、疎行列（Sparse Matrix）をまともに扱うことは、メモリ帯域という名の「ボトルネックとの戦い」を意味する。多くの若手エンジニアは、とりあえず `do` ループを回して非零要素を判定し、条件分岐で処理をスキップするようなコードを書く。だが、それは現代のCPUのパイプラインにとって「地獄」以外の何物でもない。
今回は、モダンFortranが提供する組み込み関数 `PACK` と `UNPACK` を活用し、メモリ効率を劇的に向上させつつ、コンパイラのベクトル化（SIMD）を誘発する「現場の戦術」を伝授する。
—
なぜ `do` ループによる判定が「悪」なのか
疎行列の非零要素を抽出する際、`if (A(i,j) /= 0.0) …` といった条件分岐をループ内に埋め込むと、CPUの分岐予測器は悲鳴を上げる。さらに、メモリの飛び石アクセスが頻発すれば、キャッシュラインの活用率は低下し、理論演算性能の数％も引き出せない。
これに対する回答は、「条件を満たす要素のみを連続したメモリ領域に一括で集約する」ことだ。ここで `PACK` 関数の出番となる。
PACK関数によるメモリの「圧縮」
`PACK` は、論理マスクに基づいて配列から要素を抽出し、1次元配列（ベクトル）に詰め直す。重要なのは、これが単なるデータ変換ではなく、「後続の計算におけるメモリアクセスの局所性を最大化する準備」だという点だ。
以下のコード例を見てほしい。単なる抽出ではなく、後続の演算でインデックスの追跡コストを減らすための「疎行列の圧縮形式（CSRに近い考え方）」への橋渡しである。
subroutine compress_sparse_data(matrix, mask, compressed_vec)
! 非零要素のみを抽出して連続領域に配置する
real(8), intent(in) :: matrix(:,:)
logical, intent(in) :: mask(:,:)
real(8), allocatable, intent(out) :: compressed_vec(:)
integer :: n_elements
! 1. マスクに基づいて非零要素の数をカウント
n_elements = count(mask)
allocate(compressed_vec(n_elements))
! 2. PACKによる高速な一括圧縮
! ここでメモリは連続した配列に配置され、後続のSIMD最適化が効きやすくなる
compressed_vec = pack(matrix, mask)
end subroutine compress_sparse_data
コンパイラ最適化を最大化する「設計の鉄則」
`PACK` を使用する際、単に「動けばいい」と考えてはいけない。以下の3点を意識することで、コンパイラ（Intel FortranやGNU gfortran）は驚くほど効率的なベクトル命令を吐き出すようになる。
1. 配列の列優先（Column-Major）を遵守せよ: Fortranの仕様である列優先を意識し、抽出元（`matrix`）のメモリ配置を、ループ内でインデックスが連続して変化するように設計する。
2. 一時配列の寿命管理: `PACK` の戻り値として `allocatable` を使用する場合、ヒープメモリの確保・解放がボトルネックになる可能性がある。シミュレーションのメインループ内であれば、あらかじめバッファを確保（`reallocate` を避ける）し、`intent(out)` で渡す設計にせよ。
3. コンパイラへのヒント: `-O3 -march=native -qopt-report` (ifortの場合) を付与し、最適化レポートで「Vectorized」の文字が出ているか確認すること。`PACK` が内部的にインライン展開され、ループが最適化されていることが確認できれば勝利だ。
UNPACKによる「復元」と物理モデルへの適用
計算が終われば、結果を元の格子（グリッド）へ戻さねばならない。`UNPACK` は `PACK` の逆演算であり、圧縮されたベクトルを元の行列形状へ展開する。
subroutine expand_sparse_data(compressed_vec, mask, original_shape, matrix)
real(8), intent(in) :: compressed_vec(:)
logical, intent(in) :: mask(:,:)
integer, intent(in) :: original_shape(:)
real(8), intent(out) :: matrix(original_shape(1), original_shape(2))
! 0.0で初期化した上で、マスクに従ってデータを展開
matrix = 0.0d0
matrix = unpack(compressed_vec, mask, field=matrix)
end subroutine expand_sparse_data
現場の知見：なぜこれが「セキュア」なのか
このアプローチが「堅牢」である最大の理由は、「手書きの `do` ループによる境界外参照（Out-of-bounds）を排除できる」ことにある。
境界判定をコンパイラの組み込み関数に一任することで、人間が犯しがちな「インデックスの+1/-1ミス」や「ループ回数の不整合」を言語仕様レベルで防ぐことができる。特に大規模数値計算において、デバッグに数日を溶かすようなバグの多くは、この手のループ制御で発生する。
まとめ：数理とハードウェアの調和
`PACK` / `UNPACK` を活用したデータ構造変換は、単なるコードの短縮ではない。「メモリ帯域を飽和させず、CPUのベクトルレジスタをフルに活用するための戦略的データレイアウト」である。
- 性能重視なら: `PACK` 後の連続メモリに対して、`!$omp simd` 指示文を添えてベクトル化を強制する。
- 保守性重視なら: 複雑な条件判定を `mask` 配列の生成ロジックにカプセル化し、メインの計算ロジックを極限までシンプルに保つ。
君たちが書くコードが、次の宇宙探査や次世代材料開発のシミュレーションにおいて、わずか数％の性能向上をもたらすかもしれない。その積み重ねこそが、最高峰の計算機科学を支えるエンジニアの矜持だ。現場からは以上だ。次回のコードレビューで、美しい `pack` の実装が見られることを期待している。
コメント