cblasを使ったドットプロダクトは遅い

Question

私はA^T * A（Aは2000x1000行列です）を計算したいと思います。また、私は上三角行列を解きたいだけです。内側のループでは、私は2つのベクトルの内積を解かなければならない。

ここで問題があります。 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です。あなたの操作のためにddotスタイルの内積を使用する必要がある場合、あなたはSSE2の組み込み関数を使用するために再書き込みあなたのループによりパフォーマンスが向上する可能性があります。5.

のストライドにアンロールだけループ：

#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ではなく、マトリックスサイズによって発生します。あなたの行列は非常に大きいので、キャッシュメモリに収まらない。解決方法は、3つのネストされたループを並べ替えることで、可能な限りキャッシュ内の行で処理できるようにして、キャッシュのリフレッシュを減らすことです。モデルの実装は、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; 
    }