AVX を用いた倍々精度疎行列ベクトル積の高速化 AVX Acceleration of

HPCS2013
2013/1/15
2013年ハイパフォーマンスコンピューティングと計算科学シンポジウム
High Performance Computing Symposium 2013
AVX を用いた倍々精度疎行列ベクトル積の高速化
菱沼利彰†1 藤井昭宏†1 田中輝雄†1 長谷川秀彦†2
計算性能の向上に伴い,高精度による計算が多くの場面で可能となっている.4 倍精度を効率良く実現する手法と
して,2 つの倍精度変数で 1 つの 4 倍精度変数を表現する倍々精度演算がある.本研究では,疎行列とベクトルの演
算に使われる基本演算を AVX 命令を用いて高速化し,性能を決定するパラメタについて分析を行うことにより,倍々
精度演算を AVX で行う際の効果を示した.
AVX 命令を用いた場合,同時演算数の増加,3 オペランド化によるレジスタ退避,復元処理の減尐などの効果が大
きく,SSE2 の性能と比べて,キャッシュに収まる範囲のベクトル間の演算では約 1.7 から 2.3 倍の性能となったが,
キャッシュに収まらない場合は,キャッシュアクセス,メモリアクセスが大きなボトルネックになることがわかった.
倍精度の疎行列と倍々精度のベクトルの積では,約 1.1 から 1.9 倍の性能となり,メモリアクセスはボトルネックと
ならず,疎行列の 1 行あたりの非零要素の数が性能に大きな影響を与えていることがわかった.これらの結果から,
倍々精度の疎行列ベクトル積の性能を予測する 1 つの指標を導出した.
AVX Acceleration of Sparse Matrix-Vector Multiplication in
Double-Double
Toshiaki Hishinuma†1, Akihiro Fujii†1, Teruo Tanaka†1 and Hidehiko Hasegawa†2
As computing performance is improved generation after generation, high precision computation becomes possible in many
situations. One of the efficient methods to perform quadruple precision is to use Double-Double precision which uses two double
precision variables for one quadruple precision variable. In this paper, the authors tuned basic operation kernels of sparse
matrices and vectors in Double-Double precision using AVX, and analyzed their performance.
The AVX speedup ratio of the Double-Double vector operations is from 1.7 to 2.3 when data stored in the cache. The reason of
performance acceleration is number of operations in the same time and elimination of backup and recovery values on registers by
three operands instruction, The AVX performance decreases when data not stored in the cache, because of cache hit ratio and
memory bandwidth. The AVX speedup ratio of the product of Double precision sparse matrix and Double-Double precision
vector is from 1.1 to 1.9. An average number of nonzero elements per row affects to the performance, but a memory bandwidth
does not affect to the performance. The authors define one metric to forecast the AVX performance of the product of sparse
matrix and vector in Double-Double.
1. はじめに
計算性能の向上に伴い,高精度で計算をすることが多く
の場面で可能となってきている.また,CG 法等のクリロ
フ部分空間法の収束性は丸め誤差に大きく影響されるため,
収束の改善を図るには高精度演算が有効である[1].
しかし,Fortran などに実装されている整数演算による
4 倍精度演算の実行には倍精度演算と比較してかなりの計
算時間がかかる[2].効率良く高精度計算をする手法のひと
一方,ハードウェアの進化により,Intel の Sandy Bridge
マ イ ク ロ ア ー キ テ ク チ ャ に Advanced Vector Extensions
(AVX)と呼ばれる SSE2 に代わる拡張命令が新たに導入さ
れた.
われわれは,Lisで倍々精度演算に用いられているSSE2
命令をAVX命令に置き換えることにより,ベクトル演算に
どのような性能の向上が生じるかについて,研究を進めて
きた[6].
本論文では,はじめに倍々精度演算の概要とその演算,
つに,倍精度変数を 2 つ用いて 1 つの 4 倍精度変数の値を
次にAVXの特徴と今回実験に用いる実験環境,最後に倍々
保持し,4 倍精度演算を実行する倍々精度演算という手法
精度を用いたベクトル演算と疎行列ベクトル積に対する数
がある[3].
値実験を行い,AVXによる倍々精度演算の実装の効果と有
反 復 解 法 ラ イ ブ ラ リ Lis[4] [5] で は , Intel の Single
用性について述べる.
Instruction Multiple Data (SIMD)拡張命令である Streaming
SIMD Extensions 2 (SSE2)を用いて倍々精度演算の高速化
が実装されている.
2. 倍々精度演算
倍々精度演算とは,Bailey が提案した”Double-Double”精
†1
工学院大学情報学部
Faculty of Informatics, Kogakuin University
†2
筑波大学図書館情報メディア系
Faculty of Library, Information and Media Science, University of Tsukuba
ⓒ 2013 Information Processing Society of Japan
度のアルゴリズム[3]を用い,倍精度変数 2 つを用いて 4 倍
精度変数を実装する手法である.これによって生成され
23
HPCS2013
2013/1/15
2013年ハイパフォーマンスコンピューティングと計算科学シンポジウム
High Performance Computing Symposium 2013
図 2
Figure 2
図 1
Figure 1
AVX の SIMD レジスタの構造
Architecture of AVX SIMD register
倍々精度のビット数
Bit pattern of Double-Double precision number
る”Double-Double”精度の浮動小数点変数を倍々精度と呼
ぶ.
"Double-Double"精度のアルゴリズムにおいて倍々精度
変数 a を a = a.hi + a.lo, ulp(a.hi) ≥ 2 |a.lo|(上位 a.hi と下
位 a.lo は倍精度)とする.なお,ulp(x)は x の仮数部の"unit
in the last place"を意味する[7].図 1 に倍々精度変数と
図 3
擬似 C コードによる SSE2 から AVX への置換の概要
Figure 3
Pseudo C code of SSE2 and AVX
IEEE754 準拠の 4 倍精度変数のデータ構造を示す.
倍精度の仮数部は 52bit であるため,倍々精度の仮数部
は 104bit となる.これは IEEE754 準拠の 4 倍精度の仮数部
112bit に比べて 8bit 尐なく,指数部も 11bit と 4bit 尐ない.
用し,ymm レジスタに包括されるような設計のため,SIMD
レジスタは 16 本である.
また,このような ymm レジスタの設計から,AVX 命令
しかし,精度としてはほぼ同様で,整数演算による四倍精
は,同一 ymm レジスタ上の上位と下位の 128bit 境界を越
度演算と比べ高速に実行することが可能である[2].
えての水平演算を行うことはできない.
倍々精度の四則演算は Dekker[8]と Knuth[7]の丸め誤差
AVX では,VEX プリフィックスというプリフィックス
のない倍精度の加算と乗算アルゴリズムに基づき,倍精度
方式が採用され,SSE2 が 2 オペランドの命令しか実行でき
の四則演算の組み合わせで実現する.
ないのに対し,3 または 4 オペランドの命令が実行できる
丸め誤差のない倍精度加算アルゴリズムは,加減算のみ
で成り立ち,演算数は 11flops である.丸め誤差のない倍精
度乗算アルゴリズムは,乗算 9 回と加減算 15 回によって構
成され,演算数は 24flops である.
このように,倍々精度演算は倍精度演算の組み合わせに
よって実行されるため,SIMD 命令を用いて複数の命令を
同時に発行することにより高速化が期待できる.
ようになり,レジスタ退避,復元処理の記述を省くことが
できる[9].
倍々精度演算において,SSE2 と AVX の主な違いは一命
令で倍精度変数をいくつ同時に演算できるかである.C
言語上では,アラインメントを意識する必要はあるものの,
それ以上の違いはない.
そのため,SSE2 向けのコードを AVX 向けのコードに置
き換える際には,同時処理数を変更するような改変を施す
3.
AVX の概要と実験環境
3.1 AVX の概要
AVX は SSE2 の後継となる SIMD 拡張命令である.SSE2
だけでよい.本研究では,Lis 内の SIMD 命令に対応する
ループ内でのインデックス計算,端数処理を変更し,配列
のアラインメントを合わせる作業を行った.その概要を図
3 に示す.
は浮動小数演算に使える 128bit の xmm レジスタと呼ばれ
る SIMD レジスタが 16 本使用できるが,AVX は 256bit の
3.2 実験環境
ymm レジスタと呼ばれる SIMD レジスタが 16 本使用でき
実験には,表 1 に示す計算環境を用いた.Intel Core i7
る.そのため,SSE2 は 2 つの倍精度のデータに対して SIMD
2600K は Sandy Bridge マイクロアーキテクチャであり,動
演算ができるが,AVX は 4 つの倍精度のデータに対して
作周波数は 3.4GHz,4 コア 8 スレッドの構造で, 1 コアあ
SIMD 演算ができる.厳密には AVX は今後同時処理数を増
たり乗算器と加算器を 1 つずつ持っているため,AVX を使
やせるような設計になっているが,今回の実験環境では
用し倍精度を 4 つ同時に演算した場合の倍精度演算の理論
256bit である.
ピ ー ク 性 能 は , 3.4G × 4(core) × 4(AVX) × 2(演算器) =
図 2 に AVX の SIMD レジスタの構造を示す.AVX がサ
ポートされている Sandy Bridge 以降のアーキテクチャにお
108.8GFLOPS , SSE2 を 使 用 し た 場 合 は そ の 半 分 の
54.4GFLOPS となる.
いて,xmm レジスタと ymm レジスタが 16 本ずつあるわけ
コンパイラは Intel C/C++ Compiler 12.0.3 を用い,最適化
でなく,xmm レジスタは ymm レジスタの下位 128bit を使
を行う"-O3"オプション,並列計算を行うため"-openmp"オ
ⓒ 2013 Information Processing Society of Japan
24
HPCS2013
2013/1/15
2013年ハイパフォーマンスコンピューティングと計算科学シンポジウム
High Performance Computing Symposium 2013
表 1
Table 1
CPU
実験環境
Testing environment
Intel Core i7 2600K 3.4GHz
Intel Sandy Bridge Microarchitecture
Number of core
4
L3 Cache Size
8MB
Memory
DDR3-1333 Dual Channel 16GB
Memory bandwidth
21.2GB/s (10.6×2)
OS
Fedora 16
Compiler
Intel C/C++ Compiler 12.0.3
Compile-options(AVX)
–O3–xAVX–openmp–fp-model precise
Compile-options(SSE2)
–O3–xSSE2–openmp–fp-model precise
プション,最適化により精度に影響を及ぼさないために
"fp-model precise"オプションを用い.SSE2 向けのコードに
図 4
Figure 4
ベクトル演算の性能
Performance of vector operations
実験の結果,"scale"以外の演算で AVX と SSE2 の比は約
2.1 から 2.4 倍だった.AVX の理論値は SSE2 の 2 倍のはず
は"-xSSE2"オプションを,AVX 向けのコードには"-xAVX"
であるが,今回”scale"以外の演算において AVX は SSE2 の
オプションを用いた.
2 倍以上の倍率がでている.
また,AVX 化による計算精度に変化がないか確認するた
め,全ての演算において C ソースコード上で SSE2 の結果
この原因は,AVX は VEX プリフィックスによって発行
された 3 オペランド命令が実行できることが要因の 1 つと
して考えられる.SSE2 は 2 オペランド命令しか使用できな
と減算を行い,結果が 0 になることを確認した.
いため,演算のソースオペランドを書き換えてしまうため,
後で使う変数などは temp 変数などを用意して途中結果を
4.
倍々精度ベクトル演算
保存しなければならない.しかし,AVX では 3 オペランド
命令が実行できることにより,レジスタ退避,復元のため
4.1 対象とするベクトル演算
AVX を用いた倍々精度演算の高速化の効果を確認する
の temp 変数などに対する"move"命令が削減される.2 オペ
ため,表 2 の倍々精度ベクトル演算を対象とした実験を行
ランドから 3 オペランドになることによってレジスタ退避,
った. ここで,  及び  は倍々精度のスカラー値,
復元命令が削減される例を図 5 に示す.
”nrm2”と”scale”以外の演算では,基本ループ内において
, および  は倍々精度のベクトルである.
また,各演算のロード,ストア数と演算数を表 2 に示す.
倍々精度加算 1 回と乗算 1 回によって成り立つ FMA とい
うマクロ関数を用いている[4].FMA マクロ関数の演算量
Complexity は演算カーネル内の演算数である.
は倍々精度加算 (11flops)と倍々精度乗算(24flops)のの合計
である 35flops である.FMA マクロ関数を SSE2 で実行し
表 2 倍々精度ベクトル演算
Table 2 Double-Double precision vector operations
た場合はレジスタ退避,復元処理の回数が 13 回なのに対し,
Name of operation
Operation
Load
Store
Complexity
AVX では 3 回と,基本ループ内の処理が 10 回削減されて
axpy
 =  + 
2
1
35
いた.演算量が 35flops であることを考えると,性能に大
axpyz
 =  + 
2
1
35
きな影響を与えることがわかる.
xpay
 =  + 
2
1
35
"scale"の倍率は約 1.7 倍で,他のベクトル演算の倍率に
比べて低くなっている.これは SSE2 の"scale"の性能が良
scale
 = 
1
1
24
dot
 =  ⋅ 
2
0
35
nrm2
 = ‖‖2
1
0
31
いためである.また,"scale"で用いている倍々精度乗算を
1 回行う MUL マクロ関数は,レジスタ退避,復元処理が
4.2 ベクトル演算の性能
4 スレッド化した AVX 向けのコードと SSE2 向けのコー
ドでベクトルサイズ N を 105 で計測した結果を図 4 に示す.
なお,計測は最低で 70 回繰り返し,得られた値を回数で除
算した値を用いた.
図 5
Figure 5
ⓒ 2013 Information Processing Society of Japan
3 オペランドと 2 オペランドの違い
Difference between three and two operands
25
HPCS2013
2013/1/15
2013年ハイパフォーマンスコンピューティングと計算科学シンポジウム
High Performance Computing Symposium 2013
SSE2 では 5 回なのに対し,AVX では 3 回であり,FMA マ
クロ関数を使用している他のベクトル演算と比べて 3 オペ
ランド化の効果があまり見られないことも倍率が低い原因
であると考えられる.
なお,SSE2 のソースコードに対し,コンパイルオプショ
ン"-xAVX"を使いコンパイルすると,AVX の搭載されてい
るマシン上では VEX プリフィックスにより xmm レジスタ
を用いた 3 オペランド命令が発行されるため,コンパイル
オプション"-xSSE2"を使った場合より性能が向上する.今
回の SSE2 の実験項目に対して"-xAVX"コンパイルオプシ
ョンを付けた場合,各々20 から 40%程度性能が向上し,
SSE2 コードに対し"-xAVX"を用いて 3 オペランド化したも
図 6
Figure 6
"axpy"の性能 (1 スレッド)
Performance of "axpy" operation (1 Thread)
のは,“-xAVX”を用いた 3 オペランドの AVX に対し約 1.5
倍から 1.9 倍となった.
4.3 “axpy"の性能分析(1 スレッド)
"axpy"を 1 スレッドでベクトルサイズ N を 103 から
8. 0 ×105 まで 103 ずつ増加させたときの結果を図 6 に示す.
L3 キャッシュのサイズは 8MB で,長さ 2.6×105 の倍々精
度ベクトルを 2 本格納できる.
ベクトルサイズが 103 から 1.5×105 前後までは SSE2 と比
べ約 2.3 倍だが,サイズの増加に伴い徐々に性能が低下し
最終的に約 1.7 倍となる.
ベクトルがキャッシュに収まるとき,AVX の性能は約
16.8GFLOPS となった.1 コアあたりでの AVX を用いた場
図 7
Figure 7
メモリアクセスの相対性能 (1 スレッド)
Relative performance of memory access (1 Thread)
合の理論ピーク性能は 27.2GFLOPS なので,理論ピーク性
能の 61.8%である.
キャッシュ,メモリ性能を分析するため,"axpy"のロー
この原因の 1 つとして,倍々精度の加算と乗算アルゴリ
ド,ストア命令のみを実行し性能の測定を行った.FMA マ
ズムの演算のバランスが悪いことが考えられる.Core i7
クロ関数は load 命令 6 回,store 命令 2 回によって成り立っ
2600K は,乗算器と加算器が 1 つずつ搭載されている.し
ており,これら 8 回の命令のみを AVX 命令で行い,ベク
かし,FMA マクロ関数は,加減算命令 26 回,乗算命令 9
トルサイズを 103 から 4. 0 ×105 まで 103 ずつ増加させ測定
回によって成り立ち,加減算と乗算の回数に偏りがあるた
した.なお,コンパイラによる最適化で消去されないよう,
め,乗算器と加算器が並列して動作することを前提とした
計測ループ外で標準出力への出力処理を行った."axpy"が
理論ピーク性能がでることはない.乗算と加減算のバラン
ベクトルサイズ 105 のときの性能を基準としたキャッシュ,
スを考慮すると,FMA マクロ関数の理論ピーク性能は約
メモリアクセスの相対性能(相対実行時間の逆数)を図 7 に
18.3GFLOPS となり,前述のベクトルがキャッシュに収ま
示す.
る場合の性能 16.8GFLOPS は理論ピーク性能の 91.8%であ
実験の結果,ベクトルサイズが小さい間はキャッシュ,
る.演算以外の繰り返し構文の分岐命令などを考慮すると,
メモリアクセスの性能が高いが,ベクトルサイズが約
ベクトルがキャッシュに収まるとき,"axpy"は演算器の限
1.3×105 に達すると性能が低下する結果となった.
界まで性能がでていることがわかる.
この結果から,"axpy"の性能が低下する原因の 1 つとし
また,ベクトルがキャッシュに収まらない場合,SSE2
てキャッシュ,またはメモリアクセスがボトルネックにな
はサイズの増加に対してほとんど変化がないのに対し,
っていることがわかった,データがキャッシュに収まらな
AVX はベクトルの増加に従って徐々に性能が落ち,SSE2
いとき,データサイズと実測時間から 1 秒あたりに CPU が
の 1.7 倍程度まで低下する.この原因の 1 つとして,AVX
処理しているデータ量を求めると,AVX は約 12GB/s とな
は SSE2 と比べ同時演算数の増加によってメモリへの要求
る.今回のマシンのメモリ帯域は理論値最大 21.2GB/s であ
が高いため,キャッシュストールが発生していることが考
るため,メモリの帯域限界でなく,性能の低下はキャッシ
えられる.
ュストールの発生が原因であると考えられる.
ⓒ 2013 Information Processing Society of Japan
26
HPCS2013
2013/1/15
2013年ハイパフォーマンスコンピューティングと計算科学シンポジウム
High Performance Computing Symposium 2013
4.4 “axpy"の性能分析(4 スレッド)
次に,スレッド化の効果を確認するため, 4 スレッドで
ベクトルサイズ N を 103 から 8.0×105 まで 103 ずつ増加させ
たときの"axpy"の結果を図 8 に示す.
実験の結果から,ベクトルがキャッシュに収まるとき,
AVX は SSE2 の比は約 2.3 倍となり,1 スレッドのときと
ほぼ同様の性能向上が見られた.また,4 スレッドにした
ことで 1 スレッドと比べ約 3.7 倍の向上となり,スレッド
化による効果があることがわかった.しかし,ベクトルが
キャッシュサイズに収まらない場合は,サイズの増加に従
い AVX,SSE2 どちらにおいても性能が低下し,最終的に,
AVX は SSE2 とほぼ同様の約 12GFLOPS となる.
図 9
Figure 9
"axpy"の AVX/SSE2 の比
AVX speedup ratio (Performance of AVX/SSE2)
性能が低下する原因として,キャッシュにデータが収ま
らなくなったことによるキャッシュストールの発生と,4
キャッシュ,メモリアクセスがボトルネックとなっている
スレッド化によるメモリへの要求の増加によるメモリネッ
と考えられる.また,スレッド数が多いほど早くメモリボ
クが考えられる.1 秒あたりに CPU がメモリに要求し実行
トルネックに達する傾向が見られた.
しているデータ量を計算すると,AVX においてデータがキ
ャッシュに収まるときは約 55GB/s であった.今回のマシ
4.6
“dot"の性能分析
ンのメモリ帯域は理論値最大 21.2GB/s であり,キャッシュ
1 スレッドで"dot"をベクトルサイズ N を 103 から 8.0×105
に収まる範囲においてはキャッシュアクセスが有効なため
まで 103 ずつ増加させたときの結果を図 10 に示す."dot"
メモリ性能は制約となりにくい.しかし,データがキャッ
は"axpy"とほぼ同様の傾向が見られたが,"axpy"より性能
シュに収まらなくなるとキャッシュミスが発生しやすくな
が高く,ベクトルがキャッシュに収まるとき,AVX では約
ることにより,計算性能はメモリ性能の制約うけ,1 秒あ
18.0GFLOPS, SSE2 では約 7.5GFLOPS となっていた
たりに CPU が処理しているデータ量は約 12GB/s となった.
"dot"の内部は,"axpy"と同様の FMA マクロ関数を用い
ている.”axpy"は yi = a * xi + yi を計算するのに対し,"dot"
は val = xi * yi + val を計算する.違いはマクロ関数に与える
引数だけである.
コンパイラにより生成されたアセンブリコードを解析
すると,AVX で性能が異なる原因はコンパイラにより,ス
カラー値である val のストアが基本ループの外に移動され,
Store が基本ループ内で行われないようになっていた.
FMA マクロ関数の理論ピーク性能である約 18.3GFLOPS
と比較すると,"dot"は AVX において約 98.5%の性能がでて
おり,演算器の性能が限界まで引き出されている.スレッ
図 8
Figure 8
"axpy"の性能 (4 スレッド)
ド化の影響は"axpy"とほぼ同様であった.
Performance of "axpy" operation (4 Threads)
4.5 “axpy"の性能分析(マルチスレッド化)
スレッド数を増加させたときのキャッシュストール,メ
モリネックによる性能の変化を確認するため,"axpy"にお
いてベクトルサイズ N を 103 から 8.0×105 まで 103 ずつ増加
に示す.
すべてのスレッド数において,ベクトルがキャッシュに
収まる場合 SSE2 の 2 倍以上の性能がでているが,キャッ
シュに収まらない場合,1 スレッド以外は最終的に SSE2
と同性能になった.原因は,同時演算数の増加やマルチス
レッド化によりキャッシュ,メモリへのデータ要求が増え,
図 10
Figure 10
ⓒ 2013 Information Processing Society of Japan
"dot"の性能 (1 スレッド)
Performance of "dot" operation (1 Thread)
27
HPCS2013
2013/1/15
2013年ハイパフォーマンスコンピューティングと計算科学シンポジウム
High Performance Computing Symposium 2013
5.
疎行列ベクトル積
本研究では,倍精度の疎行列と倍々精度のベクトルの積
に対する測定を行った.実際の反復法ライブラリにおいて,
入力として倍精度の行列とベクトルが与えられることから,
計算の現実的な安定化の観点では,係数行列は倍精度のま
までよいと考えた.
分析には The University of Florida Sparse Matrix Collection
(フロリダコレクション)[10]から入手した 15 種の疎行列と
"random"を用いた."random"は非零要素の分布による影響
を調べるため,行列サイズ N が 105,1 行あたり 64 個の非
零要素がランダムに分布している疎行列である.結果は行
列サイズごとに 500 回反復計測した相加平均を用いた.
疎行列の格納形式は,Compressed Row Storage (CRS)形式
図 13
Figure 13
疎行列ベクトル積の性能 (1 スレッド)
Performance of sparse matrix vector product
(1 Thread)
[11]である.CRS 形式とは,非零要素数を nnz とすると,
①サイズ N×N の正方行列 A の非零要素の値を行方向に沿
って格納する長さ nnz の倍精度配列 value,②配列 value に
格納された非零要素の列番号を格納する長さ nnz の整数配
列 index,③配列 value と index の各行の開始位置を格納す
る長さ N+1 の整数配列 ptr からなる.図 11 に CRS 形式の
データ構造を,図 12 に CRS 形式の倍精度の疎行列 A と倍々
精度ベクトルxDD の積を示す.倍々精度数間の FMA マクロ
関数と比べ,乗算と加算が 1 回ずつ尐なく,演算数は 33flops
である.疎行列ベクトル積の演算数は 33*nnz とした.
図 14
Figure 13
疎行列ベクトル積の性能 (4 スレッド)
Performance of sparse matrix vector product
(4 Threads)
図 11
Figure 11
CRS 形式のデータ構造
Data structure of CRS format
要素数の数が多い疎行列は性能向上が大きい傾向がみられ
た."random"は,平均非零要素数がほぼ同一の"olafu"と比
べ性能が低く,非零要素の分布も性能を決定する要因であ
ると考えられる.疎行列の格納形式は CRS 形式なので連続
したメモリアクセスとなるため,この原因は倍々精度ベク
トルxDD に対するキャッシュミスであると考えられる.帯
行列やフロリダコレクションの疎行列のような非零要素の
図 12
Figure 12
CRS 形式の疎行列とベクトルの積
Product of sparse matrix of CRS format and vector
5.1 疎行列ベクトル積の性能
平均非零要素数の順に並べた 1 スレッドにおける結果を
図 13 に,4 スレッドにおける結果を図 14 に示す.括弧内
は平均非零要素数を示す.
実験の結果,SSE2 と比較して AVX は約 1.1 倍から 1.9
倍となった.また,"random"以外の疎行列では,平均非零
列番号が連続している場合が多いとき,xDD は L3 キャッシ
ュに残り続ける.しかし,"random"のような規則性のない
ものではxDD に対するアクセスが不連続になり,xDD に対す
るキャッシュミスが発生していると考えられる.
4 スレッドの計測結果において,SSE2 と比較して AVX
は約 1.1 倍から 2.0 倍となった.ここでも平均非零要素数
が多い疎行列ほど向上が大きい傾向が見られた.
4 スレッドと 1 スレッドを比べると,AVX は約 2.7 から
3.7 倍,SSE2 は約 2.5 から 3.7 倍となった.マルチスレッド
化による高速化は有効である.
ⓒ 2013 Information Processing Society of Japan
28
HPCS2013
2013/1/15
2013年ハイパフォーマンスコンピューティングと計算科学シンポジウム
High Performance Computing Symposium 2013
5.2 疎行列のサイズについての分析
疎行列のサイズ,平均非零要素数,非零要素の分布が性
能に与える影響を分析するため,帯幅を m とした

if( 0 ≤ j - i ≤ m ) ai j=value

else ai j = 0
を満たすテスト用帯行列を作成し,CRS 形式で格納した.
疎行列サイズが性能に与える影響を分析するため,帯幅
m = 32 のテスト用帯行列に対し,1 スレッドで行列サイズ
N を 103 から 4.0×105 まで 103 ずつ増加させた結果を図 15
に示す. L3 キャッシュには CRS 形式の倍精度の行列と
倍々精度のベクトルが,約 1.9×104 まで格納できる.テス
図 16
ト用帯行列では,キャッシュにデータが収まるとき,AVX
Figure 16
テスト用帯行列の性能 (4 スレッド)
Performance of band matrix (4 Threads)
は SSE2 の約 1.9 倍,キャッシュに収まらないとき約 1.8 倍
となった.
5.3 平均非零要素数についての分析
次に,4 スレッドで,行列サイズ N を 103 から 4.0×105
3
まで 10 ずつ増加させた計測結果を図 16 に示す.キャッシ
ュに収まるとき,AVX は SSE2 の約 1.8 倍,キャッシュに
平均非零要素数の影響を分析するため,テスト用帯行列
の行列サイズを固定し,帯幅 m を変化させた.
1 スレッドにおいて行列サイズを 105 に固定し,帯幅を 1
収まらないとき約 1.7 倍となった.また,1 スレッドと 4
から 100 まで変化させた性能を図 17 に示す.参考のため,
スレッドを比較すると,AVX はデータがキャッシュに収ま
平均非零要素数を帯幅とみなしてフロリダコレクションと
るとき約 3.7 倍,SSE2 は約 3.8 倍となり,キャッシュに収
"random"も同一グラフ上に示した.
まらないとき AVX は約 3.5 倍,SSE2 は約 3.6 倍となった.
実験の結果,帯幅の増加に従って性能が階段状に増加す
疎行列ベクトル積では,ベクトル演算と違い AVX の性
るが,帯幅 m = 40 前後から性能の増加の傾きが緩やかにな
能が SSE2 の性能まで落ちることはない.これは,係数行
り,帯幅 m = 65 以上ではほぼ一定の性能になった.帯幅 m
列を倍精度としたためメモリへの要求が低下したと考えら
= 4 のとき,1 スレッドにおいて,AVX は 5.9GFLOPS と理
れる.CPU が処理するデータ量は 4 スレッドでキャッシュ
論ピーク性能の約 21%,SSE2 は 4.8GFLOPS と理論ピーク
に 収 ま る と き 18GB/s, キ ャ ッ シ ュ に 収 ま ら な い と き
性能の約 18%となり,AVX は SSE2 の約 1.1 倍であるのに
15GB/s であり,いずれもメモリ帯域の理論値最大である
対 し , 帯 幅 40 で は 1 ス レ ッ ド に お い て AVX が 約
21.2GB/s に達しておらず,メモリはボトルネックになって
11.2GFLOPS と 理 論 ピ ー ク 性 能 の 約 41%, SSE2 が 約
いないと考えられる.ただし,ベクトル x はキャッシュに
6.5GFLOPS と理論ピーク性能の約 48%となり,AVX は
ヒットすると仮定し,1 度だけメモリから読み込むことと
SSE2 の約 1.7 倍の向上となった.倍精度行列と倍々精度ベ
仮定している.
クトルの積は,乗算と加減算のバランスが悪く,乗算と加
これらの結果から,倍精度疎行列と倍々精度ベクトルの
減算のバランスを考慮した 1 コアの理論ピーク性能は
積では行列サイズ N による性能への影響は尐ないことがわ
18GFLOPS となるため,AVX では理論ピーク性能の 62.3%,
かる.
SSE2 は 73.2%の性能がでている.
図 15
Figure 15
テスト用帯行列の性能 (1 スレッド)
Performance of band matrix (1 Thread)
図 17
Figure 17
帯幅と性能の関係 (帯行列, 1 スレッド)
Performance of sparse matrix vector product
with bandwidth (Band matrix, 1 Thread)
ⓒ 2013 Information Processing Society of Japan
29
HPCS2013
2013/1/15
2013年ハイパフォーマンスコンピューティングと計算科学シンポジウム
High Performance Computing Symposium 2013
スが大きなボトルネックになり性能が低下する.マルチス
レッド化を行った場合,メモリへの要求が高まり,メモリ
性能による制約を受けることで,AVX は SSE2 と同様の性
能となった.
CRS 形式で格納した倍精度の疎行列と倍々精度のベク
トルの積では,AVX は SSE2 と比較して約 1.1 倍から 1.9
倍となった.テスト用に作成した帯行列において,4 スレ
ッドでデータがキャッシュに収まるとき約 1.8 倍,キャッ
図 18
Figure 18
疎行列"Dubcova2"
Sparse matrix "Dubcova2"
シュに収まらない場合 1.7 倍となり,疎行列ベクトル積に
おいてはキャッシュにデータが収まらない場合においても
メモリアクセスはネックとならなかった.
性能が階段状に変化しているのは,SSE2 では 1 回の演算
CRS 形式で格納した帯行列のサイズを固定し,平均非零
で倍精度 2 つ,AVX では倍精度 4 つのデータを処理するた
要素数を変化させたところ,帯幅の増加に伴い性能が増加
め,端数処理を行う必要があるためと考えられる.AVX は,
するが,帯幅が約 40 を超えると性能の増加の傾きが緩やか
同時演算数が増えたことにより端数処理が発生しやすく,
になり,帯幅 65 以上になるとほぼ一定の性能だった.計測
テスト用帯行列とフロリダコレクションの性能差が大きい.
結果を総合すると,疎行列ベクトル積において,性能は平
フロリダコレクションの疎行列の平均非零要素数と同
均非零要素数に依存していると考えられる,
じ帯幅のテスト用疎行列は, "random"を除きほぼ同様の性
疎行列と平均非零要素数と同じ帯幅の帯行列と比較し
能だった.”random”は,AVX では約 7GFLOPS,帯幅 64 の
た性能は約 81.1%以上となり,乱数によって非零要素を配
テスト用疎行列と比べて 57.8%,SSE2 では約 4.25GFLOPS
置した疎行列においても約 57.8%であった.
(64.1%)であった.
これらのことから,SandyBridge マイクロアーキテクチャ
フロリダコレクションの疎行列の中で性能が低かった
で AVX を用いた倍々精度演算は,キャッシュに収まるベ
平均非零要素数 15.8 の疎行列"Dubcova2"の構造を図 18 に
クトル演算では,理論ピーク性能の約 60%,キャッシュに
示す.この疎行列は行列サイズ N が 65025,非零要素数が
収まらない場合はメモリ帯域の制約を受け,メモリ性能を
1030225 の対称行列で,テスト用帯行列と比べ複雑な構造
上限とした性能になる.また,倍精度疎行列と倍々精度ベ
をしている.このような疎行列においても,テスト用帯行
クトルの積において,平均非零要素数が 40 以上の一般的な
列と比較して AVX では約 81.1%,SSE2 では 91.7%となっ
疎行列では,理論ピーク性能の約 40%,性能がでにくいラ
た.このことから,"random"のような特殊ケースを除き,
ンダムな疎行列においても約 20%以上の性能になると予測
性能が平均非零要素数に大きく依存していることがわかる.
できる.
実際のアプリケーションではある程度データが連続した疎
今後の課題として,他のマシンやより多くの疎行列に対
行列を扱うことが多いと考えられるため,平均非零要素数
して検証,分析を行うことでより詳細な予測ができると考
から性能の予測が可能であると考えられる.
えている.また,倍々精度演算の加減算と乗算のバランス
や,プリフェッチなどによるメモリアクセスの改善を行う
6.
まとめ
本研究では,倍々精度演算による疎行列とベクトルの基
本演算を AVX 命令を用いて高速化する際の効果と問題点
ことで性能向上が期待できると考えられる.
謝辞
HPCS2013 の査読者から有益なコメントを頂きま
した.ここに感謝の意を表します.
を論じた.
キャッシュに収まる範囲のベクトル間の演算では 1 スレ
ッドにおいて AVX は SSE2 と比べて性能比が約 1.7 倍から
参考文献
[1] Hasegawa, H.: Utilizing the Quadruple-Precision
2.3 倍となった.AVX は 3 オペランド命令を実行できるた
floating-Point Arithmetic Operation for the Krylov Subspace
め,SSE2 と比べて 2 倍以上の性能がでたが,倍々精度演算
Methods, The 8th SIAM Conference on Applied Linear Algebra
の加減算と乗算のバランスが悪いことにより,"axpy"演算
(2003).
において,AVX は 16.8GFLOPS,SSE2 は 7.3GFLOPS とな
[2] Bailey, D ,H.: High-Precision Floating-Point Arithmetic in
り,理論ピーク性能の 21.7GFLOPS と比べ約 61.8%しか性
Scientific Computation, computing in Science and Engineering,
能がでることはなかった.さらなる性能向上のためには,
pp. 54–61 (2005).
加減算と乗算のバランスを考慮して改善する必要がある.
[3] Bailey, D, H.:A fortran-90 double-double library.
データがキャッシュに収まらない場合は,メモリアクセ
ⓒ 2013 Information Processing Society of Japan
http://crd-legacy.lbl.gov/~dhbailey/mpdist/
30
2013年ハイパフォーマンスコンピューティングと計算科学シンポジウム
High Performance Computing Symposium 2013
HPCS2013
2013/1/15
[4] 反復解法ライブラリ Lis, http://www.ssisc.org/lis/
[5] 小武守恒,藤井昭宏,長谷川秀彦,西田晃:反復法ラ
イブラリ向け 4 倍精度演算の実装と SSE2 を用いた高速化,
情報処理学会論文誌コンピューティングシステム,Vol. 1,
No. 1, pp. 73-84 (June 2008)
[6] 菱沼利彰,浅川圭介,藤井昭宏,田中輝雄,長谷川
秀彦:反復法ライブラリ向け倍々精度演算の AVX を用いた
高速化, 情報処理学会研究報告, Vol.2012-HPC-135, No.16,
pp.1-6 (2012, 8)
[7] Knuth, D, E. : The Art of Computer Programming:
Seminumerical Algorithms,Vol. 2, Addison-Wesley (1969).
[8] Dekker, T.: A floating-point technique for extending the
available precision, Numerische Mathematik, Vol. 18,
pp. 224–242 (1971).
[9] インテル Advanced Vector Extensions プログラミン
グ・リファレンス,
http://download.intel.com/jp/software/AVE/319433-006JA.pdf
[10] The University of Florida Sparse Matrix Collection,
http://www.cise.ufl.edu/research/sparse/matrices/
[11] Barrett, R., et al.: Templates for the Solution of Linear
Systems: Building Blocks for Iterative Methods, SIAM pp. 57–
65 (1994)
ⓒ 2013 Information Processing Society of Japan
31