一维卷积的SASS实现

CUDA

对应代码和工程，见GitHub。

1、一维卷积计算过程

线性卷积:$\quad \quad y(n)=x(n)*h(n)$

设x(n)长度为N1，h(n)长度为N2 ， 则y(n)长为 N=N1+N2-1
为方便表示，在序列x(n)后添N2-1个0 ，使x(n)的长度变为 N

卷积公式 ： $y(n)=\sum\limits_{i=0}^{N-1} x(i) \times h(n-i) , \quad 0 \le n \le N-1$

可用矩阵表示为：

一个简单的计算过程示意：

x(n)=[4,3,2,1]
h(x)=[3,2,1]
y(n)=x(n)*h(n)

x  |        4  3  2  1        | 
   |                          | 
y0 |  1  2  3                 | = 4*3              = 12
y1 |     1  2  3              | = 4*2 + 3*3        = 17
y2 |        1  2  3           | = 4*1 + 3*2 + 2*3  = 16
y3 |           1  2  3        | = 3*1 + 2*2 + 1*3  = 10
y4 |              1  2  3     | = 2*1 + 1*2        = 4
y5 |                 1  2  3  | = 1*1              = 1

为保证局部性：

在实现时，卷积核H的数据先放到片上
（H点数较小时直接放到寄存器，点数稍大可放在共享内存，过大时分批从片外内存读取）

实际的原始数据X，点数较大，分块（blocking）处理。

Y的计算不按照上例中的一次性计算 $y2 = 4*1 + 3*2 + 2*3 = 16$ 。
而是 多次乘加 得到结果 ：

//Y初始化为0
regX = shareadX[0];      //共享内存读取x
Y0  += regX*H1;
Y1  += regX*H2;
Y2  += regX*H3;
regX = shareadX[1];      
Y1  += regX*H1;   
Y2  += regX*H2;
Y3  += regX*H3;
...
//计算完成后把Y 在共享内存reshape 然后放回片外

上例每次读一个X。实现时，X批量读入，计算掩盖访存延迟。

2、两个实现：16点卷积与1024点卷积

算法分析

$算法计算强度 = \frac{计算量}{数据量}=\frac{2*N*N2}{2N}=N2$
(即一维卷积计算算法强度与卷积核长度直接相关： 16点卷积计算强度为16， 1024点卷积计算强度为1024)

$硬件计算强度 = \frac{浮点计算能力}{浮点带宽}=\frac{4160 GFLOPS}{35G Float Data}=118.8$

所以16点卷积计算 ：受限于带宽 $\quad$（16 < 116.8 访存密集）
1024点卷积计算 ： 受限于核心计算能力（1024 > 116.8 计算密集）

所有算法都可以从这一角度，划分为上述两个类型之一。
计算密集的算法，可以针对硬件架构进行优化。

16点卷积

一维卷积核H：16 float （直接拷贝到寄存器）
输入数据量X：1024000 float
输出数据量Y：1024000 + 16 float （最后16个Y为了方便，没有进行计算）
线程结构 : 32个thread一个blcok，每个thread计算16点Y，一共需要2000个block。
片上存储 ： 分配544 float * 4 bytes = 2176 bytes 大小的共享内存空间用于存储。

存储结构

Shared存储

每个block负责32*16=512点Y计算
32个threads 对显存合并访问 读取X
每个thread进行4次LDG.128访问 + 一次LDG访问（32bit，末尾的32个X，只用到1/2的数据）

Shared访存
整体来看，一个线程需要访问8次共享内存（128bit形式）
eg：
第一次，所有线程访问第一行，第二次第二行，第三次第三行，第四次第四行
第五次时，需要向左移一列访问，接下来再访问三行

（注：图看不清的话，可以点击放大）

以thread 0为例 ，蓝色框是其访存范围，红色线条是访存顺序。
黑框是 thread 1 的访存范围，前后两个线程共享16个 float 数据

这样的LD/ST方法，利用了完整的共享内存带宽。

算法结构

计算示意（结合 1、一维卷积计算过程 理解）

完全避免寄存器冲突的办法

SASS代码：
conv16_final_32blockOcup.sass

1024点卷积

计算的核心循环与16点卷积一致，32点X与16点H参与计算。
区别在于下一次循环还要使用32个当前计算X中的的后16个；下一次计算使用的H是新的16个。

理论分析：由于一个小循环算完16点Y，数据都在片上，所以再套上大循环，对X、Y做双缓冲不会有速度提升
实际处理：H加双缓冲无效，寄存器双缓冲预取共享内存中的X，计算掩盖延迟。

一维卷积核H：1024 float
输入数据量X：2097152 float
输出数据量Y：2097152+1024 float （最后1024个Y为了方便，没有进行计算）
线程结构：128个thread一个blcok，每个thread计算16点Y，一共需要1024个block。

存储结构

X数据存储与16点类似，一个block存128线程*(4float*6次加载)=3072个单精度浮点。
还需要1024*4byes存放H。
共分配4096*4bytes的共享内存空间

计算结构

block内计算示意

block间计算示意

注:Y存回片外时通过共享内存reshape，有bank conflict，只用到了共享内存1/4的完整带宽。
但实测对kernel执行时间基本没有影响（占比太小）。

SASS代码：
conv1024_nobuffer.sass

3、实验结果

16点卷积

每个线程使用60个寄存器（优化后），恰好达到1个SM上32个Block的占用上限。

实测，只要一个SM上达到6个Block的占用，就可占满带宽。
6占用kernel执行时间和32个占用的时间一致。
占用多了，带宽跟不上，还是需要等待。

数据量(32*16*2000)*2 = 2*1024000FP32 = 2*1024000*4bytes/1024/1024 = 7.8125MB

下图可知，kernel平均计算时间为56us

带宽占用： 7.8125MB / 0.056ms = 139.51GB/s
与理论分析一致

1024点卷积

Visual Profiler （nvvp）分析：

带宽未占满

计算指令占比接近90%

SM占用6

实测占用 6、4、2个block的情况，kernel执行时间一致。
由于计算能力限制，占用多了也算不过来。

kernel执行时间1.4ms

数据生成及结果验证

原始数据X、卷积核H和结果Y的长度分别为N、M、P。
卷积核与原始数据用伪随机数填充。
CPU计算的卷积结果放在数据块T中。
GPU计算的卷积结果拷贝到数据块Y中。

打印前1024个T与Y的误差
Y与T差的绝对值大于1则打印出来（计算精度所在位置随数据大小变化）

4、小结

常用的优化方法：
调整算法结构（匹配硬件特性），kernel间整合（减少访存），kernel优化（数值计算过程与汇编）

16点卷积和1024点卷积恰好可划分在访存密集和计算密集这两个类型中。
是比较好的练手项目。

访存密集型的算法不需要用汇编优化，保证基本的合并访存即可。
但一些算法结构复杂，需要仔细分析，利用好片上和片外内存的带宽。

1024点卷积kernel应该达到90%+的峰值算力，未找到受限原因。
理论上1024点卷积不需要嵌套循环，不用对X做双缓冲预取，因为当前计算的Y与下一大组X没有数据依赖关系。

按现在的理解，在SM上block占用低于一个阈值的情况下，双缓冲可以带来一定的计算能力提升。
在占用足够的情况下，线程级并行、线程间的切换，其实和双缓冲是等效的。
（2017补充）

实际工程中常用FFT在频域进行快速卷积。

5、后续

二维卷积
FFT