rigayaの日記兼メモ帳

---

更新内容

---

結論から言うと、BroadwellとSkylakeで、useExp=オンの時の処理を高速化したというのがPMD_MT 高速化版+5です。

---

測定環境

---

速度の測定をした環境は以下の通り。

OS	Win10 x64	Win10 x64
CPU	i7 4770K (4C/8T)	i7 6700K (4C/8T)
CPU世代	Haswell	Skylake
Core	4.0GHz	4.3GHz
UnCore	4.0GHz	4.1GHz
キャッシュ	L3=8MB	L3=8MB
メモリ	DDR3-1600, 2ch	DDR4-2933, 2ch
CL	8-8-8-24-2	16-18-18-34-2

---

useExpの時のpmdの計算

---

pmd_mtでノイズ除去するとき、データ構造と重み関数を以下のように定義して、

typedef struct {
    short    y;  //    画素(輝度    )データ (     0 ～ 4096 )
    short    cb; //    画素(色差(青))データ ( -2048 ～ 2048 )
    short    cr; //    画素(色差(赤))データ ( -2048 ～ 2048 )
} PIXEL_YC;

PIXEL_YC *dst;  //出力フレーム
PIXEL_YC *src;  //入力フレーム
PIXEL_YC *gref; //srcにガウスぼかしを適用したフレーム

int func_weight(int x) {
  const double range = 4.0;
  const double strength2 = strength/100.0;
  const double inv_threshold2 = 1.0 / pow(2.0,threshold*0.1);
  return int((1<<16)/range * strength2
    * exp(-x*x * inv_threshold2));
}

あるピクセルは以下のように計算する。(cb, crも同様)

dst->y = src->y + (short)(( 
  ((src-max_w)->y - src->y)*func_weight((gref-max_w)->y - gref->y) + //上
  ((src-1    )->y - src->y)*func_weight((gref-1)->y     - gref->y) + //左
  ((src+1    )->y - src->y)*func_weight((gref+1)->y     - gref->y) + //右
  ((src+max_w)->y - src->y)*func_weight((gref+max_w)->y - gref->y)   //下
  ) >>16 );

というわけで、1ピクセル計算するのにfunc_weightは4回x3 (y, cb,cr)=12回呼ばれることになり、結構重い。さらにuseExpの場合、その名の通りfunc_weightの中にexpが使用されている。expを12回もやっていると、まあ、だいぶ重い計算になる。

もともとpmd_mtはこの重い重み計算を回避するため、最初にとりうるすべてのxについてあらかじめ重み計算を行ってテーブルを作っておき、あとはテーブル参照するだけというコードになっていた。

int *pmdp; //重みの計算結果を格納したテーブル

dst->y = src->y + (short)(( 
  ((src-max_w)->y - src->y)*pmdp[(gref-max_w)->y - gref->y] + //上
  ((src-1    )->y - src->y)*pmdp[(gref-1)->y     - gref->y] + //左
  ((src+1    )->y - src->y)*pmdp[(gref+1)->y     - gref->y] + //右
  ((src+max_w)->y - src->y)*pmdp[(gref+max_w)->y - gref->y]   //下
  ) >>16 );

こうすることで、重み計算を毎回行わなくて済み、だいぶ高速化できる。

---

useExpのSIMD化と表引き

---

PMD_MT 高速化版では、このPMDの計算をSIMD化することで高速化した。ただ、表引きのところは問題で、AVXまでは対応する命令がなく、結局その部分は非SIMDで行っていた。

ループ構造を取り外して、AVX2のサンプルコードはこんな感じ。

PIXEL_YC *src; //入力フレーム
PIXEL_YC *gau; //srcにガウスぼかしを適用したフレーム

__declspec(align(32)) int16_t diffBuf[64]; //一時バッファ
__declspec(align(32)) int expBuf[64];      //一時バッファ

__m256i ySrcUpperDiff, ySrcLowerDiff, ySrcLeftDiff, ySrcRightDiff;
__m256i yGauUpperDiff, yGauLowerDiff, yGauLeftDiff, yGauRightDiff;

//上下左右の差分をSIMD計算
getDiff(src, max_w, ySrcUpperDiff, ySrcLowerDiff, ySrcLeftDiff, ySrcRightDiff);
getDiff(gau, max_w, yGauUpperDiff, yGauLowerDiff, yGauLeftDiff, yGauRightDiff);

//差分の計算結果を一度バッファに格納
_mm256_store_si256((__m256i *)(diffBuf +  0), yGauUpperDiff);
_mm256_store_si256((__m256i *)(diffBuf + 16), yGauLowerDiff);
_mm256_store_si256((__m256i *)(diffBuf + 32), yGauLeftDiff);
_mm256_store_si256((__m256i *)(diffBuf + 48), yGauRightDiff);

//表引きを非SIMD実行
for (int i = 0; i < _countof(expBuf); i += 16) {
    expBuf[i+ 0] = pmdp[diffBuf[i+ 0]];
    expBuf[i+ 1] = pmdp[diffBuf[i+ 1]];
    expBuf[i+ 2] = pmdp[diffBuf[i+ 2]];
    expBuf[i+ 3] = pmdp[diffBuf[i+ 3]];
    expBuf[i+ 4] = pmdp[diffBuf[i+ 8]];
    expBuf[i+ 5] = pmdp[diffBuf[i+ 9]];
    expBuf[i+ 6] = pmdp[diffBuf[i+10]];
    expBuf[i+ 7] = pmdp[diffBuf[i+11]];
    expBuf[i+ 8] = pmdp[diffBuf[i+ 4]];
    expBuf[i+ 9] = pmdp[diffBuf[i+ 5]];
    expBuf[i+10] = pmdp[diffBuf[i+ 6]];
    expBuf[i+11] = pmdp[diffBuf[i+ 7]];
    expBuf[i+12] = pmdp[diffBuf[i+12]];
    expBuf[i+13] = pmdp[diffBuf[i+13]];
    expBuf[i+14] = pmdp[diffBuf[i+14]];
    expBuf[i+15] = pmdp[diffBuf[i+15]];
}

//表引きした重みをロード
__m256i yEUpperlo = _mm256_load_si256((__m256i *)(expBuf +  0));
__m256i yEUpperhi = _mm256_load_si256((__m256i *)(expBuf +  8));
__m256i yELowerlo = _mm256_load_si256((__m256i *)(expBuf + 16));
__m256i yELowerhi = _mm256_load_si256((__m256i *)(expBuf + 24));
__m256i yELeftlo  = _mm256_load_si256((__m256i *)(expBuf + 32));
__m256i yELefthi  = _mm256_load_si256((__m256i *)(expBuf + 40));
__m256i yERightlo = _mm256_load_si256((__m256i *)(expBuf + 48));
__m256i yERighthi = _mm256_load_si256((__m256i *)(expBuf + 56));

//取り出した重みを使用してSIMD計算を続行

まあ、あまり、速くない。

---

gather命令 (vpgatherdd)

---

さて、AVX2にはこの表引き(というか間接参照)をやる専用命令があって、

int *pDst, *pTable, *pIndex;
for (int i = 0; i < 8; i++) {
    pDst[i] = pTable[pIndex[i]];
}

これを、

int *pTable;
__m256i yDst, yIndex;
yDst = _mm256_i32gather_epi32(pTable, yIndex, 4);

とやるとvpgatherdd命令1発でやってくれるらしい(注1)。

これがあると、先ほどの表引きコードがかなり簡単になって、

PIXEL_YC *src; //入力フレーム
PIXEL_YC *gau; //srcにガウスぼかしを適用したフレーム

__m256i ySrcUpperDiff, ySrcLowerDiff, ySrcLeftDiff, ySrcRightDiff;
__m256i yGauUpperDiff, yGauLowerDiff, yGauLeftDiff, yGauRightDiff;

//上下左右の差分をSIMD計算
getDiff(src, max_w, ySrcUpperDiff, ySrcLowerDiff, ySrcLeftDiff, ySrcRightDiff);
getDiff(gau, max_w, yGauUpperDiff, yGauLowerDiff, yGauLeftDiff, yGauRightDiff);

//差分を16bit->32bitに変換
__m256i yGauUpperDiffLo = cvtlo256_epi16_epi32(yGauUpperDiff);
__m256i yGauUpperDiffHi = cvthi256_epi16_epi32(yGauUpperDiff);
__m256i yGauLowerDiffLo = cvtlo256_epi16_epi32(yGauLowerDiff);
__m256i yGauLowerDiffHi = cvthi256_epi16_epi32(yGauLowerDiff);
__m256i yGauLeftDiffLo  = cvtlo256_epi16_epi32(yGauLeftDiff);
__m256i yGauLeftDiffHi  = cvthi256_epi16_epi32(yGauLeftDiff);
__m256i yGauRightDiffLo = cvtlo256_epi16_epi32(yGauRightDiff);
__m256i yGauRightDiffHi = cvthi256_epi16_epi32(yGauRightDiff);

//重みの表引き
__m256i yEUpperlo = _mm256_i32gather_epi32(pmdp, yGauUpperDiffLo, 4);
__m256i yEUpperhi = _mm256_i32gather_epi32(pmdp, yGauUpperDiffHi, 4);
__m256i yELowerlo = _mm256_i32gather_epi32(pmdp, yGauLowerDiffLo, 4);
__m256i yELowerhi = _mm256_i32gather_epi32(pmdp, yGauLowerDiffHi, 4);
__m256i yELeftlo  = _mm256_i32gather_epi32(pmdp, yGauLeftDiffLo,  4);
__m256i yELefthi  = _mm256_i32gather_epi32(pmdp, yGauLeftDiffHi,  4);
__m256i yERightlo = _mm256_i32gather_epi32(pmdp, yGauRightDiffLo, 4);
__m256i yERighthi = _mm256_i32gather_epi32(pmdp, yGauRightDiffHi, 4);

//取り出した重みを使用してSIMD計算を続行

おお、なんかすごそうじゃないか！

ところが、この速度をHaswellで測ると…


処理対象: 1920x1080, 1024フレーム

あれー。gather使うとなんか遅くなる。

まあ、これはIntel自ら



Haswellのgatherは遅くて、次頑張るからBroadwellも買ってね的なことを言ってるのでしょーがない。

Intel Architectures Optimization Manualにも「11.16.4 Considerations for Gather Instructions」で"Example 11-43 shows some of those situations where, use of VGATHER on Intel microarchitecture code name Haswell is unlikely to provide performance benefit"(Haswellだとこういう条件ではgather使うと遅いよ)とあって、なんかいろいろ書いてあるので、まあ、基本Haswellのgather命令は遅いということなのだろう。

Haswellのgatherは要らない子?

というわけで、昔、gatherを使うコードを書いたのだけど、無効化してました。

---

それじゃSkylakeでは?

---

例のIntel Architectures Optimization Manualには、「11.16.4 Considerations for Gather Instructions」でBroadwellとSkylakeではgather命令がこんだけ速くなったアピールがされている。

実際計ってみると


処理対象: 1920x1080, 1024フレーム

と、ちゃんと速くなる。素晴らしい。

---

というわけで、どうもBroadwellも速いらしいので、BroadwellとSkylakeではgather命令を使うようにして、useExp=オンの時の処理を高速化したというのがPMD_MT 高速化版+5です。

---

ダウンロード>>
ダウンロード (ミラー) >>

ソースはこちら

---

(注1) 実際には、maskを初期化するためvpcmpeqb ymm0, ymm0, ymm0とかが入るので、2命令になったりする

---