并行程序设计大作业——稀疏神经网络推理

进击的博仔

已于 2022-09-05 11:37:24 修改

阅读量1.1k

点赞数 3

分类专栏： CUDA入门文章标签：神经网络矩阵机器学习线性代数

于 2022-06-10 01:05:15 首次发布

本文链接：https://blog.csdn.net/weixin_45773137/article/details/125200838

版权

本文详述了一次并行程序设计大作业，任务是优化稀疏神经网络的推理过程，特别是矩阵乘法和bias+relu操作。通过对cublasSgemm函数的参数配置和矩阵存储方式的优化，以及对bias和relu操作的并行化处理，成功将计算时间从16.5s减少到6.7s，显著提升了效率。

摘要由CSDN通过智能技术生成

作业背景与要求

本次作业是 2019 的一个并行竞赛题，也不记得是那个竞赛了，反正题目大致要求下述所示：将一个 60000 * 1024 的矩阵 A 当做左乘矩阵，120 个 1024 * 1024 的稀疏矩阵 B [i], $0\le i\le 119$ 作为右乘矩阵。A 与 B[i] 相乘得到矩阵 C，C 的维度为 60000 * 1024，对 C 矩阵加偏置项 bias，并用激活函数 relu 激活，激活后的 C 矩阵作为下一次相乘的矩阵 A 与 B[i+1] 相乘。即神经网络前向传播一次。运算时间通过 gettimeofday(&time,NULL); 获取。

原始代码是通过传统的矩阵乘法，利用三层循环，每个内部循环计算所求矩阵的一个元素，时间复杂度为 O(n³)，计算一次矩阵相乘的时间需要 10 分钟左右。120个计算下来不得了。所以本次作业的目的就是尽量优化代码，缩短矩阵相乘和 bias+relu 的时间。老师给上一届优化最快的一个小组奖励了樱花键盘，我们这一届也不知道会不会有，反正我不是最快的。

优化思路来源

本学期一共有 5 次家庭作业，4 次都是关于优化矩阵乘法的，我就是筛选出性能最好的一种算法应用于大作业。

第一次作业是两个方阵相乘，方阵大小在 100 ~ 2000 之间，并100 为步长共 20 组数据测试。改变循环顺序和使用 openmp 还有 openblas。

上图为改变循环顺序的结果。

上图为改变循环顺序后加openmp的结果，有两个顺序加openmp语句后输出结果不正确，就没有画出来。

在这里插入图片描述
上图为 openblas 的结果。

void gemm_OpenBlas(double *A, double *B, double *C, int m, int k, int n)
{
   
    enum CBLAS_ORDER order = CblasRowMajor;
    enum CBLAS_TRANSPOSE transposeA = CblasNoTrans;
    enum CBLAS_TRANSPOSE transposeB = CblasNoTrans;
    double alpha = 1;
    double beta = 1;
    cblas_dgemm(order,transposeA,transposeB,m,n,k,alpha,A,k,B,n,beta,C,n);
}
gemm_OpenBlas(A, B, C_yours, m, k, n);

很快啊！ $2000\times2000$ 的两个矩阵相乘只需要 0.25s。

作业三是使用 CUDA 函数继续作业一。

//cublasSgemm
cublasHandle_t s;
cublasCreate(&s);
cublasSgemm('N', 'N', m, n, k, 1.0f, d_A, m, d_B, k, 0, d_C, m);
cublasDestroy(s);
//cublasSgemm_v2
cublasHandle_t s;
cublasCreate_v2(&s);
cublasSgemm_v2(s,CUBLAS_OP_T,CUBLAS_OP_T,m,n,k,&al,d_A,k,d_B,n,&ve,d_C,n);
cublasDestroy_v2(s);
//cblas_sgemm
void gemm_OpenBlas(float *A, float *B, float *C, int m, int k, int n) {
   
    enum CBLAS_ORDER order = CblasRowMajor;
    enum CBLAS_TRANSPOSE transposeA = CblasNoTrans;
    enum CBLAS_TRANSPOSE transposeB = CblasNoTrans;
    float alpha = 1;
    float beta = 1;
    cblas_sgemm(order,transposeA,transposeB,m,n,k,alpha,A,k,B,n,beta,C,n);
}
//cuda_shared
typedef struct {
   
    int width;
    int height;
    int stride;
    float *elements;
} Matrix;

__device__ float GetElement(const Matrix A, int row, int col) {
   
    return A.elements[row * A.stride + col];
}

__device__ void SetElement(Matrix A, int row, int col,
                           float value) {
   
    A.elements[row * A.stride + col] = value;
}

__device__ Matrix GetSubMatrix(Matrix A, int row, int col) {
   
    Matrix Asub;
    Asub.width
            = BLOCK_SIZE;
    Asub.height
            = BLOCK_SIZE;
    Asub.stride
            = A.stride;
    Asub.elements = &A.elements[A.stride * BLOCK_SIZE * row
                                + BLOCK_SIZE * col];
    return Asub;
}

#define ifin(r, c, rb, cb, mat_A) \
((c + cb*BLOCK_SIZE < mat_A.width) && (r+rb*BLOCK_SIZE < mat_A.height))

__global__ void MatMulKernel_SharedMemory(Matrix A, Matrix B, Matrix C) {
   
    int blockRow = blockIdx.y;
    int blockCol = blockIdx.x;
    int row = threadIdx.y;
    int col = threadIdx.x;
    __shared__ float As[BLOCK_SIZE][BLOCK_SIZE];
    __shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE];

    Matrix subC = GetSubMatrix(C, blockRow, blockCol);
    int tot = A.width / BLOCK_SIZE + (A.width % BLOCK_SIZE != 0);
    float CValue = 0;
    for (int i = 0; i < tot; ++i) {
   
        Matrix Asub = GetSubMatrix(A, blockRow, i);

        Matrix Bsub = GetSubMatrix(B, i, blockCol);
        //if (i * BLOCK_SIZE + col < A.width)
        if(ifin(row,col,blockRow,i,A)) {
   
            As[row][col] = GetElement(Asub, row, col);
        } else As[row][col] = 0;

        if(ifin(row,col,i,blockCol,B)) {
   
            Bs[row][col] = GetElement(Bsub, row, col);
        }else Bs[row][col] = 0;
        __syncthreads();

        for (int e = 0; e < BLOCK_SIZE; ++e)
            CValue += As[row][e] * Bs[e][col];

        __syncthreads();
        if (ifin(row, col, blockRow, blockCol, C))
            SetElement(subC, row, col, CValue);
    }

}

void gemm_cuda_shared(float *A, float *B, float *C, int m, int k, int n,double *time_value) {
   
    Matrix d_A;
    d_A.width = d_A.stride = k;
    d_A.height = m;
    size_t size = k * m * sizeof(float);
    cudaMalloc(&d_A.elements, size);
    cudaMemcpy(d_A.elements, A, size, cudaMemcpyHostToDevice);
    Matrix d_B;
    d_B.width = d_B.stride = n;
    d_B.height = k;
    size = n * k * sizeof(float);
    cudaMalloc(&d_B.elements, size);
    cudaMemcpy(d_B.elements, B, size, cudaMemcpyHostToDevice);

    Matrix d_C;
    d_C.width = d_C.stride = n;
    d_C.height = m;
    size = n * m * sizeof(float);
    cudaMalloc(&d_C.elements, size);

    dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);
    dim3 dimGrid((n / dimBlock.x) + (n % dimBlock.x != 0),  (m / dimBlock.y) + (m % dimBlock.y != 0));
                 
    timeval t1,t2;
    *time_value = 0;
    cudaDeviceSynchronize();

    gettimeofday(&t1, nullptr);
    MatMulKernel_SharedMemory<<<dimGrid, dimBlock>>>(d_A, d_B, d_C);
    cudaDeviceSynchronize();
    gettimeofday(&t2, nullptr);
    *time_value+=(t2.tv_sec - t1.tv_sec) + (t2.tv_usec - t1.tv_usec) / 1000000.0;
    
    cudaMemcpy(C, d_C.elements, size, cudaMemcpyDeviceToHost);

    cudaFree(d_A.elements);
    cudaFree(d_B.elements);
    cudaFree(d_C.elements);
}
//cuda_global
__global__ void MatrixMulKernel(float *A, float *B, float *C, int m, int k, int n) {
   
// Each thread computes one element of C
// by accumulating results into Cvalue
    float Cvalue = 0;
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    if (row < m && col < n) {
   
        for (int e = 0; e < k; ++e)
            Cvalue += A[row * k + e]
                      * B[e * n + col];
        C[row * n + col] = Cvalue;
    }
}
void gemm_cuda(float *A, float *B, float *C, int m, int k, int n,double *time_value) {
   
    float *d_A, *d_B, *d_C;
    size_t size = m * k * sizeof(float);
    cudaMalloc(&d_A, size);
    cudaMemcpy(d_A, A, size, cudaMemcpyHostToDevice);

    size = k * n * sizeof(float);
    cudaMalloc(&d_B, size);
    cudaMemcpy(d_B, B, size, cudaMemcpyHostToDevice);


    size = m * n * sizeof(float);
    cudaMalloc(&d_C, size);

    dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);
    dim3 dimGrid((n / dimBlock.x) + (n % dimBlock.x != 0), (m / dimBlock.y) + (m % dimBlock.y != 0));

    timeval t1,t2;
    *time_value = 0;
    cudaDeviceSynchronize();
    
    gettimeofday(&t1, nullptr);
    MatrixMulKernel<<<dimGrid, dimBlock>>>(d_A, d_B, d_C, m, k, n);
    cudaDeviceSynchronize();
    gettimeofday(&t2, nullptr);
    *time_value += (t2.tv_sec - t1.tv_sec) + (t2.tv_usec - t1.tv_usec) / 1000000.0;
    
    cudaMemcpy(C, d_C, size, cudaMemcpyDeviceToHost);

    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_C);

}

最快的是 cublas_sgemm $2000 \times 2000$ 的矩阵只用了大约 0.0075s。其次就是 cublas_sgemm_v2 耗时 0.0175s 左右。

作业四，是稀疏方阵乘瘦矩阵，就是将方阵转化为 csr 的存储格式。然后再按照csr的方式将矩阵相乘。虽然使用csr方式进行计算耗费时间较少，但是需要的前期准备太多，而且在多个矩阵相乘过程中需要多次转换矩阵形式与分配内存，每次参与运算矩阵的 csr 格式表示都不相同，csr 的三个数组长度也不同，所以每次都需重新分配内存、释放内存，会带来较多的时间消耗。

虽然利用 csr 格式的稀疏矩阵相乘有以上的缺点，但还是对此方法进行了实验，下面是在网上拼凑的代码，原理不是很明白，只是粘取了函数部分，然后在调用函数部分做了改动。将矩阵转化为 csr 存储格式后利用cusparseSpMM：