ECE408 Lecture 8 Tiled Convolution

这次回顾ECE408 Lecture 8，这次介绍了分块卷积（Tiled Convolution）。

课程主页：

https://wiki.illinois.edu/wiki/display/ECE408

搬运视频：

https://www.youtube.com/playlist?list=PL6RdenZrxrw-UKfRL5smPfFFpeqwN3Dsz

本讲目标

了解分块卷积算法
- 分块算法的一些复杂方面；
- 输出分块与输入分块；
- 三种不同风格的输入分块加载；

我们受到内存限制吗

对于1D情况，每个输出元素需要2*MASK_WIDTH负载（每个M和N）和2*MASK_WIDTH浮点运算，受到内存限制；
对于2D情况，每个输出元素需要2*MASK_WIDTH^2加载和2*MASK_WIDTH^2浮点运算，受到内存限制；

分块1D卷积

基本思想

思路是将每个输入按照某个方式划分为分块，然后分别计算，注意需要填充一些元素：

有多少重用？考虑分块1

分析：

元素2用于计算输出4 (1×)；
元素3用于计算输出4、5 (2×)；
元素4用于计算输出4、5、6 (3×)；
元素5用于计算输出4、5、6、7 (4×)；
元素6用于计算输出4、5、6、7 (4×)；
元素7用于计算输出5、6、7 (3×)；
元素8用于计算输出6、7 (2×)；
元素9用于计算输出7 (1×)；

三种分块策略

策略1：

说明：

1.block大小覆盖输出tile；
2.使用多个步骤加载输入tile；

策略2：

说明：

block大小覆盖输入tile；
一步加载输入tile；
计算输出时关闭一些线程；

策略3：

说明：

block大小覆盖输出tile；
只加载输入tile的“核心”；
从全局内存访问halo单元；

补充：

在后续代码中TILE_WIDTH = blockIdx.x。

策略1

加载左边的halo

代码：

__shared__ float N_ds[TILE_SIZE + MASK_WIDTH - 1];
int radius = Mask_Width / 2;
int halo_index_left = (blockIdx.x - 1) * blockDim.x + threadIdx.x;

if (threadIdx.x >= (blockDim.x - radius)) {
    if (halo_index_left < 0)
    	N_ds[threadIdx.x - (blockDim.x - radius)] = 0;
    else
    	N_ds[threadIdx.x - (blockDim.x - radius)] = N[halo_index_left];
}

例子：

a[b cd e]f

cd是tile中的元素；
b是左边的halo；
e是右边的halo；
TILE_SIZE = 2, blockDim.x = 2, radius = 1, MASK_WIDTH = 3；

代码说明：

我们需要载入的是b；
threadIdx.x >= (blockDim.x - radius)是指当前index处于d对应的位置，那么减去blockDim.x，即为b对应的位置；
halo_index_left是b对应的位置；
threadIdx.x - (blockDim.x - radius)给出N_ds对应的位置，例如此处d对应的为1-(2-1)=0；

加载内部的halo

代码：

if ((blockIdx.x * blockDim.x + threadIdx.x) < Width)
	N_ds[radius + threadIdx.x] = N[blockIdx.x * blockDim.x + threadIdx.x];
else
	N_ds[radius + threadIdx.x] = 0.0f;

代码说明：

radius + threadIdx.x给出元素在N_ds中的位置；

加载右边的halo

代码：

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

例子：

a[b cd e]f

cd是tile中的元素；
b是左边的halo；
e是右边的halo；
TILE_SIZE = 2, blockDim.x = 2, radius = 1, MASK_WIDTH = 3；

代码说明：

我们需要载入的是e；
threadIdx.x < radius是指当前index处于c对应的位置，那么加上blockDim.x，即为e对应的位置；
halo_index_right是e对应的位置；
radius + blockDim.x + threadIdx.x给出N_ds对应的位置，例如此处e对应的为1+2+0=3；

完整代码

__global__ void convolution_1D(float *N, float *P, float *M, int Width) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    __shared__ float N_ds[TILE_SIZE + MASK_WIDTH - 1];
    
    int radius = Mask_Width / 2;
    
    int halo_index_left = (blockIdx.x - 1) * blockDim.x + threadIdx.x;
    if (threadIdx.x >= (blockDim.x - radius)) {
        if (halo_index_left < 0)
        	N_ds[threadIdx.x - (blockDim.x - radius)] = 0.0f;
        else
        	N_ds[threadIdx.x - (blockDim.x - radius)] = N[halo_index_left];
    }
    
    if ((blockIdx.x * blockDim.x + threadIdx.x) < Width)
    	N_ds[radius + threadIdx.x] = N[blockIdx.x * blockDim.x + threadIdx.x];
    else
    	N_ds[radius + threadIdx.x] = 0.0f;
    
    int halo_index_right = (blockIdx.x + 1) * blockDim.x + threadIdx.x;
    if (threadIdx.x < radius) {
        if (halo_index_right >= Width)
        	N_ds[radius + blockDim.x + threadIdx.x] = 0;
        else
        	N_ds[radius + blockDim.x + threadIdx.x] = N[halo_index_right];
    }
    __syncthreads();
    
    float Pvalue = 0;
    for(int j = 0; j < MASK_WIDTH; j++) {
    	Pvalue += N_ds[threadIdx.x + j]*M[j];
    }
    P[i] = Pvalue;
}

策略1的另一种实现

加载输入数据——步骤1

代码：

if (0 <= start && Width > start) { // all threads
	N_ds[threadIdx.x] = N[start];
else
	N_ds[threadIdx.x] = 0.0f;

说明：

这一步读取N_ds中前TILE_SIZE个元素；

加载输入数据——步骤2

代码：

if (threadIdx.x < (MASK_WIDTH – 1)) { // some threads
    start += TILE_SIZE;
    if (Width > start) {
    	N_ds[threadIdx.x + TILE_SIZE] = N[start];
    } else {
    	N_ds[threadIdx.x + TILE_SIZE] = 0.0f;
    }
}

说明：

这一步读取N_ds中后MASK_WIDTH - 1个元素；

完整代码

__global__ void convolution_1D(float *N, float *P, float *M, int Width) {

    int i = blockIdx.x * blockDim.x + threadIdx.x;
    __shared__ float N_ds[TILE_SIZE + MASK_WIDTH - 1];
    int radius = MASK_WIDTH / 2;
    int start = i – radius;
    
    if (0 <= start && Width > start) { // all threads
    	N_ds[threadIdx.x] = N[start];
    else
    	N_ds[threadIdx.x] = 0.0f;
    	
    if (threadIdx.x < (MASK_WIDTH – 1)) { // some threads
        start += TILE_SIZE;
        if (Width > start) {
        	N_ds[threadIdx.x + TILE_SIZE] = N[start];
        else
        	N_ds[threadIdx.x + TILE_SIZE] = 0.0f;
    }
    
    __syncthreads();
    
    float Pvalue = 0.0f;
    for (int j = 0; MASK_WIDTH > j; j++) {
    	Pvalue += N_ds[threadIdx.x + j] * Mc[j];
    }
    P[i] = Pvalue;
}

策略3

说明：

和策略1类似，区别在于不加载左右两边的填充，这部分从全局内容中读取；
其思路为，对于输出的位置j，首先计算出输入需要考虑的起始位置N_start_point = i - radius；
然后遍历MASK_WIDTH个元素，考虑N_index = N_start_point + j；
如果N_index介于This_tile_start_point和Next_tile_start_point；
- 则从N_ds中读取，索引为threadIdx.x-radius+j；
否则从全局中读取；

完整代码

__global__ void convolution_1D(float *N, float *P, float *M, int Width) {

    int i = blockIdx.x * blockDim.x + threadIdx.x;
    __shared__ float N_ds[TILE_WIDTH];
    
    N_ds[threadIdx.x] = N[i];
    
    __syncthreads();
    
    int radius = MASK_WIDTH / 2;
    int This_tile_start_point = blockIdx.x * blockDim.x;
    int Next_tile_start_point = (blockIdx.x + 1) * blockDim.x;
    int N_start_point = i - radius;
    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-radius+j] * M[j];
            else
            	Pvalue += N[N_index] * M[j];
        }
    }
    P[i] = Pvalue;
}

策略2

2D情形

所有线程都参与将N加载到共享内存中;
然后线程子集计算P;

图示：

2D卷积的线程索引的输出索引

计算公式：

col_o = blockIdx.x * TILE_WIDTH + threadIdx.x;
row_o = blockIdx.y * TILE_WIDTH + threadIdx.y;

图示：

输入块需要大于输出块

设置网格/块尺寸

dim3 dimGrid(ceil(P.width/(1.0*TILE_WIDTH)), ceil(P.height/(1.0*TILE_WIDTH)), 1);
dim3 dimBlock(TILE_WIDTH + MASK_WIDTH - 1, TILE_WIDTH + MASK_WIDTH - 1, 1);

说明：

需要有足够的块来生成所有P元素；
需要有足够的线程来加载整个输入块；

从输出坐标转移到输入坐标

代码：

int tx = threadIdx.x;
int ty = threadIdx.y;

int row_o = blockIdx.y * TILE_WIDTH + ty;
int col_o = blockIdx.x * TILE_WIDTH + tx;

int row_i = row_o – (MASK_WIDTH / 2);
int col_i = col_o – (MASK_WIDTH / 2);

在N之外加载halo的线程应返回0.0

考虑边界

float Pvalue = 0.0f;
if((row_i >= 0) && (row_i < Width) &&
	(col_i >= 0) && (col_i < Width))
	tile[ty][tx] = N[row_i*Width + col_i];
else
	tile[ty][tx] = 0.0f;
__syncthreads (); // wait for tile

并非所有线程都计算输出

if(ty < TILE_WIDTH && tx <TILE_WIDTH){
    for(i = 0; i < MASK_WIDTH; i++)
    	for(j = 0; j < MASK_WIDTH; j++)
    		Pvalue += Mc[i][j] * tile[i+ty][j+tx];
    		
    if(row_o < Width && col_o < Width)
    	P[row_o * Width + col_o] = Pvalue;
    }
}

可选项

可以将1D策略3分块卷积扩展为2D策略3分块卷积；
- 每个输入块与其对应的输出块相匹配；
- 所有halo元素将从全局内存中加载；
- 注意控制内积计算中的分歧；