pragma omp simdでvcompressが使われない


## 概要

`#pragma omp simd`を書くとIntelコンパイラの最適化でAVX-512のvcompress命令が使われず、結果速いコードが生成されないことがある。

## はじめに

私は当たり判定の処理が大好きです。少年時代のゲーム作りから現在まで、たびたび当たり判定の速い実装はないものかと試行錯誤することがあります。

ここでは、N粒子の当たり判定の結果、ヒットしたものだけをバッファに詰める、という処理を考えます。ここでバッファに詰める際には、隙間なく並べることを条件とします。

主要な当たり判定部分のコードは次のようなものです。
```
int nh = 0;
for(int k = 0; k < N;++k){
    hit_head[k]=nh;

    //#pragma omp simd   //問題の指示句
    #pragma ivdep
    #pragma vector
    for(int i = 0; i < N;++i){
        const double dx = x[k] - x[i];
        const double dy = y[k] - y[i];
        const double r2 = dx*dx + dy*dy;
        if(r2 < CUTOFF_SQ){
            hit_id[nh] = i;
            hit_dr[nh] = std::sqrt(r2);
            ++nh;
        }
    }            
}
hit_head[N]=nh;
```
粒子`k`と`i`の距離を測って、その距離(の2乗)が`CUTOFF_SQ`以下であれば、ヒットしたものとみなしてバッファ`hit_id`および`hit_dr`に詰めます。隙間なく詰めるために、ヒットした時だけ`nh`をインクリメントしています。

ここで、SIMDによるベクトル化を期待したいのは、内側の`i`ループです。しかし、ヒットした時だけ`nh`をインクリメントして詰める、という処理がなかなか難しく、私の経験ではAVX2までは、少なくともコンパイラの最適化ではベクトル化できませんでした。

一方で、AVX-512からはベクトル化できるようになりました。Xeon-PhiのAVX-512からSkylake以降のAVX-512まで共通してベクトル化できます。これは、`vcompress`というSIMD命令がAVX-512から実装されたからです。`vcompress`はベクトルレジスタうち、マスクビットの立っている要素だけを詰めるというものです。まさに今回の当たり判定にピッタリです。

## 問題

ところが、上記のコードで`//問題の指示句`と記した`#pragma omp simd`に問題がります。このコメント`//`を外して、指示句を有効化すると、コンパイラの最適化でなぜか`vcompress`が使われなくなります。

コンパイラのレポートで違いが見て取れます。コンパイルオプションは次のようにしました。
```
icpc -std=c++17 -O3 -ip -no-prec-div -xCORE-AVX512 -qopt-report=5 -qopt-zmm-usage=high -static-libstdc++
```

vcompressが使われる場合のレポート（`#pragma omp simd`をコメントアウトで無効化）
```
LOOP BEGIN at test_vcmp_vs_omp_simd.cpp(56,13)
 remark #15389: vectorization support: reference x[i] has unaligned access   [ test_vcmp_vs_omp_simd.cpp(57,42) ]
 remark #15389: vectorization support: reference y[i] has unaligned access   [ test_vcmp_vs_omp_simd.cpp(58,42) ]
 remark #15381: vectorization support: unaligned access used inside loop body
 remark #15305: vectorization support: vector length 8
 remark #15309: vectorization support: normalized vectorization overhead 0.365
 remark #15300: LOOP WAS VECTORIZED
 remark #15450: unmasked unaligned unit stride loads: 2
 remark #15457: masked unaligned unit stride stores: 2
 remark #15475: --- begin vector cost summary ---
 remark #15476: scalar cost: 56
 remark #15477: vector cost: 9.250
 remark #15478: estimated potential speedup: 5.620
 remark #15488: --- end vector cost summary ---
 remark #15497: vector compress: 2
LOOP END
```

vcompressが使われない場合のレポート（`#pragma omp simd`を有効化）
```
LOOP BEGIN at test_vcmp_vs_omp_simd.cpp(56,13)
 remark #15389: vectorization support: reference x[i] has unaligned access   [ test_vcmp_vs_omp_simd.cpp(57,42) ]
 remark #15389: vectorization support: reference y[i] has unaligned access   [ test_vcmp_vs_omp_simd.cpp(58,42) ]
 remark #15381: vectorization support: unaligned access used inside loop body
 remark #15416: vectorization support: irregularly indexed store was generated for the variable <hit_id[nh]>, masked, part of index is private   [ test_vcmp_vs_omp_simd.cpp(61,21) ]
 remark #15416: vectorization support: irregularly indexed store was generated for the variable <hit_dr[nh]>, masked, part of index is private   [ test_vcmp_vs_omp_simd.cpp(62,21) ]
 remark #15305: vectorization support: vector length 8
 remark #15309: vectorization support: normalized vectorization overhead 0.212
 remark #15301: SIMD LOOP WAS VECTORIZED
 remark #15450: unmasked unaligned unit stride loads: 2
 remark #15459: masked indexed (or scatter) stores: 2
 remark #15475: --- begin vector cost summary ---
 remark #15476: scalar cost: 56
 remark #15477: vector cost: 14.750
 remark #15478: estimated potential speedup: 3.540
 remark #15488: --- end vector cost summary ---
LOOP END
```

これらのレポートを比較すると、前者の`#pragma omp simd`を無効化したコードでは、ベクトル化による高速化の期待値が`estimated potential speedup: 5.620`となっています。

一方で、後者の`#pragma omp simd`を有効化したコードでは、`masked indexed (or scatter) stores`にもあるscatter store命令が使われていることが分かります。vcompress命令の代わりに、マスクとシャッフルで何とか格納先のインデックスを準備して、最後にscatter storeしてるのでしょう。高速化の期待値は`estimated potential speedup: 3.540`と低くなっています。

## ベンチマーク

実際に上記の当たり判定でベンチマークを行ってみました。

条件：`N=10000`, 繰り返し回数`1000`、`x[i],y[i]`の要素は`[-10.0, 10.0]`の範囲で一様乱数で生成。`CUTOFF_SQ=1.0`

CPU：Intel(R) Xeon(R) Silver 4208 CPU @ 2.10GHz

結果：

|pragma omp simd有効／無効|使用される命令|計算時間[s]|
|---|---|---|
|無効(コメントアウト)|vcompress|50.002970|
|有効|scatter store|109.771185|

store以外の計算負荷が小さいこともありますが、倍も差が出ていますね。やはりvcompressはとても有意義な命令ですね。

## まとめ

`#pragma omp simd`を使うと、現時点のIntelコンパイラの最適化では`vcompress`が使われない。`#pragma omp simd`は不必要に付けないようにしましょう。