部分ピボット LU が 8T 以上で伸びない: 現行構造のボトルネックと並列向き手法の文献調査 #123

背景

メンテナ指示 (2026-10-29): 「ピボット LU の並列性能が出ないのは既知の問題。並列化に向いた新しい手法が提案されているか research させて」。

#109 の broccoli 計測 (05:26Z コメント; EPYC 7702P、Zen2 64 コア) では、lu 1024 f64 が次のとおりで、8T 以上ほぼ伸びません。OpenBLAS も 16T 以上で平坦です (41.3 → 54.0 ms)。

この issue は、コードを読んだ結果と文献の整理だけをまとめたものです。計測と実装はしていません。

1. 現行 LU の構造とボトルネック (origin/main ad23987, faer/src/linalg/lu/partial_pivoting/factor.rs)

構造 (コードで確認済み)

Toledo / Gustavson 型の再帰 right-looking LU です。

• lu_in_place_recursion (:158) は列を半分に割ります。ブロック幅は 16 の倍数に丸めます (:174-176)。
• 処理の順序は次のとおりです。
  1. 左半分を再帰で分解する
  2. trailing update (A01 の TRSM と A11 -= A10 A01 の GEMM) を with_width の幅で実行する (:204-221)
  3. 右半分を再帰で分解する
  4. swap stage で、遅延した行交換をこのレベルの外側の列に適用する (:231-311)
• 葉 (recursion_threshold = 16, :345) は lu_in_place_unblocked (:19-80) です。m 行 × 最大 16 列の panel を完全に直列で処理します。各列について次を行います。
  • m 行の最大値探索
  • panel 内の行交換
  • スケーリング
  • matmul(.., Par::Seq) による rank-1 更新
• 並列になるのは、各再帰ノードの trailing update と swap stage だけです。幅は trailing_update_width / swap_stage_width (:108, :138) が cost ルールで決めます。sections はそれぞれ 2 と 4 (片側だけのときは 1) です。
• lookahead はありません。右半分の分解は、そのノードの trailing update がすべて終わってから始まります (ノードごとに fork-join)。

ボトルネック (推測。下の最小実験で確かめる)

• クリティカルパスは panel の連鎖です。 葉 n/16 個は直列で、それぞれ m 行にわたる level-2 演算 (メモリ律速) です。さらに、葉に近いレベルの trailing update (m × 16 × 16、m × 32 × 32 など) は小さいため、幅 1 か狭い幅で走ります。このレベル数は p に依存しません。
• par: every region at the widest width that pays, on a sub-team of that width #109 の数値を Amdahl でみると、非並列部分は約半分です。 1T 43.1 ms、8T 24.1 ms から逆算すると、非並列部分は 1T 時間の約 49% (約 21 ms) です。64T の 25.5 ms はこの値と整合します。つまり、残っているのは並列部分の効率ではなく、直列部分です。
• 同期回数は主因ではないと見ています。 再帰ノードは約 2n/16 = 128 個で、region は約 250〜400 回です。dispatch 1 回が数 µs なら合計 1 ms 前後です。
• 行交換は列方向に並列化されていますが、1 列の中の n 回の swap は逐次で、メモリ律速です。葉の中の行交換は直列です。
• OpenBLAS (記憶による。コードは未確認): getrf_parallel は lookahead を持ちますが、panel (getrf_single の再帰) は 1 スレッドで処理します。16T 以上で平坦なのは同じ理由だと思われます。

2. 並列向き手法の比較

2020〜2026 年の新提案について: 共有メモリのメニーコア CPU で GEPP の panel をそのままスケールさせる新手法は、今回の調査では見つかりませんでした (検索は数回に限っています)。最近の研究は SLATE 系 (分散 / GPU) の pivot 回避や緩和 (F, G, H) と、静的 lookahead の OpenMP 実装 [9b] が中心です。

3. 推奨

• lu() は GEPP の P を守ったまま、B (lookahead) と A (並列 panel) を入れる方向を勧めます。 tensor4all は LAPACK 互換の置換を前提にしているので、P が変わる D / E / F は lu() には入れません。入れるとしても opt-in のパラメータにします。
• どちらを先にするかは、下の最小実験の結果で決めます。
  • 葉と、幅 1 で走る update が 64T 時間の大半を占める場合: A を先にする。葉の 16 列 panel を team で並列にし、葉の閾値も見直す。
  • 幅 > 1 の update の待ちが大きい場合: B を先にする。
  • 64 コア級で効かせるには、どのみち両方が要るという見立てです ([4] でも panel 単独は 8〜12 コアで頭打ち)。
• threadteam との相性は良いと見ています。sub-team への分割 (par: every region at the widest width that pays, on a sub-team of that width #109 で入った「その幅の sub-team で実行」) が、B の panel と update の分業、および [9] の malleable / early termination にそのまま対応します。A に必要な team 内の低遅延 barrier は、threadteam 側の追加として検討が要ります。
• G / H / I は将来の別 API 候補です。 pivot を使わない solve() の高速経路として、lu() とは分けて検討します。I (混合精度) は panel 律速を解消しないので、優先度は低くします。

4. まず測るべき最小実験 (1 つ)

64T でのクリティカルパスの内訳を測ります。

• 方法: feature で囲んだタイマー (par-trace 相当) で、lu_in_place_recursion の時間を次の 4 つに分けて積算する。
  • (i) lu_in_place_unblocked (葉)
  • (ii) 幅 1 で走った trailing update
  • (iii) 幅 > 1 の trailing update (region の壁時計時間)
  • (iv) swap stage
• 条件: f64、1024 (参考として 2048)、C / F order、1 / 8 / 64T。
• 判定:
  • (i) + (ii) が 64T の時間 (約 25 ms) の大半を占める → A を先に
  • (iii) が支配的で、その中で待ちが大きい → B を先に
• この実験は、par: every region at the widest width that pays, on a sub-team of that width #109 の Amdahl 逆算 (非並列部分が約 49%) を裏付けるものでもあります。
• 実施時期: broccoli で走っている計測ジョブが終わってから。8T の内訳は evo で先に取れます。

出典

確認済み = DOI の Crossref レコードを開いて書誌を照合したもの。未確認 = 検索結果や引用でしか見ていないもの。

1. S. Toledo, "Locality of reference in LU decomposition with partial pivoting", SIAM J. Matrix Anal. Appl. 18(4), 1997. https://doi.org/10.1137/S0895479896297744 (確認済み)
2. F. G. Gustavson, "Recursion leads to automatic variable blocking for dense linear-algebra algorithms", IBM J. Res. Dev. 41(6), 1997. https://doi.org/10.1147/rd.416.0737 (確認済み)
3. A. M. Castaldo, R. C. Whaley, "Scaling LAPACK panel operations using parallel cache assignment", PPoPP 2010. https://doi.org/10.1145/1693453.1693484 ; journal 版 Castaldo, Whaley, Samuel, ACM TOMS 39(4), 2013. https://doi.org/10.1145/2491491.2491492 (確認済み)
4. J. Dongarra, M. Faverge, H. Ltaief, P. Luszczek, "Achieving numerical accuracy and high performance using recursive tile LU factorization with partial pivoting", Concurrency Computat.: Pract. Exper. 26(7), 2014. https://doi.org/10.1002/cpe.3110 (確認済み)
5. A. Buttari, J. Langou, J. Kurzak, J. Dongarra, "A class of parallel tiled linear algebra algorithms for multicore architectures", Parallel Computing 35(1), 2009. https://doi.org/10.1016/j.parco.2008.10.002 (確認済み)
6. G. Quintana-Ortí, E. S. Quintana-Ortí, R. A. van de Geijn, F. G. Van Zee, E. Chan, "Programming matrix algorithms-by-blocks for thread-level parallelism", ACM TOMS 36(3), 2009. https://doi.org/10.1145/1527286.1527288 (書誌は確認済み。incremental pivoting の記述は未確認)
7. L. Grigori, J. Demmel, H. Xiang, "CALU: A communication optimal LU factorization algorithm", SIAM J. Matrix Anal. Appl. 32(4), 2011. https://doi.org/10.1137/100788926 (確認済み)
8. S. Donfack, L. Grigori, A. K. Gupta, "Adapting communication-avoiding LU and QR factorizations to multicore architectures", IPDPS 2010. https://doi.org/10.1109/IPDPS.2010.5470348 (書誌のみ確認済み)
9. S. Catalán, J. R. Herrero, E. S. Quintana-Ortí, R. Rodríguez-Sánchez, R. van de Geijn, "A case for malleable thread-level linear algebra libraries: The LU factorization with partial pivoting", IEEE Access 7, 2019. https://doi.org/10.1109/ACCESS.2019.2895541 (確認済み)
   • 9b. S. Catalán, A. Castelló, F. D. Igual, R. Rodríguez-Sánchez, E. S. Quintana-Ortí, "Programming parallel dense matrix factorizations with look-ahead and OpenMP", Cluster Computing 23(1), 2020. https://doi.org/10.1007/s10586-019-02927-z (確認済み)
10. J. Kurzak, P. Luszczek, M. Faverge, J. Dongarra, "LU factorization with partial pivoting for a multicore system with accelerators", IEEE TPDS 24(8), 2013. https://doi.org/10.1109/TPDS.2012.242 (確認済み)
11. M. Baboulin, J. Dongarra, J. Herrmann, S. Tomov, "Accelerating linear system solutions using randomization techniques", ACM TOMS 39(2), 2013. https://doi.org/10.1145/2427023.2427025 (確認済み)。D. S. Parker, random butterfly transformation の UCLA 技術報告 (1995) は未確認 (後続論文の引用でのみ確認)。
12. J. Peca-Medlin, T. Trogdon, "Growth factors of random butterfly matrices and the stability of avoiding pivoting", SIAM J. Matrix Anal. Appl. 44(3), 2023. https://doi.org/10.1137/22M148762X (確認済み)
13. N. Lindquist ら:

• a. N. Lindquist, P. Luszczek, J. Dongarra, "Replacing pivoting in distributed Gaussian elimination with randomized techniques", ScalA 2020. https://doi.org/10.1109/ScalA51936.2020.00010
• b. N. Lindquist, M. Gates, P. Luszczek, J. Dongarra, "Threshold pivoting for dense LU factorization", ScalAH 2022. https://doi.org/10.1109/ScalAH56622.2022.00010
• c. N. Lindquist, P. Luszczek, J. Dongarra, "Using additive modifications in LU factorization instead of pivoting", ICS 2023. https://doi.org/10.1145/3577193.3593731
• d. 同著者, "Generalizing random butterfly transforms to arbitrary matrix sizes", ACM TOMS 50(4), 2024. https://doi.org/10.1145/3699714
• e. 同著者, "The stability of block eliminations and additive modifications", arXiv:2509.07305, 2025. https://arxiv.org/abs/2509.07305
• いずれも確認済み。

14. M. Gates, J. Kurzak, A. Charara, A. YarKhan, J. Dongarra, "SLATE: design of a modern distributed and accelerated linear algebra library", SC19. https://doi.org/10.1145/3295500.3356223 (確認済み)。現行の SLATE は LU の方式に PartialPiv / CALU / NoPiv / RBT / BEAM を持つ (GitHub icl-utk-edu/slate の enums で確認)。panel の多スレッド化は未確認。
15. J. Dongarra ら, "PLASMA: Parallel Linear Algebra Software for Multicore Using OpenMP", ACM TOMS 45(2), 2019. https://doi.org/10.1145/3264491 (確認済み)
16. 混合精度:

• a. E. Carson, N. J. Higham, "Accelerating the solution of linear systems by iterative refinement in three precisions", SIAM J. Sci. Comput. 40(2), 2018. https://doi.org/10.1137/17M1140819
• b. A. Haidar, S. Tomov, J. Dongarra, N. J. Higham, "Harnessing GPU tensor cores for fast FP16 arithmetic to speed up mixed-precision iterative refinement solvers", SC18. https://doi.org/10.1109/SC.2018.00050
• c. N. J. Higham, T. Mary, "Mixed precision algorithms in numerical linear algebra", Acta Numerica 31, 2022. https://doi.org/10.1017/S0962492922000022
• いずれも確認済み。

17. その他:

• iv. A. Khabou, J. Demmel, L. Grigori, M. Gu, "LU factorization with panel rank revealing pivoting and its communication avoiding version", SIAM J. Matrix Anal. Appl., 2013. https://doi.org/10.1137/120863691 (未確認)
• v. RecursiveFactorization.jl の告知 (開発者自身の主張で、査読なし). https://discourse.julialang.org/t/ann-recursivefactorization-jl/39213 (未確認)

親: #95

🤖 Generated with Claude Code