cblas gemm 时间取决于输入矩阵值 - Ubuntu 14.04
cblas gemm time dependent on input matrix values - Ubuntu 14.04
这是我的 的 扩展 ，但我单独提问是因为我真的很沮丧，所以请不要投反对票！
问题：与相同的 [=56] 相比，cblas_sgemm 调用具有 大量零 的矩阵花费更少时间的原因可能是什么=] 需要密集矩阵吗？
我知道 gemv 是为矩阵-向量乘法设计的，但为什么我不能使用 gemm 进行向量-矩阵乘法，如果它花费更少的时间，特别是对于稀疏矩阵
下面给出一个简短的代表代码。它要求输入一个值，然后用该值填充一个向量。然后它用它的索引替换每第 32 个值。因此，如果我们输入“0”，则会得到一个稀疏向量，但对于其他值，我们会得到一个密集向量。
#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;
}
我 运行 它在 Ubuntu 14.04 和 blas 3.0
中如下
erisp@ubuntu:~/uas/Whosebug$ g++ gemmcomp.cpp -o gemmcomp.o -lblas
erisp@ubuntu:~/uas/Whosebug$ ./gemmcomp.o
Enter a value to populate the vector (0 for sparse) 5
Duration = 0.0291558
erisp@ubuntu:~/uas/Whosebug$ ./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/Whosebug$ g++ gemmcomp.cpp -o gemmcomp.o -lblas
erisp@ubuntu:~/uas/Whosebug$ ./gemmcomp.o
Enter a value to populate the vector (0 for sparse) 5
Duration = 0.0301581
erisp@ubuntu:~/uas/Whosebug$ ./gemmcomp.o
Enter a value to populate the vector (0 for sparse) 0
Duration = 0.0299282
但问题是我正在尝试使用 cublas 优化其他人的代码，他成功了，并且 高效地 使用 gemm 执行此向量矩阵乘法。因此，我必须对其进行测试或明确证明此调用是 不正确的
编辑
今天我什至使用
更新了我的 blas 库
sudo apt-get install libblas-dev liblapack-dev
编辑：按照不同贡献者的建议执行了以下命令
erisp@ubuntu:~/uas/Whosebug$ 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/Whosebug$ 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)
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 为 Ubuntu 14.04（可能还有其他 Ubuntu 分布）提供的 BLAS 实现包括对某些矩阵元素为零的情况的优化。
对于 Ubuntu 14.04，可以从 here.
下载 BLAS（和 cblas）implementation/package 的源代码
解压缩该存档后，我们有一个 cblas/src 目录，其中包含 cblas API，我们还有另一个 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 call/code 中变成了 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/column 乘法运算在我添加注释 ***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. Newer published versions of this fortran routine that I could find do not contain this optimization, such as here.
这是我的
问题：与相同的 [=56] 相比，cblas_sgemm 调用具有 大量零 的矩阵花费更少时间的原因可能是什么=] 需要密集矩阵吗？
我知道 gemv 是为矩阵-向量乘法设计的，但为什么我不能使用 gemm 进行向量-矩阵乘法，如果它花费更少的时间，特别是对于稀疏矩阵
下面给出一个简短的代表代码。它要求输入一个值，然后用该值填充一个向量。然后它用它的索引替换每第 32 个值。因此，如果我们输入“0”，则会得到一个稀疏向量，但对于其他值，我们会得到一个密集向量。
#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;
}
我 运行 它在 Ubuntu 14.04 和 blas 3.0
中如下erisp@ubuntu:~/uas/Whosebug$ g++ gemmcomp.cpp -o gemmcomp.o -lblas
erisp@ubuntu:~/uas/Whosebug$ ./gemmcomp.o
Enter a value to populate the vector (0 for sparse) 5
Duration = 0.0291558
erisp@ubuntu:~/uas/Whosebug$ ./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/Whosebug$ g++ gemmcomp.cpp -o gemmcomp.o -lblas
erisp@ubuntu:~/uas/Whosebug$ ./gemmcomp.o
Enter a value to populate the vector (0 for sparse) 5
Duration = 0.0301581
erisp@ubuntu:~/uas/Whosebug$ ./gemmcomp.o
Enter a value to populate the vector (0 for sparse) 0
Duration = 0.0299282
但问题是我正在尝试使用 cublas 优化其他人的代码，他成功了，并且 高效地 使用 gemm 执行此向量矩阵乘法。因此，我必须对其进行测试或明确证明此调用是 不正确的
编辑
今天我什至使用
更新了我的 blas 库sudo apt-get install libblas-dev liblapack-dev
编辑：按照不同贡献者的建议执行了以下命令
erisp@ubuntu:~/uas/Whosebug$ 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/Whosebug$ 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)
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 为 Ubuntu 14.04（可能还有其他 Ubuntu 分布）提供的 BLAS 实现包括对某些矩阵元素为零的情况的优化。
对于 Ubuntu 14.04，可以从 here.
下载 BLAS（和 cblas）implementation/package 的源代码解压缩该存档后，我们有一个 cblas/src 目录，其中包含 cblas API，我们还有另一个 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 call/code 中变成了 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/column 乘法运算在我添加注释 ***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. Newer published versions of this fortran routine that I could find do not contain this optimization, such as here.