blas - 使用 cblas 的点积很慢

Question

我想计算乘积 A^T*A （ A 是 2000x1000 矩阵）。另外我只想解决上三角矩阵。在内部循环中，我必须解决两个向量的点积。

现在，问题来了。使用 cblas ddot() 并不比使用循环计算点积快。这怎么可能？（使用 Intel Core (TM)i7 CPU M620 @2,67GHz，1,92GB RAM）

Answer 1

该问题本质上是由矩阵大小引起的，而不是由 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;
    }

Answer 2

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 3

如果您不使用 SSE2 内在函数或使用可能无法提高性能的数据类型，您可以尝试转置矩阵，以便轻松提高使用cblas_?dot. 分块执行矩阵乘法也有帮助。

void matMulDotProduct(int n, float *A, float* B, int a_size, int b_size, int a_row, int a_col, int b_row, int b_col, float *C) {
    int i, j, k;
    MKL_INT incx, incy;

    incx = 1;
    incy = b_size;

    //copy out multiplying matrix from larger matrix
    float *temp = (float*) malloc(n * n * sizeof(float));
    for (i = 0; i < n; ++i) {
        cblas_scopy(n, &B[(b_row * b_size) + b_col + i], incy, &temp[i * n], 1);
    }

    //transpose
    mkl_simatcopy('R', 'T', n, n, 1.0, temp, 1, 1);

    for (i = 0; i < n; i+= BLOCK_SIZE) {
        for (j = 0; j < n; j++) {
            for (k = 0; k < BLOCK_SIZE; ++k) {
                C[((i + k) * n) + j] = cblas_sdot(n, &A[(a_row + i + k) * a_size + a_col], incx, &temp[n * j], 1);
            }
        }
    }
    free(temp);
}

在我的机器上，对于单精度浮点数和 2K x 2K 矩阵，此代码比 3 循环代码快大约 1 个数量级（但也比 cblas_?gemm 调用慢 1 个数量级）。（我正在使用英特尔 MKL）。