行列積とループ・ブロッキング。

インテル日本語資料の「IA-32 インテルアーキテクチャー最適化」にも記されている、ループ・ブロッキングの有効性を検証してみた。
ループ・ブロッキングとは、キャッシュサイズを考慮してループをブロックすることである。
反復のたびにおきてしまうキャッシュミスを軽減することができ、キャッシュとメモリーのパフォーマンスを向上させることができる。

実験環境は Xeon 3.0GHz であるが、今回はスレッド数は１に固定した。ちなみにキャッシュサイズは 4096 KB となっている。
今回は 512 x 512 と 1024 x 1024 の２種類のみであるが、次の４パターンについてプロファイルしてみた。
(0)　GotoBLAS （あくまでも参考）
(1)　通常の積演算
(2)　通常の積演算にループ・ブロッキングを適用
(3)　事前に行列Ｂに対して転置処理を行い、アクセス方向を修正した積演算
(4)　事前に行列Ｂに対して転置処理を行い、アクセス方向を修正した積演算に対してループ・ブロッキングを適用


ループの大きさによってかなりのバラツキがあることにおどろいた。
通常のループであれば、ブロックサイズは２～８程度と小さめで、
事前にアクセス方向を修正しておけば、２５６と大きめがよいようだ。


512 x 512
　　　　1　　　2　　4　　 8　　　16　　32　　64　 　128　 256 　　512
(0)　　35
(1)　1700
(2)　1500　765　825　1000　1327　1400　1420　1450　1450　1500
(3)　 190
(4) 　920　365　257　　228　　211　200　　196　　202　198　　196


1024 x 1024
　　　　1　　　　2　　　　4　　　8　　　16　　　32　　　64　 　128　 　256 　　512　　1024
(0) 235
(1) 60000
(2) 65000　23300　10700　9600　11200　11500　11500　11800　11700　17500　54500
(3) 2550
(4) 7600　　3050　　2120　1800　　1710　1625　　1580　　1626　1610　　2050 　2420


以下はソースコードの一部である。


　/* (1)通常の積演算 */
　for (m = 0; m ＜ M; m++) {
　　for (n = 0; n ＜ N; n++) {
　　　for (k = 0; k ＜ K; k++) {
　　　　C[m * N + n] += A[m * K + k] * B[k * N + n];
　　　}
　　}
　}


　/* (2)通常の積演算にループ・ブロッキングを適用 */
　for (m = 0; m ＜ M; m+=block) {
　　for (n = 0; n ＜ N; n+=block) {
　　　for (k = 0; k ＜ K; k+=block) {
　　　　
　　　　for (mm = m; mm ＜ m+block; mm++) {
　　　　　for (nn = n; nn ＜ n+block; nn++) {
　　　　　　for (kk = k; kk ＜ k+block; kk++) {
　　　　　　　C[mm * N + nn] += A[mm * K + kk] * B[kk * N + nn];
　　　　　　}
　　　　　}
　　　　}
　　　　
　　　}
　　}
　}

　/* (3)事前に行列Ｂに対して転置処理を行い、アクセス方向を修正した積演算 */
　for (m = 0; m ＜ M; m++) {
　　for (n = 0; n ＜ N; n++) {
　　　for (k = 0; k ＜ K; k++) {
　　　　C[m * N + n] += A[m * K + k] * B[n * K + k];
　　　}
　　}
　}

　/* (4)事前に行列Ｂに対して転置処理を行い、
　　　アクセス方向を修正した積演算に対してループ・ブロッキングを適用 */
　for (m = 0; m ＜ M; m+=block) {
　　for (n = 0; n ＜ N; n+=block) {
　　　for (k = 0; k ＜ K; k+=block) {
　　　　
　　　　for (mm = m; mm ＜ m+block; mm++) {
　　　　　for (nn = n; nn ＜ n+block; nn++) {
　　　　　　for (kk = k; kk ＜ k+block; kk++) {
　　　　　　　C[mm * N + nn] += A[mm * K + kk] * B[nn * K + kk];
　　　　　　}
　　　　　}
　　　　}
　　　　
　　　}
　　}
　}

DMA転送と処理のバランス２。

前回に続いて、SPE においてのデータ転送、処理サイズのバランスについて。

DGEMM(Double)では、２行ずつ処理していくのが効率的であると述べたが、
今回は SGEMM(Float)について、プロファイルをしてその結果から考察していく。


まずは、サイズ、使用 SPE ごとの実行時間である。

SPENUM　512x512　1024x1024　2048x2048　[msec]
　　１　　　360　　　 2600　　　22000
　　２　　　200　　　 1500　　　12000
　　３　　　150　　　 1060　　　　8200
　　４　　　125　　　 830　　　　 6400
　　５　　　113　　　 710　　　　 5400
　　６　　　105　　　 630　　　　 4700


倍精度演算と比べれば確かに早くなってはいるが、あまり期待どおりとはいえない。
今回も前回同様、処理とデータ転送のバランスを見ていく。


[データ転送 / 積和演算]
SPENUM　　　 512x512　　1024x1024　　2048x2048　[msec]
　　１　　0.17 / 0.4　　0.3 / 0.8　　0.6 / 1.7
　　～　
　　６　　0.2　/ 0.4　　0.4 / 0.8　　0.8 / 1.7

結果からは、思ったよりも積和演算が足を引っ張っているのかもしれないと考えられる。



[データ転送数]
512x512 では、512 x 4(２行、２行) x 4(float) = 8Kbyte
1024 x 1024 では 16Kbyte、2048 x 2048 では 32Kbyte となる。

[積和演算数]
512x512 では、2 x 2 x 512 = 2048 回の積和演算、
1024 x 1024 では 4096 回、2048 x 2048 では 8192回の積和演算することになる。

8Kbyte の２倍と 2048 回の積和演算、16Kbyte の２倍と 4096 回の積和演算、
32Kbyte の２倍と 8192 回の積和演算がそれぞれ同じくらいとすると、
やはり「データ：積和演算＝８：１」であるようだ。

処理を早くすることは難しいので、
それに見合うように「４行ずつ DMA 転送する」などの改良してみようと思う。
もう少し調べてみよう。

DMA転送と処理のバランス。

Cell プログラミングは SPE を上手に使うことがキーとなっている。
基本的に、PPEプログラムが SPE プログラムを立ち上げ、
必要なデータを SPE へ DMA 転送し、処理した後、DMA 転送によって PPE に返す。

SPE スレッドの立ち上げや、前処理、後処理にはある程度オーバーヘッドがかかってしまうが、
SPE はそのオーバーヘッド以上を取り戻すほどの性能を持っている。
特にSPE 内で１度に扱うデータサイズをうまく決めてあげればかなりの速度向上をすることができる。

DGEMM に関してだが、最適なデータサイズというのがあまり決まらず、DMA 転送が足を引っ張ると誤認識していた。
そのため、一度に多くのデータを持ってきてしまい、演算を繰り返すというスタイルでの処理にしようかと考えていた。
具体的には １度に DMA 転送できる最大サイズである 16Kbyte 分、DMA 転送し行数分演算する。
512 x 512 では、４行ずつ、1024 x 1024 では、２行ずつ、2048 x 2048 では１行ずついった感じだ。
512 x 512 では、4 x 4 x 512 = 8192 回の積和演算、1024 x 1024 では、2 x 2 x 1024 = 4096 回、
2048 x 2048 では 2048 回というように、DMA 転送はできるだけしない、するならたくさんとってくる。


しかし、プロファイルをとってみると予想とは違う結果になった。
なかなか良さそうなのは今までのやりかたである、「２行ずつの処理」である。
Double Buffer を用いて高速化したいのであれば、
DMA 転送と、積和演算の割合を均等にした方がより効果を得られやすい。


「２行ずつの処理」について見ていくと、簡単な流れは以下のようになっている。

STEP1　行列Ａと行列Ｃを２行ずつ DMA 転送してくる
STEP2　行列Ａ、行列Ｃにそれぞれα、βを乗じる
STEP3　行列Ｂを２行ずつ DMA 転送してくる
STEP4　行列Ａの一部と行列Ｂの一部で演算を行い、結果を行列Ｃに書き込む
STEP5　STEP3～STEP4 を行列Ｂ分繰り返す
STEP6　STEP1～STEP5 をその SPE に割りふられた分だけ繰り返す

実際に実行時間の割合が大きそうな、STEP1 と STEP4 についてプロファイルしてみた。



[DMA 転送]
512 x 512 であれば、512 x 4(行列Ａ、行列Ｃ、それぞれ２行ずつ) x 8(倍精度) = 16Kbyte
1024 x 1024 であれば 32Kbyte、2048 x 2048 であれば 64Kbyte となる。

DMA 転送に要した時間の１ループ分 （STEP1）
使用SPE数　　512　　1024　 2048　[msec]
　　　　１　　0.3　　0.6　　1.2
　　　　２　　0.3　　0.6　　1.2
　　　　３　　0.4　　0.6　　1.3
　　　　４　　0.4　　0.7　　1.4
　　　　５　　0.4　　0.8　　1.5
　　　　６　　0.4　　0.8　　1.5



[積和演算]
積和演算は１ループで４回行うことになる
(行列Ａの１行目 x 行列Ｂの１行目、行列Ａの１行目 x 行列ＣＢ２行目、……)ので、
512 x 512 では 2(行列Ａ２行分) x 2(行列Ｂ２行分) x 512 で、2048 回積和演算することになる。
1024 x 1024 では 4096 回、 2048 x 2048 では 8096 回積和演算を行うことになる。
実際には SIMD 演算（倍精度なら１度に２個ずつ）してしまうので、実際にはその半数になる。

内積に要した時間の１ループ分（STEP4）
使用SPE数　　512　　1024　 2048　[msec]
　　１～６　　0.3　　0.6　　1.3



この結果より、16Kbyte の DMA 転送 と 2048 回の積和演算、32Kbyte の DMA 転送 と 4096 回の積和演算、
64Kbyte の DMA 転送 と 8192 回の積和演算が同じくらいといえる。

数字だけみれば、「DMA 転送：積和演算＝８：１」のようである。

今回は倍精度でのみ考えているが、単精度ではまた違った組み合わせになるだろう。

Multi Core CPU の資料。

■ Intel 日本語技術資料のダウンロード ■
http://www.intel.co.jp/jp/download/index.htm


■ ITmedia ■
新約・見てわかる パソコン解体新書：
第2回　Coreマイクロアーキテクチャ [前編]
http://plusd.itmedia.co.jp/pcuser/articles/0611/01/news005.html

新約・見てわかる パソコン解体新書：
第3回　Coreマイクロアーキテクチャ に迫る[後編]
http://plusd.itmedia.co.jp/pcuser/articles/0701/24/news002.html


■ マイコミジャーナル ■
【レポート】IntelのCore Microarchitecture
http://journal.mycom.co.jp/articles/2006/03/21/intelcore/index.html

Cell での BLAS (DGEMM) - その３。

前回は、DGEMM と言いながらα、βがそれぞれ１としていた。
もちろん処理もその分少なくなってしまう。

今回はα、βの演算、さらには SIMD も用いている。まずは実行時間だが、SIMD による高速化の効果が表れている。また、行列Ａ、行列Ｂを２行ずつ扱うことで、少しだが DMA 転送に対する処理の割合を大きくすることにした。その結果、前回に比べ、実行時間が半分ほどに改善されている。

２行ずつ扱ってはいるが、Double Buffer モデルのように、DMA 転送完了待ちと、処理を並列しているわけではないので、まだまだ早くなりそうである。

gcc　　　　　　　　　　　　　　　[msec]
SPE　512X512　1024x1024　2048x2048
　1　　　370　　　 2750　　　　21500
　2　　　215　　　 1550　　　　11500
　3　　　165　　　 1150　　　　 8500
　4　　　140　　　　900　　　　 6800
　5　　　125　　　　800　　　　 6050
　6　　　120　　　　720　　　　 5450

Cell での BLAS (DGEMM) - その２。

DGEMM のバグを探すことができた。やはり Cell プログラミングでは、デバッグが容易ではない。いい加減デバッグツールを使わないとまずいのかもしれない。。

ちなみに今回のバグは調子に乗って２行ずつ演算を行っていたために起きたのだが、その部分を１行ずつにするだけで正常に実行できた。しかし、実行時間がかなり増えてしまっている。

また、xlc コンパイラも試してみた。
以下は 512x512、1024x1024、2048x2048 においての xlc/gcc の性能比較である。
もちろん Double Buffer は用いていない。これにより、xlc は gcc よりも Cell に最適化され、かなりの速度向上が見られるとした期待を裏切ることとなった。しかしまだまだなにも最適化しているわけではないので、Double Buffer 等の改善を加えた後もう一度比較してみるとよいだろう。

ちなみに Makefile に次のように追加するだけである。
　PPU_COMPILER = xlc
　SPU_COMPILER = xlc


xlc -O3 / gcc -O3　　　　　　　　　　　[msec]
SPE　　512x512　　1024x1024　　　2048x2048
　1 　1148/795　　8700/6010　　16100/47400
　2　　600/430　　4600/3170　　36300/24500
　3　　425/310　　3150/2240　　24500/16900
　4　　330/245　　2420/1730　　18800/13200
　5　　280/210　　2020/1460　　15500/11000
　6　　245/190　　1710/1270　　13200/ 9700

Cell での BLAS (DGEMM) - その１。

Cell での BLAS 、特に DGEMM（Ｃ←αＡxＢ+βＣ）を作成しているが、なかなかうまくいかない。以前は行列サイズ 2048x2048 で固定していたが、それでは使いづらいので動的にメモリ確保することにした。

PPE にてアラインされたメモリを動的に確保したいときは、次のように行う。
　posix_memalign(void **matrix_ptr, alignment, column_size * row_size * sizeof(double));

また、SPEは LS が 256Kbyte しかないため、静的に確保してしまった方がよいだろう。

行列を１次元の配列で表しているため、行ごとにはアラインされているが、列ごとには飛び飛びになってしまっている。
DMA 転送ではアラインされていないといけないので、LS に行列が乗ってしまう場合をのぞいては、
あらかじめ、Ｂ行列（Ｃ←αＡxＢ+βＣ）を転置するなどといった補正をかけることになる。
そうすれば、「行列Ａの第ｉ行と行列Ｂの第ｊ列の内積」という演算は
「行列Ａの第ｉ行と行列Ｂの第ｊ行の内積」とみることができる。

ある程度大きなサイズでないと、
SPE スレッドの起動、転置、DMA 転送などのオーバーヘッドが伴ってしまうため、
SPE に切り出すことは意味のないことになってしまう。


xlc、double-buffer 等は用いていないので、さらに高速化は可能だが、うまく動いているかものすごく怪しい。
以下はその実行時間である。

Play Station 3
CPU : Cell 3.0GHz
Mem : 256Mbyte
Compiler : ppu-gcc spu-gcc

2048x2048 : 1250 [msec]
1024x1024 : 190 [msec]
512x 512　 : 45 [msec]

GotoBLAS。

GotoBLAS という強力な BLAS があるので、使用してみた。
make は quickbuild.32bit や quickbuild.64bit があるので、とても簡単に行える。

環境に応じて
./quickbuild.32bit　または　./quickbuild.64bit
とすればよい。

すると libgoto.a というライブラリが生成される。

使用する際には、ライブラリのパスや -lgoto -lpthread といったオプションを付けることになるが
毎回コマンドで打つのは少々大変なので、以下のような Makefile を使用すると便利である。


# Makefile
OUTPUTFILE = sample
SOURCES = sample.c
OBJECT = $(SOURCES:.c=.o)

CC = gcc
CFLAGS = -O3 -Wall
LDFLAGS = -L/usr/local/lib -L# ここにライブラリのパスを記述する
LIBS = -lgoto -lpthread

.PHONY: all
all: $(OUTPUTFILE)

$(OUTPUTFILE): $(OBJECT)
　　$(CC) $(LDFLAGS) $(OBJECT) $(LIBS) -o $(OUTPUTFILE)

$(OBJECT): $(SOURCES)
　　$(CC) $(CFLAGS) -c $(SOURCES)

.PHONY: clean
clean:
rm -f *.o
rm -f $(OUTPUTFILE)

Cell SDK 3.0。

ついに Cell SDK のヴァージョンが 3.0 になりました。
今回からは BLAS や Fortran Compilerの追加など、だんだん環境が整ってきた感がありますが、肝心の実機が Play Station 3 か、Blade Center しかないってのは、どうにかならないものかと思います。価格差が半端じゃないですからね。。。


3.0 からは RHEL に対応するらしく、CentOS で代用を考えているのですが、PS3 では荷が重いでしょう。実行に支障がでそうです。。とりあえず PS3 は Fedora Core 6 のままにして、新たにクロスコンパイル環境を構築予定です。


今回追加された BLAS は、行列のサイズが一定の条件を満たさないと使用できなかったり、他の BLAS とは演算形式が異なるので、同じように使うことは難しいと思われます。ここでは使用頻度が多いであろう LEVEL3 の DGEMM を例にとって見ていきます。

・通常の BLAS (dgemm)
行列Ａ(M x K)、行列Ｂ(K x N)、行列Ｃ(M x N)という行列に対し、Ｃ←αＡＢ＋βＣという演算を行う。
Ａ、Ｂを転置、共役転置するかどうか指定できる。

・Cell での BLAS (dgemm_spu)
Ｍは２の倍数、Ｎは４の倍数、Ｋは２の倍数でなければならず、演算もＣ←ＡＢ＋Ｃという演算を行う。
Ａ、Ｂを転置、共役転置しての演算には対応していない。

単精度演算(sgemm_spu)は、Ｍは４の倍数、Ｎは１６の倍数、Ｋは４の倍数でなければならないと、更に厳しくなります。
これに漏れてしまう行列はわんさかあると思いますし、まだまだ Cell BLAS を作る価値ありそうです。やるしかない！