※NECさんより正確な情報を頂くことができたので、一部を修正しました。ただし、当時の試行錯誤も自身にとっての記録として重要なので、最低限の修正に留めました。

はじめに

マスクについて使えるようになったので、いよいよループのベクトル化に挑戦してみます。

ベクトル演算では、256要素の演算を一度に行いますが、実はループ長Nをベクトル長256で割り切れない場合、あまりの部分が発生します。あまりの部分をリマインダーループと呼ぶことにします。

このリマインダループをスカラー処理していては、せっかくのそれ以前の要素をベクトル演算しても割に合いません。SX-Auroraのintrinsic命令では、256未満の要素に対するベクトル演算もできるように命令セットが工夫されています。

ループのベクトル化からリマインダーループの扱いまでを順に試していきましょう。

ベクトル化する前のコード

まずは、ベクトル化する前のターゲットコードを示します。単純に配列bと配列cの足し算をする処理です。

int main(){
    constexpr int N = 1000;
    double b[N]={0.0};
    double c[N]={0.0};
    
    //initialize//
    for(int i = 0; i < N; ++i){
        b[i] = (double)i;
        c[i] = 1000.0;
    }

    /*main operation*/
    for(int i = 0; i < N; ++i){
        b[i] = b[i] + c[i];
    }
    
    for(int i = 0; i < N; ++i){        
        printf("%f\n", b[i]);
    }
    
    return 0;
}

実行すると、1000.0から1999.0までの数値が出力されます。

このコードの中で、初期化部分は今回は無視して、main operationと記載した足し算のループの部分をベクトル化していきます。

ループのベクトル化

例として、以下では一貫して配列の長さはN=1000とします。

あらためて、ベクトル演算では256要素ずつ一度に処理するので、ループのカウンタも256ずつ一気に進めます。つまり、

constexpr int VL = 256;
int i;
//(1)main part
for(i = 0; i <= N-VL; i+=VL){
//...ベクトル演算...//
}

//(2)reminder part
if(i < N){
 //この時点でのiは N - Nを256で割った余り値 (i == N - (N % 256))
 //もしくは (i == (N/256)*256 = (N >> 8) << 8 )
//...まだここで 残ったi ~ N-1の要素への処理が必要...//
}

となります。ループ後にもカウンタ変数iが参照できるように、iの定義をループ外で書いています。

メイン部分(1)のループの範囲に注意してください。

ループを抜けた直後の(2)の時点では、iの値はN=1000を256で割ったあまりをNからひいたもの、つまり、i==(N/256)*256=768になっているはずです。

メイン部分(1)のベクトル化

このリマインダー部分の232要素は、まずは置いておいて、メインのループ部分をベクトル化してみます。

結果から書いてしまうと次のようになります。

    constexpr int VL = 256;
    //(1) main part
    for(int i = 0; i <= N-VL; i+=VL){
        //load from array to vector resister//
        __vr vb;
        __builtin_ve_vld(vb, &amp;b[i], 8);
        __vr vc;
        __builtin_ve_vld(vc, &amp;c[i], 8);
        
        //vector add operation//
        __vr va;
        __builtin_ve_vfadd(va, vb,vc);

        //store to array from vector resister//
        __builtin_ve_vst(va, &amp;b[i], 8);
    }

ベクトル型の変数vbとvcを用意して、__builtin_ve_vld命令で配列bおよびcから256要素を読み込みます。

続いて、__builtin_ve_vfadd命令でvbとvcの各要素を足し算し、vaに結果を代入します。

最後に、__builtin_ve_vst命令でvaから配列bに書き込みます。

どうでしょう。普通のプログラムと比べてロードとストアが面倒ですが、メモリとレジスタ間のロード・ストアと思えば直感的でしょうか。やってるいることはシンプルだとおもいます。

リマインダー部分(2)のベクトル化

さて、リマインダー部分ですが、これをスカラー処理のループ(リマインダーループ)として次のように書いてしまうこともできます。

    //(2') here, i == (N >> 8) << 8  //
    for(   ; i < N; ++i){
        b[i] = b[i] + c[i];
    }

ただこれだとAVXのようにベクトル長の短い物なら良いですが、SX-Auroraのようにベクトル長が長い場合には、このリマインダーループでせっかくベクトル化が台無しになりそうです。

そこで、SX-Auroraでは、256未満の要素に対しても一度に処理できるよう、ベクトル命令のオプションがあります。

例えば、足し算をする__builtin_ve_vfadd命令では引数は3つから5つまで可変でした。

void __builtin_ve_vfadd(__vr dest, __vr src1, __vr src2
                              [, __vm mask][, int vl]);

ここで、第4引数は前回まで扱ったマスクでした。第5引数が実際に演算するベクトル長で、ここに256未満の数字を入れれば、それ以降の要素は無視するという事です。

これだけみると、今回の場合はN % 256 = N - i = 232を第5引数にセットすれば良さそうですが、果たしてうまくいくのか。気になるのはロード・ストアで、Segmentation faultになったりしないのでしょうか？やってみましょう。

とにかく、第5引数を使うとリマインダー部分のコードはループではなく1度の処理として次のように書けそうです。

//(2)reminder part
    if(i < N){
        const int vl = N - i;

        //all bit stands//
        __vm mask;
        __builtin_ve_eqvm(mask, mask, mask);
        
        //load from array to vector resister//
        __vr vb;
        __builtin_ve_vld(vb, &b[i], 8, mask, vl);
        __vr vc;
        __builtin_ve_vld(vc, &c[i], 8, mask, vl);
        
        //vector add operation//
        __vr va;
        __builtin_ve_vfadd(va, vb,vc, mask, vl);

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

余り部分の要素の長さをvl = N-iで定義しています。

足し算だけでなく、ベクトリ型へのロード・ストア命令でも、第5引数で要素数vlを渡しています。

ところで、問題は、第5引数を渡すためには第4引数も渡さなければならないということです。一方で今回はマスクによる演算の無効化は時に必要ないので、マスク変数として全てのbitを立てたものを用意しています。必要ないのに渡すというのもなんとも面倒な気もしますが、仕方ないですね（実はこのために、先にマスク編の記事を書きました。）

追記： マスクの指定を省略して、第4引数に要素数vlを入力するという記述も可能だそうです。上記の場合では例えば次のように書けます。

__builtin_ve_vld(vb, &b[i], 8, mask, vl);
↓
__builtin_ve_vld(vb, &b[i], 8, vl);

結果

最終的なコードは次のようになります。試したところ、一応これで上手く動きました。

#include <cstdio>

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

int main(){
    constexpr int N = 1000;
    double b[N]={0.0};
    double c[N]={0.0};

    //initialize//
    for(int i = 0; i < N; ++i){
        b[i] = (double)i;
        c[i] = 1000.0;
    }

    /*original loop
    for(int i = 0; i < N; ++i){
        b[i] = b[i] + c[i];
    }
    */

    constexpr int VL = 256;
    //(1) main part
    int i;
    for(i = 0; i <= N - VL; i+=VL){
        //load from array to vector resister//
        __vr vb;
        __builtin_ve_vld(vb, &b[i], 8);
        __vr vc;
        __builtin_ve_vld(vc, &c[i], 8);
        
        //vector add operation//
        __vr va;
        __builtin_ve_vfadd(va, vb, vc);

        //store to array from vector resister//
        __builtin_ve_vst(va, &b[i], 8);
    }
    
    //(2)reminder part
    if(i < N){
        const int vl = N - i;
        __vm mask;
        __builtin_ve_eqvm(mask, mask, mask);//all bit stands//
        //load from array to vector resister//
        __vr vb;
        __builtin_ve_vld(vb, &b[i], 8, mask, vl);
        __vr vc;
        __builtin_ve_vld(vc, &c[i], 8, mask, vl);
        
        //vector add operation//
        __vr va;
        __builtin_ve_vfadd(va, vb,vc, mask, vl);

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

    for(int i = 0; i < N; ++i){        
        printf("%f\n", b[i]);
    }
    
    return 0;
}

その他の書き方

今回はリマインダー部分を別途書きました。一方で、webのintrinsic チュートリアル等を参考にすると、次のようにリマインダー部分を分けずに1つのループで書いてしまう方法もあります。

#include <cstdio>
#include <algorithm>

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

int main(){
    constexpr int N = 1000;
    double b[N]={0.0};
    double c[N]={0.0};

    //initialize//
    for(int i = 0; i < N; ++i){
        b[i] = (double)i;
        c[i] = 1000.0;
    }

    /*original loop
    for(int i = 0; i < N; ++i){
        b[i] = b[i] + c[i];
    }
    */

    constexpr int VL = 256;
    //(1) main part
    int i;
    for(i = 0; i < N; i+=VL){
        const int vl = std::min(VL, N-i);
        __vm mask;
        __builtin_ve_eqvm(mask, mask, mask);//all bit stands//
        //load from array to vector resister//
        __vr vb;
        __builtin_ve_vld(vb, &b[i], 8, mask, vl);
        __vr vc;
        __builtin_ve_vld(vc, &c[i], 8, mask, vl);
        
        //vector add operation//
        __vr va;
        __builtin_ve_vfadd(va, vb,vc, mask, vl);

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

    for(int i = 0; i < N; ++i){        
        printf("%f\n", b[i]);
    }
    
    return 0;
}

この方法だと、リマインダー部分を別建てで書かなくて見やすいです。またループの範囲もoriginalと変わらないので間違えにくいですね。

一方で、この方法で心配なのは、毎回の処理にマスク処理が適用されていることです。これによるスループットの速度の低下があるのかないのか良く分からないので、個人的には前者の方法でリマインダー部を分けるのが良いかなと思います。

おわりに

いよいよベクトル化できるようになりました。

ここまで何度か行き詰りましたが、なんとか辿り着けてよかったです。

メリークリスマス！