cublas矩阵乘法效率测试

以前曾听说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");

}