GFD セミナー 2009 分科会 2009 年 9 月 6 日

SJPACK について
石岡圭一
E-mail: ishioka@gfd-dennou.org

まずは球面のスペクトル法についてのおさらい

偏微分方程式中の従属変数を球面調和関数で展開.
MXf (λ, µ, t) =nX

n=0 m=−nm(λ, µ)am(t)Ynn

ここに, λ: 経度, µ = sin θ, θ: 緯度, t: 時刻, M : 切断波数.

Ynm(λ, µ) = Pn|m|(µ)eimλ.
Pnm(µ) は Legendre 陪関数.

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

ルジャンドル陪関数による展開

球面スペクトル法には球面調和関数変換が必要. そのうち計算量の大部分を占めるのはルジャンドル陪関数に関する以下の変換である:

逆変換

gjm =MX
n=mm(µ )sm
Pjn n(m = 0, · · · , M ; j = 1, · · · , J)

正変換

smn =JX
j=1

wj gjmPnm(µj )(m = 0, · · · , M ; n = m, · · · , M )ここに, M : 切断波数, J: ガウス緯度の個数, Pnm(µ): ルジャンドル陪関数, µj : ガウス緯度, wj : ガウス重み.

必要な計算量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 が大のときは, この変換計算をいかに効率化するかが全体の計算スピードを決める.

ISPACK における球面調和関数変換のこれまでの実装
stpack(1995 年) 最も素朴な実装
smpack(1998 年) ベクトル化を追求
snpack(1999 年) ベクトル化を追求しつつ省メモリ化
sjpack(2009 年) スカラー計算機での高速化

変換のベクトル化 (snpack での話)
ベクトル計算機の能力を発揮させるためには, ベクトル長 (最内側 DO-LOOP の長さ) を十分に長くとる必要がある.
最も素直にプログラミングすると (stpack), 逆変換・正変換の計算はともに最内側 DO-LOOP の DO 変数は ガウス緯度の添字 j にとるのが自然.
しかし, この方針だと効率面で 2 つ問題がある:

1. ベクトル長は J2 = 34 M になるが, M が十分大でなければベクトル長が短くなって VP の能力が十分発揮されない. (M = 170 でもベクトル長= 128)
2. DO 変数を j にとると, 正変換の方の計算が内積型になり, 逆変換の方の計算に較べて計算スピードが落ちる (約半分程度).

→ breakthrough: DO 変数を m にとる.

そうすると, まず上記 2 の問題はクリアされる. また, ベクトル長は M とより長くなり, さらに複数種類の変数を同時に変換する場合には, ベクト
ル長は M ·(変数の種類) となり, 一気にベクトル長が長くとれるようになる.

ただし, m を DO 変数にとるには工夫が必要 (そのままのコーディングでは, 各 n に対するベクトル長は n と可変になってしまう).

図のように三角切断の一部を折り返すことによって m に関する DO 文のベクトル長を均一に M にする.

nMnM

MmMmスカラー計算機での高速化に向けて (sjpack)

ここまではベクトル機がまだ元気だった頃の話.

周知のように, 京速計算機プロジェクトから NEC・日立連合が降りたことに象徴されるように, もはやベクトル機には将来性が無い.ということで, スカラー計算機上での高速化を優先して設計していくべき(といっても, 最近の GPGPU 計算においてはベクトル計算機向けに書かれたコードが移植しやすいという話もあるが).

スカラー計算機は, ベクトル計算機に比べてメモリアクセスが遅いので, いかにキャッシュメモリの中で計算を閉じさせるかということが高速化の肝.

といっても, 素朴に実装しようとすると, 球面調和関数変換はキャッシュが効きにくい.

なぜなら, 変換に必要な球面調和関数のデータは倍精度ならJ
8 バイト × · (M + 1)(M + 2)/2 = 2J(M + 1)(M + 2) バイト.
2
例えば, M = 170, J = 256 とすると, これは 15MB. 最近の Xeon のL3 キャッシュにすら載らない.
ではどうするか? → 変換する度に Legendre 陪関数を毎回計算すればいい. そうすれば, 8 バイト × J/2 × 7 = 28J バイトの領域があればその中で計算を繰り返すことができる. J = 256 のときはこれは 7168 バイ
ト. これなら Pentium4 時代の Xeon でも L1D キャッシュに収まる.

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

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

ここに, ²mn =r(n2 − m2)/(4n2 − 1). この漸化式を 1 回計算するにm を事前に計算しておくとすれば, 乗算 3 回と減算 1 回.
は, ²mおよび1/²nn

これ以上減らしようが無いようにも見えるが,m(µ)
Pnm(µ) = αmpn nm を新に導入し, αm をうまく選ぶことによって,と αm,pnnnmµpm(µ) + pm (µ)pm(µ)=βnnn−1n+1とすることができる. こうすれば乗算 2 回と減算 1 回となり, 乗算を 1 回削減できる (このテクニックは snpack の段階で既に導入してあったが).

さらに, Pentium4 以降の x86 系プロセッサで高速化を目指すには SSE2(Streaming SIMD Extension 2) 命令を最大限活用する必要がある.
最近ではコンパイラもそれなりに賢くなってきてはいるが, SSE2 命令を活用しようとするならば, アセンブリ言語でのプログラミングは避けて通れない.ということで, Legendre 陪関数変換のコアの部分はアセンブリ言語で書いて最適化しておくべき.

実際コアとなる部分はどんなものか?
例: jlswg-fort.f

SUBROUTINE LJLSWG(JH,S,R,Y,QA,QB,W)
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION W(JH,2),Y(JH),S(2)DIMENSION QA(JH),QB(JH)
DO J=1,JH
W(J,1)=W(J,1)+S(1)*QA(J)W(J,2)=W(J,2)+S(2)*QA(J)QB(J)=QB(J)+R*Y(J)*QA(J)END DO
ENDと短いサブルーチンなので, これを (SSE2 命令等を使って) ガシガシに最適化すればいい (というか, このコアの部分をサブルーチンとして独立させてあるところが設計の肝の一つ).

SSE2 命令を使って最適化した例 (64 ビット版)
例: jlswg-sse64.f.globl ljlswg_
ljlswg_:movl(%rdi), %edi
movhpd (%rsi), %xmm0
movlpd (%rsi), %xmm0
movhpd 8(%rsi), %xmm6
movlpd 8(%rsi), %xmm6
movhpd (%rdx), %xmm1
movlpd (%rdx), %xmm1
movq 8(%rsp), %r10
shlq $3,%rdi
movq %rcx,%rsi
addq %rdi,%rsi
movq %r10,%r11addq %rdi, %r11.align 16.L0:
movaps (%rcx), %xmm4
movaps (%r8), %xmm2
movaps (%r9), %xmm3
movaps (%r10), %xmm5

movaps (%r11), %xmm8
movaps %xmm0, %xmm7
mulpd %xmm1,%xmm4
mulpd %xmm2,%xmm4
addpd %xmm4,%xmm3
mulpd %xmm2,%xmm7
addpd %xmm7,%xmm5
mulpd %xmm6,%xmm2
addpd %xmm2,%xmm8
movaps %xmm3,(%r9)
movaps %xmm5,(%r10)
movaps %xmm8,(%r11)
addq $16,%rcx
addq $16,%r8
addq $16,%r9
addq $16,%r10addq $16,%r11cmpq %rcx,%rsi
jne .L0retベンチマーク結果
とりあえず T170 の普通の設定 (J = 256, I = 512) での snpack との速度比較 (なお, sjpack では FFT として FFTJ を使っている).
なお, Flops は実際の計算量 (Legendre 陪関数の毎回計算を含む) ではな
く, Legendre 陪関数をもし事前に計算してメモリに格納していたとする
場合の変換そのものに絶対必要な計算量で算出している).

Intel Core2 Duo Xeon 3.0GHz (Woodcrest) 64bit mode
(L1: 32KB, L2: 4MB) においてSJPACK: 3.5 GFlops
SNPACK: 1.1 GFlops
Intel Pentium4 Xeon 3.06GHz (Prestonia)(L1: 8KB, L2: 512KB) においてSJPACK: 1.3 GFlops
SNPACK: 0.17 GFlops

FLTSS との比較
東大理・情報理工の須田礼仁氏の FLTSS

http://www.na.cse.nagoya-u.ac.jp/~reiji/fltss/との比較. ただし, FLTSS は FFT まで面倒を見ないので, Legendre
陪関数変換の部分のみで比較 (ifort -O3 でコンパイル).

T170 の普通の設定 (J = 256) での比較
Intel Core2 Duo Xeon 3.0GHz (Woodcrest) 64bit mode
(L1: 32KB, L2: 4MB) において
LJPACK: 3.1 GFlops
FLTSS: 1.3 GFlops
Intel Pentium4 Xeon 3.06GHz (Prestonia)(L1: 8KB, L2: 512KB) において
LJPACK: 1.4 GFlops
FLTSS: 0.73 GFlops

T341 の普通の設定 (J = 512) での比較
Intel Core2 Duo Xeon 3.0GHz (Woodcrest) 64bit mode
(L1: 32KB, L2: 4MB) において
LJPACK: 3.1 GFlops
FLTSS: 1.6 GFlops
Intel Pentium4 Xeon 3.06GHz (Prestonia)(L1: 8KB, L2: 512KB) において
LJPACK: 1.2 GFlops
FLTSS: 0.92 GFlops

課題・OpenMP による並列化はしてあり, Xeon ではそれなりに並列化効率が出ているが, Opteron では今一つ. さらに SPARC64VII では悲惨なよう (Thanks to 竹広さん). → 改善の余地あり?
そもそも, SPARC64VII ではアセンブリレベルの最適化の余地もある筈.
ただ, 私はアカウント持ってないので手が出ないが.

・GPGPU や Larrabee の動向は? SSE2 な最適化だけで「あと 10 年は戦える」か? Intel は将来的には Larrabee も CPU に内蔵して命令セットも AVX というのに整理していくらしい.

ネット上の情報源- IA32 インテルアーキテクチャソフトウェアデベロッパーズマニュアル(ftp://download.intel.co.jp/jp/developer/jpdoc/ )- Agner Fog 氏のページ (x86 アセンブリでの高速化についての情報源)(http://www.agner.org/optimize/)- IA32 SIMD の扉 (x86 の SIMD 命令についてよくまとまっている)(http://www.icnet.ne.jp/˜nsystem/simd tobira/)「ネットは広大だわ……」(by 草薙素子, 1989 年)

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

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

→ 無限領域のスペクトル法

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

**

計算例:

β 平面上の線形化された準地衡渦度方程式∂ξ∂ψ
+β= 0.
∂t∂x
ψ: 流線関数, ξ: ポテンシャル渦度,
ξ(x, y, t) = (4 − γ 2)ψ(x, y, t);∂2∂2
4=+ 2.
2∂x∂yγ: ロスビー変形半径の逆数, x: 東向き座標, y: 北向き座標, t: 時刻.
初期条件:



ξ = exp −x2 + y 22σ 2.

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

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