p 個の部分バッファを逐次合算する reduction (tridiag, matvec_colmajor, 疎×疎, compact-WY upgrade) の高並列スケーリング #93

問題

p 個の部分和/部分バッファを作り、呼び出しスレッドが逐次に合算する箇所があります。合算は O(p · size)、メモリも O(p · size) で、p = 32〜64 では Amdahl 項として効きます。構造は upstream にもともとありますが、fork がゲートを置き換えたため、どこから並列に入るかが変わっています。

tridiag の損益分岐 (推測、要素単位の粗い見積もり) は、p · rem と rem²/(2p) が等しくなる p* ≈ √(rem/2) です。n = 1024 では p* ≈ 22、n = 4096 では約 45 になり、p = 64 では n ≲ 4096 の eigh が p = 16〜32 より遅くなりえます。合算は p コアのキャッシュから読むので、CCD をまたぐコストも加わります。

OpenBLAS (記憶による、コードは未確認): dsytrd は flop の半分を dsyr2k (level 3) で処理し、dsymv の thread 化にも per-thread バッファの reduce があります。dgemv_n は行を分割するので reduce はありません。

確認済み / 推測

• 確認済み (コード): 上表の確保サイズ、逐次ループ、sparse の Seq 時の p 列確保、tridiag の入口が upstream と比べて下がったこと。
• 推測: 損益分岐の値、実際の逆転が起きる n。

計測計画 (faer#90 の計測が終わってから)

• tridiag_in_place / eigh 値のみ、f64/c64、n = 512, 1024, 2048, 4096、1/8/16/32/64T。perf で tridiag.rs:458 の合算ループの比率を確認します。p に対して単調に減らなければ確定です。
• matvec_colmajor: n×n (n = 512, 2048, 8192) と 100000×256、f64/c64。
• 疎×疎: n = 1e4/1e5/1e6、1 列あたり nnz 4〜16。

推奨

• sparse: 修正。workspace 幅を実際のスレッド数 (Seq なら 1) にし、各スレッドが自分の列をゼロにする。
• tridiag、matvec_colmajor: 計測で確認してから修正。案は 3 つあります。
  • 部分バッファ数を min(p, c·√rem) で cap する
  • 合算を行スライスで並列 reduce する
  • matvec は m が十分なら行分割にする
• upgrade_householder_factor: 問題なし (監視のみ)。

親: #95

🤖 Generated with Claude Code