CPS セミナー 2021 年 7 月 22 日

球面調和関数変換ライブラリの高速化について
石岡 圭一 (京大理・地惑)E-mail: ishioka@gfd-dennou.org

はじめに
地球流体の研究のために球面のスペクトル法をずっと使い続けてきた関係
で, それに必要な球面調和関数変換 (SHT = Spherical Harmonic Trans-form) のための数値ライブラリ (ISPACK) の開発を細々と趣味的に続けてきた (ispack-0.0 のリリースは 1997/04/15).
このセミナーでは, 球面のスペクトル法の基礎的な部分から導入し, SHT
ライブラリの高速化のためにどのようなノウハウが必要なのかを解説する.
また, 最近, ISPACK のバージョン 3.1.0 をリリースしたので, そこに含めた A64FX 向けチューニングの話も含める.

球面スペクトル法を用いた計算例:
回転球面上の 2 次元非圧縮流体の渦度方程式

∂ζ
∂ψ ∂ζ
∂ψ ∂ζ
∂ψ
+
−
+ 2Ω
= D[ζ]∂t∂λ ∂µ ∂µ ∂λ∂λζ ≡ ∇2ψ: 渦度,

ψ: 流線関数,Ω: 球の回転角速度,
λ: 経度,µ ≡ sin θ,

θ: 緯度,

t: 時刻,

D[ζ]: 高階粘性項.
球は単位時間あたり Ω/(2π) 回転する.
以下の数値計算の水平解像度は T682 (2048 × 1024 グリッド).

時間発展の例 1 (非回転の場合)(Ω = 0)(クリックして下さい. アニメーションの URL に飛びます)https://www.gfd-dennou.org/arch/ishioka/chikyu-no-butsuri/anim2000.mp4 (3.5MB)

渦度場 & 帯状平均帯状角運動量
赤: 正の渦度, 青: 負の渦度

時間発展の例 2 (高速回転の場合)(Ω = 2000)(クリックして下さい. アニメーションの URL に飛びます)https://www.gfd-dennou.org/arch/ishioka/chikyu-no-butsuri/anim2000.mp4 (5.7MB)

渦度場 & 帯状平均帯状角運動量
赤: 正の渦度, 青: 負の渦度

球面のスペクトルの応用例: 無限領域のスペクトル法
原点付近の運動が波動を励起し, それが遠方に伝播していくような問題は数多い. そのような問題を有限な領域で数値計算すると, 境界からの反射波の影響で原点付近の解が汚染される.
そこで, 通常は境界条件の設定を工夫したり (NRBC = Non-Reﬂecting

Boundary Condition), 波動を吸収するスポンジ層を導入したりして反射波を抑える.
分散性のある波動の場合, NRBC だけでは反射波を完全に抑えることはで
きず, またスポンジ層についてはその領域の広さや散逸の掛け方の設定が難しい.
そこで, 最初から領域が無限であることを考慮した計算ができないか?
無限領域のスペクトル法

アイディア: 立体射影により, 無限平面 (R2) を球面に射影し, 球面のスペクトル法を用いる.

*

*

計算例:

β 平面上の線形化された準地衡渦度方程式
∂ξ
∂ψ
+β
= 0.
∂t
∂x
ψ: 流線関数, ξ: ポテンシャル渦度,
ξ(x, y, t) = (△ − γ 2)ψ(x, y, t);

∂2
∂2
△=
+ 2.
2
∂x
∂y

γ: ロスビー変形半径の逆数, x: 東向き座標, y: 北向き座標, t: 時刻.
初期条件:




ξ = exp −x2 + y 2
2σ 2





.

ポテンシャル渦度場の時間発展 (Rossby 波の伝播を表す)
厳密解 (β = 1, γ = 0, σ = 0.2):

周期境界条件のフーリエスペクトル法による数値解 (領域: [4π × 4π]):無限領域のスペクトル法による数値解 (射影に使う球半径 = 2):

球面のスペクトル法のざっくりとしたレビュー
従属変数を球面調和関数で展開.
MX

f (λ, µ, t) =n
Xn=0 m=−nm(λ, µ).
am
(t)Y
n
n

λ: 経度, µ = sin θ, θ: 緯度, t: 時刻,
M : 切断波数.
Ynm(λ, µ) = Pn|m|(µ)eimλ.
Pnm(µ): Legendre 陪関数.

Pnm(µ) ≡

v
u
u
u
u
t

(n − m)! 1(2n + 1)(n + m)! 2nn!
n+md×(1 − µ2)m/2 n+m (µ2 − 1)ndµ(0 ≤ m ≤ n).

Legendre 陪関数による展開
球面調和関数展開では, 以下の Legendre 陪関数による展開に最も計算コストがかかる.
逆変換:

gjm =

M
X
n=m

m(µ )
sm
P
j
n n

(m = 0, · · · , M ; j = 1, · · · , J).

正変換:

sm
n =

J
X
j=1

wj gjmPnm(µj )(m = 0, · · · , M ; n = m, · · · , M ).

J: ガウス緯度の個数, µj : ガウス緯度, wj : ガウス重み.

必要な計算量

Legendre 陪関数の m ≥ 1 の成分に対しては実部・虚部を扱わねばならないことと, ルジャンドル陪関数の対称性から添字 j に関する計算は半分の J/2 でよいことを考慮すると, 逆変換・正変換に必要な計算量は N はル
ジャンドル陪関数をすべて事前に計算しておくとしても (乗算と加算が必要なことを考慮して)

J
N = · (M + 1)2 · 2 = J(M + 1)2.
2
さらに, 正変換時の数値計算で誤差が出ないことを保証するためには3 M としておく必要があるので,
J>2

3 3
N ∼ M .
2
M が大のときは, この変換計算をいかに効率化するかが全体の計算スピー
ドを決める.

x64 な CPU 上 (というか, ベクトル機でない普通の CPU) での最適化
今日, x64 な CPU が広く普及しているので, コードをそれに合わせてチューンしておくことは有意義.

x64 な CPU はメモリバンド幅が狭いので, キャッシュの活用が重要.
しかし, Legendre 陪関数変換を素朴に実装すると, 行列 × ベクトル型の演算になるので, キャッシュに載りにくい.

→ 解決策: Legendre 陪関数を変換と同時に計算する.

ただし, Legendre 陪関数を計算する計算量は変換そのものの計算量と同じオーダー. なので, Legendre 陪関数の計算コストを下げられる工夫が可能ならやるべき.

Legendre 陪関数は, 普通, 以下の漸化式で計算される.
m (µ) = (µP m(µ) − ϵmP m (µ))/ϵm .
Pn+1
n
n n−1
n+1
r

2 − m2)/(4n2 − 1).
ϵm
=
(n
n
m を事前に計算しておくと
この漸化式を 1 回計算するには, ϵm
および
1/ϵ
n
n

すれば, 乗算 3 回と減算 1 回.
変換の式および漸化式を工夫することによって, 漸化式部分のコストを削減した. 積和算命令にも馴染むものにした (1 段あたり実質積和算 1 回).

元々の漸化式
m (µ) + ϵmP m (µ)µPnm(µ) = ϵmPn n−1n+1 n+1

を 2 段階組み合わせると,m P m (µ)+{(ϵm )2+(ϵm)2}P m(µ)+ϵmϵm P mµ2Pnm(µ) = ϵmϵnnn n−1 n−2
n+1 n+2 n+2n+1となる. 赤道反対称の Legendre 陪関数について考えて,

n = m + 2l + 1 (l = 0, 1, . . .) とすると, 上の漸化式は,
mmmµ2Pm+2l+1(µ) =ϵmϵP
m+2l+2 m+2l+3 m+2l+3(µ)
2 + (ϵm2}P m+ {(ϵm))m+2l+2m+2l+1m+2l+1(µ)
m
m
+ ϵm
ϵ
P
m+2l+1 m+2l m+2l−1(µ),と表せる. ここで,mm(µ),
Pm+2l+1(µ) = µαm
p
l l

と pm
l を導入すれば,m
m pm (µ)m(µ) =ϵmαp
ϵµ2αmm+2l+2 m+2l+3 l+1 l+1l l
2 + (ϵm2}αmpm(µ)))+ {(ϵmm+2l+1l l
m+2l+2m
m pm (µ),+ ϵm
ϵαm+2l+1 m+2l l−1 l−1となる. この漸化式を簡単にするために, αml に対してm
m = −ϵm
m
m ,
ϵm
ϵ
α
ϵ
α
m+2l+2 m+2l+3 l+1
m+2l+1 m+2l l−1を満すことを課す. この式の両辺に αml を掛けて変形すると,mϵϵmm+2l+1 m+2lm=−
m αm ,
αm
α
α
l+1 l
l l−1
m
ϵm
ϵ
m+2l+3 m+2l+2を得る. この式を繰り返し適用すると,mϵm
ϵm+3 m+2m = (−1)lm αm ,
αm
α
α
1 0
l+1 l
m
ϵm
ϵ
m+2l+3 m+2l+2m を定めさえすれば, 順次, αm (l =従って, あとは αm
と
α
0
1
l

を得る.

2 3

)

.

m の定めかたには任意性があるが, 後々の便利のため, 以下のよう
αm
と
α
0
1

に定めることにする.

1

ϵm
m+1

m=
=
,
α
.
αm
1
0
m
m
m
ϵm+1
ϵm+2ϵm+3

こう定めると, 結局,
l
(−1)
m=
αm
.
α
l+1 l
m
m
ϵm+2l+3ϵm+2l+2
l αm を掛けて, 上の式となる. 先の pm(µ)
に対する漸化式の両辺に(−1)
llを用いると,pm
l+1(µ)


2µ2 − (−1)l (αm)2 (ϵm
2 + (ϵm2= (−1)l(αm)))
m+2l+2m+2l+1ll+ pm
l−1(µ).

となる.



pm
l (µ)

m を
さらに, amとblll (αm)2,
am
=
(−1)
l
l





2 ,
2 + (ϵm
m (ϵm
)
=
−a
)
bm
l
m+2l+1l
m+2l+2と導入すれば, 結局, pm
l (µ) に対する漸化式
mµ2 + bm)pm(µ) + pm (µ).
pm
(µ)
=
(a
l−1
l+1
l
l
l
m および µ2 を前もって計が得られる. この漸化式においては, amとbll

算しておけば, 漸化式一段あたり, 2 回の FMA(乗算 2 回, 加算 2 回) で計算できる. この漸化式は元の Pnm で考えれば, n を 2 つ増す漸化式に相当するので, n の一段あたり, 1 回の FMA で済んでいて, 元々の 3 回の FMAの必要性に比べて, 計算量が 1/3 になることになる.

ただし, この 2 段とばしの漸化式を Pnm(µ) の偶奇モード両方について計算していては全く意味が無い (むしろ以前のアルゴリズムの方が速いことになる) ので, どちらか一方についてのみ計算することになる. 実際には,µ = 0 付近の精度を考えて, 奇モードのみについて計算する.
隅モードについては,
m (µ) + ϵmP m (µ)µPnm(µ) = ϵmPn n−1n+1 n+1を用いて, 2 つの奇モードの重ね合わせ (の 1/µ) として計算できる.
このアルゴリズムを利用することの副産物として, 乗算の削減だけでなく,Legendre 陪関数に関するロードとストアが半分で済むので, その分のコスト削減にもつながる. 実際には, 例えば, 逆変換は以下のように計算さ
れる.



gjm =

M −m
2
X

l=0



m)pm(µ ) + µ
(tm
α
j
j
l l
l



M −m−12X

l=0m)pm(µ ).
(smαjm+2l+1 l
l

ここに,





M −m ] − 1)
m
m
m
m
s
ϵ
+
s
ϵ
(l
=
0,
1,
.
.
.
,
[
m+2l m+2l+1
m+2l+2 m+2l+2
2
tm
=
l

M −m ]).
m

 sm
ϵ
(l
=
[
m+2l m+2l+1
2

SIMD 命令の活用
x64 な CPU では, 演算性能向上のためには SIMD 命令 (AVX, AVX2,
AVX512) の活用が必須 (A64FX の話はまた後で)コンパイラが賢くなってはきているが, 少なくとも計算のホットスポット
はアセンブリでコーディングするのはかなり有効.

Intel の Haswell 以降では 1 サイクルあたり, 256 ビット長 (倍精度 4つ) の積和算が 2 組実行できる.

Intel の Skylake-SP 以降では 1 サイクルあたり, 512 ビット長 (倍精度8 つ) の積和算が 2 組実行できる.
ということで, 逆変換/正変換において, j によるループを SIMD 長ごとに分割して SIMD 命令を有効に使えるようにすればよい.

逆変換のホットスポットのコーディングの例 (SIMD 長が 8 の場合)
DO I=1,JB
DO J=1,8
Q1=Q(J,1,I)
Q2=Q(J,2,I)
Q0=Q(J,0,I)
Q3=Q(J,3,I)+S(1,2)*Q1+S(1,4)*Q2
Q4=Q(J,4,I)+S(1,1)*Q1+S(1,3)*Q2
Q5=Q(J,5,I)+S(2,2)*Q1+S(2,4)*Q2
Q6=Q(J,6,I)+S(2,1)*Q1+S(2,3)*Q2
Q1=Q1+(A(0)*Q0+A(1))*Q2
Q2=Q2+(A(2)*Q0+A(3))*Q1
Q(J,3,I)=Q3+S(1,6)*Q1+S(1,8)*Q2
Q(J,4,I)=Q4+S(1,5)*Q1+S(1,7)*Q2
Q(J,5,I)=Q5+S(2,6)*Q1+S(2,8)*Q2
Q(J,6,I)=Q6+S(2,5)*Q1+S(2,7)*Q2
Q1=Q1+(A(4)*Q0+A(5))*Q2
Q(J,1,I)=Q1
Q(J,2,I)=Q2+(A(6)*Q0+A(7))*Q1
END DOEND DO逆変換のホットスポットのコーディングの例 (AVX512 アセンブリ版)
L00:
vmovapd(%rax),%zmm14
vmovapd 64(%rax),%zmm12
vmovapd 128(%rax),%zmm13
vmovapd %zmm9,%zmm0
vfmadd231pd %zmm14,%zmm8,%zmm0
vfmadd213pd %zmm12,%zmm13,%zmm0
vmovapd %zmm11,%zmm1
vfmadd231pd %zmm14,%zmm10,%zmm1
vfmadd213pd %zmm13,%zmm0,%zmm1
vmovapd %zmm5,%zmm15
vfmadd231pd %zmm14,%zmm4,%zmm15
vfmadd213pd %zmm0,%zmm1,%zmm15
vfmadd213pd %zmm7,%zmm6,%zmm14
vfmadd213pd %zmm1,%zmm15,%zmm14
vmovapd %zmm15,64(%rax)
vmovapd %zmm14,128(%rax)
vmovapd 192(%rax),%zmm14
vfmadd231pd %zmm18,%zmm12,%zmm14
vfmadd231pd %zmm22,%zmm13,%zmm14
vfmadd231pd %zmm26,%zmm0,%zmm14
vfmadd231pd %zmm30,%zmm1,%zmm14

vmovapd %zmm14,192(%rax)
vmovapd 256(%rax),%zmm14
vfmadd231pd %zmm16,%zmm12,%zmm14
vfmadd231pd %zmm20,%zmm13,%zmm14
vfmadd231pd %zmm24,%zmm0,%zmm14
vfmadd231pd %zmm28,%zmm1,%zmm14
vmovapd %zmm14,256(%rax)
vmovapd 320(%rax),%zmm14
vfmadd231pd %zmm19,%zmm12,%zmm14
vfmadd231pd %zmm23,%zmm13,%zmm14
vfmadd231pd %zmm27,%zmm0,%zmm14
vfmadd231pd %zmm31,%zmm1,%zmm14
vmovapd %zmm14,320(%rax)
vmovapd 384(%rax),%zmm14
vfmadd231pd %zmm17,%zmm12,%zmm14
vfmadd231pd %zmm21,%zmm13,%zmm14
vfmadd231pd %zmm25,%zmm0,%zmm14
vfmadd231pd %zmm29,%zmm1,%zmm14
vmovapd %zmm14,384(%rax)
addq $448,%rax
cmpq %rax,%r10
jne L00

Skylake-SP では, 1 サイクルあたり, 積和算 2 回, ロード 2 回, ストア 1回ができる設計になっている. 上記のループは, 反復 1 回あたり, 積和算

24 回, ロード 7 回, ストア 6 回となっているので, ロードとストアがボトル
ネックになることなく (理論上は) 12 サイクルで実行できることになる.

論文:

Ishioka, K. (2018). A New Recurrence Formula for Eﬃcient
Computation of Spherical Harmonic Transform.

JMSJ, 96,

241–249.
※ Legendre 陪関数の漸化式にかかる計算コストを従来のものの 1/3 にする技法. この技法は, 他の開発者の球面調和関数変換ライブラリにも既に利用されている. cf. libsharp by Martin Reinecke( https://gitlab.mpcdf.mpg.de/mtr/libsharp/ ).
彼には, libsharp の中に “This update features signiﬁcant speedups

thanks to important algorithmic discoveries by Keiichi Ishioka(https://www.jstage.jst.go.jp/article/jmsj/96/2/96 2018-019/
and personal communication)”. と謝辞を書いてもらっている.

他の実装との速度比較
比較対象: SHTns
作者: Nathanaël Schaeﬀer (ISTerre, Université de Grenoble)

URL: http://users.isterre.fr/nschaeﬀ/SHTns/
論文: Schaeﬀer(2013): Eﬃcient spherical harmonic transforms

aimed at pseudospectral numerical simulations, Geochemistry,
Geophysics, Geosystems, vol. 14, pp.751–758.
開発履歴: Changelog を見ると, 2010 年に v1.0 がリリースされていたようだ. 現在は v3.4.6 (2021 年 5 月 15 日リリース)※実装の方向性は私のものに非常に近い (Legendre 陪関数を on-the-ﬂyで計算するとか, SIMD 命令の活用とか).

SHTns も, v3.2-r660 以降では, “New recurrence formula of Ishioka (2018) leading to faster transforms, especially for large
sizes. see https://doi.org/10.2151/jmsj.2018-019” として私の手法が取り入れられているので, それとの比較も入れる.

以前は, con-ﬁgure 時に –enable-ishioka というスィッチをオプションを付けてからmake する必要があったが, 今はそちらがデフォルトになっている. 従って, あえて –disable-ishioka というスィッチを有効にした場合とも比較しておく.

ベンチマーク結果 (AVX2)
プラットホーム:

CPU: Broadwell Xeon E5-2699v4 (22 コア)×2
コンパイラ: ifort 19.0.0.117 -xHost -qopenmp -align array64byte(SHTns は, gcc 6.3.0, option -O2 -march=native -ﬀast-math)
44 スレッドの場合のみを示す.

各スペクトルに (−1, 1) の一様乱数を与え, 逆変換 + 正変換を行った誤差評価 (SHTns’ は, SHTns を –disable-ishioka で make した場合).

ϵmax

ϵrms

ϵmax

ϵrms

SHTns’SHTns
ISPACK3
SHTns’SHTns
ISPACK3

M = 1023 M = 2047 M = 4095
1.2×10−10 5.9×10−11 2.0×10−9
1.2×10−12 4.4×10−12 2.0×10−11
5.1×10−13 1.2×10−12 5.8×10−12
7.7×10−13 2.6×10−13 6.7×10−12
7.4×10−14 1.7×10−13 3.9×10−13
4.3×10−14 8.9×10−14 1.9×10−13

SHTns’SNTns
ISPACK3
SHTns’SHTns
ISPACK3

M = 8191 M = 16383
4.9×10−10 2.5×10−8
4.3×10−11 5.9×10−11
1.8×10−11 4.1×10−11
1.5×10−12 4.7×10−11
7.2×10−13 1.2×10−12
4.3×10−13 7.9×10−13

スペクトル⇔格子点値 1 回あたりの変換にかかる時間 (sec)

M = 1023
M = 2047
bwd
fwd
bwd
fwd
SHTns’
0.0021 0.0022 0.014 0.015
SHTns
0.0015 0.0015 0.0097 0.0098
ISPACK3 0.0017 0.0014 0.011 0.0094
M = 4095
M = 8191
bwd fwd
bwd fwd
SHTns’
0.11 0.11 0.86 0.86
SHTns
0.070 0.071 0.53 0.52
ISPACK3 0.075 0.064 0.53 0.45

M = 16383
bwd fwd
7.4 7.3
4.5 4.3
4.0 3.4

スペクトル⇔格子点値 の変換における実効 GFlops 値 (このシステムにおける理論ピーク性能は, AVX 使用時の CPU 周波数を 1.8GHz とすると,

1.8 × 16 × 44 = 1267GFlops. Turbo が効いて, 2.6GHz まで上がるとするなら, 2.6 × 16 × 44 = 1830GFlops. MKL の DGEMM で4000 × 4000 と 4000 × 2000 の行列同士の行列積を計算させて計測すると, 1230GFlops).

M = 1023 M = 2047
bwd fwd bwd fwd
ISPACK3 651 823 813 940
M = 4095 M = 8191 M = 16383
bwd fwd bwd fwd bwd fwd
ISPACK3 933 1090 1052 1245 1111 1311

ベンチマーク結果その 2(AVX512)
プラットホーム:

CPU: Skylake-W Xeon W-2135 (3.7GHz, 6 コア)
コンパイラ: ifort 19.1.3.304 -xHost -qopenmp -align array64byte(SHTns は, gcc 8.3.0, option -O2 -march=native -ﬀast-math)
6 スレッドの場合のみを示す.
スペクトル⇔格子点値 1 回あたりの変換にかかる時間 (sec)

M = 1023
M = 2047
bwd
fwd
bwd fwd
SHTns
0.0044 0.0044 0.033 0.033
ISPACK3 0.0048 0.0045 0.029 0.027
M = 4095 M = 8191
bwd fwd bwd fwd
SHTns
0.20 0.19 1.6 1.5
ISPACK3 0.18 0.18 1.3 1.3

M = 16383
bwd fwd
11.0 10.9
9.9 9.6

スペクトル⇔格子点値 の変換における実効 GFlops 値 (このシステムにおける理論ピーク性能は, AVX512 使用時の CPU 周波数を 3.9GHz と仮定すると, 3.9 × 32 × 6 = 749GFlops. MKL の DGEMM で 8000 ×8000 の行列同士の行列積を計算させて計測すると, 660GFlops 程度).
M = 1023 M = 2047
bwd fwd bwd fwd
ISPACK3 235 253 310 326
M = 4095 M = 8191 M = 16383
bwd fwd bwd fwd
bwd fwd
ISPACK3 383 383 426 429
445 458

ベンチマーク結果その 3(AVX512, Thanks to 竹広さん)
プラットホーム (阪大 SQUID):

CPU: Icelake Xeon Platinum 8368 (2.4GHz, 38 コア)×2 76 ス
レッドの場合のみを示す.
スペクトル⇔格子点値 の変換における実効 GFlops 値 (このシステムにおける理論ピーク性能は, AVX512 使用時の CPU 周波数が仮に 2.4GHzのままととするなら, 2.4 × 32 × 76 = 5837GFlops)

M = 4095 M = 8191 M = 16383
bwd fwd bwd fwd bwd fwd
ISPACK3 2590 2493 3186 3181 3785 3569

A64FX でのベンチマーク結果その 1(Thanks to 齋藤泉君)
ispack-3.0.1 の状態で, 単純な Fortran90 コードの状態で A64FX(名大不老) 上でとりあえず 1 ノード (48 コア) でベンチマークをとってもらった
とき.

M = 8191
bwd fwd
ISPACK-3.0.1 1206 305
不老では A64FX が 2.2GHz で運用されているので, 理論ピーク性能は,2.2 × 32 × 48 = 3379GFlops) である筈なので, 逆変換の方で 36%程
度は出ているが, 正変換の方は 9%程度と悲惨. ということで少しは何とかしたかった.

A64FX のためのチューニングに向けてということで, x86 上と同様に, まず A64FX 上のアセンブリから勉強を始めて手でチューンすることを考えた. 以下, やっていったことを時系列
てきに散漫に書く.

A64FX 上で使えるアセンブリ (AArch64 + SVE) はとっても書き易い.
屋上屋を重ねた Intel x64 用のアセンブリ + avx512 とかに比べると,はるかに分かりやすい. また, SVE の movprtfx で実質的に 4 オペラン
ド命令が使えるのはとっても楽 (Intel にも欲しい..). 例えば,

movprfx z4.d, p0/z, z0.d
fmla z4.d, p0/m, z1.d, z2.dのように書けば, 実質的に

z4 = z0 + z1 * z2のような計算が 1 命令で実行できる.

ということで, 最初は A64FX 上でアセンブリを書き始めたのだが, A64FXはやはり命令体系は ARM アーキテクチャになったとはいえ, ハード的には富士通の石なので, 京 (SPARC64VIIIfx) の頃の残念な性質をかなり受け継いでしまっていて, よほど条件が整わないとピーク性能が出ない. 理由は:・命令のレイテンシが大きすぎる: 積和算のレイテンシが 9 サイクル, L1
からの vector load のレイテンシが 11 サイクルもあるっていったい何?
ちなみに Intel の Skylake では, それそれ, 4 サイクルと 3 サイクル.
・物理レジスタが少ない: A64FX では, 物理レジスタが 32+96 で 128しかない. ちなみに Intel の Skylake では 168ということで, よほど条件が整わない限り (例えば, DGEMM のように,

dependency chain が深くない計算を同時並列的にできるような場合以外), レイテンシの隠蔽がちゃんとできない.

ということで, A64FX 上でアセンブリを書く場合には, 命令のレイテン
シや物理レジスタの枯渇まで気にして命令の順番を考える必要があり, こ
れは A64FX の癖を知りつくしていないとできないので, 結構, 人間業を超えてしまう.
ところが, A64FX 上の Fortran コンパイラ, frtpx「は」かなり賢く, 元になる Fortran のプログラムが SIMD 化しやすく書いてあれば, かなり最適なアセンブリコードを吐いてくれる. それを眺めてみると, 「あー, そ
こで load を先読みするのね」とか, 「あー, そこでこの積和算やっとく方がレイテンシの関係で得なのね」とまあ解釈はできるが, なかなか人間業
では思いつかないような順番で命令を発行している.
ということで, 自分で A64FX 用のアセンブリを書くのは放棄して, 正変換のコア部分を frtpx が SIMD 化しやすいように, frtpx の気持になって書き換えてみた.

チューニングのために書き換えたところの例 (正変換の元の形)
DO I=1,JB
DO J=1,8Q0=Q(J,0,I)
Q1=Q(J,1,I)
Q2=Q(J,2,I)
Q3=Q(J,3,I)
Q4=Q(J,4,I)
S(1)=S(1)+Q4*Q1
S(2)=S(2)+Q3*Q1
S(3)=S(3)+Q4*Q2
S(4)=S(4)+Q3*Q2
Q1=Q1+(AC(0)*Q0+AC(1))*Q2
Q2=Q2+(AC(2)*Q0+AC(3))*Q1
S(5)=S(5)+Q4*Q1
S(6)=S(6)+Q3*Q1
S(7)=S(7)+Q4*Q2
S(8)=S(8)+Q3*Q2
Q1=Q1+(AC(4)*Q0+AC(5))*Q2
Q(J,1,I)=Q1
Q(J,2,I)=Q2+(AC(6)*Q0+AC(7))*Q1
END DOEND DOチューニングのために書き換えたところの例 (正変換の新しい形)
SV=0
DO I=1,JB
Q0V=Q(:,0,I)
Q1V=Q(:,1,I)
Q2V=Q(:,2,I)
Q3V=Q(:,3,I)
Q4V=Q(:,4,I)
SV(:,1)=SV(:,1)+Q4V*Q1V
SV(:,2)=SV(:,2)+Q3V*Q1V
SV(:,3)=SV(:,3)+Q4V*Q2V
SV(:,4)=SV(:,4)+Q3V*Q2V
Q1V=Q1V+(A(0)*Q0V+A(1))*Q2V
Q2V=Q2V+(A(2)*Q0V+A(3))*Q1V
SV(:,5)=SV(:,5)+Q4V*Q1V
SV(:,6)=SV(:,6)+Q3V*Q1V
SV(:,7)=SV(:,7)+Q4V*Q2V
SV(:,8)=SV(:,8)+Q3V*Q2V
Q1V=Q1V+(A(4)*Q0V+A(5))*Q2V
Q(:,1,I)=Q1V
Q(:,2,I)=Q2V+(A(6)*Q0V+A(7))*Q1V
END DO
DO K=1,8
S(K)=S(K)+SUM(SV(:,K))END DOA64FX でのベンチマーク結果その 2
ispack-3.1.0 で A64FX(名大不老) 上でとりあえず 1 ノード (48 コア) で
ベンチマークをとった.

M = 8191 M = 16383
bwd fwd bwd fwd
ISPACK-3.1.0 1205 1294 1316 1413ということで, 今回のチューニングによって, 正変換のスピードが逆変換並みになった. 正変換の方が若干速いのは, x86 上と同様. M = 16383のときの正変換で, ピークの 42%程度の実効性能は出ていることになる.

ギネスに挑戦?
私の知る限り, これまでに行われた自由度最大の球面調和関数変換は,Reinecke, M., & Seljebotn, D. S. (2013).
Libsharp–spherical harmonic transforms revisited..
Astronomy & Astrophysics, 554, A112.
に記述されている TL262143 (= 218 − 1) のケースだと思われる (ここでは, Intel E5-2670 (Sandy Bridge, 8 コア, 2.6GHz) を 4096 コア
分使って計算が行われている)ギネスに挑戦ではないが, 富岳を使ってこの記録を破ることを試みた(Thanks to 竹広さん).

なお, ISPACK-3.1.0 では, 球面調和関数変換を MPI 化したものも実装してある (sypack) のでそれを用いる.

MPI 実装は素朴なもので, Legendre 陪関数変換は経度方向波数 m の方向に分割して並列化しており, 経度方向の FFT を実行する際は, 各緯度円毎に並列化する. 両者をつなぐ際, 転置的な処理が必要になるが, そこは単純に MPI ALLTOALL を 1 回呼ぶだけの単純な実装 (特に通信時間の隠蔽とかには凝っていない). また, MPI+OpenMP のハイブリッド並列化をしている.

A64FX でのベンチマーク結果その 3(Thanks to 竹広さん)
ispack-3.1.0 で A64FX(富岳) 上で世界タイ記録およびその倍の切断波数でのベンチマークをとった (なお, 富岳では 2.0GHz のモードで計算)

TL262143 の計算 (256 ノード使用, 48 スレッド/ノード)
bwd: 64.3sec (実効 280TFlops)
fwd: 56.2sec (実効 321TFlops)
ϵmax = 4.71 × 10−9, ϵrms = 2.47 × 10−11
TL524287 の計算 (1024 ノード使用, 48 スレッド/ノード)
bwd: 123sec (実効 1172TFlops)
fwd: 107sec (実効 1350TFlops)
ϵmax = 1.71 × 10−8, ϵrms = 7.70 × 10−11

ということで, TL524287 の方の正変換では, 1 ノードあたり, 実効

1.32TFlops でていることになり, これは 2.0GHz 駆動の A64FX のピーク性能の 43%程はでていることになる.
また, このような大きな切断波数の計算では, 逆変換 + 正変換での相対誤差が 10−8 のレベルになりさすがに大きくなるが, まあ, 単精度よりはましな程度の精度は維持できていると見なせる.

A64FX でのベンチマーク結果その 4(Thanks to 竹広さん)
TL262143 の計算での不老と富岳の比較 (どちらも 2.2GHz 駆動モード
での計算, 256 ノード)
不老:

bwd: 55.1sec (実効 327TFlops)
fwd: 44.4sec (実効 406TFlops)
富岳:

bwd: 55.9sec (実効 323TFlops)
fwd: 49.1sec (実効 367TFlops)
富岳の方がちょっと遅い. 最初は MPI 通信に要する時間の差とかがあるのかと思ったが, 不老でも富岳でも, 正変換の場合に MPI 通信に要している
時間はどちらも 3 秒程度なので, この差が説明できない. この点については同じバイナリを使ってみるなどいろいろチェックしたがちょっと謎のまま.

まとめと TODO
まとめ

• Legendre 陪関数の漸化式の工夫など様々な技法をもりこんで, (少なくとも Intel の CPU 上である程度以上のサイズの変換では) 世界最速の

SHT ライブラリを構築した.
• ここで導入した Legendre 陪関数の漸化式の工夫については, 他の開発者の実装にも既に取り入れられてきている.

• A64FX のような CPU でもある程度の性能は出るようにチューニングを
行った.

TODO
• 経度方向の回転対称性を仮定した変換の実装• FFT について, COS, SIN 変換や, 2 次元, 3 次元の FFT の実装• Python インターフェースの整備