2022-08-24 13:52 中国石油大学（北京）后端发布于河南

关注

CUDA入门（六）：并行模式：卷积

@[TOC](NVIDIA CUDA 高度并行处理器编程（六）：并行模式：卷积)
介绍常数存储器和高速缓存

后面几个文章介绍一些重要的并行计算模式。这些模式是很多并行算法的基础。先介绍卷积，卷积广泛应用于信号处理、图像处理和计算机视觉。卷积通常涉及每个元素上的大量算术运算，比如清晰的图像和视频，计算量很大，不过还好每个输出元素的计算都是相互独立的，我们可以应用并行处理。另一方面，输入元素之间有一定程度的数据共享，还存在一些颇具挑战的边界条件，这使得卷积成为复杂分块方法和输入数据分段的一个重要案例。

背景

卷积是一种数组操作，每个输出元素都是周围输入元素的加权总和。权重是由一个输入掩码数组也叫卷积核（convolution kernel）确定的。下面卷积计算都默认padding。

在音频数字信号处理中，输入数据都是一维的。下面两图为一维卷积的实例：
在这里插入图片描述

在这里插入图片描述
M为卷积核，N为输入数据。

图7.1即计算： $P[2] = N[0]*M[0] + N[1]*M[1] + N[2]*M[2] + N[3]*M[3] + N[4]*M[4]=57$

图7.3即计算： $P[1] = 0 * M[0] + N[0]*M[1] + N[1]*M[2] + N[2]*M[3] + N[3]*M[4]=38$

计算p[1]时，由于N[1]左边只有一个元素，缺少的那个元素用 0 来代替。这些缺失的元素通常称为“幽灵元素”，在并行计算中还会引进其他类型的幽灵元素。这些幽灵元素对分块算法的复杂度和效率影响很大。

对于图像处理和计算机视觉来说，输入数据通常都是二维的。下面为二维图像卷积实例：

在这里插入图片描述
计算原理如下图：

带有边界条件的卷积，超出范围的元素当0计算：

一维并行卷积

由于卷积计算与顺序无关，所以可以使用并行解决，先声明kernel函数：

__global__ void convolution_1D_basic_kernel(float *N, float *M, float *P, 
    int Mask_Width, int Width);

N 为输入数组。
M 为卷积核。
P 为输出数组。
Mask_Width 为卷积核尺寸。
Width 为输入输出数组长度。

确定数据和线程的映射，由于输入数据是一维的，所以将网格设置为一维，每个线程负责一个输出元素的计算。
int i = blockId.x * blockDim.x + threadId.x;
确定 i 后，就可以确定 N 和 M 的索引。Mask_Width 一般为奇数（2*n+1），那么 N 要参与运算元素的索引就是 i-n ~ i+n，即 i-n+0 ~ i-n+Mask_Width。

下面是代码：

#include<stdio.h>
#include<stdlib.h>
#include<cuda.h>
#define BLOCKSIZE 1024
__global__ void convolution_1D_basic_kernel(float *N, float *M, float *P, int Mask_Width, int Half_Mask_Width, int Width){
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    float PValue = 0;
    for(int j = 0;j < Mask_Width;j++){
        if(i - Half_Mask_Width + j >= 0 && i - Half_Mask_Width + j < Width ){
            PValue += N[i - Half_Mask_Width + j]*M[j];
        }
    }
    P[i] = PValue;
}

void Convolution_1D_basic(float *N, float *M, float *P, int Mask_Width, int Half_Mask_Width, int Width){
    dim3 dimBlock(BLOCKSIZE, 1, 1);
    dim3 dimGrid(ceil((float)Width/BLOCKSIZE), 1, 1);
    float *d_N, *d_M, *d_P;
    cudaMalloc((void **) &d_N, sizeof(float)*Width);
    cudaMalloc((void **) &d_M, sizeof(float)*Mask_Width);
    cudaMalloc((void **) &d_P, sizeof(float)*Width);
    cudaMemcpy(d_N, N, sizeof(float)*Width, cudaMemcpyHostToDevice);
    cudaMemcpy(d_M, M, sizeof(float)*Mask_Width, cudaMemcpyHostToDevice);

    convolution_1D_basic_kernel<<<dimGrid, dimBlock>>>
                        (d_N, d_M, d_P, Mask_Width, Half_Mask_Width, Width);

    cudaMemcpy(P, d_P, sizeof(float)*Width, cudaMemcpyDeviceToHost);
    cudaFree(d_P);
    cudaFree(d_M);
    cudaFree(d_N);
}
int main(){
    int width = 7;
    int mask_width = 5;
    float *N = (float*)malloc(sizeof(float)*width);
    float *M = (float*)malloc(sizeof(float)*mask_width);
    float *P = (float*)malloc(sizeof(float)*width);
    for(int i = 0;i < width;++i){
        N[i] = i+1;
    }
    M[0] = 3;    M[1] = 4;    M[2] = 5;    M[3] = 4;    M[4] = 3;

    Convolution_1D_basic(N, M, P, mask_width, mask_width/2, width);

    for(int i = 0;i < width;++i){
        printf("%f ", P[i]);
    }
    return 0;
}

在这里插入图片描述
上面的 kernel 函数有两个不足：

会出现控制流多样性，靠近输入矩阵两端的线程要处理幽灵元素，即一个warp中的线程处理了不同的语句。
kernel 函数中浮点算术运算和全局存储器访问的比率仅为1.0。正如在矩阵乘法全局内存版中一样，效率很低。

常数存储器和高速缓存

先看看M的几个特征：
尺寸小，大多数卷积核某维度长度不大于10，就算是三维的也不多与于1000个数。
内容不变。
所有线程都以相同顺序访问M，从M[0]到M[Mask_Width-1]。

这几个特性都使卷积核能利用常数存储器和高速缓存。

常数存储器与其他类型存储器关系如下图所示：
在这里插入图片描述
常数存储器与全局存储器类似，所有线程都能访问其中的变量，但其中的变量不能被修改。不同设备的常数存储器大小不同，MX250通过下面程序查出常数存储器大小为65536字节。

#include<cuda.h>
#include<stdio.h>
int main(){
    int dev_count;
    cudaGetDeviceCount(&dev_count); 
    printf("number of GPU: %d\n", dev_count);
    cudaDeviceProp dev_prop;
    for(int i = 0; i < dev_count;++i){
        cudaGetDeviceProperties(&dev_prop, i);
        printf("volume of const memory of GPU %d: %d\n", i + 1, dev_prop.totalConstMem)
    }
    return 0;
}

在主机代码中使用如下语句分配和复制全局常数存储器变量：

#define MAX_MASK_WIDTH 10
__constant__ float M[MAX_MASK_WIDTH];
cudaMemcpyToSymbol(M, h_M, Mask_Width*sizeof(float));

kernel 函数访问常数存储器的方式和访问全局存储器的一样。这里不需要将M传入kernel，kernel函数通过主机代码声明的全局变量来访问。虽然常数存储器的实现也是DRAM，但是CUDA运行时系统知道常数存储器变量不会改变，所以会将其直接放入高速缓存中。

__global__ void convolution_1D_ba sic_kernel(float *N, float *P, int Mask_Width,int Width) {
    int i = blockIdx.x*blockDim.x + threadIdx.x;
    float Pvalue = 0;
    int N_start_point = i - (Mask_Width/2);
    for (int j = 0; j < Mask_Width; j++) {
        if (N_start_point  +j>=0 && N_start_point  +j<Width) {
        Pvalue += N[N_start_point + j]*M[j];
        }
    }
    P[i] = Pvalue;
}

高速缓存的示意图：
在这里插入图片描述
与GPU上的其他存储器不同的是，高速缓存一般是透明的，处理器会自动保留程序最近或经常访问的变量到高速缓存中，并记录他们的原始DRAM地址，当保留的变量再次被使用时，硬件会从他们的地址判断已经在高速缓存保留了他们的值，并直接访问高速缓存，从而消除对DRAM的访问。

在并行处理器中使用高速缓存的一个重要设计问题是缓存一致性，当一个或多个核修改缓存中的数据时容易出现。由于高速缓存通常只连接到一个处理器核，其修改的内容不易被其他存储器察觉，当一个变量被多个处理器核心上的线程所共享时，这种修改就会出问题。

但常数存储器的值不会改变，所以将其中的变量放入高速缓存不会出现缓存一致性的问题。

使用边缘元素的分块一维卷积

假定输入元素有16个，M大小为 5，分块大小为4，则卷积示例如下：
在这里插入图片描述
每个块加载需要的元素到共享存储器中，块内边缘的元素称为“边缘元素”（skirt element）或“光环元素”（halo element），比如块 0 中的N[4]，N[5]，块 1 中的N[2]，N[3]，N[8]，N[9]。块中间部分的元素称为内部元素。

声明共享存储器数组：

__shared__ float N_ds[TILE_SIZE + MASK_WIDTH - 1];
n = MASK_WIDTH/2;

加载左侧光环元素：
使用块内后 n 个线程加载

int halo_index_left = (blockIdx.x - 1) * blockDim.x + threadIdx.x;
if(threadId.x >= blockDim.x - n){
    N_ds[threadId.x - blockDim.x + n] = 
        (halo_index_left < 0)?0:N[halo_index_left];

halo_index_left 将线程索引映射到元素索引上。if 是选取后 n 个线程，N_ds索引为 0 ~ n-1 共n个数。后面的条件表达式是从全局存储器加载元素，如果是幽灵元素，就为0，不是则加载。

加载右侧光环元素：
使用块内前 n 个线程加载

int halo_index_right = (blockId.x + 1)*blockDim.x + threadIdx.x;
if(threadId.x < n){
    N_ds[threadId.x + blockDim.x + n] = 
        (halo_index_right > WIDTH)?0:N[halo_index_right];
}

上面的索引确定还是比较难的，对着上面的图应该好理解些。

加载内部元素：

N_ds[n+threadIdx.x] = N[blockIdx.x*blockDim.x + threadIdx.x];

计算：

float Pvalue = 0;
for(int j = 0; j < MASK_WIDTH;j++){
    Pvalue += M[j]*N_ds[j+threadIdx.x];
}
p[i] = Pvalue;

每个线程利用不同区域，例如线程 0 利用N_ds[0] ~ N_ds[MASK_WIDTH-1]，线程 1 利用N_ds[1] ~ N_ds[MASK_WIDTH]，以此类推。

完整代码：

#include<stdio.h>
#include<stdlib.h>
#include<cuda.h>
#define TILE_SIZE 4
#define MAX_MASK_WIDTH 10
__constant__ float M[MAX_MASK_WIDTH];

__global__ void convolution_1D_basic_kernel(float *N, float *P, int Mask_Width, int Width){
    int i = blockIdx.x*blockDim.x + threadIdx.x;
    __shared__ float N_ds[TILE_SIZE + MAX_MASK_WIDTH - 1];
    int n = Mask_Width/2;

    int halo_index_left = (blockIdx.x - 1) * blockDim.x + threadIdx.x;
    if( threadIdx.x >= blockDim.x - n){
        N_ds[threadIdx.x - blockDim.x + n] = (halo_index_left < 0)?0:N[halo_index_left];
    }

    int halo_index_right = (blockIdx.x + 1) * blockDim.x + threadIdx.x;
    if(threadIdx.x < n){
        N_ds[ threadIdx.x + blockDim.x + n] = (halo_index_right > Width)?0:N[halo_index_right];
    }

    N_ds[n + threadIdx.x] = N[i];

    __syncthreads();
    float Pvalue = 0;
    for(int j = 0; j < Mask_Width;j++){
        Pvalue += M[j] * N_ds[j+threadIdx.x];
    }
    P[i] = Pvalue;

}
void Convolution_1D_basic(float *N, float *h_M, float *P, int Mask_Width, int Width){
    cudaMemcpyToSymbol(M, h_M, Mask_Width*sizeof(float));
    dim3 dimBlock(TILE_SIZE, 1, 1);
    dim3 dimGrid(ceil((float)Width/TILE_SIZE), 1, 1);

    float *d_N,*d_P;
    cudaMalloc((void **) &d_N, sizeof(float)*Width);
    cudaMalloc((void **) &d_P, sizeof(float)*Width);
    cudaMemcpy(d_N, N, sizeof(float)*Width, cudaMemcpyHostToDevice);

    convolution_1D_basic_kernel<<<dimGrid, dimBlock>>>(d_N, d_P, Mask_Width, Width);

    cudaMemcpy(P, d_P, sizeof(float)*Width, cudaMemcpyDeviceToHost);
    cudaFree(d_P);
    cudaFree(d_N);
}
int main(){
    int width = 16;
    int mask_width = 5;
    float *N = (float*)malloc(sizeof(float)*width);
    float *M = (float*)malloc(sizeof(float)*mask_width);
    float *P = (float*)malloc(sizeof(float)*width);
    for(int i = 0;i < width;++i){
        N[i] = i;
    }
    for(int i = 0;i < mask_width;++i){
        M[i] = 1;
    }

    Convolution_1D_basic(N, M, P, mask_width, width);

    for(int i = 0;i < width;++i){
        printf("%f ", P[i]);
    }
    return 0;
}

普通卷积与分块卷积的访存对比：

普通卷积：

内部线程块的访存次数为：blockDimMask_Width 或 blockDim*（2n+1），边缘线程块除去幽灵元素： blockDim(2n+1) - n(n+1)/2。
分块卷积：
内部线程块：blockDim.x + 2n，边缘线程块：blockDim.x + n。

访存次数对比：

所以对于内部线程块，基本算法与分快算法的访存次数比值为：
( blockDim.x*(2n+1) )/(blockDim.x+2n)
对于边缘线程块，比值为：
( (blockDim.x*(2n+1) - n(n+1)/2 )/(blockDim.x+n)

如果blockDim.x 的值比n大的多，那么就近似为：2n+1=Mask_Width。

利用通用高速缓存的分块一维卷积

GPU上有L1和L2高速缓存，L1是每个SM私有的，L2是所有SM共享的，线程块访问的光环元素可能放到了L2高速缓存中。
例如分块 1 中的光环元素可能在分块 0 加载内部元素时已经放入了 L2 缓存中，但也可能没有，因为所有线程块是并行的。所以对光环元素的访问可以直接通过对访问 N 得到，这样可能直接访问 N，也有可能访问L2缓存。

声明共享存储器数组 N_ds，不需要加载光环元素，所以大小为TILE_SIZE。

__shared__ float N_ds[TILE_SIZE]

加载

N_ds[threadIdx.x] = N[i];
__syncthreads();

计算

int This_tile_start_point = blockIdx.x * blockDim.x; 
int Next_tile_start_point = (blockIdx.x + 1) * blockDim.x;
int N_start_point = i - (Mask_Width/2); 
float Pvalue = 0;
for (int j = 0; j < Mask_Width; j++) {
    int N_index = N_start_point + j;
    if (N_index >= 0 && N_index < Width) { //判断是否为幽灵元素
        if ((N_index >= This_tile_start_point)&& 
                    (N_index < Next_tile_start_point)) { //判断是否为光环元素。
             Pvalue += N_ds[threadIdx.x+j-(Mask_Width/2)]*M[j]; 
        } else {
            Pvalue += N[N_index] * M[j];
        }
    }
}
P[i] = Pvalue;

这个比上一个分块简单多了，没有复杂的加载操作，同时还可能用到L2高速缓存，但只是可能，也有可能从全局存储器中加载光环元素，这是随机的。