bugfix> c++ > 投稿

私はnumpy.triu_indices(a、1)(2番目の引数が1であることに注意してください)をavx組み込み関数でc ++に実装したいと思います。以下のコードスニペットは、私が思いついたコードのベクトル化されていないバージョンです。ここで、aは長さ(int)で、最初と2番目は2つの出力配列です

void triu_indices(int a, int* first, int* second)
{
    int index = 0;
    for (int i=0; i < a; i++)
    {
        for(int j = i+1; j < a; j++)
        {
            first[index] = i;
            second[index] = j;
            index++;
        }
    }
}

出力例として、a = 4を指定すると、

first = [0,0,0,1,1,2]
second = [1,2,3,2,3,3]

次に、これをAVX2で完全に実装します(ベクトル化された方法です)。最終的に、関数はintの配列全体で実行され、変数 a が提供されます関数、および出力配列 first および second 2つの親配列に格納されます。

明示的なAVX2組み込み関数を使用してこの関数をベクトル化する方法(つまり、コンパイラの自動ベクトル化に依存しない)として、役立つヒント(またはコードスニペット)を教えてください。 AVXを最近学習し始めたので、これがnoobの質問であれば申し訳ありません。

回答 1 件
  • まず、あなたを確認してください本当に これを行う必要があり、実際にインデックスの配列を最終結果として求めます。三角行列のデータを追跡することの一部としてではありません。 AVX2にはギャザーサポートがあり、AVX512には分散サポートがありますが、インデックスの配列を導入するとSIMDがさらに悪化します。

    三角行列のループ、およびi、jから線形へのループについては、アセンブリを使用した三角行列メモリのアドレス指定アルゴリズムを参照してください。 (各行が32バイトの境界で始まるようにインデックス付けをパディングしたい場合があります。つまり、各行の長さをAVXベクトル全体の8つの浮動小数点要素の倍数に切り上げます。これにより、 AVXベクトルを含む行列:行の最後のベクトルに次の行の先頭からの要素を含める代わりに、行の最後のパディングにガベージを格納できます。

    線形の場合-> i、j、閉形式の公式には sqrt が含まれます  (これもC ++バージョンです)、配列ルックアップが役立つ可能性がありますが、実際にはこれを行わないでください。 (たとえば、パック形式の三角行列をループする場合は、探している要素を見つけたときにルックアップが必要ないように、i、jと線形のどこにいるかを追跡してください)


    これが必要な場合:

    大きな配列の場合、ベクトル全体にかなりうまく分割され、 行の最後でのみトリッキーになります。

    4または8 int の同じベクトル内に複数の三角形の行がある最後のコーナーの特別な場合には、事前定義されたベクトル定数を使用できます。  要素。

    first = [0,0,0,1,1,2]
    
    

    大きな三角形の場合、同じ数の大きなランを保存しています( memset など) )、次の番号のわずかに短い実行など、つまり 0 の行全体を保存する sは簡単です。最後の数行を除くすべての行について、これらの実行は1つのベクトル要素よりも大きくなります。

    second = [1,2,3,2,3,3]
    
    

    繰り返しますが、単一の行内では、ベクトル化する単純なパターンです。  増加するシーケンスを保存するには、 {1,2,3,4} のベクトルから始めます  SIMDでインクリメントし、 {4,4,4,4} で追加します 、つまり _mm_set1_epi32(1) 。 256ビットAVX2ベクトルの場合、 _mm256_set1_epi32(8)  8要素のベクトルを8ずつインクリメントします。

    したがって、最も内側のループ内では、1つの不変ベクトルを格納し、 _mm256_add_epi32 を使用しています。  別のアレイに保存します。

    コンパイラはすでにオート-関数をかなり適切にベクトル化する、行末処理は手作業で行うよりもはるかに複雑ですが。 Godboltコンパイラエクスプローラーのコード( __restrict を使用)  出力配列がオーバーラップしないことをコンパイラーに伝えるため、および __builtin_assume_aligned  コンパイラーに整列されていることを伝えるため)、次のような内部ループを取得します(gccから):

    .L4:                                       # do {
        movups  XMMWORD PTR [rcx+rax], xmm0      # _mm_store_si128(&second[index], xmm0)
        paddd   xmm0, xmm2                       # _mm_add_epi32
        movups  XMMWORD PTR [r10+rax], xmm1      # _mm_store_si128(&second[index], counter_vec)
        add     rax, 16                          # index += 4  (16 bytes)
        cmp     rax, r9
        jne     .L4                            # }while(index != end_row)
    
    

    時間がある場合は、行の終わりの処理を改善するなど、これをより詳細に記述できます。例えば行の終わりで終了する部分的に重複するストアは、多くの場合良好です。または、外側のループを展開して、内側のループがクリーンアップの繰り返しパターンを持つようにします。

    次の外側のループの反復の開始ベクトルの計算は、次のような方法で実行できます。

    vsecond = _mm256_add_epi32(vsecond, _mm256_set1_epi32(1));
    vfirst  = _mm256_add_epi32(vfirst, _mm256_set1_epi32(1));
    
    

    つまり、 {0,0,0,0,...} を回す   {1,1,1,1,...} へ  すべてのベクトルを追加します。そして、 {1,2,3,4,5,6,7,8} をオンにします   {2,3,4,5,6,7,8,9} へ  すべてのベクトルを追加します。

あなたの答え