CUDA 3D convolution

overview

这是ECE408的一个作业，目标是实现3d卷积.

测试的时候使用link这个脚本对测试数据测试

课程给的测试环境是GTX1080.我用自己的RTX2070会出bug.而实验室服务器的titan xp是可以的

这个问题分为两种写法，目前只实现了一种相对好理解但效率低的写法。我认为效率低的原因是每一个线程块启动了过多的线程，而一个块内的有些线程只是为了写共享内存，而没有进入实际的卷积计算，这可能会使线程块内出现线程束分化。确定性的是如果线程束分化，对性能肯定是有影响，不过线程块内的线程束之间如果出现分化是否有这个影响不确定，应该进行量化测速。

第一种写法

线程块、网格配置

首先，我们一个线程块内的线程做的事情是：将tileCache内的所有元素写入到相应的共享内存内，我们先不做数据规约，让每一个线程写共享内存的一个具体位置。因此线程块的维度可以知道就是（输入数据划分的tile尺寸+2*卷积核半径），此时确定了线程块的维度以后，一般情况会直接使用线程块维度来计算网格维度。但这里不同！因为真正在输出相应位置写数据的线程维度是TILE尺寸的，其他的线程只是为了写共享内存。因此，网格维度需要按照TILE分块来计算。

  dim3 block(CACHE_SIZE, CACHE_SIZE, CACHE_SIZE);
  dim3 grid(ceil(x_size/(1.0*TILE_SIZE)), ceil(y_size/(1.0*TILE_SIZE)), ceil(z_size/(1.0*TILE_SIZE)));

卷积核存储

因为CNN权重共享，卷积核在一段计算过程内是不变的，因此适合写入到常量内存中，注意常量内存是全生命周期的，任何核函数都可以访问。

__constant__ float kernel[KERNEL_SIZE][KERNEL_SIZE][KERNEL_SIZE];

核函数

接下来是核函数的计算
这一块分为两个部分，第一步是将输入数据中tileCache尺寸的数据写入到共享内存中，第二步是进行卷积计算.

写共享内存

我们一个线程块内的线程执行这个tile对应的写入，而且是一对一写入。因此我们要做的是将线程的id映射到tile的坐标，同时需要考虑怎样设定tile的坐标，可以使得特定的坐标位置写入的值是0，对应超出输入数据的位置。
一种方法是，将输入数据中超出边界的位置设定为-1.这样做的目的是，可以在写入共享内存的时候对-1坐标的位置写入0，而其他位置又不超最大值边界的写入该写入的值。我们以二维的卷积为例。
对于第一个tile，线程id和写入位置的坐标之间的映射关系如下

那么推广至第i个tile，线程id和写入位置的坐标之间的映射关系如下

那么对于最大值超界的情况可以通过全局输入数据的尺寸值进行判断

    // global thread index of all blocks.
    int out_x = tx+bx; 
    int out_y = ty+by;
    int out_z = tz+bz;
    // map threads to position of input
    int cache_x = out_x -1;
    int cache_y = out_y -1;
    int cache_z = out_z -1;
	
    if (cache_x>=0 && cache_y>=0 && cache_z>=0 && cache_x < x_size && cache_y < y_size && cache_z < z_size) {
        tileCache[tz][ty][tx] = input[cache_z*(y_size*x_size)+cache_y*x_size+cache_x];
    } else {
        tileCache[tz][ty][tx] = 0.0f;
    }
    __syncthreads();

卷积计算

卷积计算这块，只需要将对应位置相乘并累计和。对于卷积核的坐标很好计算，可以通过三个循环来迭代卷积核的位置，那么在迭代卷积核的过程中，如何找到对应的输入数据应取得的坐标？

因为最终这一个block只完成一个tile的卷积计算，因此可以只取[0,tilesize)的线程来完成写操作
因为在这一个tilecache中，任何一个位置的元素都可以通过(bx+tx, by+ty, bz+tz)这个坐标获取
例如:

可以写出如下的访问方法

if (tz < TILE_SIZE && ty < TILE_SIZE && tx < TILE_SIZE) {
  for (int i =0;i<KERNEL_SIZE;i++) {
    for (int j =0; j<KERNEL_SIZE;j++) {
      for (int k=0;k<KERNEL_SIZE;k++) {
        p_value += kernel[i][j][k]*tileCache[tz+i][ty+j][tx+k];
      }
    }
  }
  //write out
}

最后再将结果写入到output.因为在output中坐标是全局的，因此需要注意，我们设定的线程数是超过整体input/output尺寸的，所以对于全局索引来说，只要计算出来超过了x_size,y_size,z_size都不用进行操作了。

write out

  if (out_z >= 0 && out_z < z_size && out_y >= 0 && out_y < y_size && out_x >= 0 && out_x < x_size) 
  //if (out_z < z_size && out_y < y_size && out_x < x_size)
      output[out_z*y_size*x_size + out_y*x_size + out_x] = p_value;
    

final

心得体会：

要注意我们如何设定我们的线程块维度和网格维度，这和线程具体要操作什么数据是相关的。

具体要将输入数据怎样进行分组，缓存、计算的时候和最终写入数据的时候分别怎样将对应的线程和内存位置对应起来。

完整代码