SIMD型拡張命令をもっと使った最適化への道のり

科学技術振興調整費「並列化コンパイラ向け共通インフラストラクチャの研究」

藤波順久(ソニー株式会社) 阿部正佳(東京大学)

2002年1月11日

第43回プログラミング・シンポジウム

内容

はじめに
SIMD命令の特徴
提案する方法
- 到達レベル
- 予備評価
- 中間表現
- 最適化の手順
Coinsコンパイラ
おわりに

1. はじめに

SIMD型のマルチメディア向け拡張命令(SIMD命令)
- VIS,MVI,MMX,SSE,3DNow!,AltiVec,…
SIMD命令を使ったコードを生成するコンパイラ
- 自動的に使ってくれることは少ない
- 拡張されたデータ型と、命令に対応する関数(intrinsic)
- 保守性と移植性の低下
SIMD命令の特徴(低並列度低コスト、特殊な演算子、データサイズの混合)を生かした最適化はまだまだある
これを自動的に行いたい

自動的にできる最適化の例

short a[1000];
short b[1000];
short c[1000];

void f() {
  int i;
  for(i=0;i<1000;i++) c[i]=a[i]+b[i];
}

Function: _f
        XOR EAX, EAX
        ALIGN 00000010, 00000008
L0001:
        MOVQ MM2, QWORD PTR _b[EAX*2]
        MOVQ MM0, QWORD PTR _a[EAX*2]
        PADDW MM0, MM2
        MOVQ MM1, QWORD PTR _a+00000008[EAX*2]
        MOVQ QWORD PTR _c[EAX*2], MM0
        MOVQ MM0, QWORD PTR _b+00000008[EAX*2]
        PADDW MM1, MM0
        MOVQ MM2, QWORD PTR _a+00000010[EAX*2]
        MOVQ QWORD PTR _c+00000008[EAX*2], MM1
        MOVQ MM1, QWORD PTR _b+00000010[EAX*2]
        PADDW MM2, MM1
        MOVQ MM0, QWORD PTR _a+00000018[EAX*2]
        MOVQ QWORD PTR _c+00000010[EAX*2], MM2
        MOVQ MM2, QWORD PTR _b+00000018[EAX*2]
        PADDW MM0, MM2
        MOVQ MM1, QWORD PTR _a+00000020[EAX*2]
        MOVQ QWORD PTR _c+00000018[EAX*2], MM0
        MOVQ MM0, QWORD PTR _b+00000020[EAX*2]
        PADDW MM1, MM0
        MOVQ MM2, QWORD PTR _a+00000028[EAX*2]
        MOVQ QWORD PTR _c+00000020[EAX*2], MM1
        MOVQ MM1, QWORD PTR _b+00000028[EAX*2]
        PADDW MM2, MM1
        MOVQ MM0, QWORD PTR _a+00000030[EAX*2]
        MOVQ QWORD PTR _c+00000028[EAX*2], MM2
        MOVQ MM2, QWORD PTR _b+00000030[EAX*2]
        PADDW MM0, MM2
        MOVQ MM1, QWORD PTR _a+00000038[EAX*2]
        MOVQ QWORD PTR _c+00000030[EAX*2], MM0
        MOVQ MM0, QWORD PTR _b+00000038[EAX*2]
        PADDW MM1, MM0
        MOVQ QWORD PTR _c+00000038[EAX*2], MM1
        ADD EAX, 00000020
        CMP EAX, 000003C8
        JL L0001
        MOV CX, WORD PTR _a+00000002[EAX*2]
        MOV DX, WORD PTR _a[EAX*2]
        ADD CX, WORD PTR _b+00000002[EAX*2]
        ADD DX, WORD PTR _b[EAX*2]
        MOV WORD PTR _c[EAX*2], DX
        MOV WORD PTR _c+00000002[EAX*2], CX
        MOV CX, WORD PTR _a+00000006[EAX*2]
        MOV DX, WORD PTR _a+00000004[EAX*2]
        ADD CX, WORD PTR _b+00000006[EAX*2]
        ADD DX, WORD PTR _b+00000004[EAX*2]
        MOV WORD PTR _c+00000004[EAX*2], DX
        MOV WORD PTR _c+00000006[EAX*2], CX
        MOV CX, WORD PTR _a+0000000A[EAX*2]
        MOV DX, WORD PTR _a+00000008[EAX*2]
        ADD CX, WORD PTR _b+0000000A[EAX*2]
        ADD DX, WORD PTR _b+00000008[EAX*2]
        MOV WORD PTR _c+00000008[EAX*2], DX
        MOV WORD PTR _c+0000000A[EAX*2], CX
        MOV CX, WORD PTR _a+0000000E[EAX*2]
        MOV DX, WORD PTR _a+0000000C[EAX*2]
        ADD DX, WORD PTR _b+0000000C[EAX*2]
        ADD CX, WORD PTR _b+0000000E[EAX*2]
        EMMS 
        MOV WORD PTR _c+0000000C[EAX*2], DX
        MOV WORD PTR _c+0000000E[EAX*2], CX
L0000:
        RETN

intrinsicの例

ADDPS命令に対して、

__m128 _mm_add_ps(__m128 a, __m128 b)

これは、aとbの4つの単精度浮動小数点値を加算する。

2. SIMD命令の特徴(1)

典型的なものは、複数のサブオペランドに同じ演算を同時実行
表1参照(参考文献に、命令セットのマニュアルのURL)
- 並列度の小さいベクトル命令のように見える
- ベクトル命令はオーバーヘッド大
  - 高い並列度を生かせる場所(多数回のループ)に使う
- SIMD命令は、オーバーヘッド小
  - 積極的に使う意味がある
特徴は、あまり大きくない並列度、低い実行コスト

SIMD命令の特徴(2)

典型的でないもの
- 飽和加減算
- 最大最小
- 飽和縮小
- 混ぜ合わせ
- 乗算の上位
- 積和・内積
- 絶対差の和
表2参照
特徴は、特殊な演算子、データサイズの混合

3.1 SIMD最適化の到達レベル

レベル0: ループのベクトル化
従来技術の適用・既存のコンパイラがほぼやってくれる
レベル1: ループ以外でも、同じ演算の並列実行が可能なものをSIMD命令にする
研究段階
レベル2: 特殊な演算子を持つ命令の利用
大幅な高速化が可能
レベル3: データサイズの混合した演算をうまく扱う
単純に同じ操作の並列実行とみなせない場合

3.2 予備評価

static unsigned char sa,sr,sg,sb;
static short da,dr,dg,db;
static short k;

static void
hLineRight(unsigned char *p,int n,
           unsigned char a,short ea,
           unsigned char r,short er,
           unsigned char g,short eg,
           unsigned char b,short eb) {
  while(n!=0) {
    *p++=b; *p++=g; *p++=r; *p++=a;
    a+=sa; r+=sr; g+=sg; b+=sb;
    if((ea+=da)>=0) { a++; ea-=k; }
    if((er+=dr)>=0) { r++; er-=k; }
    if((eg+=dg)>=0) { g++; eg-=k; }
    if((eb+=db)>=0) { b++; eb-=k; }
    --n;
  }
}

CodeWarrior 6.0の出力

mwcc -opt all -inst comp,mmx,sse -machinecodelist -c simdlp16.c

Function: _hLineRight
        CMP DWORD PTR 00000008[ESP], 00000000
        PUSH EBX
        MOV ECX, DWORD PTR 00000028[ESP]
        PUSH ESI
        MOV EDX, DWORD PTR 00000024[ESP]
        PUSH EDI
        MOV EDI, DWORD PTR 00000024[ESP]
        PUSH EBP
        MOV EBP, DWORD PTR 00000038[ESP]
        MOV ESI, DWORD PTR 00000014[ESP]
        MOV EBX, DWORD PTR 00000024[ESP]
        JE L0001
        ALIGN 00000010, 00000008
L0002:
        MOV BYTE PTR 00000000[ESI], CL
        MOV BYTE PTR 00000001[ESI], DL
        MOV BYTE PTR 00000002[ESI], BL
        MOV AL, BYTE PTR 0000001C[ESP]
        MOV BYTE PTR 00000003[ESI], AL
        MOV AL, BYTE PTR 0000001C[ESP]
        ADD AL, BYTE PTR _sa
        ADD ESI, 00000004
        MOV BYTE PTR 0000001C[ESP], AL
        MOV AX, WORD PTR 00000020[ESP]
        ADD AX, WORD PTR _da
        ADD BL, BYTE PTR _sr
        MOV WORD PTR 00000020[ESP], AX
        ADD DL, BYTE PTR _sg
        ADD CL, BYTE PTR _sb
        CMP WORD PTR 00000020[ESP], 0000
        JL L0003
        MOV AX, WORD PTR 00000020[ESP]
        INC BYTE PTR 0000001C[ESP]
        SUB AX, WORD PTR _k
        MOV WORD PTR 00000020[ESP], AX
L0003:
        ADD DI, WORD PTR _dr
        TEST DI, DI
        JL L0004
        INC BL
        SUB DI, WORD PTR _k
L0004:
        MOV AX, WORD PTR 00000030[ESP]
        ADD AX, WORD PTR _dg
        MOV WORD PTR 00000030[ESP], AX
        CMP WORD PTR 00000030[ESP], 0000
        JL L0005
        MOV AX, WORD PTR 00000030[ESP]
        INC DL
        SUB AX, WORD PTR _k
        MOV WORD PTR 00000030[ESP], AX
L0005:
        ADD BP, WORD PTR _db
        TEST BP, BP
        JL L0006
        INC CL
        SUB BP, WORD PTR _k
L0006:
        DEC DWORD PTR 00000018[ESP]
        CMP DWORD PTR 00000018[ESP], 00000000
        JNE L0002
L0001:
L0000:
        POP EBP
        POP EDI
        POP ESI
        POP EBX
        RETN

240バイト
53-54クロック(ループ1回)
3657-3658クロック(ループ100回)

目指す最適化結果

_hLineRight:
        CMP     DWORD PTR _n[ESP],0
        JZ      L2
        MOVD    MM0,DWORD PTR [ESP+_b]
        MOVD    MM1,DWORD PTR _sb
        MOVQ    MM2,QWORD PTR [ESP+_eb]
        MOVQ    MM3,QWORD PTR _db
        MOVZX   EAX,_k
        IMUL    EAX,10001H
        MOVD    MM6,EAX
        MOVQ    MM7,MM6
        PSLLQ   MM6,32
        POR     MM7,MM6
        PCMPEQB MM4,MM4 ; MM4=0FFF...FH
        PSUBB   MM1,MM4
        MOV     EAX,_p[ESP]
        MOV     ECX,_n[ESP]
L1:
        MOVD    [EAX],MM0
        ADD     EAX,4
        PADDB   MM0,MM1
        PADDW   MM2,MM3
        PXOR    MM4,MM4 ; MM4=0
        PCMPGTW MM4,MM2
        MOVQ    MM5,MM4
        PACKSSWB MM5,MM5
        PADDB   MM0,MM5
        PANDN   MM4,MM7
        PSUBW   MM2,MM4
        DEC     ECX
        JNZ     L1
L2:
        RET

110バイト
19クロック(ループ1回)
696クロック(ループ100回)
コードサイズ半分以下、実行速度5倍以上(並列度4を超える)

3.3 中間表現(1)

演算子の細かい意味を表現したい
サブオペランドのサイズや位置を表現したい
他の最適化との影響も解析したい

命令の意味を正確に表現できる、RTLベースの中間表現を採用

 r1:=r2+r3
 r4:=mem(r5+12)

中間表現(2)

仮想レジスタを仮定
parallel: 副作用の同時発生
subreg: サブオペランド

PADDW MM0,MM1
 parallel(
  subreg_I16(MM0_I64,0):=subreg_I16(MM0_I64,0)+subreg_I16(MM1_I64,0)
  subreg_I16(MM0_I64,1):=subreg_I16(MM0_I64,1)+subreg_I16(MM1_I64,1)
  subreg_I16(MM0_I64,2):=subreg_I16(MM0_I64,2)+subreg_I16(MM1_I64,2)
  subreg_I16(MM0_I64,3):=subreg_I16(MM0_I64,3)+subreg_I16(MM1_I64,3)
 )

3.4 最適化の手順

論理演算を使ったif変換
基本ブロックをDAGに変換
DAGを分解
同じ演算を複数まとめる
命令とレジスタの制約の適用
SIMDレジスタ割当

論理演算を使ったif変換

一般には:

if c then
    :
  v:=t
    :
else      →  v:=(t bitand c) bitor (f bitand bitnot c)
    :
  v:=f
    :
endif

簡略化できる場合:

if t8>=0 then      t16:=t8>=0
  t0:=t0+1     →  t0:=t0-t16
  t8:=t8-t28       t8:=t8-(t28 bitand t16)
endif

基本ブロックをDAGに変換

DAGを分解

同じ演算を複数まとめる

同じ演算、同じサイズ、異なるデータ → parallelでまとめる
判定条件が必要
- 前後関係を決める(同じ基本ブロック内なら自明)
  - 決められなければ移動できないとする
- その前後関係による経路上で依存のないことを確認
  - 移動する代入文の使用する変数の定義がない
  - 移動する代入文の定義する変数の使用がない
  - 移動する代入文の定義する変数の定義は、その代入文しかない

命令とレジスタの制約の適用(1)

ソースオペランドの一方に結果が入る場合:

parallel(
 s4:=t0+t4
 s8:=t1+t5
 s12:=t2+t6
 s16:=t3+t7
)
    ↓
parallel(
 t0:=t0+t4
 t1:=t1+t5
 t2:=t2+t6
 t3:=t3+t7
)

命令とレジスタの制約の適用(2)

比較命令が、より小さいかどうかの判定しかない場合:

parallel(
 s6:=s5<z0
 s10:=s9<z1
 s14:=s13<z2
 s18:=s17<z3
)
    ↓
parallel(
 s6:=s5>z0
 s10:=s9>z1
 s14:=s13>z2
 s18:=s17>z3
)

SIMDレジスタ割り当て

parallel式を順に処理
- SIMD命令一つで処理できる個数より多くの代入があるなら、その数だけ取り出す
- 少ないならダミーを追加
- オペランドをSIMDレジスタに対するsubreg式に置き換える
- 必要なら、あふれ処理、変数並び順の制約の生成、サブオペランド位置の変換命令の挿入を行う
残ったものは、通常の命令を生成する処理に委ねる

4. Coinsコンパイラ

Coinsとは、

2000年度より5年プロジェクト
文部科学省科学技術振興調整費による「並列化コンパイラ向け共通インフラストラクチャの研究」
コンパイラの研究開発のための基盤
特に並列化に重点を置く

コンパイラの構成

特徴

二つの中間表現
- 高水準中間表現 HIR (High level Intermediate Representation)
  抽象構文木に近い表現であり、ソース言語の構成要素と一対一に近い対応関係がある。HIRの文法は、ソース言語に依存しない。
- 低水準中間表現 LIR (Low level Intermediate Representation)
  レジスタやメモリの値の変更などを表す式の並びであり、マシン命令の意味を表現可能である。LIRの文法は、マシン非依存である。
実装はJava言語
最初はSPARCプロセッサ用のCコンパイラだが、最終的には多言語、多機種に対応

LIRでの抽象化

コード生成より前の段階では、実際のマシン命令に対応しない式が許され、次のようなものが使える。

仮想レジスタ
- 物理レジスタ割当で解消
FRAME式(スタックフレームにとられる変数のアドレス)
- フレーム割当で解消
STATIC式(固定番地に確保される変数のアドレス)

フェーズ分け

ソースプログラムをHIR-Cに変換
C言語に依存したプログラム変換
HIR-CをHIR-baseに変換
高水準最適化
HIR-baseをLIRに変換
仮想レジスタ割付
低水準最適化
フレーム割当
パターンマッチング
物理レジスタ割当
マシン語レベルでの最適化
アセンブラコード生成

SIMD命令を使った最適化をどこに入れるか

全て低水準最適化の最後に入れる

SIMD命令を使った最適化で呼び出すコード解析部が、parallel式に対応
フレーム割当以降のフェーズが、SIMD物理レジスタとsubreg式に対応
マシン命令によくある偽の依存関係に対応(排他的論理和を使ってゼロにするときは使用とはみなさないなど)
生成したSIMD命令を壊さない
変数の並びの制約を番地固定の変数の並びに反映する

現状

Coinsコンパイラはまだ完成度が低い
SIMD命令を使った最適化単独での予備的な実装
目標とする最適化をほぼ達成(他から必要な情報を得られることが前提)
残された課題
- 命令とレジスタの制約の適用を簡潔に書く方法
- SIMD命令にできなかった処理の対応など
- SIMD命令を使ったために、別の自明でない最適化が必要になる場合

5. おわりに(1)

SIMD命令の特徴
- あまり大きくない並列度
- 低い実行コスト
- 特殊な演算子
- データサイズの混合
SIMD最適化の到達レベル
- レベル0: ループのベクトル化
- レベル1: ループ以外でも、同じ演算の並列実行が可能なものをSIMD命令にする
- レベル2: 特殊な演算子を持つ命令の利用
- レベル3: データサイズの混合した演算をうまく扱う

おわりに(2)

parallelとsubregを含むRTL上で6フェーズ
- 論理演算を使ったif変換
- 基本ブロックをDAGに変換
- DAGを分解
- 同じ演算を複数まとめる
- 命令とレジスタの制約の適用
- SIMDレジスタ割当
予備実験では並列度の4を超える5倍以上の高速化