高性能计算-SGEMV矩阵向量乘(30)

1. 介绍

• 矩阵向量乘法: A * X = Y, A(M,K) X(K,1) Y(M,1);
• 实现多种并行算法及优化方法和 cublas 库 sgemv 的效率对比。

2. gpu 并行算法介绍

• 并行算法一：一个线程计算一个结果元素；
• 并行算法二：使用合并访存,需要将输入数据转置；
• 并行算法三: 合并访存+常量内存,需要将输入数据转置；
• 并行算法四: 合并访存 + 共享内存,需要将输入数据转置；

3. 源码

• common.h

点击查看代码

#include <stdlib.h>
#include <sys/time.h>
#include <stdio.h>
#define CUDA_ERROR_CHECK
#define PRINTF_LINE printf("%d final\n",__LINE__);\
fflush(stdout)

//宏定义检查API调用是否出错
#define CudaSafeCall(err) __cudaSafeCall(err,__FILE__,__LINE__)
inline void __cudaSafeCall(cudaError_t err,const char* file,const int line)
{
    #ifdef CUDA_ERROR_CHECK
    if(err!=cudaSuccess)
    {
        fprintf(stderr,"cudaSafeCall failed at %s:%d :(%d) %s\n",file,line,err,cudaGetErrorString(err));
        exit(-1);
    }
    #endif
}

//宏定义检查获取流中的执行错误，主要是对核函数
#define CudaCheckError() _cudaCheckError(__FILE__,__LINE__)
inline void _cudaCheckError(const char * file,const int line)
{
    #ifdef CUDA_ERROR_CHECK
    cudaError_t err = cudaGetLastError();
    if(err != cudaSuccess)
    {
        fprintf(stderr,"cudaCheckError failed at %s:%d :(%d) %s\n",file,line,err,cudaGetErrorString(err));
        exit(-1);
    }
    #endif
}

//CPU时间，用于串行实现计时 s
double cpuSecond()
{
    struct timeval tp;
    gettimeofday(&tp, NULL);
    return ((double)tp.tv_sec + (double)tp.tv_usec * 1e-6);
}

void CheckResult(float *arr,float *arrD,int len)
{
    for(int i=0;i<len;i++)
    {
        if(fabs(arr[i] - arrD[i])>1e-6)
        {
            printf("第%d个数值不相等,arr:%.3f arrD:%.3f\n",i,arr[i],arrD[i]);
            exit(0);
        }
    }
}

• sgemv.cu

/*
矩阵向量乘法: A * X = Y, A(M,K) X(K,1) Y(M,1)
备注:编译参数 -lcubasl -arch=sm_XX(根据显卡实际参数填写)
执行参数: 
1: mode (0:kernel_native 1:kernel_Coalesced 2:kernel_CoalescedConst 3:kernel_CoalescedShared)
2: M 1<< 的位移参数
3: K 1<< 的位移参数
*/

#include "common.h"
#include <cublas_v2.h>
#include <cuda_runtime.h>

void serialCpu(float *A,int M,int K,float *X,float *Y)
{
    for(int i=0;i<M;i++)
    {
        float tempY= 0.0;
        for(int j=0;j<K;j++)
        {
            tempY += A[i*K+j] * X[j];
            #if 0
            if(i==25 && j>1200 && j<1280)
                printf("%.f \n",tempY);
            #endif
        }
        Y[i] = tempY;
    }
}

// 并行算法一: 一个线程计算一个结果元素
__global__ void kernel_native(float *A,int M,int K,float *X,float *Y)
{
    int id = blockIdx.x * blockDim.x + threadIdx.x;   //全局id
    if(id<M)
    {
        float tempY = 0.0;
        for(int i=0;i<K;i++)
        {
            tempY += A[id*K+i] * X[i];
            // if(id==0 && i<20)
            //     printf("%.f %d %.f %d %.f\n",tempY,id*K+i,A[id*K+i],i,X[i]);
        }
        Y[id] = tempY;
    }
}

//并行算法二:使用合并访存,需要将输入矩阵A转置
__global__ void kernel_Coalesced(float *AT,int M,int K,float *X,float *Y)
{
    int id = blockIdx.x * blockDim.x + threadIdx.x;   //全局id
    if(id<M)
    {
        float tempY = 0.0;
        for(int i=0;i<K;i++)
            tempY += AT[i*M+id] * X[i];
        Y[id] = tempY;
    }
}

//并行算法三:合并访存+常量内存,需要将输入矩阵A转置
__constant__ float xConst[1<<14];//32KB,为了不影响程序,设置小一点
__global__ void kernel_CoalescedConst(float *AT,int M,int K,float *X,float *Y)
{
    int id = blockIdx.x * blockDim.x + threadIdx.x;   //全局id
    if(id<M)
    {
        float tempY = 0.0;
        for(int i=0;i<K;i++)
            tempY += AT[i*M+id] * xConst[i];
        Y[id] = tempY;
    }
}

//并行算法四:一个线程计算一个元素,合并访存 + 共享内存,需要将输入矩阵A转置；
//如果K过大,对向量X在K方向划分block,总线程数量=M
template <size_t blcokSize>
__global__ void kernel_CoalescedShared(float *AT,int M,int K,float *X,float *Y)
{
    int tid = threadIdx.x;
    int id = blockIdx.x * blockDim.x + threadIdx.x;   //全局id
    if(id<M)
    {
        __shared__ float XS[blcokSize];
        float tempY = 0.0;
        for(int k=0;k<K;k+=blcokSize)
        {
            XS[tid] = X[k+tid];
            __syncthreads();
            for(int i=0;i<blcokSize;i++)
            {
                tempY += AT[(k+i)*M + id] * XS[i];
                // if(id==25 && k<128)
                //     printf("%.f %d %.f %d %.f\n",tempY,(k+i)*M + id,AT[(k+i)*M + id],i,XS[i]);
            }
            __syncthreads();
        }
        Y[id] = tempY;
    }
}

int main(int argc,char** argv)
{
    int mode = 0;
    int M,K;
    M=K=1<<6;    //输入数据为32倍数
    if(argc>1)
        mode = atoi(argv[1]);
    if(argc==4)
    {
        M = 1 << atoi(argv[2]);
        K = 1 << atoi(argv[3]);
    }
    float *A,*AT,*X,*Y_blas,*Y_serial,*Y_parallel;//AT :存放A转置数据
    A = (float*)calloc(M*K,sizeof(float));
    AT = (float*)calloc(K*M,sizeof(float));
    X = (float*)calloc(K,sizeof(float));
    Y_blas = (float*)calloc(M,sizeof(float));
    Y_serial = (float*)calloc(M,sizeof(float));
    Y_parallel = (float*)calloc(M,sizeof(float));   
        
    //矩阵初始化
    #if 1
    srand(0);
    for(int i=0;i<M*K;i++)
        AT[(i%K)*M + i/K] = A[i] = rand()%20;    //A[i/K][i%K] -> AT[i%K][i/K]
    for(int i=0;i<K;i++)
        X[i] = rand()%10;
    #else   //FIXME 待处理:使用一下初始化代码,串行、并行、cublas 在数据长度K较大比如 12800 时计算结果不一致
    for(int i=0;i<M*K;i++)
        AT[(i%K)*M + i/K] = A[i] = i/10;
    for(int i=0;i<K;i++)
        X[i] = i%10;
    #endif
    #if 0
    for(int i=0;i<M;i++)
    {
        for(int j=0;j<K;j++)
            printf("%.f ",A[i*K+j]);
        printf("\n");
    }
    printf("============\n");
    for(int j=0;j<K;j++)
        printf("%.f ",X[j]);
    printf("\n============\n");
    #endif
    //cpu
    double cpuTime = 0.0;
    #if 0
    double start = cpuSecond();
    serialCpu(A,M,K,X,Y_serial);
    cpuTime = (cpuSecond()-start)*1000;
    #endif

    //gpu
    float *AD,*ATD,*XD,*YD; //ATD 存放AD转置数据
    cudaMalloc((float**)&AD,M*K*sizeof(float));
    cudaMalloc((float**)&XD,K*sizeof(float));
    cudaMalloc((float**)&YD,M*sizeof(float));
    float gpuTime_Blas = 0.0;
    float gpuTime_parallel = 0.0;
    
    //cublas base
    cublasHandle_t cublas;
    cudaEvent_t startD,endD;
    cudaEventCreate(&startD);
    cudaEventCreate(&endD);
    cublasCreate(&cublas);
    float alpha = 1.0,beta = 0.0;
    cudaMemcpy(AD,A,M*K*sizeof(float),cudaMemcpyHostToDevice);
    cudaMemcpy(XD,X,M*sizeof(float),cudaMemcpyHostToDevice);

    cudaEventRecord(startD);
    //内存是按列优先存储
    //Y = A * X ,结果 Y 按列存储,可以直接求 Y^T = (AX)^T=X^T * A^T,结果Y^T是正确结果的按行存储
    cublasSgemv(cublas,CUBLAS_OP_T,K,M,&alpha,AD,K,XD,1,&beta,YD,1);
    cudaEventRecord(endD);
    cudaEventSynchronize(endD);
    cudaEventElapsedTime(&gpuTime_Blas,startD,endD);
    cudaMemcpy(Y_blas,YD,M*sizeof(float),cudaMemcpyDeviceToHost);
    cublasDestroy(cublas);

    #if 0
    printf("serial\n");
    for(int j=0;j<M;j++)
        printf("%.f ",Y_serial[j]);
    printf("============\n");
    printf("blas\n");
    for(int j=0;j<M;j++)
        printf("%.f ",Y_blas[j]);
    printf("============\n");
    #endif
    // CheckResult(Y_blas,Y_serial,M);
    printf("A size: %d * %d\n",M,K);
    printf("cpu串行 耗时: %.2f ms\n",cpuTime);
    printf("cublas 耗时: %.2f ms\n",gpuTime_Blas);

    if(mode==0) //并行算法2: 一个线程计算一个结果元素
    {
        cudaMemset(YD,0,M*sizeof(float));
        int blockSize = 128;
        int gridSize = (M-1)/128 +1;
        // printf("g %d ,b %d\n",gridSize,blockSize);
        cudaEventRecord(startD);
        kernel_native<<<gridSize,blockSize>>>(AD,M,K,XD,YD);
        cudaEventRecord(endD);
        cudaEventSynchronize(endD);
        cudaEventElapsedTime(&gpuTime_parallel,startD,endD);
        cudaMemcpy(Y_parallel,YD,M*sizeof(float),cudaMemcpyDeviceToHost);
        #if 0
        for(int i=12800-10;i<12800;i++)
            printf("%.f ",Y_serial[i]);
        printf("\n");
        for(int i=12800-10;i<12800;i++)
            printf("%.f ",Y_parallel[i]);
        printf("\n");
        #endif
        printf("kernel_native  耗时: %.2f ms\n",gpuTime_parallel);
        CheckResult(Y_blas,Y_parallel,M);
    }
    else if(mode ==1)//并行算法三,合并访存
    {    
        cudaFree(AD);AD=NULL;
        cudaMemset(YD,0,M*sizeof(float));
        cudaMemset(Y_parallel,0,M*sizeof(float));
        //需要矩阵A转置
        cudaMalloc((void**)&ATD,M*K*sizeof(float));
        cudaMemcpy(ATD,AT,M*K*sizeof(float),cudaMemcpyHostToDevice);
        int blockSize = 128;
        int gridSize = (M-1)/128 +1;
        cudaEventRecord(startD);
        kernel_Coalesced<<<gridSize,blockSize>>>(ATD,M,K,XD,YD);
        cudaEventRecord(endD);
        cudaEventSynchronize(endD);
        cudaEventElapsedTime(&gpuTime_parallel,startD,endD);
        cudaMemcpy(Y_parallel,YD,M*sizeof(float),cudaMemcpyDeviceToHost);
        printf("kernel_Coalesced  耗时: %.2f ms\n",gpuTime_parallel);
        CheckResult(Y_blas,Y_parallel,M);
        cudaFree(ATD);
    }
    else if(mode ==2)   //合并访存 + 常量内存
    {
        //获取设备常量内存大小
        cudaDeviceProp deviceProp;
        cudaGetDeviceProperties(&deviceProp,0);
        size_t nBytes = deviceProp.totalConstMem;
        // printf("constMemSize: %u KB\n",nBytes/1024);    //k80 64KB

        if(M <= nBytes/4)//16KB 为了不影响其他程序,设置小一点
        {
            cudaFree(AD);AD=NULL;
            cudaMemset(YD,0,M*sizeof(float));
            cudaMemset(Y_parallel,0,M*sizeof(float));
            //需要矩阵A转置
            cudaMalloc((void**)&ATD,M*K*sizeof(float));
            cudaMemcpy(ATD,AT,M*K*sizeof(float),cudaMemcpyHostToDevice);
            cudaMemcpyToSymbol(xConst,X,M*sizeof(float));
            int blockSize = 128;
            int gridSize = (M-1)/128 +1;
            cudaEventRecord(startD);
            kernel_CoalescedConst<<<gridSize,blockSize>>>(ATD,M,K,XD,YD);
            cudaEventRecord(endD);
            cudaEventSynchronize(endD);
            cudaEventElapsedTime(&gpuTime_parallel,startD,endD);
            cudaMemcpy(Y_parallel,YD,M*sizeof(float),cudaMemcpyDeviceToHost);
            printf("kernel_CoalescedConst  耗时: %.2f ms\n",gpuTime_parallel);
            CheckResult(Y_blas,Y_parallel,M);
            cudaFree(ATD);
            // cudaFree(xConst);
        }
        else
            printf("数据量 %d 过大,超出常量内存\n",M);
    }
    else if(mode ==3)   //合并访存 + 共享内存
    {
        cudaFree(AD);AD=NULL;
        cudaMemset(YD,0,M*sizeof(float));
        cudaMemset(Y_parallel,0,M*sizeof(float));
        //需要矩阵A转置
        cudaMalloc((void**)&ATD,M*K*sizeof(float));
        cudaMemcpy(ATD,AT,M*K*sizeof(float),cudaMemcpyHostToDevice);
        const int blockSize = 128;
        int gridSize = (M-1)/128 +1;
        cudaEventRecord(startD);
        kernel_CoalescedShared<blockSize><<<gridSize,blockSize>>>(ATD,M,K,XD,YD);
        cudaEventRecord(endD);
        cudaEventSynchronize(endD);
        cudaEventElapsedTime(&gpuTime_parallel,startD,endD);
        cudaMemcpy(Y_parallel,YD,M*sizeof(float),cudaMemcpyDeviceToHost);
        printf("kernel_CoalescedShared  耗时: %.2f ms\n",gpuTime_parallel);
        #if 0
        for(int i=25;i<35;i++)
            printf("%.f ",Y_blas[i]);
        printf("\n");
        for(int i=25;i<35;i++)
            printf("%.f ",Y_parallel[i]);
        printf("\n");
        #endif
        CheckResult(Y_blas,Y_parallel,M);
        cudaFree(ATD);
    }

    cudaEventDestroy(startD);
    cudaEventDestroy(endD);
    free(A);
    free(X);
    free(Y_blas);
    free(Y_serial);
    free(Y_parallel);
    if(AD!=NULL)
        cudaFree(AD);
    cudaFree(XD);
    cudaFree(YD);
    return 0;
}

4. 测试数据

• 当矩阵 M= K = 2^14 数据量时：

测试项	耗时ms	与 cublasSgemv 浮点计算能力比较
kernel_native	87.9	9.70%
kernel_Coalesced(合并访存)	7.54	113%
kernel_CoalescedConst(合并访存+常量内存)	7.19	119%
kernel_CoalescedShared(合并访存+共享内存)	7.15	119.50%
cublasSgemv	8.55	100%

• 测试不同矩阵大小规模下，与 cublas 库函数计算能力比较

矩阵大小(M=K)	2^12	2^13	2^14	2^15
cublasSgemv(ms)	0.59	2.18	8.55	34.24
kernel_CoalescedShared(ms)	0.78	2.05	7.15	33.2
与 cublasSgemv 浮点计算能力比较	75.60%	106%	119.50%	103%

5. 结果分析：

• 合并访存对效率的提升特别明显；
• 随着矩阵规模的增大，自实现版本效率比较高，比 cublasSgemv 有略微优势。