1.3.4 二维单通道图像离散卷积_并行编程方法与优化实践-QQ阅读男生中文历史网

书名：并行编程方法与优化实践
作者名：刘文志
本章字数：1793字
更新时间：2025-02-21 21:55:47

1.3.4　二维单通道图像离散卷积

在现在流行的深度神经网络做图像识别的应用中，有一种关键的运算是计算二维图像卷积。

这里我们只考虑单通道和单精度浮点类型的二维离散卷积，输入是一个由浮点的二维数组保存的图像，还有一个m×n的模板kernel，m和n一般相等，常用的大小有3、5、7…等。设图像大小是M×N，则浮点计算次数是2×（M-m+1）×（N-n+1）×m×n。如果m和n比较大，则这是一个典型的计算访存比很高的计算过程。经过精心调优，应该可以充分利用CPU的浮点运算性能进行优化。

这里主要的难点是如何利用SIMD数据并行计算，并做好指令流水调度。对于4×4或者8×8这样的kernel，可以直接对kernel做SIMD向量化，但这种方法限制太多，对于非4倍大小的kernel会有冗余和截断，指令流水长度也会受限，所以需要换一个并行的维度。

一般情况下，输入图像和输出图像的大小不会太小，可以在输出图像上直接做寄存器分块。SSE指令的向量寄存器长度是4个float，AVX或FMA指令是8个float。前面章节讨论过，SSE和AVX的乘加延迟总和都是8个周期，FMA是10个周期。综上，对于SSE指令，可以构造8×4的寄存器分块；AVX指令可以构造8×8的分块；FMA指令可以构造10×8的分块。

下面我们就以AVX指令为例来说明如何计算一个寄存器分块。这里以3×3的kernel为例。

X64指令集有16个向量寄存器，使用8个向量寄存器保存当前计算的8×8个输出结果，首先初始化为0；从左到右每次从kernel中的取出一列值，这里一次取3个，各广播到一个向量寄存器中；从输入矩阵中每次读取一个向量，和前面广播的3个kernel值向量做乘法，并累加到输出矩阵当前的寄存器分块中；再从下一行同样的偏移位置读一个向量，做相同的操作，直到第8+3-1=10行为止。图1-6展示了计算到kernel第二列的值的乘累加位置，其中输入矩阵里深灰色的元素表示3×3的kernel需要额外读入的2行2列的数据（3-1=2）。

图1-6　单通道二维卷积的向量寄存器分块算法

一个8×8的寄存器分块计算完毕，就可以写回到输出矩阵，并可以开始计算右边紧挨着的一个新的寄存器分块。不断推进这个过程，直到输出矩阵所有位置被计算完毕。这个计算过程对cache的重用已经非常好了，不需要再做专门的局部化处理。对于边界不满足8×8的分块，可以用宏或者代码生成器生成小于8×8的寄存器分块计算过程。

下面的代码片段展示了一个8×8寄存器分块和3×3大小的kernel计算一列kernel的过程。

ymm13 = _mm256_broadcast_ss(flt + flt_w * 0 + n); // 从kernel中
ymm14 = _mm256_broadcast_ss(flt + flt_w * 1 + n); // 读入一列元素
ymm15 = _mm256_broadcast_ss(flt + flt_w * 2 + n); // 广播到3个向量
ymm10 = _mm256_loadu_ps(src + src_w * 0);  // 读取输入数组第0个向量
ymm11 = _mm256_mul_ps(ymm13, ymm10);       // 第0个向量只跟kernel最上面元素做乘加
ymm0 = _mm256_add_ps(ymm11, ymm0);          
ymm10 = _mm256_loadu_ps(src + src_w * 1); 　// 读取输入数组第1个向量
ymm11 = _mm256_mul_ps(ymm13, ymm10);       // 第1个向量跟kernel前两个元素做乘加
ymm1 = _mm256_add_ps(ymm11, ymm1);         
ymm11 = _mm256_mul_ps(ymm14, ymm10);
ymm0 = _mm256_add_ps(ymm11, ymm0);
ymm10 = _mm256_loadu_ps(src + src_w * 2);  // 读取输入数组第2个向量
ymm11 = _mm256_mul_ps(ymm13, ymm10);       // 从第2个向量开始跟全部3个kernel元素做乘加
ymm2 = _mm256_add_ps(ymm11, ymm2);         
ymm11 = _mm256_mul_ps(ymm14, ymm10);
ymm1 = _mm256_add_ps(ymm11, ymm1);
ymm11 = _mm256_mul_ps(ymm15, ymm10);
ymm0 = _mm256_add_ps(ymm11, ymm0);
ymm10 = _mm256_loadu_ps(src + src_w * 3);
ymm11 = _mm256_mul_ps(ymm13, ymm10);
ymm3 = _mm256_add_ps(ymm11, ymm3);
ymm11 = _mm256_mul_ps(ymm14, ymm10);
ymm2 = _mm256_add_ps(ymm11, ymm2);
ymm11 = _mm256_mul_ps(ymm15, ymm10);
ymm1 = _mm256_add_ps(ymm11, ymm1);
ymm10 = _mm256_loadu_ps(src + src_w * 4);
ymm11 = _mm256_mul_ps(ymm13, ymm10);
ymm4 = _mm256_add_ps(ymm11, ymm4);
ymm11 = _mm256_mul_ps(ymm14, ymm10);
ymm3 = _mm256_add_ps(ymm11, ymm3);
ymm11 = _mm256_mul_ps(ymm15, ymm10);
ymm2 = _mm256_add_ps(ymm11, ymm2);
ymm10 = _mm256_loadu_ps(src + src_w * 5);
ymm11 = _mm256_mul_ps(ymm13, ymm10);
ymm5 = _mm256_add_ps(ymm11, ymm5);
ymm11 = _mm256_mul_ps(ymm14, ymm10);
ymm4 = _mm256_add_ps(ymm11, ymm4);
ymm11 = _mm256_mul_ps(ymm15, ymm10);
ymm3 = _mm256_add_ps(ymm11, ymm3);
ymm10 = _mm256_loadu_ps(src + src_w * 6);
ymm11 = _mm256_mul_ps(ymm13, ymm10);
ymm6 = _mm256_add_ps(ymm11, ymm6);
ymm11 = _mm256_mul_ps(ymm14, ymm10);
ymm5 = _mm256_add_ps(ymm11, ymm5);
ymm11 = _mm256_mul_ps(ymm15, ymm10);
ymm4 = _mm256_add_ps(ymm11, ymm4);
ymm10 = _mm256_loadu_ps(src + src_w * 7);
ymm11 = _mm256_mul_ps(ymm13, ymm10);
ymm7 = _mm256_add_ps(ymm11, ymm7);
ymm11 = _mm256_mul_ps(ymm14, ymm10);
ymm6 = _mm256_add_ps(ymm11, ymm6);
ymm11 = _mm256_mul_ps(ymm15, ymm10);
ymm5 = _mm256_add_ps(ymm11, ymm5);
ymm10 = _mm256_loadu_ps(src + src_w * 8);
ymm11 = _mm256_mul_ps(ymm14, ymm10);      // 倒数第2个向量只跟kernel最后两个元素做乘加
ymm7 = _mm256_add_ps(ymm11, ymm7);         
ymm11 = _mm256_mul_ps(ymm15, ymm10);
ymm6 = _mm256_add_ps(ymm11, ymm6);
ymm10 = _mm256_loadu_ps(src + src_w * 9); // 读取输入数组最后一个向量
ymm11 = _mm256_mul_ps(ymm15, ymm10);      // 最后一个向量只跟kernel最后一个元素做乘加
ymm7 = _mm256_add_ps(ymm11, ymm7);

下面来分析一下这个算法的cache效率。对于一个8×8的分块，总的浮点计算量是8×8×3×3×2（乘加各算一次浮点操作）=1152；访问cache的数量（这里只以读取float个数计）是（8×（8+3-1）+3）×3=243；经过测试，支持AVX指令的SNB架构CPU的cache带宽是一个周期8个float，同时SNB架构一个周期可以发射8个单精度浮点乘法和单精度浮点加法，即一个周期内cache和浮点计算的吞吐比例是1∶2；但因为大部分cache访问并非对齐到向量长度，我们保守估计cache带宽会损失一半，即降到1∶4，但仍然高于我们前面计算过的243∶1152。综上，cache延迟可以被浮点计算掩盖，带宽足够满足浮点计算的峰值性能。

但由于cache访问量的减少会导致流水线有部分依赖，这个算法仍不能达到浮点峰值性能，表1-10展示了不同尺寸的kernel实测的浮点性能（在Haswell架构下使用FMA指令），以及与浮点峰值的比较。可以看出，kernel尺寸越大，计算访存比越高，峰值性能也越高，最高可以达到大约66%的浮点峰值性能。

表1-10　二维离散单通道卷积的性能测试