数値計算と性能,精度に関するまとめ
1.はじめに数値計算の性能,精度に関して以下の事をまとめました。
(1)SR16000のsmpとmpiの効果を対称行列のノンゼロ要素のみを使用するCG法でその性能を検証しました。
(2)ループ積分を用いて,アクセラレータ及びサーバーの性能を比較しました。
(3)サーバーにおいてRump’s例題でDD形式とIEEE形式の性能を比較しました。
(4)サーバーにおいて3次元ループ積分でDD形式とIEEE形式の性能を比較しました。

尚,チューニングの方法や精度に関する注意事項は適宜記載しています。

2 SR16000でのsmpとmpiに関して
SR16000では、1coreに関して2組の乗加算器があり一般的にsimd命令が適用される場合はsmt=off,適用されない場合はsmt=ｏｎの性能がよくなります。1ノードでの実行では以下の様になります。
Simd命令が適用されない
setenv HF_BINDPROC_STRIDE 1./cg_omp -F'prunst(threadnum(64),BIND(0))'

Simd命令が適用される
setenv HF_BINDPROC_STRIDE 2./cg_omp -F'prunst(threadnum(32),BIND(1))'
またSR16000の構成では最大8ノードで,1ノード32coreのため,1ノード当たりのコア数がノード数より多いのでハイブリッド並列の効果に特徴があります。参考までにBG/Lは1ノード2core,BG/Qは1ノード16coreとノード数が1ノード当たりのコア数より多くなります。

1ノードでのmpi数とsmp数
Simd命令が適用されるのでmpi数*smp数=32
node =1
total_tasks = 16
poe mpibind -thread 2 -hpc ./cg_mpi5
SR16000

cg法

N  2016000
スレッド数  int(32 / mpi数)
1mpi,1スレッド 1123.0154秒

実行時間(秒)
1008060

実行時間(秒)

4020

mpi数

0

1 2 3 4 5 6 7 8 16
注 mpi数=1 smp並列

1,4,8ノードでの性能
SR16000 cg法
N  10000000
1node,1mpi,1スレッド 11547.3647秒
case1 1node,2mpi,16smp
case2 1node,4mpi,8smp
case3 1node,8mpi,4smp
case4 1node,16mpi,2smp
case5 1node,32mpi,1smp
case6 4node,4mpi,32smp
case7 8node,8mpi,32smp
実行時間(秒)
1200
1000
800

実行時間(秒)

600
400
200

case

0
1 2 3 4 5 6 7

SR16000 4ノード
N=100000000
Mpi数*smp数=128

node =4
total_tasks =16
poe mpibind -thread 8 -hpc ./cg_mpi5

実行時間(秒)7500.0000
7000.0000

6500.0000

実行時間(秒)6000.0000
5500.0000
5000.0000

48 16 32 64
mpi数

SR16000 8ノード
N=100000000
Mpi数*smp数=256

node =8
total_tasks =32
poe mpibind -thread 8 -hpc ./cg_mpi5

実行時間(秒)
7000.0000
6000.0000
5000.0000
4000.0000
3000.0000
2000.0000
1000.0000
0.0000

実行時間(秒)

81632

mpi数

64

サイズNに対する反復回数
サイスNが大きくなると反復回数は増加する。
サイズNが大きくなると反復回数/Nは小さくなる。
log(反復回数)
1412
10
8
6
4
20

log(反復回数)

反復回数/N
0.1400000
0.1200000
0.1000000
0.0800000
0.0600000
0.0400000
0.0200000
0.0000000
10.82
11.5113.82
14.52
16.12
18.42

反復回数/N高エネルギー加速器研究機構のSR16000の構成は1ノードあたり32coreとなっていますが、32coreでメモリを蜜結合で共有するのではなく,4つの8coreで蜜結合が疎結合でメモリを共有する形となっています。このため、1ノード32coreではなく,4ノード8coreと考えた方が良く1ノードでは32smp並列より,4mpi,8smpのハイブリッド並列の方が性能がでます。4ノード,8ノードの場合も同様に8core16ノード,32ノード考えて実行すると効果がでます。
またsimd命令が適用できないプログラムではSmt=onとして,2smp/coreとして実行すると性能がでます。

3. ループ積分を用いて,アクセラレータ及びサーバーの性能比較

次の３つのループ積分を用いてアクセラレータ及びサーバーの性能比較をしています。

4次元積分

S221
1 1 x 1 x  y 1 x  y  z

S (s; m1 , m 2 , m 3 , m 4 , m 5 )   221

2222

20 0

 00

1dudzdydx
DCC  ( x  y  z  u )(1  x  y  z  u )  ( x  y )( z  u )E  (1  x  y  z  u )( x  z )( y  u )  ( x  y ) zu  ( z  u ) xy
2222M 2  xm1  ym2  zm3  um4  (1  x  y  z  u )m5
D   sE  M 2Cs 221(1 : 100,100,0,0,100)
解析近似解 : 0.0380004438127

26次元積分 Laporta I1 1 x1 1 x1  x2 1 x1  x2  x3 1 x1  x2  x3  x4 1 x1  x2  x3  x4  x5

I  0 0


0


0


0


0

C
dx6 dx5 dx4 dx3dx2 dx1
3
Dx7  1  x1  x2  x3  x4  x5  x6
C=(x1+x2+x3+x4+x5)*(x1+x2+x3+x6+x7)-(x1+x2+x3)**2

cc=x1*m12+x2*m22+x3*m32+x4*m42+x5*m52+x6*m62+x7*m72
D=-c*cc.+s*(x1*x2*(x4+x5+x6+x7)+x1*x5*x6+x2*x4*x7-x3*x4*x6).+t*x3*(-x4*x6+x5*x7).+p12*(x1*x3*(x4+x5+x6+x7)+x3*x4*(x6+x7)).+p22*(x2*x3*(x4+x5+x6+x7)+x3*x6*(x4+x5)).+p32*(x4*x5*(x1+x2+x3+x6+x7)+x4*x6*(x2+x3)+x1*x5*x7).+p42*(x6*x7*(x1+x2+x3+x4+x5)+x4*x6*(x1+x3)+x2*x5*x7)
m12=1.0d0, m22=1.0d0, m32=1.0d0,m42=1.0d0
m52=1.0d0, m62=1.0d0, m72=1.0d0
p12=1.0d0, p22=1.0d0, p32=1.0d0,p42=1.0d0s =1.0d0
t =1.0d0
解析近似解=0.0853513981538

8次元積分M441d
C(DiC)
I  {( x 1 , x 2 , x 3 , x 4 , x 5 , x 6 , x 7 , x 8 , x 9 ) |x 1  0, x 2  0, x 3  0, x 4  0,x 5  0, x 6  0, x 7  0, x 8  0, x 9  0,x 1  x 2  x 3  x 4  x 5  x 6  x 7  x 8  x 9  1}
x i , j, k ,....  x i  x j  x k  .....
2C  x 167 x 278 x 349 x 589  x 349 x 589 x 7  x 167 x 349 x 822 x 167 x 278 x 9  x 7 x 922222D   x 167 x 278 x 3 x 589  x 167 x 2 x 349 x 589  x 1 x 278 x 349 x 5892 x 123 x 167 x 278 x 349 x 589  2 x 1 x 2 x 349 x 589 x 7  x 3 x 589 x 72222

 x 123 x 349 x 589 x 7  x 167 x 3 x 8  x 1 x 349 x 822 x 123 x 167 x 349 x 8  2 x 167 x 2 x 3 x 8 x 9  2 x 1 x 3 x 7 x 8 x 92222

 x 167 x 2 x 9  x 1 x 278 x 9  x 123 x 167 x 278 x 922

 2 x 1 x 2 x 7 x 9  x 123 x 7 x 9
解析解 2441 7
 55.58525391567 84822

8次元積分 M 44実行結果一覧表
N  16, 演算量  850.4GFLOP

実行時間(秒),性能(GFLOPs)一覧表
CPUsmp数
実行時間
性能
備考
Pezy
8192 3.222556
264 HPC152
8192 2.438170
349 H27末
8192 1.298701
655 最終Phi5110P
240 1.838314
463 ーＯ２E5-2670
16 6.325043
134 ーＯ２HD58703200 1.542792
551 標準gmp
16 9.078043
94 ーＯ１
16 7.021060
121 ーＯ２
g9mpx
16 10.650697
80 ーＯ１
16 8.946464
95 ーＯ２
g9mpx340 4.800980
177 ーＯ１40 3.311994
257 ーＯ２
Pezy1500GFLOPs
Phi5110P 1011GFLOPs
HD5870 1088GFLOPs
gmp
g9mpx
g9mpx3

E5-2690(SB) 2.90GHz
E5-2687W v2 3.40GHz
E5-2687W v3 3.10GHz

E5-2670

2.60GHz

性能緒元一覧

アクセラレータによる性能比較
ループ積分性能比較 実行時間(秒)一覧表
プログラム Ns221
s221
laporta ilaporta i
m44m44

演算量 Pezy-SC
Phi5110P
HD5870(GFLOP) 1500GFLOPs 1011GFLOPs 1088GFLOPs
512
3711
6.911009
14.119193
6.374162
102459374 110.521333 194.452293 135.027731
6410308 16.988779
20.838768
14.426863
128 659707 930.784595 1515.707021 904.091546
16850
1.298701
1.838314
1.542792
32 217703 331.900691 298.712631 330.427198

発生した問題点
1.m44 HD5870 N=32 でゼロ割発生
2. m44 Pezy-SC 除算の性能が悪い。
（＜＝今後改善される模様)

対策
HD5870 N=16,Pezy-SC では除算が遅い
が結果は正しい計算をしている事に着目=>二分法の考え方を使用する。

HD5870のN=32での対策
1Wd  の CD
CD
計算をしている。二重 指数関数型積分ではM 44の計算を模式化すると 0  CD  1では0  W  CD  1となっている。
Wの計算で桁落ちなどに よりゼロ割が発生する。こ のため,  0を定め,    0
  CDとして, なら    0として計算する。最終 次元での式は(log(x )) 61 x 26dx

x(log(x))(log()  201 xx1 x 2となり, x  e 7 , WM   x (log(x )) 6 (log()  2x
 8231 0  CD  x 2とすると , 0  10 40 , x  10 20 , W  4.40  10 9 0  10 50 , x  10 25 , W  2.10  10 13より,  0  10 50 とすれば倍精度演算の 範囲では最終結果に影響がなく , 値は  0の減少関数で、  0   0で桁落ちで となるので二分法で  0の値を定める事が出来 る。

Pezy-SCの除算の高速化（今後はこの修正が必要なくなる)
1
しかx高速に実行されない。( x  2 149  1.4  10  45 )現時点では Pezy  SCでは単精度のこのため  0  10  40 とする。
  1.0  10  40 に対し,W Wh 0となる共通の h 0
0を求める。 h 0  1ならW W, h 0  1とすれば
0

W Wh 0W Wh 0となるので 2分法で ,00となる h 0 を求める。コーデイン グではh  1.0e00;

f , g, h : double : floatif (g   0 ){  g; }else {   0 ; h  h 0 }f  rsqrt ();W
 W * f * f * h;
fの精度不足分は h 0に反映させて積分値全 体の精度低下を抑える。

M44 修正結果と結果の精度への影響
解析解=55.5852539156784
HD5870
N=1655.5852539156735
N=3255.5852539156776
N=32で修正した結果は以前と同じ精度

Pezy-SC
N=1655.5852539156731
単精度演算を使用しても必要な精度を保っています。

4． Rump’s例題を用いたサーバーでのDD形式とIEEE形式の性能比較

rump' s例題とは , 四則演算からなる簡単 な計算で , 有効ビット数が121ビット以上を要求する計算のため , 精度、性能の検証によく使用されるもので す。
a  77617.0, b  33096.0
f  333.75b 6  a 2 (11a 2 b 2  b 6  121b 4  2)  5.5b 8
a2ba 2  5.5b 2  1  6024398689


a 2 (11a 2 b 2  b 6  121b 4  2)  5.5b 2  333.75b 6  2a54767
より, f  2  0.827396059946
2b66192

DD形式
変数の値を複数個の倍精度変数の和で表す方式。
演算順序の変更がないように注意が必要。
サーバー系では-fp-model preciseが必要
n個の倍精度変数の和を2n倍精度という。
nが大きくなると,演算量は乗算がnの3乗,加減算がnの2乗で増える。
またnが3以上となると除算で途中の結果で絶対値が大きい順に並べる注意が必要となる。
IEEE形式IEEE754-2008での4倍精度の仮数部を拡張する。
符号部1ビット,指数部15ビット,仮数部32*p-16ビットをp倍精度という。演算がコール形式となるためpが小さい場合、オーバーヘッドが大きくなる。コーデイング
では乗算、除算は比較的容易で、加減算での桁落ち数の計算で誤りが生じ易い。
サーバー系では80ビットの内部レジスタを有するものが多く拡張倍精度や、複数個の拡張倍精度の和の演算が有効になる。特にオプション-fp-model extended は性能に効果的な場合が多い。

Rump’s 例題 N=500000000
E5-2687W v3 3.10GHz 40smp
拡張倍精度+拡張倍精度 有効ビット数 130>121
10倍精度から12倍精度で実行時間が逆転

実行時間(秒)250.0000
200.0000
150.0000
100.0000

DD形式
IEEE形式50.0000

16倍精度

14倍精度

12倍精度

10倍精度

8倍精度

6倍精度

2拡張倍精度

0.0000Rump’s 例題 N=500000000
E5-2687W v2 3.40GHz 16smp
拡張倍精度+拡張倍精度 有効ビット数 130>121
10倍精度から12倍精度で実行時間が逆転

実行時間(秒)500.0000
450.0000
400.0000
350.0000
300.0000
250.0000
200.0000
150.0000
100.0000
50.0000
0.0000

DD形式

16倍精度

14倍精度

12倍精度

10倍精度

8倍精度

6倍精度

2拡張倍精度IEEE形式

5. 3次元ループ積分でのサーバーでのDD形式とIEEE形式の性能比較使用した3次元積分は以下のinfra boxと言われるものでこの問題では10進6桁正しい結果を得るには

x  1  で  10 32.8  2 108 のxが取れる
必要がある。
1 1 x 1 x  yI(s)   0 0

01dzdydx
2D(D   xys  tz(1  x  y  z)  ( x  y)2 (1  x  y  z)(1  x  y)m e22

 z(1  x  y)m f )s  500 2 , t  150 2 , m f  150, m e  0.0005,   10 302( t  m f ) 2
s
解析近似解 ln( 2 ) ln
222 s(  t  m f )
mf me
12

( t  m f ) 26
相対誤差 710以下
2 2
sln( 2 ) m e ln ( t  m f )22
m mef

Infra box
E5-2690(SB) 2.90GHz
N=1024

16smp

DD形式4倍精度のみa,b,c
X=a+b をx=1-c (a,b,c>0) と表す
コーデイング。
10倍精度で実行時間がほぼ同じ
実行時間(秒)
900
800
700
600
500
400
300
200
100
0

DD形式

12倍精度

10倍精度

8倍精度

6倍精度

4倍精度

2拡張倍精度IEEE形式

Infra box
E5-2687W v2 3.40GHz 16smp
N=1024
DD形式4倍精度のみa,b,c
X=a+b をx=1-c (a,b,c>0) と表す
コーデイング。
8倍精度から10倍精度で実行時間が逆転

実行時間(秒)
800
600
400
200
0
12倍精度

10倍精度

8倍精度

6倍精度4倍精度

2拡張倍精度

DD形式
IEEE形式

DD形式4倍精度(a,b,c>0)
X=a+b

do i=1,nx=-1024+it=x*h
y=pi*0.5q0*sinh(t)x30(i)=exp(y)/(exp(y)+exp(-y))
gw30(i)=cosh(t)*pi/(exp(y)+exp(-y))**2
end doコーデイングでは赤字の1行変更するだけで良い。
X=1-c
do i=1,nx=-1024+it=x*h
y=pi*0.5q0*sinh(t)x30(i)=1.0q0/(1.0q0+exp(-2.0q0*y))gw30(i)=cosh(t)*pi/(exp(y)+exp(-y))**2
end do