전치 행렬을보다 효율적으로 곱하는 방법은 무엇입니까?

Question

나는 이것으로 머리를 때리고있다. 매트릭스 A에 매트릭스 B을 곱하는 SSE 기반 알고리즘이 있습니다. 또한 A, B 또는 둘 모두가 바뀌는 위치에 대한 작업을 구현해야합니다. 나는 그것의 순진한 구현을했다. 아래에 표현 된 4x4 행렬 코드 (꽤 표준적인 SSE 연산이다.)는 A*B^T 연산은 A*B의 두 배 정도 걸린다. ATLAS 구현은 A*B에 대해 비슷한 값을 반환하고, 조인에 의한 곱셈과 거의 동일한 결과를 보여 주므로 효율적인 방법이 있음을 나에게 알립니다.

MM 승산 다음 MT 승산 (전치 행렬을 곱한 행렬) 용

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 
} 

는, I는 다음 명령 b3r하는 b0r 저장 적절 루프 변수 변경 대신 :

__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 개의 곱하기 및 4 개의 곱하기 명령어와 3 개를 포함한 2 개의 덧셈을 포함하는 4 개의 중괄호로 묶은 세그먼트를 대체했지만 조금 유용했습니다 . 타이밍은 동일하게 유지 (그래, 나는 코드가 내 테스트 컴파일 변경했다고 확인하기 위해 디버그 문으로 그것을 시도 말했다 디버그 문은 물론, 프로파일 링하기 전에 제거되었습니다.) :

{ // 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))); 
    } 

업데이트 : 의 행로드를 a3r으로 중괄호로 이동하여 레지스터 스래 싱을 피하려는 시도도 실패했습니다.

Answer 1

제가 생각하기에 이것은 수평 적 추가가 유용한 경우입니다. 당신은 C = A B^T를 원하지만 B는 조바꿈으로 메모리에 저장되지 않습니다. 그것이 문제이다. SoA 대신 AoS와 같은 상점입니다. 이 경우 B의 조옮김을 취하고 수직 추가를하는 것이 수평 추가를 사용하는 것보다 느립니다. 매트릭스 벡터 Efficient 4x4 matrix vector multiplication with SSE: horizontal add and dot product - what's the point? 적어도 사실입니다. 함수 m4x4 아래 코드 m4x4_vec는 m4x4T가 C는 SSE 않고 B^T를 = 않으며 m4x4T_vec는 C가 B^T usisg은 SSE를 =하지 SSE를 사용하며, 비 SSE 4 × 4 매트릭스 제품이다. 마지막 하나는 내가 생각하기를 원하는 것입니다.

참고 : 큰 행렬에 대해서는이 방법을 사용하지 않습니다. 이 경우 트랜스퍼를 먼저 사용하고 세로 추가를 사용하는 것이 더 빠릅니다 (SSE/AVX를 사용하면 좀 더 복잡한 작업을 수행 할 수 있습니다. 스트립을 SSE/AVX 폭으로 조 변경). 왜냐하면 전치가 O (n^2)가되고 행렬 곱이 O (n^3)로 바뀌기 때문입니다. 따라서 큰 행렬의 경우 전치가 중요하지 않습니다. 그러나 4x4의 경우 전치가 중요하므로 수평 추가가 우선합니다.

편집 : 나는 당신이 원하는 것을 오해했습니다. 당신은 C = (A B)^T를 원한다. 즉 (A B)만큼 빠르게해야하며 코드는 기본적으로 그냥 다음과 같이 우리는 수학을 쓸 수 A와 B
의 역할을 바꿀 거의 동일합니다

C = A*B in Einstein notation is C_i,j = A_i,k * B_k,j. 
Since (A*B)^T = B^T*A^T we can write 
C = (A*B)^T in Einstein notation is C_i,j = B^T_i,k * A^T_k,j = A_j,k * B_k,i 

비교할 경우 변경되는 유일한 두 가지는 j와 i의 역할을 교환하는 것입니다. 이 대답의 끝 부분에이 작업을 수행 할 수있는 몇 가지 코드를 넣었습니다.

#include "stdio.h" 
#include <nmmintrin.h>  

void m4x4(const float *A, const float *B, float *C) { 
    for(int i=0; i<4; i++) { 
     for(int j=0; j<4; j++) { 
      float sum = 0.0f; 
      for(int k=0; k<4; k++) { 
       sum += A[i*4+k]*B[k*4+j]; 
      } 
      C[i*4 + j] = sum; 
     } 
    } 
} 

void m4x4T(const float *A, const float *B, float *C) { 
    for(int i=0; i<4; i++) { 
     for(int j=0; j<4; j++) { 
      float sum = 0.0f; 
      for(int k=0; k<4; k++) { 
       sum += A[i*4+k]*B[j*4+k]; 
      } 
      C[i*4 + j] = sum; 
     } 
    } 
} 

void m4x4_vec(const float *A, const float *B, float *C) { 
    __m128 Brow[4], Mrow[4]; 
    for(int i=0; i<4; i++) { 
     Brow[i] = _mm_load_ps(&B[4*i]); 
    } 

    for(int i=0; i<4; i++) { 
     Mrow[i] = _mm_set1_ps(0.0f); 
     for(int j=0; j<4; j++) { 
      __m128 a = _mm_set1_ps(A[4*i +j]); 
      Mrow[i] = _mm_add_ps(Mrow[i], _mm_mul_ps(a, Brow[j])); 
     } 
    } 
    for(int i=0; i<4; i++) { 
     _mm_store_ps(&C[4*i], Mrow[i]); 
    } 
} 

void m4x4T_vec(const float *A, const float *B, float *C) { 
    __m128 Arow[4], Brow[4], Mrow[4]; 
    for(int i=0; i<4; i++) { 
     Arow[i] = _mm_load_ps(&A[4*i]); 
     Brow[i] = _mm_load_ps(&B[4*i]); 
    } 

    for(int i=0; i<4; i++) { 
     __m128 prod[4]; 
     for(int j=0; j<4; j++) { 
      prod[j] = _mm_mul_ps(Arow[i], Brow[j]); 
     } 
     Mrow[i] = _mm_hadd_ps(_mm_hadd_ps(prod[0], prod[1]), _mm_hadd_ps(prod[2], prod[3]));  
    } 
    for(int i=0; i<4; i++) { 
     _mm_store_ps(&C[4*i], Mrow[i]); 
    } 

} 

float compare_4x4(const float* A, const float*B) { 
    float diff = 0.0f; 
    for(int i=0; i<4; i++) { 
     for(int j=0; j<4; j++) { 
      diff += A[i*4 +j] - B[i*4+j]; 
      printf("A %f, B %f\n", A[i*4 +j], B[i*4 +j]); 
     } 
    } 
    return diff;  
} 

int main() { 
    float *A = (float*)_mm_malloc(sizeof(float)*16,16); 
    float *B = (float*)_mm_malloc(sizeof(float)*16,16); 
    float *C1 = (float*)_mm_malloc(sizeof(float)*16,16); 
    float *C2 = (float*)_mm_malloc(sizeof(float)*16,16); 

    for(int i=0; i<4; i++) { 
     for(int j=0; j<4; j++) { 
      A[i*4 +j] = i*4+j; 
      B[i*4 +j] = i*4+j; 
      C1[i*4 +j] = 0.0f; 
      C2[i*4 +j] = 0.0f; 
     } 
    } 
    m4x4T(A, B, C1); 
    m4x4T_vec(A, B, C2); 
    printf("compare %f\n", compare_4x4(C1,C2)); 

} 

편집 : 여기

는 C = (A B)^T를 수행하고 SSE 스칼라 함수이다. 그것들은 A B 버전만큼이나 빠르다.

void m4x4TT(const float *A, const float *B, float *C) { 
    for(int i=0; i<4; i++) { 
     for(int j=0; j<4; j++) { 
      float sum = 0.0f; 
      for(int k=0; k<4; k++) { 
       sum += A[j*4+k]*B[k*4+i]; 
      } 
      C[i*4 + j] = sum; 
     } 
    } 
} 

void m4x4TT_vec(const float *A, const float *B, float *C) { 
    __m128 Arow[4], Crow[4]; 
    for(int i=0; i<4; i++) { 
     Arow[i] = _mm_load_ps(&A[4*i]); 
    } 

    for(int i=0; i<4; i++) { 
     Crow[i] = _mm_set1_ps(0.0f); 
     for(int j=0; j<4; j++) { 
      __m128 a = _mm_set1_ps(B[4*i +j]); 
      Crow[i] = _mm_add_ps(Crow[i], _mm_mul_ps(a, Arow[j])); 
     } 
    } 

    for(int i=0; i<4; i++) { 
     _mm_store_ps(&C[4*i], Crow[i]); 
    } 
} 

Answer 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 또는 비슷한 것을 사용하면 병목 현상이 어디 있는지 알려줍니다.