c++ - cblas gemm 时间取决于输入矩阵值 - Ubuntu 14.04

Question

这是我的 earlier question 的扩展 , 但我是分开问的，因为我真的很沮丧，所以请不要投反对票!

问题:与密集矩阵的相同 cblas_sgemm 调用相比，对于具有大量零的矩阵，cblas_sgemm 调用花费更少时间的原因可能是什么？

我知道 gemv 是为矩阵- vector 乘法设计的，但为什么我不能使用 gemm 进行 vector -矩阵乘法，如果它花费更少的时间，尤其是对于稀疏矩阵

下面给出了一个简短的代表代码。它要求输入一个值，然后用该值填充一个 vector 。然后它用它的索引替换每第 32 个值。因此，如果我们输入“0”，则会得到一个稀疏 vector ，但对于其他值，我们会得到一个密集 vector 。

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

using namespace std;

int main()
{
const int m = 5000; 

timespec blas_start, blas_end;
long totalnsec; //total nano sec
double totalsec, totaltime;
int i, j;
float *A = new float[m]; // 1 x m
float *B = new float[m*m]; // m x m
float *C = new float[m]; // 1 x m

float input;
cout << "Enter a value to populate the vector (0 for sparse) ";
cin >> input; // enter 0 for sparse

// input martix A: every 32nd element is non-zero, rest of the values = input
for(i = 0; i < m; i++)
{       
A[i] = input;
if( i % 32 == 0)    //adjust for sparsity
        A[i] = i;
}

// input matrix B: identity matrix
for(i = 0; i < m; i++)
        for(j = 0; j < m; j++)
            B[i*m + j] = (i==j);

clock_gettime(CLOCK_REALTIME, &blas_start);
cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, 1, m, m, 1.0f, A, m, B, m, 0.0f, C, m);
clock_gettime(CLOCK_REALTIME, &blas_end);

/* for(i = 0; i < m; i++)
        printf("%f ", C[i]);
printf("\n\n");    */

// Print time
totalsec = (double)blas_end.tv_sec - (double)blas_start.tv_sec;
totalnsec = blas_end.tv_nsec - blas_start.tv_nsec;
if(totalnsec < 0)
{
    totalnsec += 1e9;
    totalsec -= 1;
}
totaltime = totalsec + (double)totalnsec*1e-9;
cout<<"Duration = "<< totaltime << "\n";

return 0;
}

我在带有 blas 3.0 的 Ubuntu 14.04 中运行如下

erisp@ubuntu:~/uas/stackoverflow$ g++ gemmcomp.cpp -o gemmcomp.o -lblas
erisp@ubuntu:~/uas/stackoverflow$ ./gemmcomp.o
Enter a value to populate the vector (0 for sparse) 5
Duration = 0.0291558
erisp@ubuntu:~/uas/stackoverflow$ ./gemmcomp.o
Enter a value to populate the vector (0 for sparse) 0
Duration = 0.000959521

编辑

如果我将 gemm 调用替换为以下 gemv 调用，则矩阵稀疏性无关紧要

//cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, 1, m, m, 1.0f, A, m, B, m, 0.0f, C, m);
cblas_sgemv(CblasRowMajor, CblasNoTrans, m, m, 1.0f, B, m, A, 1, 0.0f, C, 1);

结果

erisp@ubuntu:~/uas/stackoverflow$ g++ gemmcomp.cpp -o gemmcomp.o -lblas
erisp@ubuntu:~/uas/stackoverflow$ ./gemmcomp.o
Enter a value to populate the vector (0 for sparse) 5
Duration = 0.0301581
erisp@ubuntu:~/uas/stackoverflow$ ./gemmcomp.o
Enter a value to populate the vector (0 for sparse) 0
Duration = 0.0299282

但问题是，我正在尝试使用 cublas 优化其他人的代码，而他使用 gemm 成功且高效执行此 vector 矩阵乘法。因此，我必须对其进行测试或明确证明此调用是不正确

编辑

今天我什至更新了我的 blas 库

sudo apt-get install libblas-dev liblapack-dev

编辑:按照不同贡献者的建议执行了以下命令

erisp@ubuntu:~/uas/stackoverflow$ ll -d /usr/lib/libblas* /etc/alternatives/libblas.*
lrwxrwxrwx 1 root root   26 مارچ  13  2015 /etc/alternatives/libblas.a -> /usr/lib/libblas/libblas.a
lrwxrwxrwx 1 root root   27 مارچ  13  2015 /etc/alternatives/libblas.so -> /usr/lib/libblas/libblas.so
lrwxrwxrwx 1 root root   29 مارچ  13  2015 /etc/alternatives/libblas.so.3 -> /usr/lib/libblas/libblas.so.3
lrwxrwxrwx 1 root root   29 مارچ  13  2015 /etc/alternatives/libblas.so.3gf -> /usr/lib/libblas/libblas.so.3
drwxr-xr-x 2 root root 4096 مارچ  13  2015 /usr/lib/libblas/
lrwxrwxrwx 1 root root   27 مارچ  13  2015 /usr/lib/libblas.a -> /etc/alternatives/libblas.a
lrwxrwxrwx 1 root root   28 مارچ  13  2015 /usr/lib/libblas.so -> /etc/alternatives/libblas.so
lrwxrwxrwx 1 root root   30 مارچ  13  2015 /usr/lib/libblas.so.3 -> /etc/alternatives/libblas.so.3
lrwxrwxrwx 1 root root   32 مارچ  13  2015 /usr/lib/libblas.so.3gf -> /etc/alternatives/libblas.so.3gf

erisp@ubuntu:~/uas/stackoverflow$ ldd ./gemmcomp.o
    linux-gate.so.1 =>  (0xb76f6000)
    libblas.so.3 => /usr/lib/libblas.so.3 (0xb765e000)
    libstdc++.so.6 => /usr/lib/i386-linux-gnu/libstdc++.so.6 (0xb7576000)
     libc.so.6 => /lib/i386-linux-gnu/libc.so.6 (0xb73c7000)
     libm.so.6 => /lib/i386-linux-gnu/libm.so.6 (0xb7381000)
/lib/ld-linux.so.2 (0xb76f7000)
     libgcc_s.so.1 => /lib/i386-linux-gnu/libgcc_s.so.1 (0xb7364000)

Answer 1

Question: What could be the reason behind a cblas_sgemm call taking much less time for matrices with a large number of zeros as compared to the same cblas_sgemm call for dense matrices?

貌似the default libblas-dev package提供的BLAS实现Ubuntu 14.04(可能还有其他 Ubuntu 发行版)包含针对某些矩阵元素为零的情况的优化。

对于 Ubuntu 14.04，BLAS(和 cblas)实现/包的源代码可以从 here 下载。 .

解压该存档后，我们有一个包含 cblas API 的 cblas/src 目录，我们还有另一个包含 src 的目录各种 blas 例程的 F77 实现。

在 cblas_sgemm 的情况下，当指定参数 CblasRowMajor 时，cblas/src/cblas_sgemm.c 代码将调用底层的 Fortran 例程，如下所示:

void cblas_sgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA,
             const enum CBLAS_TRANSPOSE TransB, const int M, const int N,
             const int K, const float alpha, const float  *A,
             const int lda, const float  *B, const int ldb,
             const float beta, float  *C, const int ldc)
{
...
} else if (Order == CblasRowMajor)
...
  F77_sgemm(F77_TA, F77_TB, &F77_N, &F77_M, &F77_K, &alpha, B, &F77_ldb, A, &F77_lda, &beta, C, &F77_ldc);

请注意，对于此行主要调用，A 和B 矩阵的顺序在传递给 F77_sgemm 例程时是相反的。这是明智的，但我不会在这里深入探讨原因。只需注意 A 在 fortran 调用/代码中已变为 B，而 B 已变为 A。

当我们检查 src/sgemm.f 中相应的 fortran 例程时，我们看到以下代码序列:

*
*     Start the operations.
*
      IF (NOTB) THEN
          IF (NOTA) THEN
*
*           Form  C := alpha*A*B + beta*C.
*
              DO 90 J = 1,N
                  IF (BETA.EQ.ZERO) THEN
                      DO 50 I = 1,M
                          C(I,J) = ZERO
   50                 CONTINUE
                  ELSE IF (BETA.NE.ONE) THEN
                      DO 60 I = 1,M
                          C(I,J) = BETA*C(I,J)
   60                 CONTINUE
                  END IF
                  DO 80 L = 1,K
                      IF (B(L,J).NE.ZERO) THEN        ***OPTIMIZATION
                          TEMP = ALPHA*B(L,J)
                          DO 70 I = 1,M
                              C(I,J) = C(I,J) + TEMP*A(I,L)
   70                     CONTINUE
                      END IF
   80             CONTINUE
   90         CONTINUE

以上是处理没有转置 A 和没有转置 B 的情况的代码部分(对于这个 cblas row-major 是正确的测试用例)。矩阵行/列乘法运算在我添加注释 ***OPTIMIZATION 的循环开始处处理。特别是，如果矩阵元素 B(L,J) 为零，则跳过第 70 行结束的 DO 循环。但请记住，此处的 B 对应于传递给 cblas_sgemm 例程的 A 矩阵。

跳过此 do 循环允许以这种方式实现的 sgemm 函数在传递给 cblas_sgemm< 的 A 矩阵中有大量零的情况下要快得多 当指定行优先时。

根据实验，并不是所有的 blas 实现都有这个优化。在完全相同的平台上进行测试，但使用 libopenblas-dev 而不是 libblas-dev 没有提供这样的加速，即当 A 矩阵大部分为零，而不是零。

请注意，此处的 fortran (77) 代码似乎与 sgemm.f 例程的较早发布版本相似或相同，例如 here .我能找到的此 Fortran 例程的较新发布版本不包含此优化，例如 here .