cublas矩阵乘法效率测试

[cuda2010]

以前曾听说cublas的效率不是很高, 今天写了个小程序对cublas的矩阵乘法速度进行了一个测试, 发现结果并非如此. 至少就矩阵乘法来说，cublas的效率很不错, 相对CPU有非常高的加速比.
测试程序是在sdk例子simpleCUBLAS的基础上修改而成, 测试内容是分别用cublas和CPU函数计算两个N阶矩阵A和B的乘积, 然后对结果进行校验，并计算各自的flops.(每个矩阵元素加乘看作2flops)
测试所使用的显卡为GTX 295, CPU是i7 920, 均为默认频率.
显然, 矩阵乘法的速度受问题尺度影响很大, 因此具体测试了N=256-4096的情况. 结果在N=256时, cublas的速度就已经达到了CPU函数的300余倍, 而在N=4096时, cublas的速度超过了CPU的2000倍, 此时的实测浮点速度达到了约360Gflops.
当然，这里面CPU矩阵乘法的算法过于简单，并没有做什么优化. 但是即使经过充分优化，CPU的速度也只能提高1个数量级左右，仍然和cublas差距巨大.

测试具体结果:
N= 256, GPU=0.000180s(186.121Gflops), CPU=0.061939s(0.542Gflops)
N= 512, GPU=0.000934s(287.515Gflops), CPU=0.498737s(0.538Gflops)
N=1024, GPU=0.006252s(343.471Gflops), CPU=8.553767s(0.251Gflops)
N=2048, GPU=0.048440s(354.666Gflops), CPU=85.726814s(0.200Gflops)
N=3072, GPU=0.161870s(358.201Gflops), CPU=292.890472s(0.198Gflops)
N=4096, GPU=0.383253s(358.612Gflops), CPU=799.222351s(0.172Gflops)

具体的测试程序代码:

#include "cutil_inline.h"
#include "cublas.h"

#define N 1024

void simple_sgemm(const float *A, const float *B, float *C) {
    int i, j, k;
    for(i=0; i<N; i++)
    for(j=0; j<N; j++) {
        float s=0;
        for(k=0; k<N; k++) s+=A[k*N+i]*B[j*N+k];
        C[j*N+i]=s;
    }
}

int main() {    
    float *h_A=(float*)malloc(N*N*sizeof(float));
    float *h_B=(float*)malloc(N*N*sizeof(float));
    float *h_C=(float*)malloc(N*N*sizeof(float));
    float *h_C_ref=(float*)malloc(N*N*sizeof(float));
    float *d_A, *d_B, *d_C;
    unsigned int timer1=0;
    cutCreateTimer(&timer1);
    cutStartTimer(timer1);
    printf("simpleCUBLAS test running..\n");
    cublasInit();
    for(int i=0; i<N*N; i++) {
        h_A[i]=rand()/(float)RAND_MAX;
        h_B[i]=rand()/(float)RAND_MAX;
    }
    cublasAlloc(N*N, sizeof(float), (void**)&d_A);
    cublasAlloc(N*N, sizeof(float), (void**)&d_B);
    cublasAlloc(N*N, sizeof(float), (void**)&d_C);
    cublasSetVector(N*N, sizeof(float), h_A, 1, d_A, 1);
    cublasSetVector(N*N, sizeof(float), h_B, 1, d_B, 1);
    float t0, gpu_t, cpu_t, error_norm=0, ref_norm=0;
    cudaThreadSynchronize();
    t0=cutGetTimerValue(timer1);
    cublasSgemm('n', 'n', N, N, N, 1.0f, d_A, N, d_B, N, 0.0f, d_C, N);
    cudaThreadSynchronize();
    gpu_t=(cutGetTimerValue(timer1)-t0)/1000.0f;
    cublasGetVector(N*N, sizeof(float), d_C, 1, h_C, 1);
    t0=cutGetTimerValue(timer1);
    simple_sgemm(h_A, h_B, h_C_ref);
    cpu_t=(cutGetTimerValue(timer1)-t0)/1000.0f;
    printf("N=%4d, GPU=%.6fs(%.3fGflops), CPU=%.6fs(%.3fGflops)\n", 
        N, gpu_t, 1e-9*N*N*N*2/gpu_t, cpu_t, 1e-9*N*N*N*2/cpu_t);
    for(int i=0; i<N*N; i++) {
        float diff=h_C_ref[i]-h_C[i];
        error_norm+=diff*diff;
        ref_norm+=h_C_ref[i]*h_C_ref[i];
    }
    printf("Test %s\n", (sqrtf(error_norm/ref_norm)<1E-6) ? "PASSED" : "FAILED");
}

[cuda2010]

用百度搜了下，找到一篇论文也是讨论cublas的矩阵乘法效率的。这篇文章中使用的是8800GT的显卡，矩阵乘法部分测试了两个尺寸为(2560x1536)和(1536x4096)的矩阵相乘，得出的结论是cublas和自己编程的实测性能分别为530GFlops和444GFlops。不过这个数值似乎有些偏大，根据文章中给出的程序用时303ms和435ms,按照我的理解应该是106Gflops和74Gflops，正好和文章给出的数值差5倍和6倍，也许文章是把1次矩阵元素相乘计作10flops和12flops了，我觉得还是计作2flops比较合理。
这篇论文的地址是: www.ecice06.com/qikan/manage/wenzhang/091001.pdf。

[无心人_过过小日子]

intel的sse也是simd指令.呵呵.

[cuda2010]

再补充一个双精度的测试结果，只需把测试程序中的float和cublasSgemm改成double和cublasDgemm即可。双精度cublas矩阵乘法的性能基本上是单精度的1/3-1/5。

N= 256, GPU=0.000603s(55.677Gflops)
N= 512, GPU=0.004517s(59.428Gflops)
N=1024, GPU=0.030556s(70.280Gflops)
N=2048, GPU=0.239756s(71.656Gflops)
N=3072, GPU=0.802976s(72.209Gflops)
N=4096, GPU=1.901021s(72.297Gflops)

近日查阅了一些资料，证实了矩阵乘法效率不佳的只是cublas 1.x，而cublas 2.x已是顶级性能。此外，有报道称经充分优化后目前4核的CPU矩阵乘法也能达到几十GFlops的实测性能，如属实，CUDA的加速比就不是那么强大了。

[cuda2010]

GTX295单个GPU双精度的理论峰值性能是1.242*30*2=74.52Gflops，cublas在4096x4096时的实测值72.297Gflops已经很接近峰值了，绝对效率达到了97%。

单精度的实测值离峰值还有一些距离, 据说主要原因是使用shared mem时算术指令的速度会下降。此外，应该还有一个原因就是目前cublas中的代码是针对g80优化的，并没有考虑g200的dual-issue特性，如果能充分利用这一点也许效率还会有提升。

[无心人_过过小日子]

1)dual-issue特性的话应该是x3因子,而非*2因子.
2)LZ的峰值性能计算中是否少乘了8?呵呵.

[cuda2010]

dual-issue方面似乎官方的说明不多，我只能通过网上找到一些零散的资料。从目前资料看，gt200只有单精度mad+mul可以完美dual-issue，而双精度运算dual-issue，因此只能*2，从官方给出的Tesla双精度浮点峰值性能也可以看出这一点。系数8是因为双精度单元是单精度的1/8，应该没错吧。

[cuda2010]

上文中有个笔误，"而双精度运算dual-issue，因此只能*2"应该是""而双精度运算不能dual-issue，因此只能*2""

另外，刚才写了一个小程序实际测试了一下，单精度mad; mad; mad; mad; ...指令组成的循环实测浮点性能是593GFlops，(理论峰值1.242*240*2=596GFlops). 而单精度mad; mul; mad; mul; mad; mul..组成的循环实测浮点性能是823GFlops (理论峰值1.242*240*3=894GFlops))。也就是说的确很好地dual-issue了。

而双精度的情况则不同，双精度mad; mad; mad; mad; ...组成的循环实测浮点性能是74.16GFlops(理论峰值1.242*30*2=74.52GFlops)；双精度mad; mul; mad; mul; mad; mul..组成的循环实测浮点性能是55.87GFlops，刚好等于74.52的3/4，也就是mad和mul没有dual-issue而是交替执行的，这么算起来理论速度应该是(1.242*30*2+1.242*30*1)/2=55.89，符合得非常完美。这样看起来，双精度的确完全不能dual-issue。

[cuda2010]

关于CUDA sgemm的代码优化问题最近又有新的变化。blas 2.x中的sgemm使用的是伯克利Volkov的代码，在g80时代是最优的(接近理论极限值)。但是到了g200时代有不同了，因为g200支持dual-issue而Volkov的代码未考虑。就在最近CUDA官方论坛上有人发文说通过某些技巧改进Volkov的sgemm代码以激活dual-issue，实现了约10-20%的性能提升，这个应该是目前所见最高效的sgemm了。(其中在GTX295上实现了最高432GFlops的速度，的确提升不少。)

http://forums.nvidia.com/index.php?showtopic=159033

[cuidiwhere]

请问楼主有测试过非方阵的矩阵乘法吗？如果有的话，讨论下？QQ445281810，谢谢噢~

[cuidiwhere]

请教个问题：在SDK中给出的代码，红色部分实现了将结果矩阵转秩的功能，所以两种方式计算出的 矩阵乘法的结果相同。那么你是怎么解决cublas按column-major的存储方式和C语言下row-major存储方式的转换问题的？万分感谢~
void simple_sgemm(const float *A, const float *B, float *C) {
int i, j, k;
for(i=0; i<N; i++)
for(j=0; j<N; j++) {
float s=0;
for(k=0; k<N; k++) s+=A[k*N+i]*B[j*N+k];
C[j*N+i]=s;
}
}

[xin_200]

就矩阵乘法而言，目前CPU和GPU的差距并不大。CPU上的高性能矩阵相乘非常复杂，每个CPU厂商都有相应的数学库。例如INTEL MKL的BLAS。其峰值速度接近于最好的GPU结果。
目前显卡的性能提升明显慢于CPU，再考虑到将矩阵传送到显卡的时间。就矩阵相乘而言，我更看好多核CPU+大容量内存。
Nvidia自己做些例子用来吹牛就算了，包括一些研究人员在性能对比测试上，往往采用一些未经优化的代码，以突出GPU的优势，相当不严谨。

[cuda2010]

的确, 双精度的矩阵乘法目前CPU和GPU确差距不大, 四核xeon利用imkl最高可以达到70GFlops/s左右, 和GTX 295的一个GPU实测性能相当.

不过单精度差距还是比较大的, 目前CUDA的最好实测成绩是510GFlops/GPU左右(GTX 285), 而一块GTX 295的两个GPU加起来可以达到800GFlops以上, 比四核xeon的理论峰值还高一个数量级.

但目前GPU双精度差主要是构架原因, 如果CPU的性能还是保持过去发展速度的话， 等fermi出来后CPU和GPU的双精度差距可能也会大幅拉大.

关于研究人员的对比数据, 的确存在CPU代码不够优化的例子. 但是其中至少有部分并不是研究人员不想优化, 而是他们尽量对CPU和GPU代码进行了优化，最后只能达到这样的水平 (CPU和GPU都并非最优化, 但CPU代码更差)。我想也许是CPU厂商在CPU优化方面宣传力度不够，很多科研人员不了解如何优化吧。(以矩阵乘法为例, 简单使用sse+openmp很难接近imkl的成绩, 有数量级的差距)

[cuda2010]

抱歉，我上面的方法可能还不大好，因为其中的simple_sgemm中实际上还是隐含着进行了一次转置。
实际上行先列先对应于矩阵转置。一个更简单的办法是利用以下矩阵转置公式，在调用cublas时用'n', 'n'参数同时交换A, B的顺序，这样得到的结果就不需要转置了.

C=A*B -> C'=B'*A' (这里'代表转置)

代码如下：

#include "cutil_inline.h"
#include "cublas.h"

#define N 512

void simple_sgemm(const float A[N][N], const float B[N][N], float C[N][N]) {
    int i, j, k;
    for(i=0; i<N; i++)
    for(j=0; j<N; j++) {
        float s=0;
        for(k=0; k<N; k++) s+=A[i][k]*B[k][j];
        C[i][j]=s;
    }
}

float h_A[N][N], h_B[N][N], h_C[N][N], h_C_ref[N][N];

int main() {    
    float *d_A, *d_B, *d_C;
    unsigned int timer1=0;
    cutCreateTimer(&timer1);
    cutStartTimer(timer1);
    printf("simpleCUBLAS test running..\n");
    cublasInit();
    for(int i=0; i<N; i++)
    for(int j=0; j<N; j++) {
        h_A[i][j]=rand()/(float)RAND_MAX;
        h_B[i][j]=rand()/(float)RAND_MAX;
    }
    cublasAlloc(N*N, sizeof(float), (void**)&d_A);
    cublasAlloc(N*N, sizeof(float), (void**)&d_B);
    cublasAlloc(N*N, sizeof(float), (void**)&d_C);
    cublasSetVector(N*N, sizeof(float), h_A, 1, d_A, 1);
    cublasSetVector(N*N, sizeof(float), h_B, 1, d_B, 1);
    float t0, gpu_t, cpu_t, error_norm=0, ref_norm=0;
    cudaThreadSynchronize();
    t0=cutGetTimerValue(timer1);
    cublasSgemm('n', 'n', N, N, N, 1.0f, d_B, N, d_A, N, 0.0f, d_C, N);
    cudaThreadSynchronize();
    gpu_t=(cutGetTimerValue(timer1)-t0)/1000.0f;
    cublasGetVector(N*N, sizeof(float), d_C, 1, h_C, 1);shuo
    t0=cutGetTimerValue(timer1);
    simple_sgemm(h_A, h_B, h_C_ref);
    cpu_t=(cutGetTimerValue(timer1)-t0)/1000.0f;
    printf("N=%4d, GPU=%.6fs(%.3fGflops), CPU=%.6fs(%.3fGflops)\n", 
        N, gpu_t, 1e-9*N*N*N*2/gpu_t, cpu_t, 1e-9*N*N*N*2/cpu_t);
    for(int i=0; i<N; i++)
    for(int j=0; j<N; j++) {
        float diff=h_C_ref[i][j]-h_C[i][j];
        error_norm+=diff*diff;
        ref_norm+=h_C_ref[i][j]*h_C_ref[i][j];
    }
    printf("Test %s\n", (sqrtf(error_norm/ref_norm)<1E-6) ? "PASSED" : "FAILED");
}

[cuda2010]

blas在一开始设计时就考虑到了这个问题，只需修改第1和第2个参数即可满足不同的row/column major存储方式。用'T', 'T'参数应该就可满足你的要求。修改后的代码如下：

#include "cutil_inline.h"
#include "cublas.h"

#define N 512

void simple_sgemm(const float A[N][N], const float B[N][N], float C[N][N]) {
    int i, j, k;
    for(j=0; j<N; j++)
    for(i=0; i<N; i++) {
        float s=0;
        for(k=0; k<N; k++) s+=A[i][k]*B[k][j];
        C[j][i]=s;
    }
}

float h_A[N][N], h_B[N][N], h_C[N][N], h_C_ref[N][N];

int main() {    
    float *d_A, *d_B, *d_C;
    unsigned int timer1=0;
    cutCreateTimer(&timer1);
    cutStartTimer(timer1);
    printf("simpleCUBLAS test running..\n");
    cublasInit();
	for(int j=0; j<N; j++)
    for(int i=0; i<N; i++) {
        h_A[j][i]=rand()/(float)RAND_MAX;
        h_B[j][i]=rand()/(float)RAND_MAX;
    }
    cublasAlloc(N*N, sizeof(float), (void**)&d_A);
    cublasAlloc(N*N, sizeof(float), (void**)&d_B);
    cublasAlloc(N*N, sizeof(float), (void**)&d_C);
    cublasSetVector(N*N, sizeof(float), h_A, 1, d_A, 1);
    cublasSetVector(N*N, sizeof(float), h_B, 1, d_B, 1);
    float t0, gpu_t, cpu_t, error_norm=0, ref_norm=0;
    cudaThreadSynchronize();
    t0=cutGetTimerValue(timer1);
    cublasSgemm('T', 'T', N, N, N, 1.0f, d_A, N, d_B, N, 0.0f, d_C, N);
    cudaThreadSynchronize();
    gpu_t=(cutGetTimerValue(timer1)-t0)/1000.0f;
    cublasGetVector(N*N, sizeof(float), d_C, 1, h_C, 1);
    t0=cutGetTimerValue(timer1);
    simple_sgemm(h_A, h_B, h_C_ref);
    cpu_t=(cutGetTimerValue(timer1)-t0)/1000.0f;
    printf("N=%4d, GPU=%.6fs(%.3fGflops), CPU=%.6fs(%.3fGflops)\n", 
        N, gpu_t, 1e-9*N*N*N*2/gpu_t, cpu_t, 1e-9*N*N*N*2/cpu_t);
    for(int j=0; j<N; j++)
	for(int i=0; i<N; i++) {
        float diff=h_C_ref[j][i]-h_C[j][i];
        error_norm+=diff*diff;
        ref_norm+=h_C_ref[j][i]*h_C_ref[j][i];
    }
    printf("Test %s\n", (sqrtf(error_norm/ref_norm)<1E-6) ? "PASSED" : "FAILED");
}

[weizhuoqun]

接着上面的提问，是不是显卡不支持double型的运算，就不能用Function cublasDgemm()？谢谢

[cuda2010]

cublasSetVector中的incx/incy参数表示从x/y数组中每incx/incy个元素选1个，这样就可以方便地从一个二维矩阵中选出一个列向量或一个行向量。
用cublasSetMatrix可以，并且应该是更好的。不过你对第3和第6个参数的理解不对，应该是这样：

cublasSetMatrix (N, N, sizeof(float), h_A, N, d_A, N);

cublasDgemm必须硬件支持double才行，cublas手册上是这么说的：
Double precision functions are only supported on GPUs with double precision hardware.

[weizhuoqun]

您好，看到您关于CUBLAS的测试，我很感兴趣，也开始学习。其中有个疑问，就是关于数据传递的问题，您用到的是cublasSetVector(N*N, sizeof(float), h_A, 1, d_A, 1);其中的两个1表示的参数含义可以通俗点解释一下吗？根据您的例子，在这里我可不可以用cublasSetMatrix (N, N, N*N,h_A,N,h_B,N)。谢谢

[weizhuoqun]

对对，谢谢你的解答。举个例子，关于数据传递和申请GPU空间的时候，我是不是也可以用cudaMemcpy和cudaMalloc呢？

[mathsoperator]

楼主是个有心人,内容不错,赞一个!

[毕达哥拉斯半圆]

这个要关注！