C++性能优化系列——3D高斯核卷积计算(五)2D卷积分离计算

高斯卷积核具有可分离的性质，因此可以通过以下方法计算二维高斯卷积：构造一个一维高斯卷积核，将原始二维矩阵分别以行主序与列主序，与一维卷积核做卷积计算，得到的结果就是目标二维高斯卷积的结果。本篇按照上述描述的思路实现了可分离的二维高斯卷积计算，并在此基础上对计算的过程分解与重构，挖掘实现的并行性。

基线版二维高斯卷积

为了让运行时间更加稳定，增加函数的执行次数至1000

#define CONV2DREPT 1000

这里对代码实现功能做一下简单的说明。原始矩阵维度432 * 432，维度定义为x和y。高斯卷积核维度31 * 31。基于《 性能优化系列(CPU)——3D高斯核卷积计算(三)FMA向量化计算一维卷积 》中向量化实现的一维卷积的计算函数，先以x方向，按照行主序计算一维高斯卷积；在计算y方向的一维卷积时，由于y方向数据在内存中不连续，因此无法复用之前的函数，这里按照卷积定义式，取y方向上对应的元素，计算点积并保存。此外，计算过程中需要使用到一个与原始数据尺寸相同的内存块保存临时数据。
代码实现

void Conv2D_Naive(float* pSrc, int iDim[2], float* pKernel, int iKernelSize, float* pTmp)
	{

		for (int y = 0; y < iDim[1]; y++)
		{
			float* pSrcLine = pSrc + y * iDim[0];
			float* pDstLine = pTmp + y * iDim[0];
			Conv1D_opt1(pSrcLine, iDim[0], pKernel, iKernelSize, pDstLine);
		}
		
		int iHalfKernel = iKernelSize / 2;
		memset(pSrc, 0, sizeof(float) * iDim[0] * iDim[1]);
		for (int x = 0; x < iDim[0]; x++)
		{
			for (int y = 0; y < iDim[1] - iKernelSize + 1; y++)
			{
				float fTemp = 0.0f;
				for (int kx = 0, ny = y; kx < iKernelSize; kx++, ny++)
				{
					fTemp += pTmp[ny * iDim[0] + x] * pKernel[kx];
				}
				float* pD = pSrc + (y + iHalfKernel) * iDim[0] + x;
				pSrc[(y + iHalfKernel)*iDim[0] + x] = fTemp;
			}
		}

	}

执行时间

TestConv2D(Conv2D_Naive) cost total Time(ms) 1456
TestConv2D cost Time(ms) 1.456

VTune分析基线版性能瓶颈

可以看到程序整体上内存访问与向量化都存在问题。

指令执行了90亿次

内存访问方面：L1 Cache Miss 为性能瓶颈。因为做y方向卷积时，以列主序访问不连续的内存，内存访问模式导致目标地址不在缓存中，发生miss。
VTune关于L1 Bound的描述：
This metric shows how often machine was stalled without missing the L1 data cache. The L1 cache typically has the shortest latency. However, in certain cases like loads blocked on older stores, a load might suffer a high latency even though it is being satisfied by the L1. Note that this metric value may be highlighted due to DTLB Overhead or Cycles of 1 Port Utilized issues.

因为访问不连续的内存，向量化程度也受到影响。
VTune关于Cycles of 1 Post Utilized 的描述：
This metric represents cycles fraction where the CPU executed total of 1 uop per cycle on all execution ports. This can be due to heavy data-dependency among software instructions, or oversubscribing a particular hardware resource. In some other cases with high 1_Port_Utilized and L1 Bound, this metric can point to L1 data-cache latency bottleneck that may not necessarily manifest with complete execution starvation (due to the short L1 latency e.g. walking a linked list) - looking at the assembly can be helpful. Note that this metric value may be highlighted due to L1 Bound issue.

分析具体热点

根据之前的描述，在做y方向的卷积时，点积运算是最大的热点。
查看点积运算这行的反汇编

这里ICC编译器对程序通过gather和fma指令，实现了向量化，其中fma指令执行了17亿次。本人的理解是：ICC根据代码实现逻辑，分析出内存访问模式，使用gather指令将y方向上的数据合并成一个向量寄存器ymm的长度，进行向量化计算。
关于指令vgatherdps的解释如下

优化版二维卷积

优化版代码实现逻辑：x方向上一维卷积过程不变，y方向上对计算过程分解，计算过程是将一行数据与对应卷积核的一个点计算并累加更新到目标位置。
代码实现

void Conv2D_Opt(float* pSrc, int iDim[2], float* pKernel, int iKernelSize, float* pTmp)
	{
		for (int y = 0; y < iDim[1]; y++)
		{
			float* pSrcLine = pSrc + y * iDim[0];
			float* pDstLine = pTmp + y * iDim[0];
			Conv1D_opt1(pSrcLine, iDim[0], pKernel, iKernelSize, pDstLine);
		}

		int iHalfKernel = iKernelSize / 2;
		memset(pSrc, 0, sizeof(float) * iDim[0] * iDim[1]);
		for (int y = 0; y < (iDim[1] - iKernelSize + 1); y++)
		{
			float* pDstLine = pSrc + (y + iHalfKernel) * iDim[0];
			for (int kx = 0; kx < iKernelSize; kx++)
			{
				float* pSrcLine = pTmp + (y + kx) * iDim[0];
#pragma omp simd aligned(pSrcLine, pDstLine)
				for (int i = 0; i < iDim[0]; i++)
				{
					pDstLine[i] += pSrcLine[i] * pKernel[kx];
				}
			}
		}

	}

执行速度提高了一倍多

TestConv2D(Conv2D_Opt) cost total Time(ms) 693
TestConv2D cost Time(ms) 0.693

VTune分析优化版性能瓶颈

从整体执行情况看，内存访问的问题已经不存在了。

总指令减少了20亿次，CPI从0.646降低到0.408

和基线版本对比，热点代码（y方向的计算）CPI下降到了0.391.
查看反汇编

这里vfmadd213ps指令一共只执行了不到9亿次，CPI大概在0.4左右，
对比基线版本vfmadd231ps指令执行17亿次，CPI大概1.0。可以看到改变内存访问方式后，编译器对fma的指令做出了调整，向量化指令fma对内存数据的依赖更低，在执行数量和执行速度上都得到了优化。

总结

本文基于卷积的定义式，对可分离的高斯卷积进行实现，并分析了存在的性能问题。在此基础上对计算过程进行重组，最终获得了两倍的加速倍数。