Ubuntu 14.04 中矩阵稀疏性对 cblas sgemm 的影响

Question

我最近发现，如果矩阵中有 "large" 个零，则 cblas_sgemm 调用矩阵乘法的性能会显着提高。它改进到它击败它的 cublas 表亲大约 100 倍的程度。这很可能归因于某些 稀疏性自动检测 和 cblas_sgemm 函数进行的适当格式转换。

不幸的是，它的 cuda 对应物即 cublasSgemm 没有表现出这种行为。

所以，问题是，对于 可能有 大量零的矩阵，我如何才能在 cublasSgemm 上获得相同类型的优化。

cblas_sgemm使用什么技术自动调整到稀疏矩阵？

请不要推荐 cuSparse / CUSP 等，因为

我事先不确定输入矩阵的稀疏性
我正在研究一种迭代算法，在最初的几次迭代中，矩阵 可能是 稀疏的，但随着时间的推移逐渐变得密集。

提前致谢

编辑以包含重现上述场景的代码

#include <iostream>
#include <stdio.h>
#include <time.h>
#include <cblas.h>
#include <cublas_v2.h>
using namespace std;

int main()
{
const int m = 5000; 

timespec blas_start, blas_end, cublas_start, cublas_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

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

// input matrix B: identity matrix
// col major = row major
for(i = 0; i < m; i++)
    for(j = 0; j < m; j++)
    {
        if (i==j)           
            B[j*m + i] = 1;
        else
            B[j*m + i] = 0;     
    }

clock_gettime(CLOCK_REALTIME, &blas_start);
cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, 1, m, m, 1, A, m, B, m, 0, C, m);
clock_gettime(CLOCK_REALTIME, &blas_end);
    /*
for(i = 0; i < 12; i++)
    printf("%f ", C[i]);
    */

//cublas section

cudaError_t cudaStat;   
cublasHandle_t handle;
cublasCreate(&handle);
//Declaring Device Variables
float *A_d, *B_d, *C_d;

//Allocating Memory for Device Variables
cudaStat = cudaMalloc(&A_d, sizeof(float)*m);
if(cudaStat != cudaSuccess) printf("Error Allocating Memory for A_d\n");

cudaStat = cudaMalloc(&B_d, sizeof(float)*m*m);
if(cudaStat != cudaSuccess) printf("Error Allocating Memory for B_d\n");

cudaStat = cudaMalloc(&C_d, sizeof(float)*m);
if(cudaStat != cudaSuccess) printf("Error Allocating Memory for C_d\n");

// Moving values of A, B onto Device variables
cublasSetVector(m, sizeof(float), A, 1, A_d, 1);
cublasSetMatrix(m, m, sizeof(float), B, m, B_d, m); 

// Do the actual multiplication
float alpha = 1.0f, beta = 0.0f;
cudaDeviceSynchronize();
clock_gettime(CLOCK_REALTIME, &cublas_start);

cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, 1, m, m, &alpha, A_d, 1, B_d, m, &beta, C_d, 1);

cudaDeviceSynchronize();
clock_gettime(CLOCK_REALTIME, &cublas_end);

cublasGetVector(m, sizeof(float), C, 1, C_d, 1);
    /*
for(i = 0; i < 12; i++)
    printf("%f ", C[i]);
    */

// Print times
// blas 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<<"BLAS Time = "<< totaltime << "\n";

//cublas
totalsec = (double)cublas_end.tv_sec - (double)cublas_start.tv_sec;
totalnsec = cublas_end.tv_nsec - cublas_start.tv_nsec;
if(totalnsec < 0)
{
    totalnsec += 1e9;
    totalsec -= 1;
}
totaltime = totalsec + (double)totalnsec*1e-9;
cout<<"CUBLAS Time = "<< totaltime << "\n";

return 0;
}

运行它得到如下结果

malang@ubuntu:~/uas/Whosebug$ nvcc -arch=sm_12 blascomp.cu -o blascomp.o -lblas -lcublas
malang@ubuntu:~/uas/Whosebug$ ./blascomp.o
BLAS Time = 0.000964504
CUBLAS Time = 0.0365322

编辑

在@Eric

回答后编辑

使用cublasSgemv 大大提高了GPU 的性能。但是，我仍然有这个问题，即 cblas_sgemm 对于 CPU 上的 稀疏矩阵 效率更高。可能的原因是什么？

EDIT 根据@Eric @osgx @Robert Crovella

的建议执行了以下命令

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)

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

Answer 1

您的代码有问题 - 您使用了错误的 BLAS API。您使用矩阵-矩阵-乘法例程 gemm() 进行向量-矩阵-乘法运算。

对于 vec-mat-mul 或 mat-vec-mul 你应该使用 gemv()。当然 gemm() 可以用只有 1 行的矩阵给出正确的结果。但这是 gemv() 应该处理的意外极端情况，所以你可能无法在 GPU and/or CPU.

上获得最佳性能

您可以更改为 gemv() 并再次进行基准测试。

编辑

这是我使用单线程 MKL 的基准测试结果。 A 和 B 的值与您的代码中的值相同。我无法在 CPU 上重现“0.000964504s”的结果。您可以检查结果向量的正确性。您的 cblas 库有可能存在错误。

使用gemm()

BLAS Time = 0.0169784
CUBLAS Time = 0.00356155

使用gemv()

BLAS Time = 0.0167557
CUBLAS Time = 0.0013809

EDIT2

我现在可以使用包 libblas-dev.

在 unbuntu 14.04 上重现 FAST 结果

原因在下面的问题中得到了回答。

在特定版本的 BLAS 中，有检查零元素的代码。检查成本为 O(n^2)，因此值得在成本为 O(n^3) 的矩阵-矩阵乘法上执行此操作。

对于GPU gemm()，由于计算顺序不同（逐块而不是逐行），这样的优化可能不可行。但这对于 GPU gemv() 是可行的，可以节省从全局内存加载矩阵所花费的时间。

Ubuntu 14.04 中矩阵稀疏性对 cblas sgemm 的影响

Impact of matrix sparsity on cblas sgemm in Ubuntu 14.04

cuda

blas

cublas

ubuntu-14.04