LoginSignup
2
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

@pochman

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

Posted at

はじめに

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

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

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

ターゲットとする問題

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

ベクトルする前のコードのメインループ部分は次のようなものです。

//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ループ分だけです。以下のようにしてみました。

//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なのですね。

これを使ったコードは次のようになりました。

//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で切り替えているので使う時はご注意ください。

#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によるベクトル化も、だんだんと慣れてきました。
面白くてキリがないのですが、

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

メリークリスマス！

2
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
2
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?