Gauss-Seidelのループ依存を測定し、ループ展開で修正する

前回、Gauss-Seidel法を探求しました。その結果、数学的にはGauss-Seidelが収束に半分のイテレーションしか必要としないことがわかり、これは数値テストケースでも裏付けられました。しかし、イテレーション数が半分であるにもかかわらず、Gauss-SeidelソルバーはJacobiソルバーの4〜5倍の時間がかかります。これは考えてみると奇妙な事実です。紙の上ではより優れた方法が、なぜ数値実装では劣るパフォーマンスを示すのでしょうか？ご覧いただいたように、その説明は、Jacobiとは異なり、Gauss-Seidelの更新ルールにはループ依存があり、コンパイラがコードをベクトル化するのを妨げているという点にあります。今日、私たちが答えようとする質問は、これで行き止まりなのか、それともJacobiのハードウェア効率を犠牲にすることなく、Gauss-Seidelの収束の利点を何らかの方法で回復できるのか、ということです。そして、その答えは「できる」のです！ただし、2つの異なるルートがあります。1つは、コンパイラを最大限に支援する比較的汎用的な方法（コンパイラとCPUの仕組みを少し理解する必要があります）、もう1つは、特定の2Dポアソンテストケースに特化した、数学に導かれる方法です。どちらの道も平坦ではありませんが、その過程で多くのことを学ぶことができます。最終的には両方を検討しますが、今のところは最初のルートをたどり、ループ展開とは何かを発見します。待てよ、問題は何でしたっけ？少し巻き戻しましょう。前回、同じ2Dポアソン問題で2つの更新ルールを比較しました。素朴なJacobiイテレーションと、それより少し洗練されたそのいとこであるGauss-Seidelです。 ``` do j = 2, n-1 do i = 2, n-1 u(i, j) = 0.25_dp*(dx2*b(i,j) + v(i+1,j) + v(i-1,j) & + v(i,j+1) + v(i,j-1)) end do end do ``` Jacobiは常に古いイテレート `v` から読み込み、新しい配列 `u` に書き込みます。対照的に、Gauss-Seidelはインプレースで更新します。 ``` do j = 2, n-1 do i = 2, n-1 u(i, j) = 0.25_dp*(dx2*b(i,j) + u(i+1,j) + u(i-1,j) & + u(i,j+1) + u(i,j-1)) end do end do ``` これは1行の違いで、より高速に収束する方法を手に入れることができますが、指定された許容誤差に到達するために必要なスイープ数はちょうど半分です。紙の上では、Gauss-Seidelが圧倒的に勝つはずです。それなのに、実際に数値を実行してみると、Gauss-Seidelは必要なイテレーション数が半分であるにもかかわらず、壁時計時間でJacobiの4〜5倍の時間がかかりました。その1行の違いのどこかに、私たちに多大なコストをかけているものがあり、それは収束理論のどこにも現れていません。 ループ依存性の具体化 原因は、私たちが発見したように、Gauss-Seidelの更新における `u(i-1,j)` が、それを読み取る時点ですでに上書きされていることです。それは前のスイープからの古い値ではなく、このスイープで新しく計算された隣接値です。まさにこれがGauss-Seidelがより速く収束する理由です。1回のスイープ内で情報を伝播させているので、完全なイテレーションを待つ必要がありません。しかし、それはまた、イテレーション `i` がイテレーション `i-1` が `u(i-1,j)` への書き込みを完全に終了するまで開始できないことを意味します。Jacobiにはそのような制約はありません。Jacobiスイープでのすべての読み込みは `v` から行われ、そのスイープでは何も `v` を変更しません。これはループ依存性の教科書的な定義であり、コンパイラがループをベクトル化できるか、命令を並べ替えられるか、アウトオブオーダーコアでイテレーションをオーバーラップできるかを決定する前に推論する必要があるまさにその種のものです。Jacobi：そのような依存性なし、コンパイラの自由度は高い。Gauss-Seidel：ハードな依存性、コンパイラの自由度は低い。 「拘束されている」ことの数値化 それは良い話ですが、測定できない限りはまだ話に過ぎません。そこで今回は、コンパイラの言葉を鵜呑みにするのではなく、OSACA（Open Source Architecture Code Analyzer）を使用して、両方のカーネルのコンパイル済みアセンブリに直接向けました。OSACAはマークされたアセンブリループを読み込み、理想化されたアウトオブオーダーコアを仮定して、ループ本体ごとに2つの数値を報告します。クリティカルパス（CP）は、単一イテレーション内で実行される最長の依存チェーンの長さ、およびループ依存（LCD）は、次のイテレーションからクロスする必要があるそのチェーンの一部です。LCDがここで重要な数値です。これはイテレーションあたりのサイクルのハードフロアであり、巧妙なスケジューリングではそれを下回ることはできません。以下は、当社の `-O3 -march=native` カーネルを、プレーンな `do` ループ（`do concurrent` やその他の特別なことはなし）でコンパイルした結果です。 ``` Summary of the osaca analysis of the Jacobi and standard Gauss-Seidel kernels. Both kernels have been compiled using gfortran 15.3 with options -O3 -mtune=native -march=native. The precise numbers reported may depend on the exact CPU and compiler you’re using. CP (cycles) LCD (cycles) Port-pressure floor (cycles) Governing bound Jacobi (vectorized, 4 elements/it.) 24 1 3 throughput Gauss-Seidel (scalar, 1 element/it.) 24 12 2.5 latency ``` すぐに2つのことがわかります。第一に、クリティカルパスは両方のカーネルで同一で24サイクルです。これは理にかなっています。基本的に両方のケースで浮動小数点演算の同じチェーンです。3つの加算： `tmp1 = v(i+1, j) + v(i-1, j) + v(i, j+1) + v(i, j-1)` 1つの乗算： `tmp2 = b(i, j) * dx2` 1つのFMA： `tmp3 = 0.25_dp * tmp1 + 0.25_dp*tmp2` ロードとストアとともに。算術演算は変更されていません。オーバーラップが許可されているものだけが変更されています。第二に、そしてこれがすべてを決定します：Jacobiのループ依存は、わずか1サイクル（ループカウンターのインクリメントのみで、実際の計算とは無関係）です。Gauss-SeidelのLCDは12サイクルであり、OSACAはそれがどこから来るのかを正確に示しています：`46 | 12.0 | vaddsd %xmm1, %xmm0, %xmm0 | [46, 47, 48]` それをたどると、それは西隣の項です：%xmm1は、このスイープで新しく書き込まれた `u(i,j)` を保持しており、次のループトリップで `u(i-1,j)` として直接読み込まれます。3つの連鎖した命令（西隣の加算、北隣の加算、および最終的な乗算）は、それぞれ4サイクルのレイテンシを持ち、合計でちょうど12サイクルになります。これはコンパイラが怠けているのではなく、アルゴリズム自体が要求するため、`u(i-1,j)` が存在する前に `u(i,j)` の計算を開始する方法がないことをコンパイラが正しく報告しているのです。もう1つの手がかりは、OSACAを必要とせずにアセンブリから直接読み取れる命令のニーモニック自体です。Jacobiの内部ループは、パケット化された倍精度命令（`vaddpd`、`vmulpd`、`vfmadd213pd`）で完全に構築されており、それぞれが単一の256ビットAVXレジスタで一度に4つのグリッドポイントを処理します。Gauss-Seidelの内部ループは、まったく同じ命令のスカラー形式（`vaddsd`、`vmulsd`、`vfmadd213sd`）を使用しており、一度に1つのグリッドポイントを処理します。`pd`/`sd` の接尾辞だけで、コンパイラがベクトル化できたループとできなかったループが一目でわかります。 数値は合っていますか？ 各カーネルについて、CPUの実際の持続的なイテレーションあたりのコストは、スループットフロアとLCDの両方の制約によって下限が定められており、したがって実際に支配的になるのは、2つのうち大きい方です。Gauss-Seidelの場合、それはLCD（12サイクル、2.5サイクルのスループットフロアに対して）です。Jacobiの場合、それはスループットフロア（3サイクル、無視できる1サイクルのLCDに対して）です。言い換えれば、Gauss-Seidelはレイテンシバウンドであり、Jacobiはスループットバウンドです。まったく異なる2つのボトルネックであり、比較が単に両方の同じ列を読む問題ではない理由です。各カーネルの支配的なバウンドをグリッドポイントあたりで正規化する（Jacobiの値を4で割る、なぜならベクトルイテレーションあたり4つの要素を処理するため）と、スイープあたりの予測される減速は約 ``` (12 / (3 / 4)) = 16倍 ``` となります。これは、壁時計で測定した4〜5倍よりも高い値です。しかし、壁時計は解法までの総時間を測定し、Gauss-Seidelは収束するためにJacobiの半分のスイープしか必要としないことを覚えておいてください。それを考慮に入れると、測定値から示唆される実際のスイープあたりの減速は ``` (T_GS / T_J) / (N_GS / N_J) ≈ (4–5) × 2 ≈ 8–10倍 ``` となり、OSACAの予測値の2倍以内の範囲に収まります。これは、メモリシステムへの可視性がゼロの静的モデルとしては、健全なレベルの一致です。OSACAが特定するメカニズムは、Gauss-Seidelがハードな12サイクルの壁でレイテンシバウンドになり、Jacobiがc