MKL 是否针对 *major 订单优化 cblas？

Question

我正在使用mklcblas_dgemm，目前与CblasRowMajor、CblasNoTrans、CblasNotrans一起用于我的矩阵。

我知道c 是一种以行为主的语言，而dgemm 是一种以列为主的算法。如果我链接到mkl，我很想知道切换矩阵的顺序是否会对cblas_dgemm 算法产生任何影响。 mkl 是否足够聪明，可以在幕后做一些我试图优化矩阵乘法的事情？如果不是，使用mkl 执行矩阵乘法的最佳方法是什么？

Answer 1

TL;DR：简而言之无论您使用 row-major 还是 column-major 排序来执行矩阵-矩阵乘法都无关紧要使用 MKL（和其他 BLAS 实现）。

---

我知道 c 是一种以行为主的语言，而 dgemm 是一种以列为主的算法。

DGEMM不是列主要算法，它是用于计算具有一般矩阵的矩阵-矩阵乘积的 BLAS 接口。 DGEMM（以及大部分 BLAS）的通用参考实现是Netlib's，它是用 Fortran 编写的。它假设列优先顺序的唯一原因是因为 Fortran 是一种列优先顺序语言。 DGEMM（和相应的 BLAS 3 级函数）不是专门用于列主要数据。

DGEMM 计算什么？

基础数学中的 DGEMM 执行 2D matrix-matrix multiplication。用于乘以 2D 矩阵的标准 算法要求您沿其行遍历一个矩阵，并沿其列遍历另一个矩阵。为了执行矩阵-矩阵乘法，AB = C，我们将 A 的行乘以 B 的列strong> 生成 C。因此，输入矩阵的顺序无关紧要，因为一个矩阵必须沿其行遍历，而另一个矩阵必须沿其列遍历。

使用 MKL 研究行优先和列优先 DGEMM 计算

英特尔 MKL 足够聪明，可以在后台利用这一点，并为 row-major 和 column-major 数据提供完全相同的性能。

cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, ...);

和

cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, ...);

将以相似的性能执行。我们可以用一个相对简单的程序来测试一下

#include <float.h>
#include <mkl.h>
#include <omp.h>
#include <stdio.h>

void init_matrix(double *A, int n, int m, double d);
void test_dgemm(CBLAS_LAYOUT Layout, double *A, double *B, double *C, const MKL_INT m, const MKL_INT n, const MKL_INT k, int nSamples, double *timing);
void print_summary(const MKL_INT m, const MKL_INT n, const MKL_INT k, const int nSamples, const double *timing);

int main(int argc, char **argv) {
    MKL_INT n, k, m;
    double *a, *b, *c;
    double *timing;
    int nSamples = 1;

    if (argc != 5){
        fprintf(stderr, "Error: Wrong number of arguments!\n");
        fprintf(stderr, "usage: %s mMatrix nMatrix kMatrix NSamples\n", argv[0]);
        return -1;
    }

    m = atoi(argv[1]);
    n = atoi(argv[2]);
    k = atoi(argv[3]);

    nSamples = atoi(argv[4]);

    timing = malloc(nSamples * sizeof *timing);

    a = mkl_malloc(m*k * sizeof *a, 64);
    b = mkl_malloc(k*n * sizeof *a, 64);
    c = mkl_calloc(m*n, sizeof *a, 64);

    /** ROW-MAJOR ORDERING **/
    test_dgemm(CblasRowMajor, a, b, c, m, n, k, nSamples, timing);

    /** COLUMN-MAJOR ORDERING **/
    test_dgemm(CblasColMajor, a, b, c, m, n, k, nSamples, timing);

    mkl_free(a);
    mkl_free(b);
    mkl_free(c);
    free(timing);
}

void init_matrix(double *A, int n, int m, double d) {
    int i, j;
    #pragma omp for schedule (static) private(i,j)
    for (i = 0; i < n; ++i) {
        for (j = 0; j < m; ++j) {
            A[j + i*n] = d * (double) ((i - j) / n);
        }
    }
}

void test_dgemm(CBLAS_LAYOUT Layout, double *A, double *B, double *C, const MKL_INT m, const MKL_INT n, const MKL_INT k, int nSamples, double *timing) {
    int i;
    MKL_INT lda = m, ldb = k, ldc = m;
    double alpha = 1.0, beta = 0.0;

    if (CblasRowMajor == Layout) {
        printf("\n*****ROW-MAJOR ORDERING*****\n\n");
    } else if (CblasColMajor == Layout) {
        printf("\n*****COLUMN-MAJOR ORDERING*****\n\n");
    }

    init_matrix(A, m, k, 0.5);
    init_matrix(B, k, n, 0.75);
    init_matrix(C, m, n, 0);

    // First call performs any buffer/thread initialisation
    cblas_dgemm(Layout, CblasNoTrans, CblasNoTrans, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc);

    double tmin = DBL_MAX, tmax = 0.0;
    for (i = 0; i < nSamples; ++i) {
        init_matrix(A, m, k, 0.5);
        init_matrix(B, k, n, 0.75);
        init_matrix(C, m, n, 0);

        timing[i] = dsecnd();
        cblas_dgemm(Layout, CblasNoTrans, CblasNoTrans, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc);
        timing[i] = dsecnd() - timing[i];

        if (tmin > timing[i]) tmin = timing[i];
        else if (tmax < timing[i]) tmax = timing[i];
    }

    print_summary(m, n, k, nSamples, timing);
}

void print_summary(const MKL_INT m, const MKL_INT n, const MKL_INT k, const int nSamples, const double *timing) {
    int i;

    double tavg = 0.0;
    for(i = 0; i < nSamples; i++) {
        tavg += timing[i];
    }
    tavg /= nSamples;

    printf("#Loop | Sizes  m   n   k  | Time (s)\n");
    for(i = 0; i < nSamples; i++) {
        printf("%4d %12d %3d %3d  %6.4f\n", i + 1 , m, n, k, timing[i]);
    }

    printf("Summary:\n");
    printf("Sizes  m   n   k  | Avg. Time (s)\n");
    printf(" %8d %3d %3d %12.8f\n", m, n, k, tavg);
}

在我的系统上会产生

$ ./benchmark_dgemm 1000 1000 1000 5
*****ROW-MAJOR ORDERING*****

#Loop | Sizes  m   n   k  | Time (s)
   1         1000 1000 1000  0.0589
   2         1000 1000 1000  0.0596
   3         1000 1000 1000  0.0603
   4         1000 1000 1000  0.0626
   5         1000 1000 1000  0.0584
Summary:
Sizes  m   n   k  | Avg. Time (s)
     1000 1000 1000   0.05995692

*****COLUMN-MAJOR ORDERING*****

#Loop | Sizes  m   n   k  | Time (s)
   1         1000 1000 1000  0.0597
   2         1000 1000 1000  0.0610
   3         1000 1000 1000  0.0581
   4         1000 1000 1000  0.0594
   5         1000 1000 1000  0.0596
Summary:
Sizes  m   n   k  | Avg. Time (s)
     1000 1000 1000   0.05955171

我们可以看到 column-major 排序时间和 row-major 排序时间之间几乎没有区别。 column-major 0.0595 秒，row-major 0.0599 秒。再次执行此操作可能会产生以下结果，其中行优先计算快 0.00003 秒。

$ ./benchmark_dgemm 1000 1000 1000 5
*****ROW-MAJOR ORDERING*****

#Loop | Sizes  m   n   k  | Time (s)
   1         1000 1000 1000  0.0674
   2         1000 1000 1000  0.0598
   3         1000 1000 1000  0.0595
   4         1000 1000 1000  0.0587
   5         1000 1000 1000  0.0584
Summary:
Sizes  m   n   k  | Avg. Time (s)
     1000 1000 1000   0.06075310

*****COLUMN-MAJOR ORDERING*****

#Loop | Sizes  m   n   k  | Time (s)
   1         1000 1000 1000  0.0634
   2         1000 1000 1000  0.0596
   3         1000 1000 1000  0.0582
   4         1000 1000 1000  0.0582
   5         1000 1000 1000  0.0645
Summary:
Sizes  m   n   k  | Avg. Time (s)
     1000 1000 1000   0.06078266