2012-10-03 74 views
0

參考http://blogs.msdn.com/b/xiangfan/archive/2009/04/28/optimize-your-code-matrix-multiplication.aspx瞭解SSE3矩陣乘法優化

template<> 
void SeqMatrixMult4(int size, float** m1, float** m2, float** result) 
{ 
    Transpose(size, m2); 
    for (int i = 0; i < size; i++) { 
     for (int j = 0; j < size; j++) { 
      __m128 c = _mm_setzero_ps(); 

      for (int k = 0; k < size; k += 4) { 
       c = _mm_add_ps(c, _mm_mul_ps(_mm_load_ps(&m1[i][k]), _mm_load_ps(&m2[j][k]))); 
      } 
      c = _mm_hadd_ps(c, c); 
      c = _mm_hadd_ps(c, c); 
      _mm_store_ss(&result[i][j], c); 
     } 
    } 
    Transpose(size, m2); 
} 

爲什麼還有2個_mm_hadd_ps(c, c)之後的最內圈?爲了驗證我的理解:此代碼從m1加載4個浮點數,並從m2加載另外4個浮點數,然後將它們相乘得到4個浮點數(__m128)。然後我將它們總結爲c(此時,它仍然是4個花車?)。然後在for循環I hadd這個結果兩次?這是做什麼的?


我的代碼稍微改寫產生錯誤的結果似乎

long long start, end; 
__m128 v1, v2, vMul, vRes; 
vRes = _mm_setzero_ps(); 

start = wall_clock_time(); 
transpose_matrix(m2); 
for (int i = 0; i < SIZE; i++) { 
    for (int j = 0; j < SIZE; j++) { 
     float tmp = 0; 
     for (int k = 0; k < SIZE; k+=4) { 
      v1 = _mm_load_ps(&m1[i][k]); 
      v2 = _mm_load_ps(&m2[j][k]); 
      vMul = _mm_mul_ps(v1, v2); 

      vRes = _mm_add_ps(vRes, vMul); 
     } 
     vRes = _mm_hadd_ps(vRes, vRes); 
     _mm_store_ss(&result[i][j], vRes); 
    } 
} 
end = wall_clock_time(); 
fprintf(stderr, "Optimized Matrix multiplication took %1.2f seconds\n", ((float)(end - start))/1000000000); 

// reverse the transposition 
transpose_matrix(m2); 
+0

這似乎工作在4x4矩陣?我需要它爲1024x1024到4K x 4k矩陣 –

+0

剛剛使用MKL/Goto BLAS/ATLAS。對於大型矩陣乘法,您需要進行多級緩存阻塞,但如果您以前從未做過,則需要幾周才能實施和調整。 SIMD只會給你緩存阻塞的一小部分性能提升。 –

回答

4

haddps不總結的所有四個元素在載體中。兩個haddps指令需要獲得完整的水平和。

如果我們編號向量{c0,c1,c2,c3}的元素,第一個haddps產生{c0+c1, c2+c3, c0+c1, c2+c3}。第二個產生{c0+c1+c2+c3, <same thing in the other lanes>}

+0

夠公平的,但是'_mm_store_ss'只會存儲矢量的第一個元素?當我閱讀MSDN上的文檔時沒有清楚。 –