ここでは、主に自分の中での整理、という意味で現在の書き直し版コード(名前
は lu2 ということにします)の概要をまとめておきます。
lu2 の特色は、
-
行列の格納順序を途中で転置することで、行スワップとピボットサーチの両
方を連続アクセスにする
-
パネル分解は普通の再帰。(Gustavson のアルゴリズム)
-
右上側の処理(通常なら DTRSM) に、DGEMM と再帰的 DTRSM を組み合わせた
新しいアルゴリズムを開発
というところです。最後のところは、普通の計算機では若干計算量が増えてし
まうのですが、 GRAPE-DR のような行列積だけが速くなる計算機では重要にな
ります。また、普通の計算機でも DGEMM のチューニングは DTRSM よりも容易
だったりするので、場合によってはメリットがあります。
転置の部分の関数を以下に示すます。これは 16 行の転置専用の関数です。
void transpose_rowtocol16(int n, double a[][RDIM], double at[][n],
int istart)
{
int i,j,k;
const int m=16;
double atmp[m][m];
int mend;
#ifdef TIMETEST
BEGIN_TSC;
#endif
//#pragma omp parallel for private(i,j,k,atmp)
for(i=istart;i<n;i+=m){
for(k=0;k<m;k++){
v2df * ak = (v2df*) a[i+k];
v2df * akk = (v2df*) atmp[k];
akk[0] =ak[0];
akk[1] =ak[1];
akk[2] =ak[2];
akk[3] =ak[3];
akk[4] =ak[4];
akk[5] =ak[5];
akk[6] =ak[6];
akk[7] =ak[7];
}
for(j=0;j<m;j++){
v2df * atk= (v2df*)(at[j]+i);
atk[0]=(v2df){atmp[0][j],atmp[1][j]};
atk[1]=(v2df){atmp[2][j],atmp[3][j]};
atk[2]=(v2df){atmp[4][j],atmp[5][j]};
atk[3]=(v2df){atmp[6][j],atmp[7][j]};
atk[4]=(v2df){atmp[8][j],atmp[9][j]};
atk[5]=(v2df){atmp[10][j],atmp[11][j]};
atk[6]=(v2df){atmp[12][j],atmp[13][j]};
atk[7]=(v2df){atmp[14][j],atmp[15][j]};
}
}
#ifdef TIMETEST
END_TSC(t,2);
#endif
}
これは特にこったことはしていなくて、行優先で格納された行列を16行読み出
して連続アドレス領域に一旦格納し、それを転置しながら格納します。
16行(128バイト)を単位にしているのは、 Core i7 ではキャッシュを2ライン
づつ勝手にに読んでくる機能があるので、その辺が性能的には具合がよいから
です。現在のコードでは、転置前の行列を連続アドレスの小行列にコピーして
から、いきなり転置行列に書き込んでいますが、一旦連続アドレスでの転置し
た小行列に書き込んでから最終的な縦長の行列に書くほうが具合がよいかもし
れません。
この部分の転送速度は、現在のところ Core i7 920 で 4GB read+write (合計
8GB) に 1.4 秒、 5.6GB/s というところです。転置前の行列サイズは最大で
4.25MB、平均ではずっと小さくて殆どの場合にL3 にいるはずであることを考え
るともうちょっと速くてもよいはずで、まだかなり工夫の余地があります。
L3 のバンド幅自体は 50GB/s 程度あるのに対してその 1/10 程度の速度しか
でていないからです。
もっとも、L3レイテンシが 50サイクル程度なので、16ワード、つまり128バイ
トアクセス毎に50サイクルかかるのかもしれません。それでは話にならないの
で、16x16 の行列全部を L1 までもってくるようなプリフェッチをかけること
ができれば改善しそうです。GCC の __builtin_prefetch では今一つ細かい制
御ができないので、
asm("prefetcht0 %0"::"m"(a[i+m][0]):"memory");
(prefetchnta を使うべき？)のように、実際に SSE プリフェッチを手で入れ
るべきかもしれません。 Ci7 は32KB もの巨大な L1 をもつのですが、
1度に読み込む行列は 2KB しかありません。なので、かなり先まで
プリフェッチをだしておいてもよいはずです。
逆方向の転置は以下のようにしています。
void transpose_coltorow16(int n, double a[][RDIM], double at[][n],
int istart)
{
int i,j,k;
const int m=16;
double atmp[m][m];
#ifdef TIMETEST
BEGIN_TSC;
#endif
#pragma omp parallel for private(i,j,k,atmp)
for(i=istart;i<n;i+=m){
for(k=0;k<m;k++){
v2df * ak = (v2df*) (at[k]+i);
v2df * akk = (v2df*) atmp[k];
akk[0] =ak[0];
akk[1] =ak[1];
akk[2] =ak[2];
akk[3] =ak[3];
akk[4] =ak[4];
akk[5] =ak[5];
akk[6] =ak[6];
akk[7] =ak[7];
}
for(j=0;j<m;j++){
v2df * atk= (v2df*)(a[i+j]);
atk[0]=(v2df){atmp[0][j],atmp[1][j]};
atk[1]=(v2df){atmp[2][j],atmp[3][j]};
atk[2]=(v2df){atmp[4][j],atmp[5][j]};
atk[3]=(v2df){atmp[6][j],atmp[7][j]};
atk[4]=(v2df){atmp[8][j],atmp[9][j]};
atk[5]=(v2df){atmp[10][j],atmp[11][j]};
atk[6]=(v2df){atmp[12][j],atmp[13][j]};
atk[7]=(v2df){atmp[14][j],atmp[15][j]};
}
}
#ifdef TIMETEST
END_TSC(t,3);
#endif
}
要するに、構造は同じです。構造は同じで、メモリアクセスの量も同じなので
すが、こちらのほうが2倍近く速度が速くなっています。これは、読出しアド
レスの連続性が高い、つまり、 row-to-col の場合には長さ n のベクトルか
ら 16 ワードを読む、というのを、 n 個の別のベクトルに順番に適用するの
に対して、 col-to-row では 16 本の長さ n のベクトルから 16 語ずつ読む
ので、フットプリントが小さい、というのが1つの理由でしょう。
但し、2倍速いだけで、十分に速い、というわけでは決してないので、改善の
余地がないわけではありません。
asm("prefetcht0 %0"::"m"(at[k][i+m*3]):"memory");
asm("prefetcht0 %0"::"m"(at[k][i+m*3+8]):"memory");
をいれてみた場合のバンド幅変化は以下のようでした
ctor N=8192 16384
no prefetch 0.2443 0.2613
prefetchnta 0.3177
prefetcht0 0.3180
prefetcht1 0.3235
prefetcht2 0.3239 0.3235
rtoc N=8192 16384
no prefetch 0.2327 0.2043
prefetchnta 0.2813 0.2506
prefetcht0 0.2789
prefetcht1 0.2608
prefetcht2 0.2605
builtin(0,0) 0.2778
で、 prefetch で 2-3割は上がりますが、5倍とかには到底いかないです。
試しに、rtoc のケースで小行列へのコピーと小行列からのコピーを別々に測
定して見ると、
小行列へのコピー 0.5746
小行列からのコピー 0.6643
で、小行列「から」に意外に時間がかかっています。これを、転置を小行列の
範囲でもう一度やるようにすると、最後の連続アドレスでのコピーは 0.8語/
サイクル程度まで向上しますが、転置自体に 1語/サイクル程度の時間がかかっ
て総合的にはかえって速度低下します。
通常の計算機で実行する分には、例えば 40Gflops で 34k の行列だと 700秒
程度になり、1秒を問題にする意味はあまりありません。が。 1Tflops 出すと
これが 25秒で終わってしまうので、行列乗算でないところを5秒程度まで
抑え込む必要があります。このためには、つめられるものは極限までつめる必
要があるわけです。現在のコード(2009/9/13現在)では、34k での実行時間が
行交換 1.3
行スケール 0.6
行->列転置 1.4
列->行転置 0.8
転置後の行列乗算 2.0
2列でのソルバ 1.5
合計 7.6 (秒)
という程度ですが、これをあと 2-3秒はつめたいところです。
これらの処理は殆ど L3 にあるデータに対する操作なので、全体にもう
数倍は速くできてもよいはずなので、これは不可能ではないと思います。
と、話がずれてますね。再帰のところは、実際に書いてみると
void column_decomposition_recursive(int n,
double a[n][RDIM],
int m,
double awork[][n],
int pv[],
int i)
{
int j, k;
int ip,ii;
double ainv;
if (m <= 16){
column_decomposition_with_transpose(n, a, m,awork, pv,i);
}else{
column_decomposition_recursive(n, a, m/2, awork, pv,i);
process_right_part(n,a,m/2,awork, pv,i,i+m);
column_decomposition_recursive(n, a, m/2, awork, pv+m/2,i+m/2);
for(ii=i+m/2;ii<i+m;ii++){
swaprows(n,a,pv[ii-i],ii,i,i+m/2);
}
for(ii=i+m/2;ii<i+m;ii++){
scalerow(n,a,1.0/a[ii][ii] ,ii,i,i+m/2);
}
}
}
と単純な構造です。日本語での疑似コードにして見ると
幅が 16 以下になったら
転置してパネル分解する
でなければ
左半分をパネル分解する
左半分の分解結果を使って右側半分を更新する
右側半分を分解する
右側半分の結果を使って、左半分の行交換とスケーリングをする
です。
右上側の処理(通常なら DTRSM) に、DGEMM と再帰的 DTRSM を組み合わせた
新しいアルゴリズムですが、本来するべき処理は
for(ii=i;ii<i+m;ii++)
for(j=ii+1;j<i+m;j++)
for (k=i+m;k<n;k++)
a[j][k] -= a[j][ii]*a[ii][k];
というものです。図で書くと以下のような感じです。
+-----------+----------------------------------+
|\ | | |
| \ | | |
| \ | | |
| \ | | |
| \ | | |
|-----\ | X |
| \ | |
| \ | |
| \ | |
| \ | |
+-----------+----------------------------------+
X の位置のデータから、左の下三角行列の 「---」と、 X の上の列との内積
を引くわけですが、単純な行列積ではなくて、1つ上の行の要素は変形してか
らの操作になります。わかりにくいので、右側の行列の列1行についての操作
にしてみると、
for(ii=i;ii<i+m;ii++)
for(j=ii+1;j<i+m;j++)
b[j] -= a[j][ii]*b[ii];
ですね。この形式だと b の計算はいわゆる recurrent なもので、 b[i] が求
まらないと b[i+1] が計算できないですが、数学的には、新しい b[k] は
b[0] から b[k] までの線形結合ですから、a を変形しておけば
for(j=i+1;j<i+m;j++){
c[j]=0;
for(ii=i;ii<=j;ii++)
c[j] += anew[j][ii]*b[ii];
}
という形にできます。これは、anew も上半分がない、ということを別にすれ
ば単純なベクトル行列積なので、k の方向も考えると行列積になって、
効率が半分になるのをいいことにすれば GRAPE-DR や、他の高速な行列計算ラ
イブラリで実行できます。
anew は、単位行列を上の処理と同じ計算で変形すればよい、つまり、
for(ii=0;ii<m;ii++)
for(j=ii+1;j<m;j++)
for (k=0;k<m;k++)
b[j][k] -= a[j][ii]*b[ii][k];
というような計算をすればよいわけですが、ここでも DGEMM を引っ張り出す
ことを考えてみます。書いたコードを示します。
static void solve_triangle_for_unit_mat_recursive(int n,
double a[][RDIM],
int nb,
double b[][nb],
int m)
{
int i,ii,j,k;
if (m < 16){
solve_triangle_for_unit_mat_internal(n, a, nb, b,m);
return;
}
const int mhalf = m/2;
solve_triangle_for_unit_mat_recursive(n, a, nb, b,mhalf);
dgemm( mhalf, mhalf, mhalf, -1.0, &(a[mhalf][0]), RDIM,
&(b[0][0]), nb, 1.0, &(b[mhalf][0]),nb );
double bwork[mhalf][mhalf];
double bwork2[mhalf][mhalf];
for (j=0;j<mhalf;j++)
for (k=0;k<mhalf;k++)bwork[j][k]=0.0;
for (j=0;j<mhalf;j++)bwork[j][j]=1.0;
solve_triangle_for_unit_mat_recursive(n, (double(*)[])(&a[mhalf][mhalf]),
mhalf, bwork,mhalf);
for(i=0;i<mhalf;i++)
for(j=0;j<mhalf;j++)
bwork2[i][j]=b[i+mhalf][j];
dgemm(mhalf, mhalf, mhalf, 1.0, (double*)bwork,mhalf,
(double*)bwork2, mhalf, 0.0, &(b[mhalf][0]),nb );
for (j=0;j<mhalf;j++)
for (k=0;k<j+1;k++)b[mhalf+j][mhalf+k]=bwork[j][k];
}
このコードに表現した手順は、b に単位行列がはいっているとして、
まず左上 1/4 だけを変形する
変形した結果と a の左下 1/4 との積を bの左下 1/4 から引く
右下 1/4 を変形する
変形した結果と b の左下 1/4 との積を bの左下 1/4 から引く
というものです。これにでてくる行列積に単純に DGEMM を使うと、
計算量が2倍になって若干損ですが、 GRAPE-DR 等が使えるならそれでも
そのほうがメリットがあるわけです。普通の計算機なら、3角行列と
正方行列の乗算専用のルーチンを書けば問題はありません。