mal*_*ang 2 c++ blas ubuntu-14.04
这是我之前的问题的延伸,但我是分开问的,因为我真的很沮丧,所以请不要低估它!
问题:与密集矩阵的相同cblas_sgemm调用相比,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;
}
我在使用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
编辑
如果我用以下gemv调用替换gemm调用,那么矩阵稀疏性无关紧要
//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来执行此向量矩阵乘法.因此,我必须对其进行测试或明确证明此调用不正确
编辑
我今天甚至通过使用更新了我的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)
问题:与密集矩阵的相同cblas_sgemm调用相比,cblas_sgemm调用背后的原因可能是使用大量零的矩阵花费更少的时间？
似乎由Ubuntu 14.04(以及可能的其他Ubuntu发行版)的默认libblas-dev包提供的BLAS实现包括对某些矩阵元素为零的情况的优化.
对于Ubuntu 14.04,可以从此处下载BLAS(和cblas)实现/包的源代码.
解压缩该存档后,我们有一个cblas/src包含cblasAPI 的目录,我们有另一个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已成为Bfortran通话/代码,并B已成为A.
当我们检查相应的fortran例程时src/sgemm.f,我们看到以下代码序列:
*
* 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行主试验情况).矩阵行/列乘法运算在我添加音符的循环开始处理***OPTIMIZATION.特别地,如果矩阵元素B(L,J)为零,则跳过在线70处关闭的DO循环.但请记住,B这里对应于A传递给cblas_sgemm例程的矩阵.
跳过这个do-loop允许以这种方式实现的sgemm函数对于在指定row-major时A传递给矩阵中存在大量零的情况要快得多cblas_sgemm.
实验上,并非所有blas实现都有这种优化.在完全相同的平台上进行测试但是使用libopenblas-dev而不是libblas-dev提供这样的加速,即当A矩阵大多为零时基本上没有执行时间差异,而不是基本上没有.
请注意,这里的FORTRAN(77)的代码似乎是相似或相同于旧版本的发布sgemm.f常规,如这里.我发现的这个fortran例程的较新发布版本不包含此优化,例如此处.