@[TOC](NVIDIA CUDA 高度并行处理器编程（六）：并行模式：卷积)介绍常数存储器和高速缓存 后面几个文章介绍一些重要的并行计算模式。这些模式是很多并行算法的基础。先介绍卷积，卷积广泛应用于信号处理、图像处理和计算机视觉。卷积通常涉及每个元素上的大量算术运算，比如清晰的图像和视频，计算量很大，不过还好每个输出元素的计算都是相互独立的，我们可以应用并行处理。另一方面，输入元素之间有一定程度的数据共享，还存在一些颇具挑战的边界条件，这使得卷积成为复杂分块方法和输入数据分段的一个重要案例。 背景 卷积是一种数组操作，每个输出元素都是周围输入元素的加权总和。权重是由一个输入掩码数组也叫卷积核（convolution kernel）确定的。下面卷积计算都默认padding。 在音频数字信号处理中，输入数据都是一维的。下面两图为一维卷积的实例： M为卷积核，N为输入数据。 图7.1即计算： 图7.3即计算： 计算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高速缓存，但只是可能，也有可能从全局存储器中加载光环元素，这是随机的。