格子 QCD における混合精度ソルバー
金森 逸作 (理研 R-CCS)
2023 年 12 月 8 日 (online)
第 20 回 HPC-Phys 勉強会Outline••

混合精度ソルバーのアルゴリズム
ドメインウォール型格子フェルミオンでの比較
benchmark には Bridge++ に手を加えたものを利用
開発リソース：「富岳」ra000001, Wisteria Odyssey (筑波大学際共同利用)•(近似式の検証：精度が足りずに（？）最近ハマってしまったこと)

金森逸作: 第 20 回 HPC-Phys 勉強会, 2023 年 12 月 8 日 (online)

2/15

格子 QCD での線形ソルバー
格子 QCD （格子上の場の理論）•

時空を 4 次元の格子に切って、場の量子論記述
無限自由度の系の正則化とか、対称性とか。単に計算機に載せるための近似ではない•

計算機上で、経路積分を実行（モンテカルロ積分）
RRdU dψ exp(−Sgauge [U ] − ψD[U ]ψ) = dU det D[U ] exp(−Sgauge [U ]•••

U : ゲージ場 (グルーオン)
ψ: フェルミオン（クォーク）。フェルミオン 2 つ (u,d quark) の寄与なら det D† D
D: Dirac 演算子。離散化の仕方はいろいろある。 Wilson, Domainwall, staggered etc.
格子上では巨大な疎行列 (rank: O(108 )–O(1010 ))ソルバー？ Ax = bx = A−1 bR••分布の作成（配位生成）: det D D = dφdφ exp(−φ(D† D)−1 φ)
フェルミオンの相関関数 : hψψi = D−1•D(x, y)η(y) = δ(x, x0 ) を解くと η(y) = D−1 (y, x0 ) = hψ(y)ψ(x0 )i が求まる
計算時間の 8 割程度はソルバー パラメータによる†

金森逸作: 第 20 回 HPC-Phys 勉強会, 2023 年 12 月 8 日 (online)

3/15

なぜ混合精度か•••

単精度演算のほうが倍精度演算より速い: SIMD, GPUAI の影響で、単精度や半精度の演算性能が大きく向上
メモリアクセスでも単精度は倍精度の半分
QCD の計算はメモリバンド幅・通信律速「京」は倍精度演算しかなかったが、データを単精度で保持することで B/F 比を稼いだ (LDDHMC)••

隣接通信の通信量も、単精度は倍精度の半分ソルバーの反復数（行列・ベクトル積の回数）は増える

金森逸作: 第 20 回 HPC-Phys 勉強会, 2023 年 12 月 8 日 (online)

4/15

残差反復 (preconditioned Richardson iterations)
欲しいもの: x s.t. Dx = b
近似解は手元にある: xn(approx. x)r = b − Dxn として、xn + D−1 r は解になっている:
D(xn + D−1 r) = Dxn + r = b••−1r なら計算できるものとする
D−1 r は無理でも、近似的な Dapp.
input: x0

r0 := b − Dx0
2 i = 03 while |ri | is not small enough do−1ri
4xi+1 := xi + Dapp
5ri+1 := b − Dxi+16i := i + 1

1−1に単精度（低精度）のソルバーを用いる
この Dapp
金森逸作: 第 20 回 HPC-Phys 勉強会, 2023 年 12 月 8 日 (online)

5/15

混合精度ソルバー
黒：倍精度

赤： 単精度

r0 := b − Dx0
2 i = 03 while |ri | is not small enough do(single)
4
ri:= to single(ri /|ri |)−1 (single)
ri5
δx(single) := DappKrylov 部分空間を使った反復法で解く（単精度ソルバー）
6xi+1 := xi + to double(δx(single) ) × |ri |7ri+1 := b − Dxi+1
8i := i + 11••

演算のほとんどが、単精度ソルバー残差を倍精度で計算することが重要

金森逸作: 第 20 回 HPC-Phys 勉強会, 2023 年 12 月 8 日 (online)

6/15

単精度（低精度）ソルバーの問題点
反復法: （例）CG 法で Ax = b を解く

r0 = b − Ax0
2 p0 = r03 i=04 while |ri | is not small enough do(r ,p )
内積 (a, b) の計算
5αi = (pi i,Api i )1

7xi+1 = xi + αi pi
ri+1 = ri − αi Api

8βi =9pi+1 = ri + βi pii := i + 1

610(ri+1 ,ri+1 )(ri ,ri )•••小さ過ぎる••

内積やノルムが正しく求まらない部分和を足し上げるなどの工夫が必要| 内積 | ≪ | ノルム | のときは桁落ちの恐れ

とくに BiCGStab 法•
ノルムの計算内積、ノルムの計算：要素数は O(108 )–O(1010 )
単精度の有効桁：7 桁程度 （半精度：3 桁程度）
= 108
単精度だと、108 + |{z}1•

内積やノルムの計算だけ倍精度（あるいは多倍長）にするのも選択肢
並列計算では、結合則の破れ(x + y) + z 6= x + (y + z) もある
reproducible summation (Ahrens+ 2020) の利用:
Lattice 2023, Kate Clark
https://indico.fnal.gov/event/57249/contributions/270632/

金森逸作: 第 20 回 HPC-Phys 勉強会, 2023 年 12 月 8 日 (online)

7/15倍精度、単精度 + 倍精度、 半精度 + 倍精度
フェルミオンは Möbius Domainwall fermion: 32 × 8 × 8 × 12 (行列のランクは 3.5 × 106 )半精度は任意精度ライブラリ mpfr で代用 (*) (mpreal https://github.com/advanpix/mpreal 経由)
幸谷さんの「多倍長精度数値計算」（森北出版）が参考になりました(*) 仮数部分は 11bit, 指数部分は半精度より大きい
fp32: target low prec (tol)2=1e-10

fp64

mpfr: signifcant=11bit, target low prec (tol)2=1e-4

1

1double prec.
1e-28

low prec solver
outer (in double)
1e-28

-5

1x10-15

residual norm squared

residual norm squared

1x10-10516 mults

1x10-20

1x10-5

1x10-5

1x10-10

1x10-10

1x10-15588 low prec. mults

1x10-20

1x10-15

1x10-25

1x10-25

1x10-30

1x10-30

1x10-30

0200

400

600

800

1000

1x10-35

0200

iterations

400

600
iterations

800

1000

762 low prec. mults

1x10-20

1x10-25

1x10-35•••

residual norm squared

1x10low prec solver
outer (in double)
1e-28

1x10-35

7.5e-27

0200

400

600

800

1000

iterations

収束条件：|ri |2 < 10−28
単精度：|ri |2 < 10−10
半精度：|ri |2 < 10−4

金森逸作: 第 20 回 HPC-Phys 勉強会, 2023 年 12 月 8 日 (online)

8/15

内積とノルムの精度を上げる••（左）low prec: 半精度相当 (11 bits) ＋内積とノルムは倍精度相当 (53 bits)（左）low prec: すべて半精度相当 (11 bits)
22mpfr: signifcant=11bit, target low prec (tol) =1e-4

mpfr: signifcant=11bit (norm: 53bits), target low prec (tol) =1e-4

11
low prec solver
outer (in double)
1e-28

low prec solver
outer (in double)
1e-28

-5

1x10-5-10

1x10

1x10-10-15

762 low prec. mults

1x10-20

1x10-25

1x10-15762 low prec. mults

1x10-20

1x10-251x10

7.5e-27-30

7.0e-27-30

1x10

1x10-35

1x10residual norm squared

residual norm squared

1x10-35

0200

400

600

800

1000

1x100200

iterations

400

600

800

1000

iterations

差はわずか：収束直前で、7% (7.5 × 10−27 vs. 7.0 × 10−27 )系が大きくなると差も広がると予想
金森逸作: 第 20 回 HPC-Phys 勉強会, 2023 年 12 月 8 日 (online)

9/15

より現実的な場合
Wisteria Odssey 12 nodes で実行 (筑波大学際共同利用)

格子サイズ：483 × 12 （行列ランク 1.9 × 108 ）••••収束条件：|ri |2 < 10−26
単精度：|ri |2 < 10−12
ノルム・内積が倍精度の方が多少速い
単精度はもっと早く打ち切った方がよさそう
single: (norm: single), target low prec (tol)2=1e-12

doublesingle: (norm: double), target low prec (tol)2=1e-12
1double1e-26

1x10-5

1x10-10low prec solver
outer (in double)
1e-26
1x10-5

1x10-10439 sec., 11460 mults

1x10-20

1x10-10

1x10-15

residual norm squared

1x10-15

residual norm squared

residual norm squared

low prec solver
outer (in double)
1e-26

1x10-5605 sec., 29362 low prec. mults

1x10-20

1x10-15

1x10-20

1x10-25

1x10-25

1x10-25

1x10-30

1x10-30

1x10-30

1x10-35

0

5000

10000

15000

20000

25000

30000

1x10-35

05000

iterations

倍精度

10000

15000

20000

25000

iterations

30000

585 sec., 27825 low prec. mults

1x10-35

05000

10000

15000

20000

25000

単精度＋倍精度

単精度＋倍精度ノルム・内積は単精度

ノルム・内積は倍精度

金森逸作: 第 20 回 HPC-Phys 勉強会, 2023 年 12 月 8 日 (online)

30000

iterations

10/15単精度のソルバーの収束精度依存性
single: (norm: single), target low prec (tol)2=1e-12

single: (norm: single), target low prec (tol)2=1e-10

single: (norm: single), target low prec (tol)2=1e-8

1low prec solver
outer (in double)
1e-26

low prec solver
outer (in double)
1e-26

1x10-15

residual norm squared

1x10-10605 sec., 29362 low prec. mults

1x10-20

1x10-5

1x10-5

1x10-10

1x10-10

1x10-15430 sec., 20231 low prec. mults

1x10-20

1x10-15

1x10-25

1x10-25

1x10-30

1x10-30

1x10-30

0

5000

10000

15000
iterations

20000

25000

1x10-35

30000

0

5000

single: (norm: double), target low prec (tol)2=1e-12

10000

15000
iterations

20000

25000

1x10-35

30000

1x10-10585 sec., 27825 low prec. mults

1x10-20

1x10-15

residual norm squared

1x10-10residual norm squared

1x10-10436 sec., 20195 low prec. mults

1x10-20

1x10-15

1x10-25

1x10-25

1x10-30

1x10-30

1x10-30

20000

25000

単精度：|ri |2 < 10−1220000

30000

1x10-35

25000

30000

05000

10000

15000
iterations

20000

25000

単精度：|ri |2 < 10−1030000

394 sec., 18067 low prec. mults

ノルム・内積は倍精度

1x10-20

1x10-25

15000
iterations

15000
iterations

low prec solver
outer (in double)
1e-26
1x10-5

10000

10000

low prec solver
outer (in double)
1e-26
1x10-5

5000

5000

1

1x10-500single: (norm: double), target low prec (tol)2=1e-8

1low prec solver
outer (in double)
1e-26

1x10-35

ノルム・内積は単精度

single: (norm: double), target low prec (tol)2=1e-10

11x10-15

389 sec., 18059 low prec. mults

1x10-20

1x10-25

1x10-35

residual norm squared

low prec solver
outer (in double)
1e-26

residual norm squared

1x10-5residual norm squared

11x10-35

0

5000

10000

15000
iterations

20000

25000

30000

単精度：|ri |2 < 10−8※ 上記より更に精度を落としたら実行時間が増加した
金森逸作: 第 20 回 HPC-Phys 勉強会, 2023 年 12 月 8 日 (online)

11/15

精度に関わる小ネタ

金森逸作: 第 20 回 HPC-Phys 勉強会, 2023 年 12 月 8 日 (online)

12/15

分数べき乗の近似

Pnai
分数べき乗の近似：x = a0 + i=1 x+bi†(QCD の Rational HMC で 行列 D D の分数べき乗が必要になる)
1/4

係数を求めるプログラム AlgRemez https://github.com/maddyscientist/AlgRemez
近似式の精度を確認したい （正しい係数を使っていることを確認したい）••

approx vs. pow: x1/4

relative error
-11

1.8

1x10

pow(x,0.25)
approx

1.6

double5.571223e-15

1x10-121.4
1x10-131.2
11x10-140.8
1x10-150.6
0.4

1x10-16

0.2
00

12

3

4

56

7

8

9

10

1x10-170.0001

0.001

0.01

相対誤差が AlgRemez が出した 5.571223e-15 より大きい
金森逸作: 第 20 回 HPC-Phys 勉強会, 2023 年 12 月 8 日 (online)

0.1

110

桁落ち？
a0 = 5.5803685240057019
13/15

おそらく、、、

mpfr を使った検証
abosolute error

relative error

-13

-11

1x10

1x10

doublempfr (200 bits)
5.571223e-15

1x10-14

doublempfr (200 bits)
5.571223e-15

-12

1x10

1x10-13
1x10-15

1x10-14

1x10-16

1x10-15
1x10-16

1x10-17

1x10-17
1x10-18
1x10-190.0001

1x10-18
0.001

0.01

0.1

110

1x10-190.0001

0.001

0.01

0.1

110

おそらく x < 1 は相対誤差（左）、1 < x は絶対誤差（右）になっている。ドキュメントで仕様を確認する必要あり。

金森逸作: 第 20 回 HPC-Phys 勉強会, 2023 年 12 月 8 日 (online)

14/15

まとめと議論•••

格子 QCD で用いている混合精度ソルバー紹介
単精度＋倍精度：単精度ソルバーの収束条件に注意
内積、ノルムの精度も大事：精度が落ちる前に打ち切るのが良さそう
今回は調べていないが、多倍長も有効かも••

半精度もたぶん行ける？ To do: A64FX の半精度演算
体積が大きくなるともっと大変になる
cf. メトロポリステストで ∆H = H − H0 = O(1010 ) − O(1010 ) ∼ 0.1

金森逸作: 第 20 回 HPC-Phys 勉強会, 2023 年 12 月 8 日 (online)H ∝ volume

15/15

Backup Slides

金森逸作: 第 20 回 HPC-Phys 勉強会, 2023 年 12 月 8 日 (online)

appendix-i

部分和Bridge++ での例: ノルムの計算 (「富岳」向けのコードから)
123456789101112131415161718192021222324

set threadtask(ith, nth, is, ns, Nouter);// range for each thread
svreal t yt; set vec(pg, yt, 0.0);// yt=0 (SIMD variable)for (int ex = 0; ex < Nex; ++ex) {svreal t tmp yt;
set vec(pg, tmp yt, 0.0);// tmp yt=0
for (int i = is; i < ns; ++i) {real t *v = &vp[VLEN * Ninner * (i + Nouter * ex)]; // real t: float or double
svreal t xt, vt;
set vec(pg, xt, 0.0);// xt=0for (int in = 0; in < Ninner; ++in) {load vec(pg, vt, &v[VLEN * in]);
add norm2 vec(pg, xt, vt);// xt += vt*vt (1)}add vec(pg, tmp yt, xt);// tmp yt += xt (2)}add vec(pg, yt, tmp yt);// yt += tmp yt (3)}real t a;reduce vec(pg, a, yt);// summation in SIMD (4)#pragma omp barrier// summation over threads then over MPI processes
ThreadManager::reduce sum global(a, ith, nth);
return a;

場は 3 次元配列のような構造を持つ（1 次元配列 + index で疑似 3 次元配列）
そのため自然に 3 重ループになる(5) and (6)金森逸作: 第 20 回 HPC-Phys 勉強会, 2023 年 12 月 8 日 (online)

appendix-ii