12. HPL 書き直しその3 (2009/9/14 書きかけ)
lu2 では、基本データ構造は単純な row major としましたが、実際にコードを
書いて性能を測ってみるとやはりこれは問題があるように見えます。それは、
例えば hugetlb を使うにしても、細い幅で縦にアクセスする時にどうも遅い、
ということ、もうひとつは、横にアクセスする時にも、並列化効率は低い、と
いうことです。
例えば行をスケールする時には、複数の行を並列に処理でき、この時には複数
のコアが全く違うところをアクセスするので高速にできます。ところが、行交
換では、1行づつやっていくので並列化は行内ですることになり、連続アドレス
のベクトルを分割してアクセスすることになるのであまり並列化効率が上がり
ません。
これを改善する方法ですが、以下の2つがありえるように思います。
-
転置のサイズをずっと大きくする
-
単純な row major でなく、ブロック化する
現在は幅 16 くらいで転置していて、これは基本的には転置した前後に行列が
l3 に入ることを期待しているからです。しかし、はいっているにもかかわら
ず速度が 0.3words/cycle, つまり 2.67ghz では 14gb/s 程度しかでないなら、
主記憶を使っても速度はあまり変わりません。この場合には、
幅が 16 とか 32 とか、もっと広い幅でのアクセスをしたほうがよいような気
もします。
しかし、これは実際にやってみると速度が 1/4 程度になってしまって全く話
になりません。
ブロック化は、例えば 2048 くらいの長さに行を切って、それが縦に並ぶ、と
いう形にする、というものです。この、n x 2048 のブロックを横に並べて行
列を作るわけです。 そうすると、行交換等は幅 2048 のブロック単位に並列
化すればアドレス空間は完全に独立になります。また、転置したものも同様に
ブロック化しておくと、行列乗算の処理をアドレス空間を分離して並列化でき
ます。転置についても有利になるはずです。
この場合に効率の観点から1つだけ問題なのは、DGEMM を呼ぶのが
煩雑になる、つまり、スレッド並列されている DGEMM を呼ぶためには
横長の行列に代入しなおして、結果も戻す必要があることです。
但し、コードはかなり煩雑になるので、あまりやりたくない、という気もしま
す。
行交換については、もうちょっと賢いやり方はないか？というのが問題です。
例えば、行交換とスケールで、広い領域を2度アクセスするのはいかにも非効
率です。また、行交換には相当程度並列性があるはずです。つまり、ほとんど
の交換は独立なはずです。もちろん、問題は、独立でない行交換があったらど
うするかです。
独立でない交換とは、例えば 2 行目を 8行目と交換して、次に3行目を交換し
たらそこに前の2行目がはいってきてしまった、というものです。ブロックサ
イズが大きくなると、起こる確率はもちろんあがるので、無視して処理するわ
けにはいきません。
まあ、そうはいっても割合は低いはずなので、単純に、独立でないものは分離
してシリアルに処理して、それ以外を並列に処理、というのでよいように思い
ます。
但し、 swaprows については、現在でも 15GB/s ていどの速度なので改善の余
地はそれほどないかもしれません。であれば、単純に、 swaprows で上側に代
入する時にスケールもする、というので十分でしょう。これは実装しています。
void swaprows_simple_with_scale(int n, double a[n][RDIM],
int row1, int row2,
int cstart, int cend, double scale)
{
int j;
if (row1 != row2){
int jmax = (cend+1-cstart)/2;
v2df *a1, *a2, tmp, tmp1;
v2df ss = (v2df){scale,scale};
a1 = (v2df*)(a[row1]+cstart);
a2 = (v2df*)(a[row2]+cstart);
for(j=0;j<(jmax & (0xfffffffe));j+=2){
__builtin_prefetch(a1+j+32,1,0);
__builtin_prefetch(a2+j+32,1,0);
tmp = a1[j];
a1[j]=a2[j];
a2[j]=tmp*ss;
tmp1 = a1[j+1];
a1[j+1]=a2[j+1];
a2[j+1]=tmp1*ss;
}
if (jmax &1){
tmp = a1[jmax-1];
a1[jmax-1]=a2[jmax-1];
a2[jmax-1]=tmp*ss;
}
}else{
scalerow(n,a,scale ,row2,cstart,cend);
}
}
まあ、どうして初めからそうしてなかったんだ？というような話です。
これで、 34 k で交換+スケーリングが大体 1.3秒で、スケーリングにかかっ
ていた 0.6 秒が減りました。
2列にまで再帰したところでの変形では、BLAS1 でいう DAXPY な操作がでてき
ます。これについては sse を使ってさらに4重にアンロールすると、バンド幅
として 20GB/s 程度になり、まあいいかな、という程度です。
この程度までチューニングすると、ピボットのための最大値探索が結構遅いも
のとして見えてきます。
Enter n, seed, nb:N=34816 Seed=1 NB=256
read/set mat end
copy mat end
Emax= 1.136e-07
Nswap=0 cpsec = 2784.11 wsec=696.654 40.3857 Gflops
swaprows time=3.43235e+09 ops/cycle=0.176578
scalerow time=1.56005e+07 ops/cycle=38.8497
trans rtoc time=3.14886e+09 ops/cycle=0.192475
trans ctor time=1.89551e+09 ops/cycle=0.319744
trans mmul time=5.14042e+09 ops/cycle=1.76856
trans nonrec cdec time=2.94967e+09 ops/cycle=0.205472
trans vvmul time=5.03337e+08 ops/cycle=1.20412
trans findp time=2.44507e+09 ops/cycle=0.247877
solve tri u time=3.80914e+08 ops/cycle=9.14013e-05
solve tri time=2.67092e+10 ops/cycle=11.6182
matmul nk8 time=0 ops/cycle=inf
matmul snk time=67816 ops/cycle=142993
trans mmul8 time=0 ops/cycle=inf
trans mmul4 time=1.39435e+09 ops/cycle=1.73866
trans mmul2 time=8.13527e+08 ops/cycle=1.49
DGEMM time=1.83959e+12 ops/cycle=15.4628
これは、小細工ではどうしようもないのですが、
Intel の技術メモ に詳しく議論したものがあって SSE2 を使ってオールアセ
ンブラでやれば 2-3倍速くなりそうです。これは、自力で書くものでもなくて
GotoBLAS の idamax を使うので OK でした。
Enter n, seed, nb:N=34816 Seed=1 NB=256
Emax= 1.136e-07
Nswap=0 cpsec = 2782.51 wsec=696.245 40.4095 Gflops
swaprows time=3.43022e+09 ops/cycle=0.176688
scalerow time=1.56175e+07 ops/cycle=38.8077
trans rtoc time=3.14043e+09 ops/cycle=0.192992
trans ctor time=1.89603e+09 ops/cycle=0.319656
trans mmul time=5.14364e+09 ops/cycle=1.76745
trans nonrec cdec time=1.40592e+09 ops/cycle=0.431089
trans vvmul time=4.97509e+08 ops/cycle=1.21822
trans findp time=9.07189e+08 ops/cycle=0.668083
solve tri u time=3.85358e+08 ops/cycle=9.03471e-05
solve tri time=2.66933e+10 ops/cycle=11.6251
matmul nk8 time=0 ops/cycle=inf
matmul snk time=59096 ops/cycle=164093
trans mmul8 time=0 ops/cycle=inf
trans mmul4 time=1.39377e+09 ops/cycle=1.73939
trans mmul2 time=8.11773e+08 ops/cycle=1.49322
DGEMM time=1.84005e+12 ops/cycle=15.459
現時点で、大きく残っているのは k=8 の行列乗算で、k=4 の場合に対して演
算量は2倍。メモリアクセスは半分なんだが時間は2倍でほぼ1秒かかっている。
k=8 の時には GotoBLAS を使っているが、素晴らしく速い、というわけでも
ありません。それ以前にそもそも問題なのは k=4 のケースです。k=2 の場合に
サイクル当り 1.5 演算です。これはメモリバンド幅としては、アクセス 3 語
あたりに4演算するので、9/8語/サイクルとなり、主記憶一杯一杯です。キャッ
シュにはいってるはずなのにどうして遅いのかはともかく、そんな悪い数字で
はありません。で、希望としては、k=4 になったらメモリアクセスに対する演
算数は倍になるので、3演算くらいしてくれ、と思うわけです。ここはシング
ルコア実行なので、4演算が理論限界ですが、それでもまあ3くらいはいってく
れてもよさそうなものです。現在使っている k=4 のコードは以下です。
static void matmul_for_nk4(int n1, double a[][n1],
int n2, double b[][n2],
int n3, double c[][n3],
int n)
{
register v2df b00 = (v2df){b[0][0],b[0][0]};
register v2df b01 = (v2df){b[0][1],b[0][1]};
register v2df b02 = (v2df){b[0][2],b[0][2]};
register v2df b03 = (v2df){b[0][3],b[0][3]};
register v2df b10 = (v2df){b[1][0],b[1][0]};
register v2df b11 = (v2df){b[1][1],b[1][1]};
register v2df b12 = (v2df){b[1][2],b[1][2]};
register v2df b13 = (v2df){b[1][3],b[1][3]};
register v2df b20 = (v2df){b[2][0],b[2][0]};
register v2df b21 = (v2df){b[2][1],b[2][1]};
register v2df b22 = (v2df){b[2][2],b[2][2]};
register v2df b23 = (v2df){b[2][3],b[2][3]};
register v2df b30 = (v2df){b[3][0],b[3][0]};
register v2df b31 = (v2df){b[3][1],b[3][1]};
register v2df b32 = (v2df){b[3][2],b[3][2]};
register v2df b33 = (v2df){b[3][3],b[3][3]};
v2df * a0 = (v2df*) a[0];
v2df * a1 = (v2df*) a[1];
v2df * a2 = (v2df*) a[2];
v2df * a3 = (v2df*) a[3];
v2df * c0 = (v2df*) c[0];
v2df * c1 = (v2df*) c[1];
v2df * c2 = (v2df*) c[2];
v2df * c3 = (v2df*) c[3];
int nh = n>>1;
int j;
if (nh & 1){
j=nh-1;
c0[j] -= a0[j]*b00+a1[j]*b01+a2[j]*b02+a3[j]*b03;
c1[j] -= a0[j]*b10+a1[j]*b11+a2[j]*b12+a3[j]*b13;
c2[j] -= a0[j]*b20+a1[j]*b21+a2[j]*b22+a3[j]*b23;
c3[j] -= a0[j]*b30+a1[j]*b31+a2[j]*b32+a3[j]*b33;
nh--;
}
for(j=0;j<nh;j+=2){
c0[j] -= a0[j]*b00+a1[j]*b01+a2[j]*b02+a3[j]*b03;
c0[j+1] -= a0[j+1]*b00+a1[j+1]*b01+a2[j+1]*b02+a3[j+1]*b03;
c1[j] -= a0[j]*b10+a1[j]*b11+a2[j]*b12+a3[j]*b13;
c1[j+1] -= a0[j+1]*b10+a1[j+1]*b11+a2[j+1]*b12+a3[j+1]*b13;
}
for(j=0;j<nh;j+=2){
c2[j] -= a0[j]*b20+a1[j]*b21+a2[j]*b22+a3[j]*b23;
c2[j+1] -= a0[j+1]*b20+a1[j+1]*b21+a2[j+1]*b22+a3[j+1]*b23;
c3[j] -= a0[j]*b30+a1[j]*b31+a2[j]*b32+a3[j]*b33;
c3[j+1] -= a0[j+1]*b30+a1[j+1]*b31+a2[j+1]*b32+a3[j+1]*b33;
}
}
ループを2重にアンロールするので、最初に端数の処理をします。それから、
まず c の0、1行目を計算してから独立のループで 2, 3 行目を処理します。
大原 さんの解析によると、レジスタとL1 の転送速度限界はread/writeそれぞれ 2語
/Cycle (L2, L3 はその半分、 1/4)とのこ
となので、L3 にきているだけのデータだと上の形のループでは 7+2語アクセ
スで 16演算なので、 L3 からデータがくるなら 14 サイクルで 16 演算、L1
からきたら3.5 サイクルで 16 演算になってピークがでるかも、という感じで
す。ただ、L1 に入るデータサイズにしたとしても、レジスタの本数が足りな
いような気がします。 xmm レジスタは16本なので、基本的に入らないわけで
す。そうすると、bが全部溢れているとすれば32演算に対して 22 リードで、
L1 に全部はいってスケジューリングも理想的として 3 演算です。なので、
ここはちゃんと b がレジスタに止まるようなコードを書くことが必須です。
出力アセンブラをみながら色々やってみたところでは、以下のコードは b を
完全にレジスタに置くようです (gcc 4.1.2, -O2 の場合)
static void matmul_for_nk4(int n1, double a[][n1],
int n2, double b[][n2],
int n3, double c[][n3],
int n)
{
v2df * a0 = (v2df*) a[0];
v2df * a1 = (v2df*) a[1];
v2df * a2 = (v2df*) a[2];
v2df * a3 = (v2df*) a[3];
int nh = n>>1;
int j;
int k;
v2df * cv;
v2df * cvv;
register v2df b0;
register v2df b1;
register v2df b2;
register v2df b3;
register v2df b4;
register v2df b5;
register v2df b6;
register v2df b7;
for(k=0;k<4;k+=2){
cv = (v2df*) c[k];
cvv = (v2df*) c[k+1];
b0 = (v2df){b[k][0],b[k][0]};
b1 = (v2df){b[k][1],b[k][1]};
b2 = (v2df){b[k][2],b[k][2]};
b3 = (v2df){b[k][3],b[k][3]};
b4 = (v2df){b[k+1][0],b[k+1][0]};
b5 = (v2df){b[k+1][1],b[k+1][1]};
b6 = (v2df){b[k+1][2],b[k+1][2]};
b7 = (v2df){b[k+1][3],b[k+1][3]};
for(j=0;j<nh;j++){
register v2df aa0 = a0[j];
register v2df aa1 = a1[j];
register v2df aa2 = a2[j];
register v2df aa3 = a3[j];
register v2df x = aa0*b0;
x+= aa1*b1;
x+= aa2*b2;
x+= aa3*b3;
cv[j]-=x;
x = aa0*b4;
x+= aa1*b5;
x+= aa2*b6;
x+= aa3*b7;
cvv[j] -= x;
}
}
}
この時に、性能は 2.6 演算/サイクル程度でした。アセンブラ出力を見ると
.L19:
movapd (%rax,%rsi), %xmm6
addl $1, %edi
movapd (%rax,%r12), %xmm3
movapd %xmm6, %xmm1
mulpd %xmm10, %xmm6
movapd %xmm3, %xmm0
mulpd %xmm14, %xmm1
movapd (%rax,%rbp), %xmm4
mulpd %xmm13, %xmm0
movapd (%rax,%rbx), %xmm5
mulpd %xmm9, %xmm3
movapd %xmm5, %xmm2
mulpd %xmm7, %xmm5
addpd %xmm1, %xmm0
movapd %xmm4, %xmm1
mulpd %xmm11, %xmm2
addpd %xmm6, %xmm3
mulpd %xmm12, %xmm1
mulpd %xmm8, %xmm4
addpd %xmm0, %xmm1
movapd (%rax,%r10), %xmm0
addpd %xmm3, %xmm4
addpd %xmm1, %xmm2
addpd %xmm4, %xmm5
subpd %xmm2, %xmm0
movapd %xmm0, (%rax,%r10)
movapd (%rax,%r8), %xmm0
subpd %xmm5, %xmm0
movapd %xmm0, (%rax,%r8)
addq $16, %rax
cmpl %r11d, %edi
jne .L19
という感じです。 %xmm[n] は SSE2 命令のレジスタなので、その間の演算が殆
どで、メモリからのロード (movapd (%rax,%r12), %xmm3 とか)が6、ストア
(movapd %xmm0, (%rax,%r8) とか)が2で、メモリアクセスの観点からはベスト
なコードです。ループ内の演算数はあまり多くないし、依存関係が多いので、
スケジューリングは上手くいっていないかもしれません。ループアンロールを
するともうちょっと改善する可能性はあります。
また、このループでは L1 に入る範囲に切る必要があり、さらに少なくとも
k=8 の場合には並列化して速度を稼ぐ必要もあるので、そういうチューニング
が残っています。