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

Make Gauss-Seidel great again!

---

［画像: Gauss-Seidel法の収束実験結果を示すグラフ］

---

前回は、Gauss-Seidel法を探求しました。その結果、数学的にはGauss-Seidel法は収束に必要な反復回数が半分で済むという結論に至り、数値テストケースでも裏付けられました。しかし、反復回数が半分で済むにもかかわらず、Gauss-SeidelソルバーはJacobi法よりも4〜5倍時間がかかりました。これは、紙の上ではより優れた手法が、数値実装では劣るという、一見不可解な事実です。その理由は、Jacobi法とは異なり、Gauss-Seidel法の更新規則にはループ依存性があり、コンパイラがコードをベクトル化するのを妨げていることにあります。したがって、本日回答しようとする問いは、「これは終着点なのか、それともJacobi法のハードウェア効率を犠牲にすることなく、Gauss-Seidel法の収束性の利点を何らかの方法で回復できるのか？」ということです。

そして、その答えは「はい、できます！」です。ただし、2つの異なるルートがあります。1つは、コンパイラを最大限に支援する（ただし、コンパイラとCPUの仕組みを少し理解する必要があります）比較的汎用的な方法、もう1つは、2次元ポアソンテストケースに非常に特化した、数学に導かれる方法です。どちらの道も平坦ではありませんが、その過程で多くのことを学ぶでしょう。最終的には両方を探求しますが、まずは最初のルートをたどり、「ループ展開」とは何かを発見します。

待て、問題は何だったか？

一瞬巻き戻しましょう。前回、同じ2次元ポアソン問題に対して、素朴なJacobi反復と、それよりわずかに洗練されたそのいとこであるGauss-Seidel法の2つの更新規則を比較しました。

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回の完全な反復を待つ代わりに、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なし、何も特別なことはなし）でコンパイルした結果です。

Jacobiおよび標準Gauss-SeidelカーネルのOSACA解析の概要。両方のカーネルはgfortran 15.3で -O3 -mtune=native -march=native オプションを使用してコンパイルされました。報告される正確な数値は、使用するCPUとコンパイラによって異なる場合があります。

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)の計算を開始する方法はない、とコンパイラが正しく報告しているのです。

アセンブリコードをOSACAなしで見てもわかるもう一つの手がかりは、命令のニーモニック自体です。Jacobiの内部ループは、すべてパケット化された倍精度命令 — vaddpd、vmulpd、vfmadd213pd — で構成されており、それぞれが単一の256ビットAVXレジスタで一度に4つのグリッドポイントを処理します。Gauss-Seidelの内部ループは、全く同じ命令のスカラー形式 — vaddsd、vmulsd、vfmadd213sd — を使用しており、一度に1つのグリッドポイントを処理します。pd/sdの接尾辞を見るだけで、コンパイラがベクトル化できたループとできなかったループが一目でわかります。

数値は一致するか？

各カーネルについて、CPUの実際の持続的な1反復あたりのコストは、スループットフロアとLCDの両方の制約によって下限が定められているため、実際に支配的になるのは、2つのうち大きい方です。Gauss-Seidelの場合、それはLCD（12サイクル、対するスループットフロアは2.5サイクル）です。Jacobiの場合、それはスループットフロア（3サイクル、対する無視できるほどの1サイクルのLCD）です。つまり、Gauss-Seidelはレイテンシバウンドであり、Jacobiはスループットバウンドです。まったく異なる2つのボトルネックであり、比較が単に両方の同じ列を読むだけの問題ではない理由です。

各カーネルの支配的な制約をグリッドポイントあたりに正規化する（Jacobiの場合は4で割る、なぜならベクトル反復あたり4要素を処理するため）と、スイープあたりの予測遅延はほぼ次のようになります。

\[\dfrac{12}{3 / 4} = 16\times\]

これは、壁時計で測定された4〜5倍よりも高い値です。しかし、壁時計は解法までの総時間を測定し、Gauss-Seidelは収束するためにJacobiの半分のスイープしか必要としないことを思い出してください。それを差し引くと、測定値から示唆される実際のスイープあたりの遅延は次のようになります。

\[\left(\frac{T extrm{_GS}}{T extrm{_J}}\right) \Big/ \left(\frac{N extrm{_GS}}{N extrm{_J}}\right) \approx (4\text{–}5) \times 2 \approx 8\text{–}10\times\]

これは、静的なモデルとしては良好な一致の範囲内であり、メモリシステムを全く可視化できません。OSACAが特定したメカニズム — Gauss-Seidelはレイテンシバウンドでハードな12サイクル壁、Jacobiは快適にスループットバウンド — は正確ですが、正確な数値には少し塩が必要です。

ある程度の問題サイズでは、Jacobiの実際のボトルネックは、OSACAが報告する3サイクルの計算フロアよりも、DRAM帯域幅である可能性が非常に高く、これにより真のスイープあたりのコストが膨らみ、残りのギャップが埋まるでしょう。

したがって、理論、静的解析、測定の3つすべてが一致して、診断が確認されました。Gauss-Seidelのより速い収束は、コンパイラが単独では壊せない、シリアル化された12サイクル/要素の依存性チェーンの代償として完全に購入されています。これにより、明白な疑問が生じます：その依存性チェーンはアルゴリズムに根本的なものなのか、それともループの書き方のアクシデントなのか？

では、ループ展開とは何か？

ループ展開は、最適化コンパイラのプレイブックの中で最も古いトリックの1つであり、そのアイデアは驚くほど単純です。ループ本体を反復ごとに1回実行し、毎回ブランチとカウンタ更新のオーバーヘッドを支払う代わりに、本体を2回、4回、または8回書き出し、カウンタを一度にその数だけ進めます。

簡単な例を見てみましょう。まだポアソン問題とは関係ありません。

do i = 1, n
   a(i) = 2.0_dp * a(i)
end do

これを2倍に展開すると、次のようになります。

do i = 1, n, 2
   a(i)   = 2.0_dp * a(i)
   a(i+1) = 2.0_dp * a(i+1)
end do

nが奇数の場合の簡単なブックキーピングはありますが、ここでは便宜上無視します。ブランチが少なく、カウンタのインクリメントが少なく、そして—最も重要なこととして—単一のループトリップ内でコンパイラに公開される独立した作業量が2倍になります。最後の点は本当に魅力的です。隣り合って配置された2つのステートメント間に依存性がないため、コンパイラはそれらをインターリーブしたり、単一のベクトル命令にパックしたり、ハードウェアが好むようにスケジュールしたりできます。展開はCPUの作業量を減らすわけではありません。一度に見るべきものを増やすことで、スケジューラの仕事を楽にします。したがって、当然、尋ねるべき質問は次のようになります：同じトリックはGauss-Seidelを救うのか？

うまくいかない誘惑的な近道

明白なことを試して、Gauss-Seidelの内部ループを2倍に展開し、他の何も変更しないようにしましょう。

do i = 2, n-2, 2
   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))
   u(i+1, j) = 0.25_dp*(dx2*b(i+1,j) + u(i,j)     + u(i+2,j)     &
                                     + u(i+1,j-1) + u(i+1,j+1))
end do

2番目のステートメントをよく見てください。それは、最初のステートメントがちょうど2行上で計算を終えたばかりの値、u(i,j)を読み取ります。ループ依存性をまったく取り除いていません。単にそのコピーをループ本体内に移動させただけで、それはクロス反復依存性ではなく、同じ反復内依存性になりました。コンパイラは依然としてこれらの2つのステートメントを独立に扱うことができず、一緒にベクトル化することもできず、厳密な順序で実行しなければなりません。展開のブックキーピング（より大きなコード、奇数nの残りを処理する必要性）のコストを支払い、何も得られませんでした。

これは少し考える価値があります。なぜなら、これは本当に一般的な落とし穴だからです。展開は依存性を解消するのに役立つように見えますが、ベクトル化を可能にするためにしばしば使われるため、真の依存性チェーンを壊すためには単独では全く何も効果がありません。既存の並列性を公開するには非常に優れています（私たちの些細なa(i) = 2*a(i)の例のように）が、作成しなければならない並列性に対しては完全に無力です。Gauss-Seidelはこの2番目のカテゴリにきれいに収まります。依存性はコンパイラが見逃したスケジューリングの事故ではなく、更新規則の数学に焼き付けられています。

したがって、展開だけでは行き止まりの場合、何が足りないのでしょうか？答えは、展開が効果を発揮する前に、依存性自体を書き換えて解消する必要があるということです。それは、単なるループ構造の変更ではなく、代数に手を伸ばすことを意味します。

わかりました、Gauss-Seidelのためにそれをやりましょう！

代数を正直かつ読みやすくするために、今は完全な2次元グリッドのことは忘れて、1次元ポアソン方程式、(u'' = -f)を均一なグリッドで離散化したものとして考えます。

\[\frac{u extrm{_i+1} - 2 u extrm{_i} + u extrm{_i-1}}{\Delta x^2} = -f extrm{_i}.\]

(u extrm{_i})について解くと、更新規則が得られます。

\[u extrm{i} = \dfrac{1}{2}\left(u extrm{{i-1}} + u extrm{_{i+1}} + \Delta x^2 f extrm{_i}\right),\]

したがって、驚くことではありませんが、以前の(0.25)の役割を果たしていた定数は(c = \tfrac12)です。左から右へスイープし、インプレースで更新すると、Gauss-Seidelは次のように読み取ります。

do i = 2, n-1
   u(i) = c*(u(i-1) + u(i+1) + dx2*b(i))
end do

ここでの1つの隣接要素だけが問題です：u(i-1)は前のループトリップで上書きされましたが、u(i+1)はまだ変更されていない古い値です（スイープはまだそこに到達していません）。それが依存性の全体であり、他のものを邪魔するものなしに完全に露呈しています。

再帰として記述する

スイープ自体の進行に依存しないすべてを、1つの列あたりの量にまとめます。

\[t extrm{i} = u extrm{{i+1}} + \Delta x^2 f extrm{_i},\]

これはブロックが開始する前に完全にわかっています（古いデータであり、現在計算しているものには何も触れられていません）。更新は、クリーンな1項のアフィン再帰に収束します。

\[u extrm{i} = c\left(u extrm{{i-1}} + t extrm{i}\right), \qquad u extrm{{i-1}} \equiv \texttt{um},\]

ここで、umは、このブロックに供給される最後に更新された値を短縮したものです。これは、依存性チェーンを(u)の単なる再帰として可視化するまさにその形状です — それ以上でもそれ以下でもありません。

2ステップ先まで再帰を解く

(u extrm{i})を計算し、それから (u extrm{{i+1}})を計算するのではなく、ブロックが開始する前に存在した量のみを基準に、両方について再帰を明示的に解きます。

\[u extrm{_i} = c,(\texttt{um} + t extrm{_i})\]

\[u extrm{_{i+1}} = c,(u extrm{i} + t extrm{{i+1}}) = c\Big(c,(\texttt{um}+t extrm{i}) + t extrm{{i+1}}\Big) = c^2,\texttt{um} + c^2,t extrm{i} + c,t extrm{{i+1}}\]

何が起こったか見てみましょう：(u extrm{_{i+1}})はもはや(u extrm{_i})を参照していません。両方の出力は、um、(t extrm{i})、および(t extrm{{i+1}})の純粋な関数になりました。このペアの更新が始まる前に存在したすべての量です。プレーンな展開では触れることができなかったペア内の依存性（前のセクションの「展開のみ」カーネル）は、代数的に解消されました。

コードでは、(c extrm{_2} \equiv c^2)と定義すると、次のカーネルが得られます。

do i = 2, n-2, 2
   um   = u(i-1)
   t0   = u(i+1) + dx2*b(i)
   t1   = u(i+2) + dx2*b(i+1)
   u(i)   = c *(um + t0)
   u(i+1) = c2*um + c2*t0 + c*t1 ! u(i) に依存しない！
end do

2つの独立した式、一方の出力を他方の計算で読み取らないもの。コンパイラがシリアルではなく並行してスケジュールするために必要なものです。

2Dに戻る：新しいものは何も必要ない

ここで1Dの迂回が役立ちます。実際の2Dステンシルでは、北、南、ソースの寄与は、i依存性を生成する上で全く役割を果たしません。これらは、Gauss-Seidelがすでに終了した行の値か、まだ触れていない行の値ですが、iループ自体が計算しているものではありません。したがって、それらは定義した同じ集約された(t extrm{_i})に折りたたまれます。

\[t extrm{_i} ;\longrightarrow; \texttt{tmp} extrm{_i} = \Delta x^2 b(i,j) + u(i+1,j) + u(i,j+1) + u(i,j-1),\]

そして、まったく同じ2行の置換がそのまま適用され、展開されたカーネルが得られます。

real(dp), parameter :: c  = 0.25_dp
real(dp), parameter :: c2 = 0.0625_dp
do j = 2, n-1
   do i = 2, n-2, 2
      um   = u(i-1, j)
      tmp1 = dx2*b(i, j)   + (u(i+1, j) + u(i, j+1) + u(i, j-1))
      tmp2 = dx2*b(i+1, j) + (u(i+2, j) + u(i+1, j+1) + u(i+1, j-1))
      u(i, j)   = c *(um + tmp1)
      u(i+1, j) = c2*um + (c2*tmp1 + c*tmp2)
   end do

   ! グリッドポイント数が奇数の場合を処理します。
   if (mod(n, 2) == 1) then
      tmp1 = dx2*b(n-1, j) + (u(n, j) + u(n-1, j+1) + u(n-1, j-1))
      u(n-1, j) = c*(u(n-2, j) + tmp1)
   end if
end do

末尾のifブロックは、nが奇数の場合、つまり展開されたループが必要とする残りのブックキーピングと同じものを処理します。概念的に新しいものはありません。

実際に機能するか？

さて、新しいカーネルは元のカーネルよりも少し複雑に見えます。しかし、数学的にはすべて同じものですが、本当に効率が良いのでしょうか？理論を機械で確認する時が来ました。以前と同じワークフロー：-O3 -march=native -Sでコンパイルし、内部ループをマークし、OSACAを実行します。

                                     ポート圧力 (サイクル)                                      
     |  0   - 0DV  |  1   |  2   -  2D  |  3   -  3D  |  4   |  5   |  6   |  7   ||  CP  | LCD  |
--------------------------------------------------------------------------------------------------
  74 |             |      |             |             |      |      |      |      ||      |      |   .L5:
  75 |             |      | 0.50   0.50 | 0.50   0.50 |      |      |      |      ||      |      |   vmovsd 16(%rax), %xmm1
  76 |             |      | 0.50   0.50 | 0.50   0.50 |      |      |      |      ||      |      |   vmovsd 8(%rax), %xmm0
  77 | 0.00        | 0.00 |             |             |      | 0.50 | 0.50 |      ||      |      |   addq $16, %rdx
  78 | 0.50        | 0.50 | 0.50   0.50 | 0.50   0.50 |      |      |      |      ||  4.0 |      |   vaddsd 8(%rax,%rsi,8), %xmm1, %xmm1
  79 | 0.50        | 0.50 | 0.50   0.50 | 0.50   0.50 |      |      |      |      ||      |      |   vaddsd (%rax,%rsi,8), %xmm0, %xmm2
  80 | 0.50        | 0.50 | 0.50   0.50 | 0.50   0.50 |      |      |      |      ||  4.0 |      |   vaddsd 8(%rax,%rcx,8), %xmm1, %xmm1
  81 | 0.50        | 0.50 | 0.50   0.50 | 0.50   0.50 |      |      |      |      ||      |      |   vaddsd (%rax,%rcx,8), %xmm2, %xmm2
  82 | 0.00        | 0.00 |             |             |      | 0.50 | 0.50 |      ||      |      |   addq $16, %rax
  83 | 0.50        | 0.50 | 0.50   0.50 | 0.50   0.50 |      |      |      |      ||  4.0 |      |   vfmadd231sd -8(%rdx), %xmm5, %xmm1
  84 | 0.50        | 0.50 | 0.50   0.50 | 0.50   0.50 |      |      |      |      ||      |      |   vfmadd231sd -16(%rdx), %xmm5, %xmm2
  85 | 0.50        | 0.50 | 0.50   0.50 | 0.50   0.50 |      |      |      |      ||  4.0 |      |   vmulsd %xmm4, %xmm1, %xmm1
  86 | 0.50        | 0.50 |             |             |      |      |      |      ||      |      |   vaddsd %xmm3, %xmm2, %xmm0
  87 | 0.50        | 0.50 |             |             |      |      |      |      ||      |      |   vmulsd %xmm4, %xmm0, %xmm0
  88 | 0.50        | 0.50 |             |             |      |      |      |      ||  4.0 |      |   vfmadd132sd %xmm6, %xmm1, %xmm2
  89 | 0.50        | 0.50 |             |             |      |      |      |      ||  4.0 |  4.0 |   vfmadd132sd %xmm6, %xmm2, %xmm3
  90 |             |      |             |             |      | 1.00 |      |      ||  1.0 |      |   vunpcklpd %xmm3, %xmm0, %xmm0
  91 |             |      | 0.00        | 0.00        | 1.00 |      |      | 1.00 ||  0.0 |      |   vmovupd %xmm0, -16(%rax)
  92 | 0.00        | 0.00 |             |             |      | 0.00 | 1.00 |      ||      |      |   cmpq %rax, %rdi
  93 |             |      |             |             |      |      |      |      ||      |      | * jne .L5

       5.50          5.50   4.00   4.00   4.00   4.00   1.00   2.00   2.00   1.00    29.0    4.0

最後の2つの列を見ると、2つのことが際立つはずです。クリティカルパスCPはわずかに長くなりましたが、ループ依存性（LCD）はループトリップあたり12サイクルから4サイクルに減少しました — そして各トリップは現在2つの要素をカバーしているため、要素あたり約6倍の減少です。

これが静的な予測です。実際の機械が同意するかどうかを確認する時です — そして、より興味深いことに、なぜ同意するのかを、コンパイルフラグを1つずつ見ていきます。

実験1 — ベースライン、最適化なし

まず、整合性チェックです。-O0でコンパイルしたため、スケジューリングの工夫は一切行われません。

次に進む前に確認する価値のあることが2つあります。第一に、3つのソルバーすべてが、以前の投稿と同じ反復回数で解に到達します — 展開されたカーネルが、代数を含めて、テキストブックのものとまったく同じものを計算していることを確認します。第二に、-O0では、展開されたカーネルはすでにテキストブックカーネルよりも信頼性高く高速であり、Gauss-SeidelがJacobiよりも理論的に予測される比率（2倍）にほぼ正確に一致しています — コンパイラスケジューリングが活用されなくても、ループ本体ごとに2つの独立した計算を公開することが、ループオーバーヘッドを半分にすることから純粋に何かを買うという nice confirmation です。しかし、誰も-O0コードを出荷しないので、これはベースラインであり、結果ではありません。

実験2 — クリティカルな実験

次に、実際の最適化でコンパイルしますが、ベクトル化を意図的にオフにします（-O3 -march=native -mtune=native -fno-tree-vectorize）。これにより、見えるすべての利得がSIMDに起因しないようになります。

これが重要な数値です。ベクトル化が明示的に無効になっている場合、展開されたカーネルは、テキストブックカーネルよりもスイープあたり6倍高速であり、Jacobiのスイープあたりのコストに5%以内です。このフラグではどちらのGauss-Seidelカーネルもベクトル化できないため、この利得は、グリッドポイントを1つの命令にパックすることとは関係ありません。それは、より微妙なものに由来します：ループ依存性の減少により、CPUのアウトオブオーダースケジューラは、テキストブックカーネルの12サイクル依存性チェーンがアイドル状態にしていた実行ポートを埋める、スカラー作業を連続する反復でオーバーラップさせることができます。待機に費やすサイクルが少なく、発行される命令数は同じです。これは、OSACAテーブルのLCD対スループットの話の経験的な確認であり、ベクトル化から完全に分離されています。

実験3 — ベクトル化を再度オンにする

最後に、同じ3つのカーネル、同じフラグですが、ベクトル化を許可した場合です。

両方のGauss-Seidelカーネルは、マイクロ秒単位で完全に変更されていません。これは、アセンブリを知っていれば予想されることです。gfortranはこのフラグに関係なく、どちらのGSカーネルも決してベクトル化しないため、フラグが作用するものは何もありません。ベクトル化は、実験2のスピードアップの背後にあるメカニズムではありませんでした。それは別の話であり、Jacobiにのみ適用されます。

Jacobiの場合でさえ、利得は現実ですが控えめです：166 µs → 115 µs、約1.4倍であり、4つの倍精度数を単一のAVX命令にパックすることから期待される約4倍ではありません。これは、LCDの話が予測したが、それ自体では完全に説明できなかった副産物です。ベクトル化は、カーネルの計算がボトルネックである場合にのみ効果があり、Jacobiの内部ループは、ループ依存性が全くないにもかかわらず、この問題サイズではほとんどメモリバウンドです。これは、命令レベルの分析だけでは示唆されるよりもはるかに低い達成可能な速度アップを、DRAM帯域幅が制限しています。

3つの実験を並べて見ると、結論は曖昧ではありません：ベクトル化ではなく、ループ依存性の欠如（または削減）がスケジューラにその仕事を行わせます。 ベクトル化は、利用可能な場合に良いボーナスです。Jacobiにはまだ役立ちますが、単純に推測するほどではありません。しかし、それはGauss-Seidelが遅かった理由ではなく、LCDの削減はすでに、単独で、Jacobiとのパフォーマンスギャップのほぼ全体を埋めるのに十分でした。実際、この展開されたGauss-Seidelカーネルは、Jacobiのものよりも3秒近く速くソリューションを提供します。ミッション完了！あるいは、そう見える...

しかし、ニーモニックをよく見てください

ここに、CP/LCDの数値だけを見て勝利を宣言した場合に見落としがちな詳細があります。OSACA出力の最も右側を見ると、それはカーネルの最も内側のループ（i）のアセンブリです。そのテーブルのすべての算術命令は、依然としてスカラー演算です — vaddsd、vmulsd、vfmadd231sd、vfmadd132sd。これらはすべて、パケット化されたpd命令ではなく、%xmmレジスタで個々の倍精度数を操作しています。これを、完全にvaddpd/vmulpd/vfmadd213pdで256ビット%ymmレジスタを一度に4つのグリッドポイントで処理するJacobiの内部ループと比較してください。コンパイラはこのループをベクトル化しませんでした。 それは、Jacobiが行うように、4つの要素を1つの命令にパックするのではなく、2つの独立したスカラーチェーン（u(i,j)用とu(i+1,j)用）を公開することによって、Jacobiのスループットに到達しました。2つのレーンが実際に交差するのは、最後の行、90行目です：ストアの直前に2つの独立して計算されたスカラーを単一の128ビットレジスタに結合するvunpcklpdです。これはベクトル化された計算ではありません。それは2つのスカラー結果のベクトル化された書き込みであり、実際の作業の上に乗り、スピードアップの原因ではない、小さく無料のボーナスです。

この区別は、一見思われるよりも重要であり、次の投稿のクリフハンガーの種です。コンパイラが真のSIMDではなく、ポート間の命令レベル並列性を通じてJacobiレベルのスループットを達成しているため、このカーネルは依然として無料の広いベクトル幅を受け取ることができず、スレッド間で分割することもできません。そのわずかな4サイクルペアツーペア依存性は、どれほど安価であっても、まだ逐次的な引き継ぎであり、どのdo concurrentまたはOpenMPディレクティブも安全に分割することはできませんでした。シングルコアのラップトップでは、これは完全な成功です。しかし、スケーリングが必要なものにとっては行き止まりです。そして、結局のところ、グリッドが十分に大きくなれば、数値計算の全体的なポイントはスケーリングです。Gauss-Seidelカーネルがベクトル化され、マルチスレッド化されるようにしたい場合、展開と巧妙な代数だけではそこまで到達できません。グリッドポイントを訪問する順序を完全に変更する必要があります。これはまさに、次回、赤黒順序が登場する場所です。

---

謝辞 — クレジットを適切な人に帰さなければなりません。このGauss-Seidelの展開トリックを私に示したのはIvan Pribecです。彼がosacaを私のレーダーに乗せてくれたので、ここで提示したことはすべて彼なしでは不可能でした。本当にありがとう！

私の他の記事を読みたい場合

• Gauss-SeidelはJacobiより速い。それとも？
• Jacobi法：素朴な実装から最新のFortranマルチスレッド実装へ
• Fortranは数値線形代数の基本を教えるのにPythonより優れているか？