6

私はこれで頭を殴り続けています。Amatrixを matrixで乗算するための SSE ベースのアルゴリズムがありBます。A、B、またはその両方が転置される操作も実装する必要があります。以下に示す 4x4 マトリックス コード (これはかなり標準的な SSE 操作だと思います) の素朴な実装を行いましたが、A*B^T操作には約 2 倍の時間がかかりますA*B。ATLAS 実装は、 に対して同様の値を返しA*B、転置による乗算に対してほぼ同じ結果を返します。これは、これを行う効率的な方法があることを示唆しています。

MM-乗算:

m1 = (mat1.m_>>2)<<2;
n2 = (mat2.n_>>2)<<2;
n  = (mat1.n_>>2)<<2;

for (k=0; k<n; k+=4) {
  for (i=0; i<m1; i+=4) {
    // fetch: get 4x4 matrix from mat1
    // row-major storage, so get 4 rows
    Float* a0 = mat1.el_[i]+k;
    Float* a1 = mat1.el_[i+1]+k;
    Float* a2 = mat1.el_[i+2]+k;
    Float* a3 = mat1.el_[i+3]+k;

    for (j=0; j<n2; j+=4) {
      // fetch: get 4x4 matrix from mat2
      // row-major storage, so get 4 rows
      Float* b0 = mat2.el_[k]+j;
      Float* b1 = mat2.el_[k+1]+j;
      Float* b2 = mat2.el_[k+2]+j;
      Float* b3 = mat2.el_[k+3]+j;

      __m128 b0r = _mm_loadu_ps(b0);
      __m128 b1r = _mm_loadu_ps(b1);
      __m128 b2r = _mm_loadu_ps(b2);
      __m128 b3r = _mm_loadu_ps(b3);

      {  // first row of result += first row of mat1 * 4x4 of mat2
        __m128 cX1 = _mm_add_ps(_mm_mul_ps(_mm_load_ps1(a0+0), b0r), _mm_mul_ps(_mm_load_ps1(a0+1), b1r));
        __m128 cX2 = _mm_add_ps(_mm_mul_ps(_mm_load_ps1(a0+2), b2r), _mm_mul_ps(_mm_load_ps1(a0+3), b3r));
        Float* c0 = this->el_[i]+j;
        _mm_storeu_ps(c0, _mm_add_ps(_mm_add_ps(cX1, cX2), _mm_loadu_ps(c0)));
      }

      { // second row of result += second row of mat1 * 4x4 of mat2
        __m128 cX1 = _mm_add_ps(_mm_mul_ps(_mm_load_ps1(a1+0), b0r), _mm_mul_ps(_mm_load_ps1(a1+1), b1r));
        __m128 cX2 = _mm_add_ps(_mm_mul_ps(_mm_load_ps1(a1+2), b2r), _mm_mul_ps(_mm_load_ps1(a1+3), b3r));
        Float* c1 = this->el_[i+1]+j;
        _mm_storeu_ps(c1, _mm_add_ps(_mm_add_ps(cX1, cX2), _mm_loadu_ps(c1)));
      }

      { // third row of result += third row of mat1 * 4x4 of mat2
        __m128 cX1 = _mm_add_ps(_mm_mul_ps(_mm_load_ps1(a2+0), b0r), _mm_mul_ps(_mm_load_ps1(a2+1), b1r));
        __m128 cX2 = _mm_add_ps(_mm_mul_ps(_mm_load_ps1(a2+2), b2r), _mm_mul_ps(_mm_load_ps1(a2+3), b3r));
        Float* c2 = this->el_[i+2]+j;
        _mm_storeu_ps(c2, _mm_add_ps(_mm_add_ps(cX1, cX2), _mm_loadu_ps(c2)));
      }

      { // fourth row of result += fourth row of mat1 * 4x4 of mat2
        __m128 cX1 = _mm_add_ps(_mm_mul_ps(_mm_load_ps1(a3+0), b0r), _mm_mul_ps(_mm_load_ps1(a3+1), b1r));
        __m128 cX2 = _mm_add_ps(_mm_mul_ps(_mm_load_ps1(a3+2), b2r), _mm_mul_ps(_mm_load_ps1(a3+3), b3r));
        Float* c3 = this->el_[i+3]+j;
        _mm_storeu_ps(c3, _mm_add_ps(_mm_add_ps(cX1, cX2), _mm_loadu_ps(c3)));
      }
  }
// Code omitted to handle remaining rows and columns
}

MT 乗算 (転置行列で乗算された行列) では、代わりに次のコマンドを使用して b0r から b3r に保存し、ループ変数を適切に変更しました。

__m128 b0r = _mm_set_ps(b3[0], b2[0], b1[0], b0[0]);
__m128 b1r = _mm_set_ps(b3[1], b2[1], b1[1], b0[1]);
__m128 b2r = _mm_set_ps(b3[2], b2[2], b1[2], b0[2]);
__m128 b3r = _mm_set_ps(b3[3], b2[3], b1[3], b0[3]);

スローダウンの原因の一部は、一度に行を取得することと、列を取得するために毎回 4 つの値を格納する必要があることの違いによるものだと思われますが、これについては別の方法で B の行を取得し、次に As の列を掛けると、コストが 4 列の結果を格納するようにシフトされます。

また、B の行を行として取り込んで_MM_TRANSPOSE4_PS(b0r, b1r, b2r, b3r);から、転置を実行しようとしました (そのマクロには追加の最適化があるのではないかと考えました) が、実際の改善はありません。

表面的には、これはより高速であるべきだと思います...関係するドット積は行ごとに実行されるため、本質的により効率的ですが、ドット積をまっすぐに実行しようとすると、同じことをしなければならなくなります結果を保存します。

ここで何が欠けていますか？

追加:明確にするために、行列を転置しないようにしています。それらに沿って反復したいと思います。問題は、私が知る限り、_mm_set_ps コマンドが _mm_load_ps よりもはるかに遅いことです。

また、A 行列の 4 行を格納し、1 つのロード、4 つの乗算、2 つの加算を含む 4 つの中かっこで囲まれたセグメントを 4 つの乗算命令と 3 で置き換えるバリエーションも試しましたhaddsが、ほとんど役に立ちませんでした。タイミングは同じままです (もちろん、テスト コンパイルでコードが変更されたことを確認するために debug ステートメントを使用して試しました。もちろん、プロファイリングの前にデバッグ ステートメントは削除されました)。

    {  // first row of result += first row of mat1 * 4x4 of mat2
      __m128 cX1 = _mm_hadd_ps(_mm_mul_ps(a0r, b0r), _mm_mul_ps(a0r, b1r));
      __m128 cX2 = _mm_hadd_ps(_mm_mul_ps(a0r, b2r), _mm_mul_ps(a0r, b3r));
      Float* c0 = this->el_[i]+j;
      _mm_storeu_ps(c0, _mm_add_ps(_mm_hadd_ps(cX1, cX2), _mm_loadu_ps(c0)));
    }

    { // second row of result += second row of mat1 * 4x4 of mat2
      __m128 cX1 = _mm_hadd_ps(_mm_mul_ps(a1r, b0r), _mm_mul_ps(a1r, b1r));
      __m128 cX2 = _mm_hadd_ps(_mm_mul_ps(a1r, b2r), _mm_mul_ps(a1r, b3r));
      Float* c0 = this->el_[i+1]+j;
      _mm_storeu_ps(c0, _mm_add_ps(_mm_hadd_ps(cX1, cX2), _mm_loadu_ps(c0)));
    }

    { // third row of result += third row of mat1 * 4x4 of mat2
      __m128 cX1 = _mm_hadd_ps(_mm_mul_ps(a2r, b0r), _mm_mul_ps(a2r, b1r));
      __m128 cX2 = _mm_hadd_ps(_mm_mul_ps(a2r, b2r), _mm_mul_ps(a2r, b3r));
      Float* c0 = this->el_[i+2]+j;
      _mm_storeu_ps(c0, _mm_add_ps(_mm_hadd_ps(cX1, cX2), _mm_loadu_ps(c0)));
    }

    { // fourth row of result += fourth row of mat1 * 4x4 of mat2
      __m128 cX1 = _mm_hadd_ps(_mm_mul_ps(a3r, b0r), _mm_mul_ps(a3r, b1r));
      __m128 cX2 = _mm_hadd_ps(_mm_mul_ps(a3r, b2r), _mm_mul_ps(a3r, b3r));
      Float* c0 = this->el_[i+3]+j;
      _mm_storeu_ps(c0, _mm_add_ps(_mm_hadd_ps(cX1, cX2), _mm_loadu_ps(c0)));
    }

更新：a0rそうです。レジスタのスラッシングを回避するためにto の行のロードをa3r中括弧に移動することも失敗しました。

4

2 に答える 2

2

役立つかもしれないいくつかの推奨事項：

  • アラインされていないメモリを使用しないでください (これらの _mm_loadu* は低速で​​す)。
  • メモリに順次アクセスしていないため、キャッシュが失われます。そのメモリへの実際のアクセスの前に行列を転置してみてください。これにより、CPU が可能な限りキャッシュをフェッチして使用できるようになります。このように、以下__m128 b0r = _mm_set_ps(b3[0], b2[0], b1[0], b0[0]); // and b1r, etc..は必要ありません。アイデアは、4 つのコンポーネント全体を順番に取得することです。SSE コードが呼び出される前にメモリを再編成する必要がある場合は、そうしてください。
  • これを内側のループにロードしています: _mm_load_ps1(a0+0) (a1、a2、a3 についても同じ)しかし、内側のループのすべての反復で一定です。これらの値を外部にロードして、いくつかのサイクルを保存できます。以前のイテレーションから再利用できるものに注目してください。
  • プロフィール。Intel VTune などを使用すると、ボトルネックがどこにあるかがわかります。
于 2013-05-14T20:48:29.677 に答える