Outline of Programming Massively Parallel Processors Parallel Concepts Part

写在最前

本页面内容仅对于毫无基础的初学者想要快速了解可供参考,深度学习者仍建议自行阅读Programming Massively Parallel Processors书籍原文.

在本页面的所有内容中,我会优先针对整体内容进行讲解,为读者构建基础的印象感知,同时方便温习使用.此外,由于此部分内容强依赖代码,所以会有对实现代码的讲解以及参考代码.

以下内容为基于自身理解的梗概,内容难免有错漏,欢迎各位阅读者批评指正!

再次声明,以下仅供参考

Convolution 卷积

Introduction 引入

有时我们会遇到加权求和的问题,比如针对一个一维向量的每一个元素,它的加权输出为以它为中心的五个连续元素的加权和,如下图模式:

通过一个加权规则,对输入向量的每一个元素进行依次计算,得到输出向量

类似这样的加权规则,就是针对一维向量的卷积,其中,就是卷积核(convolution kernal,filter)

在以上例子中,对于边界元素,由于其左右两侧元素量不足,针对超出边界的部分,通常默认补0处理.

2D Convolution 二维卷积

在实际运用中,我们通常遇到的是二维的卷积问题

如上所示,针对每一个输入元素,都会利用卷积核进行加权计算.图中的卷积核是卷积核.对于边界元素,依旧采取补零操作,超出部分默认为0.

针对这一实例,我们给出以下基础CUDA代码实现

__global__
void convolution_2D_basic_kernal(float* N, float* F, float* P, int r, int width, int height){
    /*
    N:输入矩阵指针
    F:卷积核指针
    P:输出矩阵指针
    r:卷积核半径
    width:矩阵宽度
    height:矩阵高度
    */
    //计算输出矩阵元素索引
    int outCol = blockIdx.x*blockDim.x + threadIdx.x;
    int outRow = blockIdx.y*blockDim.y + threadIdx.y;
    float Pvalue = 0.0f;
    //确定卷积位置以及求和累加
    for (int fRow = 0; fRow < 2*r+1; fRow++){
        for (int fCol = 0; fCol < 2*r+1; fCol++){
            inRow = outRow - r + fRow;
            inCol = outCol - r + fCol;
            if (inRow >= 0 && inRow < height && inCol >= 0 && inCol < width){
                Pvalue += F[fRow][fCol] * N[inRow*width + inCol];
            }
        }
    }
    //写入输出矩阵
    P[outRow][outCol] = Pvalue;
}

然而,对于上述核函数,有极大的性能短板.

首先,即是内存带宽限制.上述函数的计算密度极低,而且几乎读取写入操作全部都是在显存进行的,这会导致其内存利用极低,需要进行优化.

对此,由于卷积核是固定的,我们可以考虑利用常量内存(constant memory)以及缓存机制来进行优化.因为访问constant memory比访问显存更快,所以可直接将卷积核元素写入constant memory,进而省下一部分访问开销.

此外,对于并行,应当进行分块(tiling)操作,让每一个block对应一个分块,将其元素载入共享内存(shared memory)以加大内存访问效率.由于每个块内容载入shared memory时,考虑到卷积核大小,输入的元素会实际上比块更大,多出来的一圈元素叫做halo cells,如下图浅蓝色部分所示

对应代码如下:

#define IN_TILE_DIM 32
#define OUT_TILE_DIM ((IN_TILE_DIM) - 2*(FILTER_RADIUS)) //定义输入矩阵分块和输出矩阵分块大小,显然输入矩阵比输出矩阵分块大
__constant__ float F_c[2*FILTER_RADIUS + 1][2*FILTER_RADIUS + 1];
__global__
void convolution_tiled_2D_const_mem_kernal(float* N, float* P, int width, int height){
    int col = blockIdx.x*OUT_TILE_DIM + threadIdx.x - FILTER_RADIUS;
    int row = blockIdx.y*OUT_TILE_DIM + threadIdx.y - FILTER_RADIUS;
    //载入输入分块
    __shared__ N_s[IN_TILE_DIM][IN_TILE_DIM];
    if (row >= 0 && row < height && col >= 0 && col < width){
        N_s[threadIdx.y][threadIdx.x] = N[row*width + col];
    }else{
        N_s[threadIdx.y][threadIdx.x] = 0;
    }
    __syncthreads();
    //计算输出矩阵元素
    int tileCol = threadIdx.x - FILTER_RADIUS;
    int tileRow = threadIdx.y - FILTER_RADIUS;
    //边界元素分支控制,且限制并停止分块边界部分线程
    if (col >= 0 && col < width && row >= 0 && row < height){
        if (tileCol >= 0 && tileCol < OUT_TILE_DIM && tileRow >= 0 && tileRow < OUT_TILE_DIM){
            float Pvalue = 0.0f;
            for (int fRow = 0; fRow < 2*FILTER_RADIUS + 1; fRow++){
                for (int fCol = 0; fCol < 2*FILTER_RADIUS + 1; fCol++){
                    Pvalue += F[fRow][fCol]*N_s[tileRow+fRow][tileCol+fCol];
                }
            }
            P[row*width + col] = Pvalue;
        }
    }
}