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

本篇基于可分离卷积的性质，按照 X Y Z的顺序，依次计算每个维度的一维卷积。

代码实现

因为是按照X Y Z的计算顺序，因此只能够在计算X维度的卷积时，复用之前实现的一维卷积计算函数。Y维度的计算是将一个Z平面上的二维数据中每行与卷积核中一个点相乘，并将31个点的卷积核计算出的结果累加至一行，更新到中间缓存的目标位置。Z维度的计算是将一个Z平面的二维数据和卷积核中的一个点相乘，并将31个点的卷积核计算出的结果累加至一个二维平面，更新到结果的目标位置。这里对Y 和 Z维度的计算都是通过编译器ICC实现向量化。
代码实现如下：

void GaussSmoothCPU3DBase(float* pSrc, int iDim[3], float* pKernel, int kernelSize[3], float* pDst, float* pBuffer)
	{
		//结果正确
		int iSliceSize = iDim[1] * iDim[0];
		int nCenter = kernelSize[0] / 2;
		for (int z = 0; z < (iDim[2]); z++)
		{
			float* pSrcSlice = pSrc + z * iSliceSize;
			float* pBuffSlice = pBuffer + z * iSliceSize;
			memset(pBuffSlice, 0, iSliceSize * sizeof(float));
			for (int y = 0; y < iDim[1]; y++)
			{
				float* pSrcLine = pSrcSlice + y * iDim[0];
				float* pDstLine = pBuffSlice + y * iDim[0];
				Conv1D_Opt_Cmb(pSrcLine, iDim[0], pKernel, kernelSize[0], pDstLine);
			}
			for (int y = 0; y < (iDim[1] - kernelSize[0] + 1); y++)
			{
				float* pDstLine = pSrcSlice + (y + nCenter) * iDim[0];
				memset(pDstLine, 0, iDim[0] * sizeof(float));
				for (int kx = 0; kx < kernelSize[0]; kx++)
				{
					float* pSrcLine = pBuffSlice + (y + kx) * iDim[0];
#pragma omp simd
					for (int i = 0; i < iDim[0]; i++)
					{
						pDstLine[i] += pSrcLine[i] * pKernel[kx];
					}
				}
			}
		}

		for (int z = 0; z < (iDim[2] - kernelSize[0] + 1); z++)
		{
			float* pDstSlice = pDst + (z + nCenter) * iSliceSize;
			memset(pDstSlice, 0, iSliceSize * sizeof(float));
			for (int kx = 0; kx < kernelSize[0]; kx++)
			{
				float* pSrcSlice = pSrc + (z + kx) * iSliceSize;
#pragma omp simd
				for (int i = 0; i < iSliceSize; ++i)
				{
					pDstSlice[i] += pKernel[kx] * pSrcSlice[i];
				}
			}
		}
	}

执行时间

GaussSmoothCPU3DBase cost Time(ms) 1010.5

多线程并行化

计算逻辑不变，基于OpemMP实现多线程并行化。因此线程数设置8与16分别测试其执行耗时情况，并选择速度最快的版本。
代码实现如下：

void GaussSmoothCPU3DBase(float* pSrc, int iDim[3], float* pKernel, int kernelSize[3], float* pDst, float* pBuffer)
	{
		//结果正确
		int iSliceSize = iDim[1] * iDim[0];
		int nCenter = kernelSize[0] / 2;
#pragma omp parallel for num_threads(16) schedule(dynamic)
		for (int z = 0; z < (iDim[2]); z++)
		{
			float* pSrcSlice = pSrc + z * iSliceSize;
			float* pBuffSlice = pBuffer + z * iSliceSize;
			memset(pBuffSlice, 0, iSliceSize * sizeof(float));
			for (int y = 0; y < iDim[1]; y++)
			{
				float* pSrcLine = pSrcSlice + y * iDim[0];
				float* pDstLine = pBuffSlice + y * iDim[0];
				Conv1D_Opt_Cmb(pSrcLine, iDim[0], pKernel, kernelSize[0], pDstLine);
			}
			for (int y = 0; y < (iDim[1] - kernelSize[0] + 1); y++)
			{
				float* pDstLine = pSrcSlice + (y + nCenter) * iDim[0];
				memset(pDstLine, 0, iDim[0] * sizeof(float));
				for (int kx = 0; kx < kernelSize[0]; kx++)
				{
					float* pSrcLine = pBuffSlice + (y + kx) * iDim[0];
#pragma omp simd aligned(pSrcLine, pDstLine)
					for (int i = 0; i < iDim[0]; i++)
					{
						pDstLine[i] += pSrcLine[i] * pKernel[kx];
					}
				}
			}
		}

#pragma omp parallel for num_threads(8) schedule(dynamic)
		for (int z = 0; z < (iDim[2] - kernelSize[0] + 1); z++)
		{
			float* pDstSlice = pDst + (z + nCenter) * iSliceSize;
			memset(pDstSlice, 0, iSliceSize * sizeof(float));
			for (int kx = 0; kx < kernelSize[0]; kx++)
			{
				float* pSrcSlice = pSrc + (z + kx) * iSliceSize;
#pragma omp simd
				for (int i = 0; i < iSliceSize; ++i)
				{
					pDstSlice[i] += pKernel[kx] * pSrcSlice[i];
				}
			}
		}
	}

执行时间

GaussSmoothCPU3DBase cost Time(ms) 218.4

VTune分析性能问题

指令执行情况如下：

其中，为了执行结果稳定，重复调用函数30次。
线程占用率如下：

可以看到线程大部分时间还是在做有用工作的。

计算X Y维度卷积的性能状态：

整体上没有突出的性能问题。

热点语句是Y维度的FMA运算。

查看反汇编，其中broadcast指令CPI比较高。关于指令的解释如下：
broadcast指令CPI理论值为1。这里抓取的CPI为1.4，略低于理论值。这里执行这个指令的原因是将一维卷积核的一个点展开成一个向量，但是根据反汇编中broadcast指令的执行次数和fmadd是一个数量级的，推断ICC在这里应该是内层循环每次迭代都做了一次broadcast，但显然有更高效的做法：只将卷积核展开一次，并保存在寄存器中复用，效率会更高。
计算Z维度卷积的执行状态：

整体上内存访问为主要的性能问题。

在Z维度FMA运算时CPI最高。

查看反汇编指令，其中加载内存与broadcast指令CPI高。
计算X Y维度和计算Z维度过程类似，为什么CPI差距会这么大呢？在内存访问上，XY维度的计算每次迭代，所有内存读写操作都2 * 432 * 432 的内存块内进行；Z维度的计算每次迭代，内存读写跨度在332 * 432 *432 ，因此造成内存访问成为性能瓶颈。由Z维度计算的逻辑设计所限制，每次迭代必须要访问这么多的内存数据，因此目前针对内存访问的问题，还没找到很好的解决方案。

总结

本文按照 X Y Z的维度顺序，实现了3D高斯卷积的计算，同时基于OpenMP技术，实现了多线程并行化。同时分析了Z维度计算时造成内存瓶颈的原因。