NVIDIA CUDA 高度并行处理器编程(六):并行模式:卷积


介绍常数存储器和高速缓存

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

背景

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

音频数字信号处理中,输入数据都是一维的。下面两图为一维卷积的实例:
NVIDIA CUDA 高度并行处理器编程(六):并行模式:卷积

NVIDIA CUDA 高度并行处理器编程(六):并行模式:卷积
M为卷积核,N为输入数据。

图7.1即计算:
NVIDIA CUDA 高度并行处理器编程(六):并行模式:卷积

图7.3即计算:
NVIDIA CUDA 高度并行处理器编程(六):并行模式:卷积

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

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

NVIDIA CUDA 高度并行处理器编程(六):并行模式:卷积

一维并行卷积

  1. 由于卷积计算与顺序无关,所以可以使用并行解决,先声明kernel函数:
__global__ void convolution_1D_basic_kernel(float *N, float *M, float *P, 
	int Mask_Width, int Width);
  • N 为输入数组。
  • M 为卷积核。
  • P 为输出数组。
  • Mask_Width 为卷积核尺寸。
  • Width 为输入输出数组长度。
  1. 确定数据和线程的映射,由于输入数据是一维的,所以将网格设置为一维,每个线程负责一个输出元素的计算。
    int i = blockId.x * blockDim.x + threadId.x;

  2. 确定 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;
}

NVIDIA CUDA 高度并行处理器编程(六):并行模式:卷积
上面的 kernel 函数有两个不足:

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

常数存储器和高速缓存

先看看M的几个特征:

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

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

常数存储器与其他类型存储器关系如下图所示:
NVIDIA CUDA 高度并行处理器编程(六):并行模式:卷积
常数存储器与全局存储器类似,所有线程都能访问其中的变量,但其中的变量不能被修改。不同设备的常数存储器大小不同,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;
}

高速缓存的示意图:
NVIDIA CUDA 高度并行处理器编程(六):并行模式:卷积
与GPU上的其他存储器不同的是,高速缓存一般是透明的,处理器会自动保留程序最近或经常访问的变量到高速缓存中,并记录他们的原始DRAM地址,当保留的变量再次被使用时,硬件会从他们的地址判断已经在高速缓存保留了他们的值,并直接访问高速缓存,从而消除对DRAM的访问。

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

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

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

假定输入元素有16个,M大小为 5,分块大小为4,则卷积示例如下:
NVIDIA CUDA 高度并行处理器编程(六):并行模式:卷积
每个块加载需要的元素到共享存储器中,块内边缘的元素称为“边缘元素”(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;
}

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

  • 普通卷积:
    NVIDIA CUDA 高度并行处理器编程(六):并行模式:卷积
    内部线程块的访存次数为: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高速缓存,但只是可能,也有可能从全局存储器中加载光环元素,这是随机的。

文章出处登录后可见!

已经登录?立即刷新

共计人评分,平均

到目前为止还没有投票!成为第一位评论此文章。

(0)
心中带点小风骚的头像心中带点小风骚普通用户
上一篇 2022年6月13日 下午12:37
下一篇 2022年6月13日 下午12:39

相关推荐