2015年4月14日火曜日

疎行列のデータ配置

疎行列の続きデス。今回はデータ構造について直接説明するのではなく、
その構造のなかでどうデータを配置するかという話デス。

配置には大きく分けて2つの課題がありますデス。
  • ソート
  • 分散
データは必ずしも均一ではありません。しかしソートすると、 その偏りが制御できます。例えば、
■■
■■■
■■■
■■■■
■■■■
■■■■
■■
■■■■
■■
■■■
■■
■■■
なんて非0の数字が並んでいる疎行列があったとしましょう。そのまま、例えばGPGPUで4スレッド毎にまとめて計算すると(縦4つずつ1セットにして繰り返すと)、以下のようになります。□が4スレッド内の他のスレッドを待っている時間です。
■■□□
■■■□
■■■□
■■■■
■■■■
■■■■
■■□□
■■■■
■■□
■■■
■■□
■■■
待ち時間が多いですね。これを、非0の数字の数でソートすると、以下のようになります。待ち時間の□を減らすことができます(この例ではたまたま□が0個まで減らせています)。ソートすることで、仕事の不均衡を偏らせて、必要ない計算をする箇所を減らすことができたわけです。
■■■■
■■■■
■■■■
■■■■
■■■
■■■
■■■
■■■
■■
■■
■■
■■
反復法で解く場合、反復の前にソートして行を入れ替え、反復法を終えた後に元に戻せば済むため、こういったソートなどの重たい処理も、あまり性能に影響を与えることなく実施できます。

この方法は利用するデータ構造に依存せず利用できます。特に、GPGPU向けに、ELLPACK-Rと組み合わせて広く利用されています。例えば、Sparse matrix-vector multiplication on GPGPU clusters: A new storage format and a scalable implementationでは、pJDS (padded Jagged Diagonals Storage)法でソートする説明があります。

一方、OpenMPなんかで並列化する場合、static schedulingを使うことが多く、ソートした上で均一になるよう分散させたほうが都合がよくなります。例えば、4並列で処理する場合は例えば以下のようにソートした上で分散させるわけです。
■■■■
■■■
■■
■■■■
■■■
■■
■■■■
■■■
■■
■■■■
■■■
■■
例えば、地球シミュレータ向け線形ソルバーでは、pDJDS (parallel Descending-order Jagged Diagonal Storage)法でソートして均一になるよう分散して扱う説明があります。

スレッド数が少ない間ならこれでいいんですが、Xeon Phiのように244スレッドなどで処理することになると、これで長さが不十分になってしまいます。どうしたものか。ここから先は、仕事でやってる内容なので秘密です。

2015年4月12日日曜日

疎行列のデータ構造

疎行列のデータ構造のまとめの続き。
ELLPACKを中心に、ELLPACK-R形式、SELLPACK形式についてまとめます。

前回と同じく以下の配列があるとします。

1 7 0 0
0 2 8 0
5 0 3 9
0 6 0 4

前回説明したELLPACK形式。これに行中の0以外の値の個数を加えたのがELLPACK-R形式です。

column         values         row length
0 1 * 1 7 * 2
1 2 * 2 8 * 2
0 2 3 5 3 9 3
1 3 * 6 4 * 2

今回は縦方向に並べます。

column: 0 1 0 1 1 2 2 3 * * 3 *
values: 1 2 5 6 7 8 3 4 0 0 9 0
rl:     2 2 3 2

これのSpMVするプログラムはこうなります。
for (i = 0; i < M; ++i)
    for (j = 0; j < rl[i]; ++j)

        Y[i] += values[i + j*M] * X[column[i + j*M]];
この方式の利点はif文が必要なくなるということです。またGPGPUの場合4スレッドや8スレッドがセットで同じ回数繰り返す必要があります。そのため、回る回数が事前に分かると効率良くスレッド群(wrap)を割り当てることができます。更に縦方向にデータが並ぶため複数スレッドでデータを共有でき、メモリ効率も向上します。そんな次第で、ELLPACK-R形式はGPGPU向けの標準データ構造として広く使われています。

4スレッドや8スレッドでセットで処理するデータが、一箇所にまとまるよう複数行ごとに分割したのがSELLPACK形式です。Sliced ELLPACK形式の名の通り、複数行(S行)ごとに分断したデータ構造です。以下はS=2の例。

column         values
0 1 1 7
1 2 2 8
--------------------
0 2 3 5 3 9
1 3 * 6 4 *

column:    0 1 1 2 0 1 2 3 3 3
values:    1 2 7 8 5 6 3 4 9 0
row_start: 0 4 10

データ構造ではrow_startを使うことでパディングを減らします。CSR形式に近くなりました。ただし、Sで分けたグループ内では1行内の0以外の数が同じになるようパディングします。これのSpMVするプログラムはこうなります。
for (i = 0; i < M; ++i) {
    off = row_start[i/S];
    next = row_start[i/S + 1];
    for (j = 0; j < (next - off ) / S; ++j)

        Y[i] += values[i%S + j*S + off] * X[column[i%S + j*S + off]];
}
このSELLPACK形式は元々GPGPU向けに開発されましたが、実はSIMDプロセッサのように4演算毎、8演算毎に高速に演算できる場合、Sを8などに固定したSELLPACK形式が効率よいことが分かり、SIMD向けに注目されているデータ構造です。

更にSELLPACK形式のデータを例えば8つずつセットになるようパディングするなどしたSELL-P形式もありますが、ここでは省きます。

続きはまた今度

疎行列のデータ構造

最近、仕事で疎行列ばっか弄ってるのでまとめ。
COO形式、CSR形式、ELLPACK形式についてまとめます。

例えば以下の配列があるとします。

1 7 0 0
0 2 8 0
5 0 3 9
0 6 0 4

単純に行列として記録すると、上記を1次元に並べたデータ構造を使います。

values: 1 7 0 0 0 2 8 0 5 0 3 9 0 6 0 4

これとベクトルXを掛けてベクトルYを計算するプログラムはこうなります。
for (i = 0; i < M; ++i)
    for (j = 0; j < N; ++j)
        Y[i] += values[i*N+j] * X[j];

上記のデータだと0の部分が多いので、これを圧縮したのがCOO形式です。Coordinate形式の名の通り、座標を記録します。valuesに0が無くなるのが特徴です。

row:    0 0 1 1 2 2 2 3 3
column: 0 1 1 2 0 2 3 1 3
values: 1 7 2 8 5 3 9 6 4

これとベクトルXを掛けてベクトルYを計算する(Sparse Matrix-Vector multiplication、SpMVとよく呼びます)プログラムはこうなります。
for (i = 0; i < NZ; ++i)
    Y[row[i]] += values[i] * X[column[i]];
COO形式でもデータのrowの部分に同じ値が何度も出て来るので、これを圧縮したのがCSR形式です。Compressed Sparse Row形式の名前の通り、行を圧縮します。row_startの中に同じ数字が出てこなくなります。

row_start: 0 2 4 7 9
column:    0 1 1 2 0 2 3 1 3
values:    1 7 2 8 5 3 9 6 4

row_startは各行のデータが何番目から始まるかを示します。例えば、row_start[2]が4というのは、3行目(=2+1)のデータは、values[4]から始まるという意味です。

これのSpMVするプログラムはこうなります。
for (i = 0; i < M; ++i)
    for (j = row_start[i]; j < row_start[i+1]; ++j)

        Y[i] += values[j] * X[column[j]];
データは圧縮できていいんですが、計算に関節参照が増えます。ベクトル化しようとすると1行あたりのデータ数が一定じゃないのも辛い。

次は全然別のELLPACK形式。これはELLPACKというというelliptic problemsを解くためのシステムで利用されていたデータ構造ということで、この名前らしい。1行あたりのデータ数を決めて、パディングします。

column         values
0 1 * 1 7 *
1 2 * 2 8 *
0 2 3 5 3 9
1 3 * 6 4 *

*がパディングです。データは本来は縦方向に並べるのですが、今回は横方向に並べます。valuesの*は0を埋めておきます。1行あたりのデータ数が決まっているので、rowのデータは必要ありません。

column: 0 1 * 1 2 * 0 2 3 1 3 *
values: 1 7 0 2 8 0 5 3 9 6 4 0

これのSpMVするプログラムはこうなります。NPRが1行辺りのvalueの個数です。
for (i = 0; i < M; ++i)
    for (j = i*NPR; j < (i+1)*NPR; ++j)
        if (values[j] != 0)

            Y[i] += values[j] * X[column[j]];
if文で分岐せずに全て計算してしまっても構いません。そうするとこうなります。
for (i = 0; i < M; ++i)
    for (j = i*NPR; j < (i+1)*NPR; ++j)

        Y[i] += values[j] * X[column[j]];
valueの値の0を掛けるので、Xは何でもいいのです。ただし、領域外を参照しないようにcolumnの*には適当な値を設定する必要があります。0でもいいですが、キャッシュに無駄が出ないように一つ前と同じ値を埋めるのもいい方法です。同じ値を埋めると、データは以下のようになります。*の所にそれぞれ一つ前と同じ、1や2や3が入ります。

column: 0 1 1 1 2 2 0 2 3 1 3 3
values: 1 7 0 2 8 0 5 3 9 6 4 0

この方式の利点は規則的ってことです。またiとjのループを入れ替えると縦方向に処理ができます。大抵M >> NPRなので、この方がベクトル長が伸ばせて効率が良くなります。

欠点として、パディングの分データ量が増えます。1行あたりの0以外の値の数がほぼ均一なら問題はないのですが、長い行が混ざっていると無駄なデータが必要になってしまいます。

それを改善するのが、ELLPACKとCOOのハイブリッド方式です。1行につき2つのデータをELLPACK形式で保持し、残りのデータはすべてCOO形式にまとめます。COO形式が一つ一つのデータに依存関係がないため、どこからでも開始できるので、こうして他のデータを補完する形で使えるわけです。

ELLPACK部分
column: 0 1 1 2 0 2 1 3
values: 1 7 2 8 5 3 6 4

COO部分
row:    2
column: 3
values: 9

これのSpMVするプログラムはELLPACKとCOOを繋げた形になります。ここでは省略。

このハイブリッド形式は、ELLPACK部分が最適に作れて良さそうなデータに見えるのですが、1行あたりの0以外の値の数が均一である必要があるとか、要素の長い行は少なくないといけないとか、色々制限が多いので、実際の例ではあまり使われません。

続きはまた今度