「Cell Speed Challenge 2008 規定課題部門」 2008 年 6 月

Cell スピードチャレンジ 2008 で用いた高速化手法
須藤 郁弥 †

本田 耕一 †

松田 健護 †

坂内 恒介 †篠原 歩 ‡† 東北大学工学部 電気情報・物理工学科

柴田 悠希 †小林 隼人 ‡石野 明 ‡‡ 東北大学大学院 情報科学研究科

Cell スピードチャレンジ 2008 規定課題部門では，連立一次方程式求解の高速化に取りくんだ．本開発では求解アルゴリズムとして LU 分解法を用いた．また，高速化として LU 分解のブロック化，SIMD 演算の利用，ループアンローリング，DMA ダブルバッファリング，帯行列 LU 分解の計算量削減などを行った．結果
として，係数行列が密行列で与えられた場合には 7 SPE 利用時に最大 77 Gflops の性能を実現できた．これは，ピーク性能 179.2 Gflops（動作周波数 3.2 GHz）の 43% にあたる．さらに，帯行列 LU 分解では密行列と比較して 2 倍から 5 倍の高速化に成功した．

1規定課題概要

kkk+32nCell スピードチャレンジ 2008 の規定課題は連立一次方程式の求解であった．コンテストでは，与えられた係数行列 A（n × n 行列）および右辺行列 B（n × m 行列）にkk+32

nk+32

対し，AX = B を満たす解行列 X を高速に求めることがn課題であった．係数行列 A は行方向に，右辺行列 B は列n

Step2

Step1

方向に要素が格納されて与えられ，行列の各要素は floatk型であった．また，問題の制約として，入力データは利用可能であるメインメモリの容量 80 MB の半分の 40 MB

kk+32

nk+32k ← k + 32以下であり，n は 32 の倍数と定められていた．

2kn求解アルゴリズム

Step3

我々のプログラムでは求解アルゴリズムとして LU 分
図 1: 縦ブロックガウス法の処理解法を採用した．処理の流れは以下の通りである．

1. 係数行列 A の LU 分解
LU 分解により，A = LU を満たす下三角行列 L および上三角行列 U を求める．LU 分解は right-looking

3.1

LU 分解には，縦ブロックガウス法 [1] と呼ばれるブロック化アルゴリズムを用いた．縦ブロックガウス法では，データを列方向に分割し，定めたブロック幅ごとに法 [1] により行った．

2. 前進代入LU 分解の結果から，前進代入により LY = B となる行列 Y を求める．

LU 分解を行う．今回，問題サイズ n が 32 の倍数という制約から，ブロック幅は 32 とした．図 1 は (n − k) × (n − k)

3. 後退代入
後退代入により U X = Y を解き，解行列 X を得る．

3

LU 分解のブロック化

領域 (k = 0, 32, ..., n − 32) を更新する際の処理を示したものである．各ステップで行う処理を次に示す．

• Step1
Step1 では更新する領域の左 32 列を先行的に LU分解する．計算は行方向にサイクリック分割して各LU 分解の高速化
本節では，LU 分解に対して行った主な高速化についてSPE へ割り当てた．このとき，問題サイズの制約から各 SPE の担当領域はローカルストア (以下 LS)述べる．実装にあたっては [2, 3, 4] を参考にした．

1

に格納可能なサイズとなるため，Step1 におけるメ

for(i=0;i<n;i+=2){c = spu_nmsub(ai+1, bi+1, c);c = spu_nmsub(ai, bi, c);}インメモリとのデータ転送は Step1 の開始時および終了時のみでよい．また，SPE 内に列データを全て保有しているため，SPE 間の通信のみでピボッティングを実現できる．

• Step2 および Step3
Step2，Step3 では残りの領域を更新する．計算はfor(i=0;i<n;i+=2){c = spu_nmsub(ai, bi,spu_nmsub(ai+1, bi+1, c));}

列方向にブロック分割して各 SPE へ割り当て，一回あたりのデータ転送量がなるべく大きくなるようにした．ここで，Step2 の計算結果は LS に格納可

図 2: SIMD 演算の入れ子化能なサイズとなる．この計算結果は Step3 の計算に用いるため，LS に記憶しておくことでこの分のデータ転送を削減できる．また，ブロック化により*(vector float*)&a[k] =spu_nmsub(a1, b1,*(vector float*)&a[k]);*(vector float*)&a[k+4] =spu_nmsub(a2, b2,*(vector float*)&a[k+4]);

更新が 32 段ずつ進むようになるため，データ転送回数は 1/32 にできる．さらに，多段の情報を用いた更新となるため，演算が 2 重ループから 3 重ループになり，後に述べるループアンローリングの効果が大きく表れた．

3.2

vector float t1 , t2 ;t1 = spu_nmsub(a1, b1,*(vector float*)&a[k]);t2 = spu_nmsub(a2, b2,*(vector float*)&a[k+4]);*(vector float*)&a[k]= t1;*(vector float*)&a[k+4] = t2;

SIMD 演算の利用

Cell では 16 バイト長のベクタデータに対する演算を 1命令で実行する SIMD 演算が利用できる．行列の各要素が float 型（4 バイト）であるため，SIMD 演算により 1命令で 4 要素の処理を行うことができる．プログラム中の演算に SIMD 演算を用いることで，高速化ができた．

3.3

図 3: vector float 型の一時変数の利用

ループアンローリングプログラム中でループアンローリングを行うことで，高速化の効果が見られた．さらに，本開発ではループア

vector float 型の一時変数に SIMD 演算の結果を格ンローリングを行った際に，場合に応じて以下の操作を納してから，元のベクタに戻すようにすることで，行った．

展開分の演算が高速化できた．

1. SIMD 演算の入れ子化
ループ中の SIMD 演算で更新されるベクタが一定

3.4である場合，SIMD 演算を入れ子にした．この操作

DMA ダブルバッファリング

DMA 転送が必要な部分では DMA ダブルバッファリングを用い，演算中に次に必要な要素の読み込みや結果の書き込みを行うことで通信時間を隠蔽した．

により，多くの問題で解精度の向上が見られた．ここで，図 2 中の a, b, c はメモリ上のベクタデータとする．

2. 一時変数の利用
配列 a に対し，*(vector float*)& a[k] と記述する

3.5

帯行列 LU 分解の計算量削減

ことで，配列中の連続した 4 要素をベクタとして我々は帯幅の小さい帯行列について，計算する必要の扱うことができる．この場合，図 3 に示すようにない領域が現れることに注目し計算量の削減を図った．

2

kk+32+2w

k+32

for(j=0;j<i;j++)b[i] -= a[i][j]*b[j];

k
Step2
k+32
Step1

0temp = b[i];for(j=0;j<i;j++)temp -= a[i][j]*b[j];b[i] = temp

Step3
k+32+w0

図 5: スカラ変数の利用

前進代入・後退代入の高速化

4
図 4: 帯行列 LU 分解の計算量削減

前進代入・後退代入においても SIMD 演算や DMA ダブルバッファリングなどの基本的な高速化を行った．本帯行列の性質を利用した計算量削減を行うため，LU 分解を開始する前に次の手順で行列の半帯幅 w を調べた．
半帯幅 w は各対角要素の上下左右に並ぶ非零の要素数が節では，特に前進代入・後退代入に対して効果が見られた高速化手法について述べる．

w 個以下であることを保証する．
1. w = 0 で初期化する．

4.1

解ベクトルの複数本処理

2. SPE の id を用いて i = id として，第 i 行について，行の左端および右端から数えて最初に非零の要素が現れる位置をそれぞれ調べ，対角要素の左に並ぶ非度に 3 本ずつ解ベクトルを読み込んで処理を行うように零の要素数 wLi および右に並ぶ非零の要素数 wRi

した．このとき，生じるループに対して，LU 分解と同様を求める．のループアンローリングを行い，さらに高速化を図った．

係数行列の読み込み回数を削減するため，各 SPE で一

3. 各 SPE のもつ max(wLi , wRi ) のうち最大のものをw とする．閾値 σ に対し，w > σ ならば，密行列4.2として LU 分解を行う．

4. 全ての行を調べ終わったら半帯幅として w を出力する．そうでないならば i = i + 7（7 SPE 利用時）

スカラ変数への格納

前進代入および後退代入では，配列の 1 要素に繰り返しアクセスを行った．この際，図 5 のように配列の要素として，2 へ戻る．

をスカラ変数に格納しそのスカラ変数に対してアクセスを行うことで，処理時間を削減する事ができた．

w を用いると，LU 分解の各ステップで計算が必要な領域は図 4 のように表すことができる．Step 1 において，ピボット交換は第 k 行から第 (k + 32 + w) 行の間で行わ4.3れるため，ピボット交換後も第 (k + 32 + w) 行以降の要右辺行列の行交換

素はすべて 0 である．また，範囲中の最下行とのピボット交換が起こった場合でも，Step 2 の領域で非零の要素右辺行列は要素が列方向に配置されているため，行交は第 (k + 32 + 2w) 列までに収まる．さらに，Step 3 の計換を行いたい要素を取り出すためにはその周辺の要素ご算では Step 1 および Step 2 の領域を参照するため，そと取り出す必要があり，行交換の効率は大変悪い．そのれぞれで考慮した行・列についてのみ計算すれば十分でため，係数行列で行ったピボット交換の操作を保管し，前ある．

進代入・後退代入の処理時に一括して解行列のピボット交換を行った．

以上のように帯行列 LU 分解の計算量を削減できる．

3

の担当領域をブロック分割により決めることで，一度に表 1: 決勝ラウンド結果なるべく大きいサイズを読み込むようにし，データ転送

問題

nm

時間

順位行列の性質の効率化を図った．このため，問題 3 のように n が大き12768224

1896

148

3位3位非対称密行列い場合には十分な効果が得られたが，n が小さい問題で非対称密行列は分割により各 SPE の担当領域も小さくなるために効率

342,688256

120,000

186120

2位4位非対称密行列が落ちてしまったと考えられる．対策としては，行方向非対称密行列にもブロック幅を考慮して 32 × 32 のブロック単位で処5671,792
322,720

1131,0721

74 3 位
64 4 位FAILED

対称密行列理するようにし，係数行列のデータレイアウト変換を行対称密行列

うことが考えられる．

89102,240
2,656
3,200

224
100
7113
212
62

2位
3位
1位

三重対角行列帯行列

6対角優位行列
ブロック行列LU 分解法による連立一次方程式求解の高速化を行った．LU 分解のブロック化などによりデータ転送の回数を削減した上で DMA ダブルバッファリングを行い，また時間は msec

5まとめ

評価結果

演算に SIMD 命令やループアンローリングを用いることで，様々な問題に対して高速に解を出力するプログラム以上の高速化を行ったプログラムを用いて，n = 3200，を作成することができた．また，係数行列が帯行列であっm = 1 として係数行列が密行列である場合の求解を行った．この場合の性能は 77 Gflops であり，動作周波数 3.2た場合には効率よく計算量の削減を行うことができた．

GHz，7 SPE 利用時のピーク性能 179.2 Gflops の 43% の性能を実現できた．また，Cell スピードチャレンジ 2008

参考文献

規定課題部門本選の結果は表 1 の通りであった．問題 7 で[1] 規定課題問題作成 WG． マルチコアプログラミング

は，[5] を参考に実装したアルゴリズムを用いたところ誤差が大きくなってしまい正しい解を出力できなかった1 が，コンテスト Cell スピードチャレンジ 2008 規定課題参考資料，2008．

それ以外の全ての問題で 4 位以内に入ることができた．さ[2] 東芝セミコンダクター社 Cell オンラインチュートリ
アルらに，係数行列のサイズ n が大きい問題（問題 3）および，係数行列が帯行列で与えられた問題（問題 8，10）では，いずれも 2 位以内となった．

[3] Sony Computer Entertainment Inc. SPU C/C++ 言問題 8，問題 10 に対し，計算量削減を行わなかった場

語拡張 Version 2.3, 2006．

合の求解時間はそれぞれ 192 msec，284 msec であった．

[4] FIXSTARS Cell プログラミングチュートリアル

よって，計算量削減により，求解時間は問題 8 ではおよそ半分に，問題 10 では 1/5 にまで短縮できたことが分か[5] W. H. Press, S. A. Teukolsky, W. T. Vetterling, andB. P. Flannery．丹慶勝市ほか訳．Numerical Recipes
in C．技術評論社，pp．62-63，1993．

る．また，提出したプログラムでは先に述べた領域よりも余分に計算を行っていたため，実際には求解時間をさらに短縮できる．具体的には，問題 8 で 103 msec，問題10 で 42 msec に短縮できることが確認できた．ちなみに，
問題 7（三重対角行列）に対し同方法を用いた場合の求解時間は 28 msec であり，さらにこの場合は十分な解精度が得られた．計算量削減を行わない場合の求解時間は186 msec であったため，この場合も計算量を大幅に削減できていることが分かる．
本開発では，LU 分解の Step 2，Step 3 において各 SPE
1 練習問題，予選問題の三重対角行列に対しては正しい解を出力した

4