参考链接： C++ fma()
本篇基于可分离卷积的性质，按照 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维度计算时造成内存瓶颈的原因。
