使用cblas的點積很慢

Question

我想計算產品A^T * A（A是2000x1000矩陣）。另外我只想解決上三角矩陣。在內部循環中，我必須解決兩個向量的點積。

現在，這是問題所在。使用cblas ddot（）並不比使用循環計算點積快。這怎麼可能？ （使用Intel Core TM i7 CPU M620 @ 2,67GHz，1,92GB RAM）

Answer 1

CBLAS點積實際上僅僅是一個稍微展開的循環中的計算。 netlib Fortran就是這樣：

 DO I = MP1,N,5 
     DTEMP = DTEMP + DX(I)*DY(I) + DX(I+1)*DY(I+1) + 
$   DX(I+2)*DY(I+2) + DX(I+3)*DY(I+3) + DX(I+4)*DY(I+4) 
    END DO 

ie。剛剛展開至5

一個箭步如果你必須使用一個ddot風格的點積爲您的操作循環，你可能會通過重新編寫循環性能提升使用SSE2內部函數：

#include <emmintrin.h> 

double ddotsse2(const double *x, const double *y, const int n) 
{ 
    double result[2]; 
    int n2 = 2 * (n/2); 
    __m128d dtemp; 

    if ((n % 2) == 0) { 
     dtemp = _mm_setzero_pd(); 
    } else { 
     dtemp = _mm_set_sd(x[n] * y[n]); 
    } 

    for(int i=0; i<n2; i+=2) { 
     __m128d x1 = _mm_loadr_pd(x+i); 
     __m128d y1 = _mm_loadr_pd(y+i); 
     __m128d xy = _mm_mul_pd(x1, y1); 
     dtemp = _mm_add_pd(dtemp, xy); 
    } 

    _mm_store_pd(&result[0],dtemp); 

    return result[0] + result[1]; 
} 

（未經測試，從未編譯過，請注意購買者）。

這可能比也可能比標準的BLAS實現更快。您可能還想調查是否進一步循環展開可以提高性能。

Answer 2

該問題主要由矩陣大小造成，而不是由ddot引起。您的矩陣非常大，以至於它們不適合緩存內存。解決方案是重新排列三個嵌套循環，以便儘可能使用高速緩存中的一行來完成，以減少高速緩存刷新。對於ddot和daxpy方法都遵循模型實現。在我的電腦上，耗時約爲15：1。換句話說：永遠不會，永遠也不會從我們在學校學到的「行時間列」方案中編程矩陣乘法。

/* 
     Matrix product of A^T * A by two methods. 
     1) "Row times column" as we learned in school. 
     2) With rearranged loops such that need for cash refreshes is reduced 
     (this can be improved even more). 

     Compile: gcc -o aT_a aT_a.c -lgslcblas -lblas -lm 
    */ 

    #include <stdio.h> 
    #include <stdlib.h> 
    #include <time.h> 
    #include <cblas.h> 

    #define ROWS 2000 
    #define COLS 1000 

    static double a[ROWS][COLS]; 
    static double c[COLS][COLS]; 

    static void dot() { 
     int i, j; 
     double *ai, *bj; 
     ai = a[0]; 
     for (i=0; i<COLS; i++) { 
     bj = a[0]; 
     for (j=0; j<COLS; j++) { 
      c[i][j] = cblas_ddot(ROWS,ai,COLS,bj,COLS); 
      bj += 1; 
     } 
     ai += 1; 
     } 
    } 

    static void axpy() { 
     int i, j; 
     double *ci, *bj, aij; 
     for (i=0; i<COLS; i++) { 
     ci = c[i]; 
     for (j=0; j<COLS; j++) ci[j] = 0.; 
     for (j=0; j<ROWS; j++) { 
      aij = a[j][i]; 
      bj = a[j]; 
      cblas_daxpy(COLS,aij,bj,1,ci,1); 
     } 
     } 
    } 

    int main(int argc, char** argv) { 
     clock_t t0, t1; 
     int i, j; 

     for (i=0; i<ROWS; ++i) 
     for (j=0; j<COLS; ++j) 
      a[i][j] = i+j; 

     t0 = clock(); 
     dot(); 
     t0 = clock(); 
     printf("Time for DOT : %f sec.\n",(double)t0/CLOCKS_PER_SEC); 
     axpy(); 
     t1 = clock(); 
     printf("Time for AXPY: %f sec.\n",(double)(t1-t0)/CLOCKS_PER_SEC); 

     return 0; 
    }