Document

なんでもセミナー 3/4
ベクトル化のススメ
吉野寿宏 (コンピュータ科学・米澤研)
<[email protected]>
Introduction
ベクトル処理とは？

ベクトル = 要素を並べたもの
しかも各要素は互いに独立
 演算は要素ごとに行われる


データ並列性(DLP)が出てくる
ベクトルの演算は通常ループで書かれる
 独立なので要素間での演算順は不同、
または並列にも処理可能

たとえば

配列 a と b の和を c に計算する
for(int i = 0; i < 100; i++)
c[i] = a[i] + b[i];

もし要素が小さい型だったら…

たとえば 2 要素いっぺんに計算できたら
命令数はほぼ半分で済む
• そうできるための条件は？
SIMD

Single Instruction Multiple Data の略
一命令で複数データを演算する
 演算単位幅 > 要素データ幅な場合


プロセッサがそのための命令を持って
いることもある

例: MMX
• 浮動小数点レジスタ(64bit)を使いまわして、
• 8bit 整数×8 / 16bit 整数×4 /
32bit 整数×2 / 64bit 整数として扱える
ベクトル処理の長所・短所

長所


データ並列性を活かして高速化を図れる
(特にデータ幅が小さい場合には有効)
短所

周囲との整合性をとるために前準備が
必要になるなど、コードが複雑になる
• おかげで「読めねーよ！」と言われたことも
あるけど…
実例 1:
文字列スキャナ
文字列スキャナ

文字列とパターンが与えられ、文字列中
にパターンが何回出現するかを求める


DNA の処理とかでありそうな感じですね
例: 文字列=“aaaaaa”、パターン=“aa”
⇒5回

aaaaaa, aaaaaa, aaaaaa, aaaaaa,
aaaaaa
文字列スキャナ

ほぼ自明な実装
unsigned count(const char *str, const char *pat)
{
unsigned result = 0;
size_t str_len = strlen(str), pat_len = strlen(pat);
for(int idx = 0; idx <= str_len - pat_len; idx++) {
for(int j = 0, k = idx; j < pat_len; j++, k++)
if(str[k] != pat[j]) goto next_round;
}
result++;
next_round: ;
}
return result;
文字列スキャナ

ほぼ自明な実装
unsigned count(const char *str, const char *pat)
{
unsigned result = 0;
size_t str_len = strlen(str), pat_len = strlen(pat);
for(int idx = 0; idx <= str_len - pat_len; idx++) {
for(int j = 0, k = idx; j < pat_len; j++, k++)
if(str[k] != pat[j]) goto next_round;
}
result++;
next_round: ;
}
return result;
パターン比較
に注目！
パターン比較の処理 in depth

文字同士を比較して等しいか検査


s = s1s2s3…, p = p1p2p3… として、
s1=p1 かつ s2=p2 かつ s3=p3 かつ…
文字単位の比較操作は互いに独立

たとえば s1=p1 の結果は、
• s1, p1 のみに依存
• 他の文字の比較には影響しない
パターン比較の処理 in depth

1 文字 = 8bit
たいていの CPU はレジスタ長 32bit
 バイトアクセスは遅いこともある

• アラインメント調整のためシフタが使われる
• パイプラインストールを引き起こすことも
• 文字列比較くらいでは起こりませんが

いっそベクトル式に比較してみたら？
パターン比較系ベクトル式

ある程度条件を緩める必要がある

範囲外アクセスしても例外が起きない
• 勿論そんなに突拍子もないところではなくて、
せいぜい終端の隣のワードくらい

文字列のアラインメントがワード単位で
とれている
• つまりはアドレスがワードサイズの倍数に
なっていること
• Solaris/SPARC などでは、アラインメントが
とれていないアクセスで Bus Error になる
パターン比較系ベクトル式
末尾以外はワード単位で単純に比較
 末尾ではパターンの存在する部分のみ
比較

a a a a a a \0
比較
比較
a a a a a \0
一致
一
致
⇒ 一致
パターン比較系ベクトル式

C で書くとこんな感じ (リトルエンディアンの場合)

一致していれば 1、不一致なら 0 を返す
int match(char *str, char *pat) // str, pat must be aligned
{
unsigned long *s = (unsigned long *) str;
unsigned long *p = (unsigned long *) pat;
size_t pat_len = strlen(pat);
int i;
if(strlen(str) < pat_len) return 0;
}
for(i = 0; i < (pat_len >> 2); i++)
if(s[i] ^ p[i]) return 0;
if(pat_len & 3)
if((s[i] ^ p[i]) & ( (1UL << 8 * (pat_len & 3)) - 1 ))
return 0;
return 1;
文字列スキャナ・ベクトル版

長くなったので別紙を見てください

最外周に 4 回のループ
• 内周にしたほうがキャッシュ効率よい？


パターン・マスクをシフトしながら比較
大体 10%～20% くらい速くなった

gcc 3.3.3 -O0 で比較
• 実は -O3 ではナイーブ版に負けるんですが

10M 文字の文字列・数文字程度の
パターンを探す
練習問題

パターン比較のコードに「リトルエンディ
アン」と書いてありました。
じゃあビッグエンディアン用のコードは
どういう風に変わるの？

ヒント:
最後にあまりがある場合、マスクしてる
とこが変わりますね
問題の条件設定を利用しよう

もし「パターンが 8 文字以下」という制限
がついていたら

32bit 値の配列でたかだか 2 要素
• アラインメントあわせて比較するため、ずらす
分で 1 要素加えても高々 3 要素


比較のループを手動で展開できる
以上の変更を行うと、だいたいナイーブ
なコードの倍くらいの速度が出る

-O3 でも 10% 程度速い
論理演算と数値演算
を使い分ける
論理演算いろいろ

NOT (否定, ~)


2 度適用すると戻る
~ (~A) = A
NOT
0
1
1
0
XOR
0
1
0
0
1
1
1
0
XOR (排他的論理和, ^)


交換則 A ^ B = B ^ A
結合則
(A ^ B) ^ C = A ^ (B ^ C)
論理演算いろいろ

AND (論理積, &), OR (論理和, |)
AND
0
1

0
0
0
1
0
1
OR
0
1
0
0
1
1
1
1
互いに双対な関係
•
•
•
•
交換則
結合則
分配則
ド・モルガン則
A & B = B & A, A | B = B | A
(A & B) & C = A & (B & C)
A & (B | C) = (A & B) | (B & C)
~A | ~B = ~(A & B)
マスクを駆使しよう

たとえば
「整数 x より大きい最小の 4 の倍数」
を求めるには？

いろいろ書き方はあります
((x + 3) / 4) * 4 は分岐しないので賢い
 x % 4 == 0 ? x / 4 : x / 4 + 1 とか

• SETcc 使えば分岐なくていいけどね
マスクを駆使しよう

((x + 3) / 4) * 4 という式がどういう処理
をするのか考えよう
/4 = 2bit 右シフト、*4 = 左シフト
 左シフト時下位は 0 で埋められる
⇒ /4 * 4 は最下位 2bit を 0 にすること


マスクを使うこともできる

(x + 3) & (~3) でも OK
XOR は便利なのだ

特定ビットの反転に使う



異なっている部分をビット
単位で探す


1 と XOR をとると NOT
0 とならそのまま
XOR をとって 0 なら一致
とも言える
情報量が減らない


真理値表に 0 と 1 が同数
Cf. AND, OR
XOR
0
1
0
0
1
1
1
0
論理演算と数値演算の差

論理演算は繰り上がり・下がりがない
あるビットの演算は周囲によらない
⇒ ある 1 ビットのみに注目すればよい


繰り上がり・下がりは周囲に影響する

複数要素を同時処理する場合には面倒
• 要素間を跨がないように注意が必要

数値演算は、影響範囲を制限しないと
使いにくいことも
• マスクなどをかけて演算
1 であるビット数を数える

ワード中で 1 になっているビット数を数える
処理

たとえば 0x12345678 ⇒ 13
• 0001 0010 0011 0100 0101 0110 0111 1000

ループで回すのはよくやりますね
unsigned result = 0;
for(int i = 32; i--; x >>= 1)
if(x & 1) result++;
1 であるビット数を数える

しかしこんなんでも OK
unsigned BX_(unsigned x) {
x -= (x >> 1) & 0x77777777;
x -= (x >> 2) & 0x33333333
return x - (x >> 3) & 0x11111111;
}
unsigned COUNT(unsigned x) {
x = BX_(x) + BX_(x >> 4);
return (x & 0x0F0F0F0F) % 255;
}
1 であるビット数を数える
いわゆる分割統治法
 4bit だけで考えてみましょう

BX_(x) {
return x - (x >> 1) & 7 – (x >> 2) & 3
– (x >> 3) & 1;
}
x=
b3
b2
b1
b0
b3
b2
b1
b0
b3
b2
b1
b0
b3
b2
b1
b0
1 であるビット数を数える:
数学的な正当性
BX_ は、4bit ずつの区画ごとに、1 で
あるビット数を数える
 2 進数 x = b3b2b1b0 が表現する数は
23b3 + 22b2 + 21b1 + 20b0

(x >> 1) & 0x7 = 22b3 + 21b2 + 20b1
 (x >> 2) & 0x3 = 21b3 + 20b2
 (x >> 3) & 0x1 = 20b3
 したがって、BX_(x) = b3 + b2 + b1 + b0

1 であるビット数を数える:
数学的な正当性

BX_ によって 4bit ごとにビット数が出る


たとえば 0x10223204 のように
4bit シフトして足すと 8bit 分になる
はみ出した分はマスク処理
 ⇒ 0x01040504


最後に全部和をとる
28 mod 255 = 1 を利用
 ⇒ 1 + 4 + 5 + 4 = 14

高速化のテクニック
高速実行しやすいコードとは

CPU の気持ちになって考えてみるに…

メモリのアクセス局所性がある
• データキャッシュがきく

メモリアクセスのパターンが単純
• ストライドアクセス(線形なアクセス)とかで
データプリフェッチがうまく動くこと

条件分岐ミスが起こらない
• 分岐予測ミスのペナルティは結構大きい
• そもそも条件分岐がなかったり、静的予測で
ほとんど当たる場合は嬉しいかも
条件分岐

最近のプロセッサは投機実行を行う
予測ミスがあるとパイプラインがクリア
されるので、ペナルティ大
 そもそも分岐を減らすのが効果的


条件分岐を行うのはどういう場合か？
for, while, do-while などループの終了
条件
 三項演算子による分岐
 if 文による分岐

IA-32 プロセッサのおはなし

静的予測エンジンは分岐方向で判断




負方向 ⇒ ループ条件と考え成立と予測
正方向 ⇒ 脱出と考え不成立と予測
これだけを考えれば、ループでは最後の一回
のみが予測ミス
ということは


ループはこの法則に従えばよく当たる
if 文などでほぼ均等に両方向分岐するような
場合は当たりにくいかもしれない
IA-32 プロセッサのおはなし

SETcc 命令

分岐するかわりに、その条件が成立する
かどうかをレジスタに代入する命令
• 成立するなら 1、不成立なら 0 になる
• 三項演算子(pred ? expr : expr)に有効

分岐するわけではないのでパイプライン
はクリアされない
IA-32 プロセッサのおはなし

SETcc 命令を使って三項演算を行う

L1:
L2:
L3:
eax = ebx > 0 ? a1 : a2
cmp
jg
mov
jmp
mov
…
ebx, 0
L2
eax, a2
L3
eax, a1
条件分岐が
なくなった！
cmp
setg
dec
and
add
…
ebx, 0
eax
eax
eax, a2 – a1
eax, a1
実例 2:
αブレンディング
αブレンディング

透過処理
各ピクセル、各色チャネルについて:
 結果 = A ×α + B ×(1 - α)

• 色空間での内分点
• 比率αは「α値」と呼ばれる

画像データ: 1 ピクセル = 16～24bit
各チャネルは 5or6bit ～ 8bit 程度の幅
 ベクトル処理向きのサイズ

24bit/pixel ビットマップ
R(0,0) G(0,0) B(0,0) R(1,0) G(1,0) B(1,0)
… R(w-1,0) G(w-1,0) B(w-1,0)
R(0,1) G(0,1) B(0,1)
padding
それぞれのチャネルは 8bit 幅
 各行の末尾は 4-byte aligned

…
…
アルファブレンディング
ナイーブな例
BYTE *out = <output bitmap>;
BYTE *in1 = <src1>, *in2 = <src2>;
DWORD dwSize = <bitmap size>;
DWORD dwSizeOfPad = <padding size>;
for(int y = 0; y < height; y++) {
for(int x = 0; x < width; x++) {
DWORD red =
in1[0] * α + in2[0] * (256 - α);
// the same code for green and blue channel
out[0] = red; // and green, blue
in1 += 3, in2 += 3, out += 3;
}
in1 = dwSizeOfPad; // and in2, out
}
ループをまとめる

ナイーブなコードでは 2 重ループに
なっていた
縦方向のループと横方向のループ
 色チャンネルでもループすると 3 重に


しかしメモリが連続配置されているなら
一つにまとめてよい
ビットマップでは確かに連続した領域
 行間に padding は入っているが、バイト
単位処理には影響しない

ループをまとめてみた
アルファブレンディング
BYTE *out = <output bitmap>;
BYTE *in1 = <src1>, *in2 = <src2>;
DWORD dwSize = <bitmap size>;
while(dwSize--) {
DWORD temp =
*in1++ * α + *in2++ * (256 - α);
*out++ = temp >> 8;
}

ループ一つだけになりました
 簡潔ですね
ループをまとめてみた
アルファブレンディング

R, G, B 全て同じ計算なのでループは
一つで済む


padding はイメージに影響しないので、
計算してもしなくてもよい
⇒ ループにするコストを考えてフラットに
しかし、まだ要素(バイト)単位で計算

せっかくレジスタは 32bit あるのに…
複数要素の同時計算

α値の乗算はそのままでは隣にはみ出す

α値は 0～256 (固定小数点と考え、後で 256
で割る)、各チャネルは 0～255
⇒ 計算結果の最大値は 0xFF00 (16bit)
R
B
G
×
R
α
R×α
G×α
B×α
R×α
赤枠部分のみが欲しい
複数要素の同時計算

2回にわけて計算すればよい

マスクをとってあとで合成すればよし
B
R
α
×
R
G
×
α
R×α
G×α
B×α
R×α
影響範囲が
重ならない
ここはレジスタを
はみ出す
複数要素の同時計算

α値を固定小数点で考えているので、最後に
除算が必要


ここでは 8bit 右シフト
ビットシフトを前に持ってきてもよい
R
G
×
α
G×α
R×α
これでレジスタに
収まる
複数要素の同時処理

0, 2 バイト目
B
R×α
α値の乗算
 8ビットシフト
B×α
8ビットシフト
R
• 下 8bit は必ず 0


α値の乗算
あとはマスクかけて
合成
R×α
B×α
1, 3 バイト目

α
×


R
G
R
×
G
α
G×α
R×α
複数要素の同時処理

こんなんでました
DWORD *out = <output>;
DWORD *in1 = <src1>, *in2 = <src2>;
DWORD dwSize = <bitmap size>;
while(dwSize--) {
DWORD temp;
temp = (*in1 & 0x00FF00FF) * α
+ (*in2 & 0x00FF00FF) * (256 – α);
*out = (temp >> 8) & 0x00FF00FF;
}
temp = ((*in1++ & 0xFF00FF00) >> 8) * α
+ ((*in2++ & 0xFF00FF00) >> 8) * (256 – α);
*out++ |= temp & 0xFF00FF00;
動作速度の計測

Pentium IIIm 800MHz の WindowsXP マシン


gcc 3.3.3 -O3 (Cygwin)
64×64×24bpp のビットマップ


キャッシュに載るくらいのサイズってことで
それぞれ 100,000 回ずつ
ナイーブ
ループ
併合
同時
処理
時間 (秒)
8.067
11.076
6.435
speedup
1
0.638
1.254
動作速度の計測

1024×768×24bpp のビットマップだったら？


キャッシュからはみ出すのでミスが起こる
それぞれ 200 回ずつ
ナイーブ
時間 (秒)

5.891
ループ
併合
5.849
あまり有意な差は認められなくなった
同時
処理
5.769
実験結果の考察

ループをまとめただけでは嬉しくない


ナイーブなコードより遅かったし
同時処理で 25% 程度高速に
4 要素処理するのにナイーブなループ
(= 3 要素分)と同等量の命令数
 単純計算で ¾ だから高々～33%


キャッシュミスのペナルティは結構大きい

キャッシュからはみ出すような大きな
ビットマップではそのペナルティが優位
ライフゲームを
ベクトル化する
ライフゲーム

J. H. Conway 考案によるゲーム
2次元の、区画(セル)に分けられた世界
 各セルには生物が存在するかしないか
の状態


セルオートマトン


次世代での状態は、現在の状態と周囲
8 方向のセルの状態により決定する
ライフゲームは万能チューリング機械
ライフゲームのルール

誕生するパターン


周囲 3 セルが生存
生命が維持するパターン

周囲 2 セルまたは 3 セルが生存
＝生命なし
＝生命あり
ライフゲームのルール

死滅するパターン
＝生命なし
前述のもの以外
＝生命あり

過密状態
過疎状態
ライフゲームを作るには

自分と周囲のセル数を数える
たかだか隣接セルしか見ない
⇒局所的
 隣の要素に影響する
⇒境界部分で特別処理が必要


自分の状態と周囲のセル数から次世代
の状態を決める

各セルごとに独立な計算
ライフゲームを作るには

自分と周囲のセル数を数える
たかだか隣接セルしか見ない
⇒局所的
 隣の要素に影響する
⇒境界部分で特別処理が必要


自分の状態と周囲のセル数から次世代
の状態を決める

各セルごとに独立な計算
次世代の状態を求める


死 = 0, 生 = 1 とする
周囲のセル数はたかだか 8 セルしかないので、1 セル
あたり 4bit で領域は十分


32bit レジスタ 1 ワードあたり 8 セル分を保持できる
まず、わかりやすいように表を書いてみる
周囲
0
1
2
3
4
5
6
7
8
0→
0
0
0
1
0
0
0
0
0
1→
0
0
1
1
0
0
0
0
0
状態
次世代の状態を求める
周囲
0
1
2
3
4
5
6
7
8
0→
0
0
0
1
0
0
0
0
0
1→
0
0
1
1
0
0
0
0
0
状態

この表を論理圧縮してみよう



0:0011 + 1:0010 + 1:0011 = 0:0011 + 1:001*
周囲の生存セル数 OR 現在状態 == 3
となる
3 なら 1、それ以外で 0 となる論理を組む

上 2bit を反転して 4bit 全体を AND
次世代の状態を求める

C で書くとこんな感じ



state: 現在の状態, count: 周囲の生存セル数
next: 次世代状態
もちろんそれぞれ 4bit/セルである
next
next
next
next
= (count | state) ^ 0xCCCCCCCC;
&= next >> 1;
&= next >> 2;
&= 0x11111111;
ライフゲームを作るには
(Revisited)

自分と周囲のセル数を数える
たかだか隣接セルしか見ない
⇒局所的
 隣の要素に影響する
⇒境界部分で特別処理が必要


自分の状態と周囲のセル数から次世代
の状態を決める

各セルごとに独立な計算
周囲のセルを数える

隣接セルに影響を受ける

したがってひとつのマップ内のみで完結
させることは難しい
• 実際にはできます (3数の和なので、和を
とるたびにひとつ前に書き込めばいい)

後段の処理はひとつのマップでよい
• マップを 2 つ用意して交互に使用する


境界部の処理も必要
周囲 = (自分+周囲) - 自分
周囲のセルを数える

縦方向の和と横方向の和に分割して
考える
縦 3 行の和、の 3 列分の和をとる
 または横 3 列の和、の 3 行分でもよい

または
周囲のセルを数える
縦方向の和:
オフセットを加えてロード、足し算
 横方向の和:
シフトして足し算


たかだか 9 までにしかならないので、
セル境界ははみ出さない

↑こういう保証が重要
境界の処理

ビットシフトによって横方向の境界を
超えてしまう
前の行
現在の行
次の行
通常隣り合って(間隔をあけず)配置
端の処理を別ループにする？
 番兵(緩衝地帯)を設けて後で調整する？

境界の処理

とりあえず番兵にしてみました
行と行の間に 1 セル以上の隙間を開け
 ループは一つで計算しておいて
 最後に隙間部分をマスクしてクリア

完成品

結構長いのでソースコードはどこかに
おいておきます

てきとー GUI フロントエンドをつけて
アプリケーションにしたててみました

XGA の画面で最大化表示して 230fps
ちょいくらい
参考資料

今回作成したソースコードなど


IA-32 Architecture Software Developer’s
Manual



http://www.yl.is.s.u-tokyo.ac.jp/~tossy-2/nandemo/
あたりにおいておく予定です
http://developer.intel.com/design/pentium4/manuals/in
dex_new.htm
命令セットリファレンスなどのドキュメントのほか、
パフォーマンス最適化のためのマニュアルなどもあります
Wikipedia 「ライフゲーム」

http://ja.wikipedia.org/wiki/%E3%83%A9%E3%82%A
4%E3%83%95%E3%82%B2%E3%83%BC%E3%83
%A0

Download Report