SX-Aurora TSUBASAでvector intrinsicを使ってみる５：ベクトル化とベンチマーク


## はじめに

全開で簡単ではありますがループのベクトル化もできるようになりました。

とはいえ本当に高速化できたかどうかはわかりません。特にintrinsicで無理やりベクトル化しても速くならない、ということはよくあります。（これは`#pragma`で無理やりベクトル化した場合もそうですね。）

ということで今回は、実際に時間計測をしてみて、ベクトル化によるスピードアップを体感してみます。

## ターゲットとする問題

一次元のグリッドデータに対して、二階微分を二次精度の数値差分で求める問題を取り上げます。

ベクトルする前のコードのメインループ部分は次のようなものです。
```c++
//original loop
const double coef = 1.0/(2.0 * xi_delta*xi_delta);
c[0] += (b[1] - 2.0 * b[0] + 0.0) * coef;
for(int i = 1; i < N-1; ++i){
    c[i] += (b[i+1] - 2.0*b[i] + b[i-1]) * coef;
}
c[N-1] += (0.0 - 2.0 * b[N-1] + b[N-2]) * coef;
```
配列`b`にグリッドデータとしての関数の値が入っており、その二次精度の差分の結果を配列`c`にいれています。インデックス`i`に対して、その前後1グリッドの`b[i-1], b[i+1]`にアクセスするところをどのようにベクトル化するかがポイントになりそうです。

ベンチマークにおいてはループ長`N=10000`の上記のループを、さらに`NK=1000`回繰り返して、掛かった時間を計測しました。

ソースファイル全文は最後に載せます。

## ベクトル化しないときの計算時間

まずは、ベクトル化せずに、実測してみます。ベクトル化しない場合のコンパイルオプションは次のようにしました。
```
nc++ -std=c++17 -O2 -mno-vector
```
ここで、`-mno-vector`オプションでベクトル化を禁止しています。

この場合の計算時間は
```
time  =0.593170 [s]
```
でした。これが規準タイムです。

## コンパイラの自動ベクトル化を利用したとき計算時間

今度は、上記と同様のコードに対しコンパイラの自動ベクトル化を行うことでベクトル化します。コンパイルオプションは次のようにしました。
```
nc++ -std=c++17 -O2
```

無事にベクトル化できた場合は次のようなメッセージが出るでしょう。
```
nc++: vec( 101): source.cpp, line 49: Vectorized loop.
```

この場合の計算時間は
```
time  =0.011747 [s]
```
でした。ベクトル化しなかったときに比べて約50倍速くなってますね。

## Intrinsicによるベクトル化したときの計算時間（3回ロード版）

今度はIntrinsicでベクトル化してみます。intrinsicで書き換えたのは、上記のmainループ分だけです。以下のようにしてみました。
```c++
//intrinsic (1) bからのloadを3回する場合//

constexpr int VL = 256;

//original loop
const double coef = 1.0/(2.0 * xi_delta*xi_delta);

c[0] += (b[1] - 2.0 * b[0] + 0.0) * coef;
    
//(1) main part
int i;
for(i = 1; i <= N-1 - VL; i+=VL){            
    __vr vb;
    __builtin_ve_vld(vb, &b[i], 8);
    __vr vc;
    __builtin_ve_vld(vc, &c[i], 8);

    __vr vd;
    __builtin_ve_vld(vd, &b[i-1], 8); 
    
    __vr ve;
    __builtin_ve_vld(ve, &b[i+1], 8); 
    
    __builtin_ve_vfadd(ve, ve, vd);   
    __builtin_ve_vfnmsub(vd, vb, 2.0, ve);
    __builtin_ve_vfmadd(vc, vd, coef, vc);
    
    //store to array from vector resister//
    __builtin_ve_vst(vc, &c[i], 8);
}   

if(i < N-1){                
    const int vl = N-1 - i;
    __vm mask;
    __builtin_ve_eqvm(mask, mask, mask);//all bit stands// 

    __vr vb;
    __builtin_ve_vld(vb, &b[i], 8, mask, vl);
    __vr vc;
    __builtin_ve_vld(vc, &c[i], 8, mask, vl);

    __vr vd;
    __builtin_ve_vld(vd, &b[i-1], 8, mask, vl);

    __vr ve;
    __builtin_ve_vld(ve, &b[i+1], 8, mask, vl);

    __builtin_ve_vfadd(ve, ve, vd, mask, vl);   
    __builtin_ve_vfnmsub(vd, vb, 2.0, ve, mask, vl);
    __builtin_ve_vfmadd(vc, vd, coef, vc, mask, vl);

    //store to array from vector resister//
    __builtin_ve_vst(vc, &c[i], 8, mask, vl);
}

c[N-1] += (0.0 - 2.0 * b[N-1] + b[N-2]) * coef;
```

この場合は、`b[i]`を先頭に256要素をベクトルレジスタ`vb`にロードするのに加え、それぞれ前後に1つずらした`b[i-1]`を先頭とした要素を`vd`に、`b[i+1]`を先頭とした要素を`ve`にロードします。その後、add演算やFMA演算、FMAの一種であるfnmsub演算を使って差分を計算しています。

実は先の自動ベクトル化によって出力されるアセンブラも、ほぼこれに近いコードが出力されていました。

さて、このソースを自動ベクトル化と同じオプションでコンパイルすると、次のように該当ループの部分がベクトル化できなかったと旨のメッセージが出ます。ただ、アセンブラを見てみるとベクトル命令が使われているので、おそらく大丈夫と思います。
```
nc++: vec( 103): source.cpp, line 65: Unvectorized loop.
```

この場合の計算時間は
```
time  =0.012196 [s]
```
でした。ベクトル化しなかったときに比べて約48倍速くなってますね。また、自動ベクトル化と遜色ない速度です。

## Intrinsicによるベクトル化したときの計算時間（1回ロード＋roration）

自動ベクトル化も上記のintrinsicの実装も、配列`b`からのロードを開始点をずらして3度も行っていました。キャッシュやメモリの速度にも依りますが、なんだかもったいない気もします。

そこで、ベクトルレジスタ上で要素ずらす(roration)ことで、ロードの回数を1回にしてみました。ベクトルレジスタのrotationには`__builtin_ve_vmv`を使います。

（もしかすると、rotationをしない単純なレジスタのコピーでもこの命令を使うのかもしれません。仕様書がないのでわかりませんが、命令リストと睨めっこする限りはそんな感じです。）

rotationは基本的に要素を前進方向にだけずらしますが、レジスタの幅256をはみ出した部分はインデックスの0側から回り込んで書き込まれます。よって、シフトではなくrotationなのですね。

これを使ったコードは次のようになりました。
```c++
//intrinsic (2) bからのloadの回数を1回にしてrotateで代用する場合//

constexpr int VL = 256;

//original loop
const double coef = 1.0/(2.0 * xi_delta*xi_delta);

c[0] += (b[1] - 2.0 * b[0] + 0.0) * coef;
    
//(1) main part
int i;
for(i = 1; i <= N-1 - VL; i+=VL){            
    __vr vb;
    __builtin_ve_vld(vb, &b[i], 8);
    __vr vc;
    __builtin_ve_vld(vc, &c[i], 8);

    __vr vd;
    __builtin_ve_vmv(vd, vb, 1);  //vd[0:255] = (vd[255],vd[0:254]) //rotate 1//
    __builtin_ve_lsv(vd, 0, b[i-1]);  //vd[0] = b[i-1] //

    __vr ve;
    __builtin_ve_vmv(ve, vb, 255);  //vd[0:255] = (vd[1:255],vd[0]) //rotate 255(-1)//
    __builtin_ve_lsv(ve, 255, b[i+VL]);  //vd[255] = b[i+VL] //

    __builtin_ve_vfadd(ve, ve, vd);   
    __builtin_ve_vfnmsub(vd, vb, 2.0, ve);
    __builtin_ve_vfmadd(vc, vd, coef, vc);
    
    //store to array from vector resister//
    __builtin_ve_vst(vc, &c[i], 8);
}   

if(i < N-1){                
    const int vl = N-1 - i;
    __vm mask;
    __builtin_ve_eqvm(mask, mask, mask);//all bit stands// 

    __vr vb;
    __builtin_ve_vld(vb, &b[i], 8, mask, vl);
    __vr vc;
    __builtin_ve_vld(vc, &c[i], 8, mask, vl);

    __vr vd;
    __builtin_ve_vmv(vd, vb, 1);  //vd[0:255] = (vd[255],vd[0:254]) //rotate 1//
    __builtin_ve_lsv(vd, 0, b[i-1]);  //vd[0] = b[i-1] //

    __vr ve;
    __builtin_ve_vmv(ve, vb, 255);  //vd[0:255] = (vd[1:255],vd[0]) //rotate 255(-1)//
    __builtin_ve_lsv(ve, vl-1, b[i+vl]);  //vd[255] = b[i+vl] (=b[N-1]) //

    __builtin_ve_vfadd(ve, ve, vd, mask, vl);   
    __builtin_ve_vfnmsub(vd, vb, 2.0, ve, mask, vl);
    __builtin_ve_vfmadd(vc, vd, coef, vc, mask, vl);

    //store to array from vector resister//
    __builtin_ve_vst(vc, &c[i], 8, mask, vl);
}

c[N-1] += (0.0 - 2.0 * b[N-1] + b[N-2]) * coef;
```

今回はrotationを使って、1要素ズレた`vd`および`ve`を作りました。ただし、結局回り込んでしまった端の1要素に関しては、`b[i-1]`や`b[i+VL]`を追加で書き込んでいます。その命令は`__builtin_ve_lsv`です。この部分がちょっと損しているでしょうね。

とはいえ、256要素のロードは1度で済みました。

さて、この場合の計算時間は次のようになりました。
```
time  =0.033916 [s]
```
あれ？　ベクトル化しない場合よりは17倍も速いですが、自動ベクトル化や3回ロードする場合に比べると1/3の速度です。う～ん。難しい。

## 全ソースコード

最後に全ソースコードを載せます。`#if`で切り替えているので使う時はご注意ください。
```c++
#include <cstdio>
#include <cmath>
#include <chrono>
#include <algorithm>

//the definition of vector type. 2048 means 8 byte x 256 //
typedef double __vr __attribute__((__vector_size__(2048)));  


template <typename DURATION>
double DoubleSec(DURATION d) {
    return 1.0e-9 * (double)std::chrono::duration_cast<std::chrono::nanoseconds>(d).count();
}


int main(){
    constexpr int NK = 1000;
    constexpr int N = 100000;
    double* b = new double[N];   //function
    double* c = new double[N];   //result of second order derivative
    const double xi_min = -8.0;
    const double xi_max = 4.0;
    const double xi_delta = (xi_max - xi_min) / (double)N;

    //initialize//
    for(int i = 0; i < N; ++i){
        const double xi = xi_min + xi_delta * (double)i;
        const double r= exp(xi);
        b[i] = r*exp(-r);
        c[i] = 0.0;
    }

    auto start_tm = std::chrono::high_resolution_clock::now();


    for(int k=0; k < NK; ++k){
        
#if 0
        //original loop
        const double coef = 1.0/(2.0 * xi_delta*xi_delta);
        c[0] += (b[1] - 2.0 * b[0] + 0.0) * coef;
        for(int i = 1; i < N-1; ++i){
            c[i] += (b[i+1] - 2.0*b[i] + b[i-1]) * coef;
        }
        c[N-1] += (0.0 - 2.0 * b[N-1] + b[N-2]) * coef;
#elif 1
        //intrinsic (1) bからのloadを3回する場合//

        constexpr int VL = 256;

        //original loop
        const double coef = 1.0/(2.0 * xi_delta*xi_delta);

        c[0] += (b[1] - 2.0 * b[0] + 0.0) * coef;
            
        //(1) main part
        int i;
        for(i = 1; i <= N-1 - VL; i+=VL){            
            __vr vb;
            __builtin_ve_vld(vb, &b[i], 8);
            __vr vc;
            __builtin_ve_vld(vc, &c[i], 8);

            __vr vd;
            __builtin_ve_vld(vd, &b[i-1], 8); 
            
            __vr ve;
            __builtin_ve_vld(ve, &b[i+1], 8); 
            
            __builtin_ve_vfadd(ve, ve, vd);   
            __builtin_ve_vfnmsub(vd, vb, 2.0, ve);
            __builtin_ve_vfmadd(vc, vd, coef, vc);
            
            //store to array from vector resister//
            __builtin_ve_vst(vc, &c[i], 8);
        }   
        
        if(i < N-1){                
            const int vl = N-1 - i;
            __vm mask;
            __builtin_ve_eqvm(mask, mask, mask);//all bit stands// 

            __vr vb;
            __builtin_ve_vld(vb, &b[i], 8, mask, vl);
            __vr vc;
            __builtin_ve_vld(vc, &c[i], 8, mask, vl);

            __vr vd;
            __builtin_ve_vld(vd, &b[i-1], 8, mask, vl);

            __vr ve;
            __builtin_ve_vld(ve, &b[i+1], 8, mask, vl);

            __builtin_ve_vfadd(ve, ve, vd, mask, vl);   
            __builtin_ve_vfnmsub(vd, vb, 2.0, ve, mask, vl);
            __builtin_ve_vfmadd(vc, vd, coef, vc, mask, vl);

            //store to array from vector resister//
            __builtin_ve_vst(vc, &c[i], 8, mask, vl);
        }
        
        c[N-1] += (0.0 - 2.0 * b[N-1] + b[N-2]) * coef;
#else
        //intrinsic (2) bからのloadの回数を1回にしてrotateで代用する場合//

        constexpr int VL = 256;

        //original loop
        const double coef = 1.0/(2.0 * xi_delta*xi_delta);

        c[0] += (b[1] - 2.0 * b[0] + 0.0) * coef;
            
        //(1) main part
        int i;
        for(i = 1; i <= N-1 - VL; i+=VL){            
            __vr vb;
            __builtin_ve_vld(vb, &b[i], 8);
            __vr vc;
            __builtin_ve_vld(vc, &c[i], 8);

            __vr vd;
            __builtin_ve_vmv(vd, vb, 1);  //vd[0:255] = (vd[255],vd[0:254]) //rotate 1//
            __builtin_ve_lsv(vd, 0, b[i-1]);  //vd[0] = b[i-1] //

            __vr ve;
            __builtin_ve_vmv(ve, vb, 255);  //vd[0:255] = (vd[1:255],vd[0]) //rotate 255(-1)//
            __builtin_ve_lsv(ve, 255, b[i+VL]);  //vd[255] = b[i+VL] //

            __builtin_ve_vfadd(ve, ve, vd);   
            __builtin_ve_vfnmsub(vd, vb, 2.0, ve);
            __builtin_ve_vfmadd(vc, vd, coef, vc);
            
            //store to array from vector resister//
            __builtin_ve_vst(vc, &c[i], 8);
        }   
        
        if(i < N-1){                
            const int vl = N-1 - i;
            __vm mask;
            __builtin_ve_eqvm(mask, mask, mask);//all bit stands// 

            __vr vb;
            __builtin_ve_vld(vb, &b[i], 8, mask, vl);
            __vr vc;
            __builtin_ve_vld(vc, &c[i], 8, mask, vl);

            __vr vd;
            __builtin_ve_vmv(vd, vb, 1);  //vd[0:255] = (vd[255],vd[0:254]) //rotate 1//
            __builtin_ve_lsv(vd, 0, b[i-1]);  //vd[0] = b[i-1] //

            __vr ve;
            __builtin_ve_vmv(ve, vb, 255);  //vd[0:255] = (vd[1:255],vd[0]) //rotate 255(-1)//
            __builtin_ve_lsv(ve, vl-1, b[i+vl]);  //vd[255] = b[i+vl] (=b[N-1]) //

            __builtin_ve_vfadd(ve, ve, vd, mask, vl);   
            __builtin_ve_vfnmsub(vd, vb, 2.0, ve, mask, vl);
            __builtin_ve_vfmadd(vc, vd, coef, vc, mask, vl);

            //store to array from vector resister//
            __builtin_ve_vst(vc, &c[i], 8, mask, vl);
        }
        
        c[N-1] += (0.0 - 2.0 * b[N-1] + b[N-2]) * coef;
#endif

    }
    auto end_tm = std::chrono::high_resolution_clock::now();

    printf("c[N/2]=%f\n", c[N/2]);
    printf("time  =%f [s]\n", DoubleSec(end_tm -start_tm));
    
    delete[] b;
    delete[] c;
    return 0;
}
```

## 結果

あらためて四ケースのタイムです
```
time  =0.593170 [s]   no-vectorization
time  =0.011747 [s]   auto vectorization
time  =0.012196 [s]   intrinsic: 3 load
time  =0.033916 [s]   intrinsic: 1 load + rotation
```
自動ベクトル化が一番速いのは納得です（ベンダーのコンパイラ部門の方の努力に簡単に勝てるわけないので）が、3回ロードした方が1回ロード＋rotationよりも3倍も速いんですね。単にキャッシュに乗っているからなのか？

もう一点分からないのは、ロード・ストアにおいて、メモリ側のアライメント境界から読み書きすることがIntel AVX系では重要みたいですが、SX-Auroraの場合もアライメント境界か否かでコストが違うのだろうか？

## おわりに

Intrinsicによるベクトル化も、だんだんと慣れてきました。
面白くてキリがないのですが、

おっと、そろそろサンタさん来る時間なので寝ます。

メリークリスマス！