2

行列とベクトルおよび行列と行列の乗算関数を作成する必要がありますが、SSE コマンドについて頭を悩ませることができません。

行列とベクトルの次元は常に 4 の倍数です。

次のようなベクトル間乗算関数を作成できました。

void vector_multiplication_SSE(float* m, float* n, float* result, unsigned const int size)
{
    int i;

    __declspec(align(16))__m128 *p_m = (__m128*)m;
    __declspec(align(16))__m128 *p_n = (__m128*)n;
    __declspec(align(16))__m128 *p_result = (__m128*)result;

    for (i = 0; i < size / 4; ++i)
        p_result[i] = _mm_mul_ps(p_m[i], p_n[i]);

    // print the result
    for (int i = 0; i < size; ++i)
    {
        if (i % 4 == 0) cout << endl;
        cout << result[i] << '\t';
    }
}

そして今、行列とベクトルの乗算を実装しようとしています。

これが私がこれまでに持っているものです：

void multiply_matrix_by_vector_SSE(float* m, float* v, float* result, unsigned const int vector_dims)
{
    int i, j;

    __declspec(align(16))__m128 *p_m = (__m128*)m;
    __declspec(align(16))__m128 *p_v = (__m128*)v;
    __declspec(align(16))__m128 *p_result = (__m128*)result;

    for (i = 0; i < vector_dims; i += 4)
    {
        __m128 tmp = _mm_load_ps(&result[i]);
        __m128 p_m_tmp = _mm_load_ps(&m[i]);

        tmp = _mm_add_ps(tmp, _mm_mul_ps(tmp, p_m_tmp));
        _mm_store_ps(&result[i], tmp);

        // another for loop here? 
    }

    // print the result
    for (int i = 0; i < vector_dims; ++i)
    {
        if (i % 4 == 0) cout << endl;
        cout << result[i] << '\t';
    }
}

この関数は完全に間違っているように見えます。つまり、正しく機能しないだけでなく、間違った方向に進んでいるようにも見えます。

ベクトル行列と行列行列の乗算の実装を手伝ってくれる人はいますか? いくつかのサンプルコードと非常に詳細な説明をいただければ幸いです

アップデート

ここに私の試み番号2があります：

例外で失敗しAccess reading violationますが、それでも近いと感じます

void multiply_matrix_by_vector_SSE(float* m, float* v, float* result, unsigned const int vector_dims)
{
    int i, j;

    __declspec(align(16))__m128 *p_m = (__m128*)m;
    __declspec(align(16))__m128 *p_v = (__m128*)v;
    __declspec(align(16))__m128 *p_result = (__m128*)result;

    for (i = 0; i < vector_dims; ++i)
    {
        p_result[i] = _mm_mul_ps(_mm_load_ps(&m[i]), _mm_load_ps1(&v[i]));
    }

    // print the result
    for (int i = 0; i < vector_dims; ++i)
    {
        if (i % 4 == 0) cout << endl;
        cout << result[i] << '\t';
    }
}

更新 2

void multiply_matrix_by_vector_SSE(float* m, float* v, float* result, unsigned const int vector_dims)
{
    int i, j;
    __declspec(align(16))__m128 *p_m = (__m128*)m;
    __declspec(align(16))__m128 *p_v = (__m128*)v;
    __declspec(align(16))__m128 *p_result = (__m128*)result;

    for (i = 0; i < vector_dims; ++i)
    {
        for (j = 0; j < vector_dims * vector_dims / 4; ++j)
        {
            p_result[i] = _mm_mul_ps(p_v[i], p_m[j]);
        }
    }

    for (int i = 0; i < vector_dims; ++i)
    {
        if (i % 4 == 0) cout << endl;
        cout << result[i] << '\t';
    }
    cout << endl;
}
4

1 に答える 1

8

トリックなど何もなければ、行列とベクトルの乗算は、ベクトルと行列の行の間の内積の集まりにすぎません。あなたのコードは実際にはその構造を持っていません。実際に内積として書きます（テストされていません）：

for (int row = 0; row < nrows; ++row) {
    __m128 acc = _mm_setzero_ps();
    // I'm just going to assume the number of columns is a multiple of 4
    for (int col = 0; col < ncols; col += 4) {
        __m128 vec = _mm_load_ps(&v[col]);
        // don't forget it's a matrix, do 2d addressing
        __m128 mat = _mm_load_ps(&m[col + ncols * row]);
        acc = _mm_add_ps(acc, _mm_mul_ps(mat, vec));
    }
    // now we have 4 floats in acc and they have to be summed
    // can use two horizontal adds for this, they kind of suck but this
    // isn't the inner loop anyway.
    acc = _mm_hadd_ps(acc, acc);
    acc = _mm_hadd_ps(acc, acc);
    // store result, which is a single float
    _mm_store_ss(&result[row], acc);
}

一度に複数の行を処理する、ベクターからのロードを再利用する、いくつかの独立した依存チェーンを作成してスループットをより有効に活用するなど、いくつかの明らかなトリックがあります (以下を参照)。また、mul/add コンボに FMA を使用するという非常に単純なトリックもありますが、サポートはまだそれほど普及していません (2015 年にはサポートされていませんでしたが、2020 年にはかなり普及しています)。

これから行列 - 行列の乗算を作成できますが (結果の場所を変更した場合)、最適ではありません (以下を参照)。

一度に 4 つの行を取得する (テストされていません):

for (int row = 0; row < nrows; row += 4) {
    __m128 acc0 = _mm_setzero_ps();
    __m128 acc1 = _mm_setzero_ps();
    __m128 acc2 = _mm_setzero_ps();
    __m128 acc3 = _mm_setzero_ps();
    for (int col = 0; col < ncols; col += 4) {
        __m128 vec = _mm_load_ps(&v[col]);
        __m128 mat0 = _mm_load_ps(&m[col + ncols * row]);
        __m128 mat1 = _mm_load_ps(&m[col + ncols * (row + 1)]);
        __m128 mat2 = _mm_load_ps(&m[col + ncols * (row + 2)]);
        __m128 mat3 = _mm_load_ps(&m[col + ncols * (row + 3)]);
        acc0 = _mm_add_ps(acc0, _mm_mul_ps(mat0, vec));
        acc1 = _mm_add_ps(acc1, _mm_mul_ps(mat1, vec));
        acc2 = _mm_add_ps(acc2, _mm_mul_ps(mat2, vec));
        acc3 = _mm_add_ps(acc3, _mm_mul_ps(mat3, vec));
    }
    acc0 = _mm_hadd_ps(acc0, acc1);
    acc2 = _mm_hadd_ps(acc2, acc3);
    acc0 = _mm_hadd_ps(acc0, acc2);
    _mm_store_ps(&result[row], acc0);
}

現在、4 つの FMA ごとに 5 つのロードしかありませんが、行が展開されていないバージョンでは 1 つの FMA ごとに 2 つのロードがありました。また、4 つの独立した FMA、または FMA 収縮のない add/mul ペアがあり、いずれにしても、パイプライン化/同時実行の可能性が高まります。たとえば、Skylake はサイクルごとに 2 つの独立した FMA を開始でき、完了するまでに 4 サイクルかかるため、両方の FMA ユニットを完全に占有するには、8 つの独立した FMA が必要です。おまけとして、これらの 3 つの水平方向の追加は、最終的には水平方向の合計に対して比較的うまく機能します。

異なるデータ レイアウトは最初は欠点のように思えます。行列とベクトルの両方から単純にベクトル ロードを実行し、それらを乗算することはもはや不可能です (最初の行列の小さな行ベクトルを、最初の行列の小さな行ベクトルで乗算することになります)。 2 番目の行列、これは間違っています)。しかし、完全な行列 - 行列乗算は、本質的に行列に多くの独立したベクトルを乗算しているという事実を利用できます。これは、実行する必要がある独立した作業でいっぱいです。水平合計も簡単に回避できます。実際には、行列とベクトルの乗算よりもさらに便利です。

重要なのは、行列 A から小さな列ベクトルを取り、行列 B から小さな行ベクトルを取り、それらを小さな行列に乗算することです。これは、慣れ親しんだ方法とは逆に聞こえるかもしれませんが、計算が常に独立しており、水平方向の操作がないため、SIMD ではこの方法がうまく機能します。

例（テストされていません。行列が展開係数で割り切れる次元を持っていると仮定し、x64 が必要です。それ以外の場合はレジスタが不足します）

for (size_t i = 0; i < mat1rows; i += 4) {
    for (size_t j = 0; j < mat2cols; j += 8) {
        float* mat1ptr = &mat1[i * mat1cols];
        __m256 sumA_1, sumB_1, sumA_2, sumB_2, sumA_3, sumB_3, sumA_4, sumB_4;
        sumA_1 = _mm_setzero_ps();
        sumB_1 = _mm_setzero_ps();
        sumA_2 = _mm_setzero_ps();
        sumB_2 = _mm_setzero_ps();
        sumA_3 = _mm_setzero_ps();
        sumB_3 = _mm_setzero_ps();
        sumA_4 = _mm_setzero_ps();
        sumB_4 = _mm_setzero_ps();

        for (size_t k = 0; k < mat2rows; ++k) {
            auto bc_mat1_1 = _mm_set1_ps(mat1ptr[0]);
            auto vecA_mat2 = _mm_load_ps(mat2 + m2idx);
            auto vecB_mat2 = _mm_load_ps(mat2 + m2idx + 4);
            sumA_1 = _mm_add_ps(_mm_mul_ps(bc_mat1_1, vecA_mat2), sumA_1);
            sumB_1 = _mm_add_ps(_mm_mul_ps(bc_mat1_1, vecB_mat2), sumB_1);
            auto bc_mat1_2 = _mm_set1_ps(mat1ptr[N]);
            sumA_2 = _mm_add_ps(_mm_mul_ps(bc_mat1_2, vecA_mat2), sumA_2);
            sumB_2 = _mm_add_ps(_mm_mul_ps(bc_mat1_2, vecB_mat2), sumB_2);
            auto bc_mat1_3 = _mm_set1_ps(mat1ptr[N * 2]);
            sumA_3 = _mm_add_ps(_mm_mul_ps(bc_mat1_3, vecA_mat2), sumA_3);
            sumB_3 = _mm_add_ps(_mm_mul_ps(bc_mat1_3, vecB_mat2), sumB_3);
            auto bc_mat1_4 = _mm_set1_ps(mat1ptr[N * 3]);
            sumA_4 = _mm_add_ps(_mm_mul_ps(bc_mat1_4, vecA_mat2), sumA_4);
            sumB_4 = _mm_add_ps(_mm_mul_ps(bc_mat1_4, vecB_mat2), sumB_4);
            m2idx += 8;
            mat1ptr++;
        }
        _mm_store_ps(&result[i * mat2cols + j], sumA_1);
        _mm_store_ps(&result[i * mat2cols + j + 4], sumB_1);
        _mm_store_ps(&result[(i + 1) * mat2cols + j], sumA_2);
        _mm_store_ps(&result[(i + 1) * mat2cols + j + 4], sumB_2);
        _mm_store_ps(&result[(i + 2) * mat2cols + j], sumA_3);
        _mm_store_ps(&result[(i + 2) * mat2cols + j + 4], sumB_3);
        _mm_store_ps(&result[(i + 3) * mat2cols + j], sumA_4);
        _mm_store_ps(&result[(i + 3) * mat2cols + j + 4], sumB_4);
    }
}

このコードのポイントは、計算を非常に SIMD に適したものにするのが簡単で、浮動小数点ユニットを飽和させる独立した演算を多数使用し、同時に比較的少ない負荷を使用することです (そうしないとボトルネックになる可能性があります)。数が多すぎるだけで、L1 キャッシュを見逃す可能性があることはさておき)。

このコードを使用することもできますが、インテル® MKL と競合するものではありません。特に、タイリングが非常に重要な中規模または大規模なマトリックスの場合。これを AVX にアップグレードするのは簡単です。小さな行列にはまったく適していません。たとえば、2 つの 4x4 行列を乗算するには、効率的な 4x4 行列の乗算を参照してください。

于 2015-11-19T19:54:09.723 に答える