第 28 回数値流体力学シンポジウム
E06-2

係数行列のビット圧縮による疎行列反復解法の高速化
Performance improvement of iterative solver using bit-compression for a sparse matrix⃝小野謙二，理研 計算科学研究機構, 神戸市中央区港島南町 7-1-26, E-mail：keno@riken.jp
Kenji Ono, RIKEN AICS, 7-1-26 Minatojima-minami-cho, Chuo-ku, Kobe, Japan

A nobel Bit-representation/compression technique is proposed to enhance the performance of iterative methods fora large-size sparse matrix. This technique is applied to the implementation of iterative kernels with Dirichlet andNeumann boundary conditions. The first advantage of this approach is that it reduces memory traﬃc from mainmemory to Core and eﬀectively utilizes SIMD units with cache. Secondly, the proposed implementation can replace
if-branch statements with mask operations using the bit expression. This promotes the optimization of code duringcompilation and run-time. The Red-Black SOR and BiCGstab algorithms are employed to investigate the proposedimplementation. Consequently, the proposed approach achieves 3.5 times faster than a naı̈ve implementation on
both Intel and Fujitsu Sparc architectures.

はじめに非圧縮性流体の数値解法には圧力の Poisson 方程式が現れ，離散化により疎行列を係数行列にもつ連立一次方程式が導かれる．従来より，大規模疎行列の連立一次方程式の解法として様々な反復計算方法が提案されているが (2) (3) (4) (5) ，近年の計算機資源の増大に伴い扱える問題規模も大きくなり，より高速な解法や実装が期待されている．
演算性能は，対象とするアーキテクチャ上で利用するアルゴリズムとその実装によって大きく変化する．近年のスーパーコンピュータ（スパコン）は，超並列プロセス，マルチコア，SIMD 演算器などの特徴を備えているものが主流となっている．また，今後のスパコンの開発は，一定の消費電力内での高性能化の方向に向かっており，コアあたりのメモリ実装量が少なくなる可能性もあり，省メモリなプログラミングも再考すべき点である．
あるアーキテクチャ上でのアルゴリズムの性能予測に
Roofline モデル (1) が提案され，性能の概要を説明できる簡便な方法として利用されている．Roofline モデルでは，Operational Intensity が大きければそのアーキテクチャの性能を十分に引き出せ，理論性能に近い演算性能となる．言い換えると，プログラムのループ中の演算数と必要な変数のロードの比と，アーキテクチャの演算処理性能とメモリバンド幅の比によって，実現可能な性能が予測できる．このモデルは疎行列計算の場合には，メモリからコアへのデータ転送がボトルネックとなり，演算性能が律速となることをうまく説明している．また，性能を改善するためには，Operational Intensity を大きくすることが必要であること述べている．このためには，演算数を多くする（高コストであるが高収束の反復アルゴリズムを用いる）か，メモリからのデータ転送量を削減する指針が得られる．本論文は，後者のアプローチにより，反復計算の性能を改善するものである．

1.

基礎方程式と係数行列のビット表現
三次元の非圧縮性 NavierStokes 方程式の分離解法で用いられる Poisson 方程式 (1) を構造格子上で離散化して考える．

2.

∇ (∇p) = div (∂u) ≡ ϕ,∂tは，Neumann 条件か Dirichlet 条件として扱えるので，この二種の境界条件を計算中に自動的に取り込めるように，スキームを構築する (7) ．
まず，準備として，境界条件を識別する離散 Heaviside
関数 (2) を導入する．
0 (Boundary Condition)
H =(2)
1 (F luid)
これにより，計算式中の圧力 p は，離散 Heaviside 関数
により通常の圧力と境界条件としての圧力 pBC を式中に取り込む．

p = p H + (1 − H) p BC

したがって，式 (1) の半離散式は，式 (2) と式 (3) により，次式のようになる．

∑(∇p H N)

nl = hϕ −ll∑()1 − H lN ∇p BCnlll(4)ここで h，H N ，n l は，それぞれ，格子幅，Neumann
条件を識別する離散 Heaviside 関数，セルの l 方向の法線ベクトルを表す．同様の実装は，プログラム中においてif 文を使っても記述できるが，式 (4) の実装の方がよりSIMD ユニットの最適化の促進が期待できる．Dirichlet
条件についても，同様に識別離散 Heaviside 関数 H D を導入することにより，スキーム中で取り込める．例えば，Fig. 1 において，セルの e 面の圧力勾配は，次式のように表現できる．

(∇p H N)e(1)

ここで p，u，ϕ はそれぞれ，圧力，速度，Poisson 方程式のソース項を表す．この方程式を直交等間隔格子上で離散化して得られる代数方程式から導かれる連立一次方程式について考える．ステンシル計算であれば，直交格子に限らず，同様の方法が適用可能である．式 (1) に対して二次精度の差分近似を行うと，7 点ステンシルのスキームが得られる．
実用的なシミュレータでは，計算領域内の任意の部分に配置された物体を扱う必要があるため，壁面境界条件を効率的に扱うことが要求される．多くの圧力境界条件(3)=1{p i+1 H eDh (} N)+ 1 − H eD p BCi+1 − p i H e(5)最終的に，連立一次方程式 Ax = b の形式は式 (6) となる．
∑∑()
H lN
pH D H N l − pl )l∑(nl
= h2 ϕ − h
1 − H lN ∇p BCl(6)l∑()H lN n l−1 − H lD p BCll

1c 2014 by JSFM
Copyright ⃝第 28 回数値流体力学シンポジウム
E06-2∇pw = 0∇pe =n!

w
i-1+ c_t * p(i ,j ,k+1) + c_b * p(i ,jdp = ( (ss + b(i,j,k) ) / dd - pp ) * omg
p(i,j,k) = pp + dp
res = res + dble(dp*dp)* dble( ibits(bp(i,j,k), Active, 1) )end doend doend doend dop BCi+1 − pihp BCi+1eii+1s,k-1)[Bit-reps code]

Fig. 1: Neumann and Dirichlet boundary conditions forcell i in two dimensions. A Neumann BC is applied at
the west cell face, which is solidly shaded. A Dirichlet
BC is employed at the east cell face, where the boundary
value is given by the pressure pBC
i+1 .

全ての境界条件の寄与は，式 (4) の右辺にまとめられている．左辺の係数行列に現れる係数の要素は，例えば，対角項では 0 ∼ 6 である．0 はセルの 6 面全てが壁面で特殊な場合であり，計算しない．非対角項の係数は 0 または1 である．つまり，非対角項の係数は 1 ビットで表現でき，対角項も 3 ビットあればよい．したがって，あるセルの反復式に必要な係数は合計 9 ビットあれば十分である．実際の実装においては，Fig. 2 に示すように，4 バイトの整数配列にその他の計算に必要なフラグ情報と合わせて格納している．BC Diag と BC Ndag x が対角項と非対角項の係数を表している（x は方向を示す）．BC N x と BC D xは，セルの各面に対する Neumann 条件と Dirichlet 条件を示す．
これらのビット情報は，計算が始まる前の初期化段階で計算し，配列のエンコードして保持しておく．エンコー
ドとデコードは，単純にシフト演算を用いて，次のようにコーディングしている．

Encoding;
inline int onBit (int idx, const int s) {return ( idx | (0x1<<s) );}Decoding;#define BIT_SHIFT(a,b) ( (a >> b) & 0x1 )
ビット圧縮を用いた反復法の実装
Red-Black SOR 反復法について，ビット表現による実装を示す．
2 つの RB-SOR の実装を示す．ひとつは，ナイーブな実装で係数行列を浮動小数点の配列 pn(i,j,k,n) に保持する方法である．もうひとつは，本提案であるビット表現による実装方法である．下記に，Fortran コードを示す．

3.

[Naive code]
do color=0,1
do k=1,kx
do j=1,jx
do i=1+mod(k+j+color,2), ix, 2
c_w = pn(i,j,k,1)
c_e = pn(i,j,k,2)
c_s = pn(i,j,k,3)
c_n = pn(i,j,k,4)
c_b = pn(i,j,k,5)
c_t = pn(i,j,k,6)
dd = pn(i,j,k,7)
pp = p(i,j,k)
ss = c_e * p(i+1,j ,k ) + c_w * p(i-1,j ,k+ c_n * p(i ,j+1,k ) + c_s * p(i ,j-1,kdo color=0,1
do k=1,kx
do j=1,jx
do i=1+mod(k+j+color,2), ix, 2
idx = bp(i,j,k)
c_e = real( ibits(idx, BC_Ndag_E, 1) )c_w = real( ibits(idx, BC_Ndag_W, 1) )c_n = real( ibits(idx, BC_Ndag_N, 1) )c_s = real( ibits(idx, BC_Ndag_S, 1) )c_t = real( ibits(idx, BC_Ndag_T, 1) )c_b = real( ibits(idx, BC_Ndag_B, 1) )d0 = real( ibits(idx, BC_Diag+0, 1) )d1 = real( ibits(idx, BC_Diag+1, 1) )d2 = real( ibits(idx, BC_Diag+2, 1) )dd = d2*4.0 + d1*2.0 + d0
pp = p(i,j,k)
ss = c_e * p(i+1,j ,k ) + c_w * p(i-1,j ,k )+ c_n * p(i ,j+1,k ) + c_s * p(i ,j-1,k )+ c_t * p(i ,j ,k+1) + c_b * p(i ,j ,k-1)
dp = ( (ss + b(i,j,k) ) / dd - pp ) * omg
p(i,j,k) = pp + dp
res = res + dble(dp*dp)* dble( ibits(idx, Active, 1) )end doend doend doend do

コード中の配列 p, b, bp, pn は，それぞれ，圧力，式(6) の右辺，ビット情報配列（Fig. 2），係数行列を示す．
ibits(bp(i,j,k), Active, 1) は，計算セルの残差計算が有効（流体）か無効（壁面内）であるかを識別するフラグである．
両実装を比較する．まず，必要な配列（メモリ要求量）については，ナイーブ実装では係数行列 pn とソース項 b，ビット情報 bp，圧力 p の 10 単位である．一方，ビット表現実装では，p, b, bp の 3 単位である．したがって，ビット表現実装はナイーブ実装に比べて約 1/3 のメモリ要求量となる．
次に，Operational Intensity（Flop/Byte）について調べる．ナイーブ実装では, pn, b, bp で 9 つのロードが発生する．加えて，圧力 p については，変数のアドレス値を調べると，一度のロードで i-1, i, i+1, j-1, j+1 のアクセスがキャッシュ内に入ることがわかった．k-1, k+1のアクセスについては別途ロードが発生するため，都合12 回のロードと 1 回のストアがループ中で生じる．一方，ビット表現実装では，b, bp の 2 つのロード，圧力 p は一度のロードで i-1, i, i+1, j-1, j+1 のアクセスまでキャッシュ内に入っていたので都合 5 回のロードと 1 回のストアがループ中で生じる．
浮動小数点の演算量は，ナイーブ実装は 26，ビット表
現実装では 30 である 1 ．これらのデータから Flop/Byte
値が，Table 1 のように計算できる．
1
ここで，足し算，引き算，かけ算は 1flop，割り算は 8flops とした（Sparcアーキテクチャの演算カウントに基づく）．

))

2c 2014 by JSFM
Copyright ⃝第 28 回数値流体力学シンポジウム
E06-20911 122329 30 31
ActiveState

BC_D_T
BC_D_B
BC_D_N
BC_D_S
BC_D_E

BC_D_W
BC_N_T
BC_N_B
BC_N_N
BC_N_S
BC_N_E

BC_N_W
BC_Ndag_T
BC_Ndag_B
BC_Ndag_N
BC_Ndag_S

BC_Ndag_E
BC_Ndag_W

BC_Diag(0~6)

17 18Fig. 2: Bit representation. Several bits required for the bit-representation are encoded into this array. This example
includes diagonal(BC Diag), non-diagonal (BC Ndag x), Neumann boundary(BC N x), Dirichlet boundary (BC D x),
cell state (State), and activeness (Active) of a cell. Other bits are used for more complicated processes.

Tab. 2: Specification of evaluation machines. TRIAD scores are measured by the STREAM benchmark(11) .
Architecture

ClockCore

CPU(GHz)
Xeon X5650
Xeon E5-2670
Xeon E5-1680
Sparc VIIIfx
Sparc IXfx

PeakCache

MemoryTheoreticalTRIAD(GFlops/s)

(MB)

(GB)

BW (GB/s)(GB/s)

2.66
2.6
3.0
2.068882211127.7
166.4
96.0
128.0

122025
6

1664
64
1664
102—64225949
36

1.85

16

1236.512328550

Tab. 1: Comparison of characteristic for two types ofimplementation.

Memory Requirement
Load & StoreArithmetic
F/B

Naı̈veBit-Reps.

10 unit
12+1
26
0.5

3 unit5+130
1.25

IXfx（FX10）の性能の違いは，コンパイラの最適化の相違である．textitFujitsu Venus VIIIfx のコンパイラには，ibits(a, b, 1) のシフト演算を SIMD を利用する最適化が施されているが，Fujitsu Venus IXfx のコンパイラにはその実装がされていないためである 2 ．以上により，ビット表現により，使用メモリ量を抑えることでキャッシュの有効活用ができ，性能劣化のタイミングが遅れることがわかった．また，異なるアーキテクチャで提案手法が有効に働くことがわかった．
次に，Fig. 6 にスレッド並列時の性能の測定結果を示す．問題のサイズはキャッシュサイズより大きな 256 を選んだ．Intel アーキテクチャでは，ナイーブ実装で 4 スレッド以上のときメモリバウンドの傾向が見られる．ビット表現実装は，より性能が高いがコア数が増えると，やはりメモリバウンドの傾向は見える．このとき，最大 3倍以上の性能改善を果たしている．
Fujitsu Venus では，VIIIfx と IXfx で異なる傾向を示す．VIIIfx では，ナイーブ実装で Intel の場合と同様に 6スレッド以上でメモリバウンドである．ビット表現実装
では，8 スレッドまでリニアに性能向上し，最大 2.5 倍の性能改善を果たしている．IXfx では 16 スレッドまでスケールしているが，スレッド数が増えたときにビット表現実装の方がより性能が高い．
Table 1 に示す反復カーネルの分析と実測値の結果をFig. 3 に示す．図中，実線，破線，一点鎖線は各アーキテクチャの理論性能を示す．ナイーブ実装は F/B=0.5 程度で，各アーキテクチャの予測性能から少し低い実測性能となっているが，SIMD の効果はある．ビット表現実
装では傾向は異なり，Westmere(X5650) と Sparc VIIIfx
は F/B=1.25 の値に応じて性能改善が見られる．一方，Sparc IXfx は SIMD 最適化が効いていないので性能改善代は小さい．
以上のように，係数行列の係数をビット表現で圧縮保持する方法は，メモリロード量を抑制するとともに省メモリ効果により演算性能を 3 倍以上改善できることが明らかとなった．

性能評価
前節に示したナイーブ実装とビット表現実装について，Table 2 に示す計算機で性能を評価した．
評価する問題は，式 (1) と同じ性質をもつ三次元の熱伝導問題 ∇2 φ = 0 を Dirichlet/Neumann 境界条件で解く問題 (2) ，および三次元のキャビティフロー問題を扱った．計算は単精度計算である．演算性能の測定には，Performance Monitor library (PMlib) (8) , (10) を用いた．
PMlib はユーザが宣告した演算数と処理前後のタイミン
グ測定データから演算速度を推定する．この方法はアルゴリズムレベルの演算数の推定を元にしているため，実行時の精確な浮動小数点演算数ではないが，アーキテクチャ間の性能比較の点では有効な方法と考えられる．また，PAPI によってハードウェアカウンタの値を取得し，性能を推定するモードも提供する．
Fig. 5 に逐次実行時の測定結果を示す．Fujitsu Venus
IXfx を除き，ビット表現実装の方がナイーブ実装と同等か高い性能を示すことがわかる．特に，問題サイズが64 より大きく（L3 キャッシュサイズ以上）なったところ
で，ナイーブ実装は性能劣化が見られるが，ビット表現実装ではその性能劣化が遅れている．これは，ビット表現実装はアルゴリズムを省メモリで実装できているので，キャッシュを有効に使えているためと考えられる．また，
Fujitsu Venus VIIIfx（京コンピュータ）と Fujitsu Venus

2
ちなみに，Intel コンパイラは ibits(a, b, x) も最適化するが，京のコンパイラは現時点で x=1 のモードのみである．

3c 2014 by JSFM
Copyright ⃝4.

第 28 回数値流体力学シンポジウム
E06-2

参考文献

Attainable Performance (GFlops/s)

最後に，Fig. 4 にビット表現実装を用いた三次元非圧
縮熱流体コード FFV-C(9) の京コンピュータでの測定結果を示す．京コンピュータは 10PFlops であるが，FFV-Cコードは約 12%の実効性能を達成している．

(1) S. Williams, S., Waterman, A. and Patterson, D.:
Roofline; An Insightful Visual Performance Modelfor Multicore Arch. Commun. ACM, Vol.52 No.4
(2009) 65–76

Westmere
Sparc VIIIfx
Sparc IXfxWestmere bit-repsWestmere naive
Sparc VIIIfx bit-reps
Sparc VIIIfx naive
Sparc IXfx bit-reps
Sparc IXfx naive100(2) Yokokawa, M. : Vector-Parallel Processing of theSuccessive Overrelaxation Method. Japan Atomic
Energy Research Institute JAERI-M Report No. 88017 (1988) in Japanese(3) Willcock, J. and Lumsdaine, A.: Accelerating
sparse matrix computations via data compression.
Proc. 20th Annual ICS ’06 (2006) 307–316

10(4) Tang, W. T., et al.:Accelerating Sparse
Matrix-vector Multiplication on GPUs Using Bitrepresentation-optimized Schemes. Proc. of SC1326 (2013) 1–12
1!"!#$!

!"!#!

!

(5) Van der Vorst, H. A. : Bi-CGSTAB: A Fast andSmoothly Converging Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems. SIAM J.
Sci. and Stat. Comput. Vol.bf 13 No.2 (1992) 631644

Operational Intensity (Flops/Byte)

Fig. 3: Performance analysis of Roofline model for naı̈ve
and bit-reps implementation.

(6) Ono, K. and Kawashima, Y. : Multicolor SOR
Method with Consecutive Memory Access Implementation in a Shared and Distributed Memory
Parallel Environment. Lecture Notes in Computational Science and Engineering, Vol.74 (2010) 183–191

Ideal
FFV-C

1x106

(7) Ono, K., Chiba, S., Inoue, S., and Minami, K. :Performance Improvement of Iterative Methods using a Bit-Representation Technique for Coeﬃcient
Matrices. Vecpar 2014, (2014)

GFLOPS

1x105

1x104

(8) Ono, K., Kawashima, Y. and Kawanabe, T.:
Data Centric Framework for Large-scale Highperformance Parallel Computation. Procedia Computer Science, Vol.29 (2014) 2336–2350

1x103

1x102

(9) http://avr-aics-riken.github.io/ffvc\_package/

1x101
1

1x101

1x102

1x103

1x104

1x105

(10) http://avr-aics-riken.github.io/PMlib/

Number of Processes(11) http://www.cs.virginia.edu/stream
Fig. 4: Measured performance of FFV-C code on the K
computer with 82,944 nodes. Each node has 8 cores.

まとめ
構造格子系の圧力 Poisson 方程式を離散化して得られる連立一次方程式に対して，その係数をビット表現で圧縮する手法について，その方法と効果を示した．提案方法は，主記憶からコアまでのデータロードのコストを減少させることができ，かつ省メモリな点でどのアーキテクチャにも適用可能で，効果がある．提案手法を Intel お
よび Sparc のアーキテクチャで実験し，ナイーブな実装に比べて最大 3.5 倍の性能向上を観測した．本実装方法
は，SIMD を有効に活用でき，かつメモリ要求量が少ないため，将来のアーキテクチャにおける高速化にも非常に有望な方法である．

5.

謝辞 本研究の一部は，京コンピュータを利用して得られた．

4c 2014 by JSFM
Copyright ⃝第 28 回数値流体力学シンポジウム
E06-2
44bit
naive

bit
naive
3GFlops/s

GFlops/s

321021163264

1280256

1632Size(a) Intel Xeon X5650.

2564bit
naive

bit
naive

33GFlops/s

GFlops/s

128(b) Intel Xeon E5-2670.

42

1064

Size

21163264

1280256

Size

163264128256Size(c) Fujitsu Sparc Venus VIIIfx.

(d) Fujitsu Sparc Venus IXfx.

Fig. 5: Comparison of serial performance of each machine. The problem size varies ranging from 163 to 2563 .

5c 2014 by JSFM
Copyright ⃝第 28 回数値流体力学シンポジウム
E06-2303025bit
naive

20GFlops/s

GFlops/s

25151050bit
naive

2015105

12468100121246Core(a) Intel Xeon X5650.

12141630

2525bitnaive

20

GFlops/s

GFlops/s

10(b) Intel Xeon E5-2670.

3015151055

124608Corebit
naive

20

1008Core

1

2

4

68

10121416Core(c) Fujitsu Sparc Venus VIIIfx.

(d) Fujitsu Sparc Venus IXfx.

Fig. 6: Comparison of thread parallel performance of each machine. The problem size is chosen to 2563 so that thedata resides in main memory.

6c 2014 by JSFM
Copyright ⃝