性能¶
このページの数値はすべて、同じ作業負荷 — 内殻イオン化テーブルの生成 — を Windows 上の 16 コア Ryzen 9 9950X で測ったものです。絶対時間より比のほうが他の環境へ持ち越し やすいので、絶対時間を書いてある箇所は上記のマシンでの値だと読んでください。
このページで使う用語
- 行とはテーブルの 1 行、つまり 1 つのチャネル (1 つの元素と 1 つの副殻 — Fe K や Au L3) の、1 つのビームエネルギー \(E_0\) における \(F(s)\) のことです。 「200 keV の Fe K」が 1 行で、v4/v5 データセットは 525 チャネル・14,796 行を 持ちます。
HIGHは本番用の求積プリセット (src/l0_numerics.jlのHIGH_SETTINGS)、-t 4は 1 プロセスに Julia スレッド 4 本という意味です。- ε ノードとは、放出電子のエネルギー ε についての求積のノード (節点) 1 つの
ことです。どの行もこのノードにわたって積分し (
HIGHでは 96 個)、スレッドが 分担するのはこのループです。 - 部分波とは、放出電子の連続状態波の角運動量成分 \(l'\) (Dirac の場合は κ) の 1 つで、それぞれが別々の動径方程式になります。200 keV の Fe K の行では、最上位の ε ノードで \(l_{\max}\) が 42 に達し、κ の値にしておよそ 85 個 (\(2 \cdot 42 + 1\)) です。下位の ε ノードはもっと少なくて済みます。
- Miller 漸化式は球ベッセル関数 \(j_\lambda(qr)\) の評価法です。必要な最大の \(\lambda\) よりずっと高い次数から下向きの漸化を始め、最後に規格化します。開始 次数から \(\lambda_{\max}\) までの段は、何も保存しない助走です。
- LPT (longest processing time first) は「重い仕事から先に配る」という スケジューリング規則です。
- ビット同一とは出力のバイト列が変わらないことで、ここに載せる最適化はすべて これを合格条件にしています。再現性の規律 を参照してください。
時間はどこへ行くのか¶
本番の 1 行 (200 keV の Fe K、HIGH 求積) を v4 処方で、シングルスレッド、後述の
最適化を施した後に、v3/v4 データセットの 161 点 s グリッド上で測ったものです。v5 は
321 点を出荷しており、s グリッドとともに増えるのは角度積分だけなので
(docs/notes/tail_contract_2026-08-09.md)、321 点ではその割合がここに示すより大きくなります:
| 領域 | 割合 | 性格 |
|---|---|---|
| 球ベッセル漸化式 | 33 % | \(\lambda\) について逐次的、動径点をまたいでは完全に独立。Miller の開始次数から \(\lambda_{\max}\) までの段のうち約 60 % は何も保存しない助走 |
| 角度積分 | 29 % | 部分波和の \((l', \lambda)\) 項ごとに Legendre 漸化式と PCHIP 評価が 1 回ずつ |
| κ 分解 Dirac 連続状態 (RK4) | 24 % | κ ごとに連立 1 階方程式 2 本 |
| その他すべて | 14 % |
これより前のプロファイル — ベッセル 58 % / 動径積分の内側ループ 24 % — は、8 レーンの ベッセルカーネルと角度融合パスの前に取ったものです。大きなものを 1 つ取り除くたびに 全体像は変わるので、毎回、着手前にプロファイルを取り直してください。
何が得られたか¶
以下の各段階はすべてビット同一です。ここには総和の順序を変えるものは 1 つも無く、
muladd・fma・@simd を縮約 (reduction) に使うものもありません。
v3 データセットを生成したコードに対して¶
| 段階 | 利得 |
|---|---|
| 8 レーン SIMD 球ベッセル + 角度融合パス | 4.3× |
動径積分テーブルのループ入れ替え、Core.Box の修正、q レーンの作業 |
2.72× — フリート A/B (48 ジョブ × 交互 2 パス) での実測 |
合わせると v3 の生成器に対しておよそ 11.7× です。8 レーンというのは AVX-512 レジスタ
1 本の幅で、8 つの動径点が 1 命令で Miller 漸化式を 1 段進みます。Core.Box の修正は
Julia のクロージャ捕捉の落とし穴を取り除いたものです — ntuple のクロージャに
捕まえられたループ変数は box 化され、8 点のグループごとにそのための動的ロードと
割り当てを払っていました。これは q レーンの作業中に、配備済みの 8 レーンベッセル
カーネル (sph_jl_tile!) へのディスパッチで見つかりました。
この 48 ジョブのうち 47 はパス間でビット同一でした。残る 1 つの例外は 再現性の規律 に書いた種類の一過性のもので、コードの性質では ありません。
v4 処方について (2026-08-08)¶
当初、v4 処方 (κ 分解 Dirac 連続状態) は 161 点グリッド上で v3 より 1 行あたり
2.2–2.9× のコストがかかりました。記録によれば、この差の大半は物理ではなく、Dirac 経路が
一度も受けていなかった最適化によるものです (後述の段階の後では v3 が平均 5.1 s/行、
v4 が 6.3 s/行)。そこで本番実行の前にもう一度プロファイルを取りました。「同じ値を二度
計算しない」(あるいは同じデータを二度流さない) という種類のビット同一の段階をさらに
5 つ重ねて、本番の 1 行で 3.9× を得ました — -t 4 で 24.7 s から 6.3 s へ、本番
5 行の平均です (200 keV の Fe K 自体は 31.6 s → 7.9 s):
| 段階 | 利得 | 内容 |
|---|---|---|
| RK4 のポテンシャル標本を κ をまたいで共有 | 1.72× | 実行時間の 37 % が 1 つのスプライン評価で、κ に依存しない点でおよそ 85 本の部分波ごとにやり直していた |
| \(R(Q)\) の \(Q_+\) 側を ε ノードごとに 1 回へ巻き上げ | 1.21× | \(Q_+^2 = k_i^2 + k_f^2 - 2 k_f k_i \cos\theta\) は方位角にも \(K\) にも依存しないのに、161 × 48 回も再補間していた |
| Legendre 漸化式を格子点 8 個でインターリーブ | 1.15× | この漸化式は除算 1 回のレイテンシ律速 (実測約 16 サイクル)。格子点どうしは独立 |
ε ループを :greedy + 降順 (LPT) に |
1.25× | コストは ε とともに増えるので、連続チャンクの分配ではスレッドが遊んでいた |
| q レーン SIMD 累積を Dirac 動径テーブルへ移植 | 1.21× | 非相対論経路にしか書かれていなかったが、その経路は v4 で出荷経路ではなくなった |
最後の 1 つが教訓です。あるコード経路が「比較専用」から「これが出荷される」へ昇格する とき、その最適化は一緒には付いてきません。その経路が一度も受けていない最適化を 探してください。
測って却下したものも含む完全な記録は docs/notes/speedup_v4_2026-08-08.md です。
並列化: スレッドよりプロセス¶
スレッドレベルの並列化は ε ノードにわたって走り、8 スレッドよりずっと手前で飽和します。 同じコア数を、スレッド数の少ないプロセスに多く分けるほうが劇的に良くなります:
| 構成 | 結果 |
|---|---|
| 4 プロセス × 8 スレッド | 基準、CPU 使用率 68 % |
| 8 プロセス × 4 スレッド | 2.26×、CPU 使用率 94 % |
| 16 プロセス × 2 スレッド | CPU 使用率 100 %、ただし 8 プロセスの場合に対して 1.16× だけ |
重要なのは最後の行です。CPU 使用率は上がったのにスループットはほとんど動いていません。 増えた CPU 時間はメモリ待ちに費やされたからです。この作業負荷は相当程度メモリ帯域 律速なので、プロセスを足すやり方はすでに頭打ちの手前 (曲線の膝) まで来ています。
これは、コアを増やすより SIMD を選ぶ根拠でもあります。1 命令で 8 要素を処理すれば、 動かすバイトあたりの仕事が増えます。
本番ではこれが 8 プロセス (フリートのスクリプトではレーンと呼びます) × 4 スレッドの 配置です。v4 の測定では 43.7 行/分で、14,796 行の生成は約 5.6 h と見積もられました (v4 の実行そのものは 5 h 16 min、321 点の v5 の実行は 6 h 47 min かかりました)。
角度融合パスの後の再測定 (2026-08-05)
上の表は元の測定で、角度融合パスが割り当てを 96 % 削る前に取ったものです。
最適化の監査では、2.26× の交絡要因として GC 圧力の可能性が指摘されていました。
tools/bench_e1/ のハーネス (48 ジョブ、HIGH、161 点グリッド、v3 処方、構成
ごとに 2 パス。プロトコルは同ディレクトリの README に書いてあり、結果ファイルは
ローカルに保管して追跡していません) で、そのパス後のコードで比較をやり直し
ました: 8 × 4 が依然として最速の配置で、1 × 32 はそのスループットの 0.72–0.75、
16 × 2 は 0.86–0.90 でした。順位が変わらなかったので 8 × 4 が本番配置のままです。
--gcthreads=1 は約 3 % のコストで、本番では常に渡しています
(src/gen_production.jl を参照)。これは Windows の GC クラッシュを無くすものでは
なく、さらされる機会を減らすだけです。
なぜ C++ や C# へ移植しないのか¶
仮定ではなく実測です:
- Julia は、この形のスカラー数値カーネルでは C++ と本質的に同じ速さです。
- GC が総実行時間に占めるのは 4.3 % です。
- したがって素直な移植で買えるのは 1.0–1.2× です。
- 監査した 43 件の最適化案のうち、言語を変えることで可能になるものはゼロでした。
律速しているのは FP64 除算のスループット、FP 加算のレイテンシ連鎖、そしてこの プロジェクト自身のビット同一の規律の 3 つで、どれも実装言語を変えても変わりません。 決定は Julia に留まることです。見直すのは、上流の GC 修正が入った後もランタイムが 長時間の本番実行を壊し続ける場合に限り、そのときでも最初の一手は小さな FFI カーネルを 切り出すこと (C で書いたホットループ 1 つを Julia の外部関数インタフェース経由で呼ぶ) であって、全部を書き直すことではありません。
なぜ GPU ではないのか¶
| FP64 スループット | |
|---|---|
| RTX 4060 (コンシューマ向け Ada、FP64 = FP32/64) | 0.24 TFLOPS |
| AVX-512 付き Ryzen 9950X | ≈ 2.2 TFLOPS |
実際に手元にあるハードウェアでは、FP64 では CPU のほうが 9× (9 倍) 上です。FP32 へ 落とす道筋は示せていません。いま書かれている Miller 漸化式はおよそ 10⁻²⁵⁰ から 10²⁵⁰ の 間を振れ、FP32 では担えない再スケーリングの段を必要としますし、FP32 で安全な定式化を 私たちは持っていません。ホットスポットは逐次的な依存連鎖とメモリ律速の内積で、この コードではどちらも GPU にうまく載ることが示されていません。ここでの結論は手元の ハードウェアに限ったものです。
もし見直すことがあれば、指標は FLOPS ではなく帯域です — A100 の 2 TB/s、MI300X の 5.3 TB/s に対して DDR5 は 80–90 GB/s です。再検討の条件は次の 3 つがすべて揃うことです: 作業負荷が 20–50 倍に増えること、アルゴリズム側の案が尽きること、精度ゲートがもっと 緩い \(|\Delta F / F|\) に移ること。
うまくいかなかったもの¶
もっともらしく見え、否定するのに時間がかかったので記録しておきます:
| 案 | 実測 |
|---|---|
| 漸化式で除算の代わりに逆数を使う | 1.03× — 連鎖はレイテンシ律速 (しかも c*(1/X) は c/X と丸めが違うので、どのみちビット同一にならない) |
--heap-size-hint |
効果ゼロ。GC 回数 151 → 151 |
| BigInt 階乗の代わりに参照表 | 割り当て 0.777 → 0.775 GB |
動径ループの転置 + @simd |
局所では 8.3×、全体では 1.27×、しかもビット同一ではない → 却下 |
まだ残っている手¶
- 動径積分テーブルを E₀ の行をまたいで共有する。\(R\) は E₀ に依存しないのに、1 つの チャネルが持つ ~30 本の E₀ 行それぞれで、その 82 % を計算し直しています。残っている 構造的な利得ではこれが最大ですが、まだ実測していません。監査の見積りは、段階版 (E₀ に依存しない第 1 区間の ε ノードをキャッシュする。~30 MB) で全体 1.15–1.25 倍、 完全版 (ε のマスタ格子へ転置する) で全体 2.0–2.5 倍、82 % がコストゼロで消えた場合の Amdahl 上限が約 5.6 倍です。ビット同一は前提にできません — 完全版は ε の標本化と 行チェックポイントの単位を変えるので、リファクタではなくデータセット生成の話になります。
- 対数動径グリッドを再利用する。メッシュは内側の領域で対数なので、\(x = qr\) の 相異なる値は約 6.2 M から約 1.2 M へ潰れます — 全体で約 1.9× です。ビット同一 ではない (q グリッドが変わる) ので、次のデータセット生成まで保留しています。
- λ と Miller の開始次数を動径点ごとに打ち切る。連続状態波は \(r^{l+1}\) が \(e^{-60}\) に 達するところから種を置くので、小さい半径では高い λ を持つ部分波の項は一切寄与しません — \(l = 42\) では種がグリッドの 92 % より先にあります。そこで \(\lambda_{\max}\) に上限を 付ければ Miller 漸化式もそれに合わせて短くなり、全体でおよそ 1.2× の価値が あります。ビット同一ではない (開始次数が値を相対 10⁻¹³ ほど動かす) ので、 求積をどのみち変える世代まで保留しています。
- さらに 6 件の監査済み候補がまだ判定待ちで、
docs/notes/speedup_audit_2026-08-05.jsonに一覧してあります (同じ台帳に、かつて判定待ちだった項目のうち v4 で採用したもの・ 2026-08-05 の角度側書き換えで対象が消えたもの・上の λ 打ち切りと一緒に保留したものを、 2026-08-19 付で記録してあります)。
測って却下したもの: 寄与の小さい部分波を、組み立てる前に事前に除外することです。 部分波の 95–99 % は有意性フィルタを通過するので、飛ばせるものがありません。
ベンチマークの規則¶
- BenchmarkTools.jl は別の環境から使ってください
(
julia --project=/some/scratch/benchenv)。リポジトリは依存ゼロのままにします。 - 統計量は
minを報告します。2–3 % 未満の差はノイズとして扱ってください。 - 全体の数値とマイクロベンチマークは日常的に食い違います — 8.3× のループ高速化が 全体では 1.27× でした。両方を書いてください。
tools/bench_e1/のベンチマークドライバは全コアを飽和させ、PowerShell 7+ を 必要とします。