移动端GPU矩阵乘优化
移动端GPU目前主要有3家供应商， Qualcomm的Adreno系列，Arm的mali系列和Imagination的PowerVR GPU。主流开发语言包括OpenCL、OpenGL以及Vulkan，本文不对各个语言的应用进行讨论，仅以OpenCL为例。不同设备的体系结构差异很大，即使同一供应商的设备，也存在多个系列，因此优化策略也有不同。本文仅介绍纹理内存在Adreno和Mali设备上所带来的性能提升。
在GPU上的内存一般分为两种，一种是普通内存，OpenCL中叫做buffer内存，一种是纹理内存(Texture内存)， OpenCL中叫做Image内存。纹理内存和普通的buffer内存是通过不同的硬件单元来加载和写入的。除此之外，移动端设备上，不同的GPU架构下，对纹理内存和buffer内存的访问都存在差异；例如高通设备上，纹理内存的读可以使用L1 Cache，Mali设备上虽然无此差异，但是Mali确在最近几代GPU架构的迭代中不断的增强纹理内存的访存能能力。
本文主要从以下几个方面展开：
测试环境介绍
基础优化版本
Adreno设备的Texture方案
Mali设备的Texture方案
Mali(ValHall)的FMA方案
其他优化方案简介
测试环境及指标介绍
本文测试设备使用Qualcomm 865芯片和MTK的天玑1000芯片，对应GPU为Adreno 650 及Mali的G77 MP9，峰值数据是实际测试乘加计算的吞吐，非理论峰值.
| GPU | FP16峰值(FP32峰值) |
|---|---|
| Adreno 650 | 1407(764) |
| Mali G77(MP9) | 881(445) |
测试数据
矩阵维度为：A的维度为M x K, B的维度为K x N, C维度为M x N， 其中(M=N=K=1024); 测试数据采用float16 随机数进行测试。
指标计算
评价指标采用GFLOPS , 计算方式为(M * N * K) * 2 / 1024 / 1024 / 1024 / computeTime(s);
使用OpenCL的event机制对计算kernel计时，计时之前会循环调用10次该kernel进行warm up；随后对该kernel循环调用20次，取平均值作为执行时间。
基础优化版本
直接实现版本
首先按照矩阵乘法的计算公式，实现最简单版本作为base，如下如图所示，A矩阵的第一行乘以B矩阵的第一列得到C矩阵对应行列的一个元素。
代码实现如下：
// global_work_size = {N, M}#pragma OPENCL EXTENSION cl_khr_fp16 : enable__kernel void gemm_opt(__global half* A, __global half* B, __global half* C, int M, int N, int K){int idx = get_global_id(0); // 0--(N-1)int idy = get_global_id(1); // 0--(M-1)if(idx > N || idy > M) return ;int a_index = idy * K;int b_index = idx;half cval = 0;for(int i = 0; i < K; i++){cval += A[a_index + i] * B[b_index + i * N];}int c_index = idy * N + idx;C[c_index] = cval;}
该实现版本性能如下：
该版本可以看出，对于矩阵A的访问步长为 K * sizeof(float), 显然不满足GPU访存合并的原则。其次，计算过程中存在大量的数据重复加载，例如A矩阵的第一行数据，会在计算第一行每一列数据的时候被反复加载。
合并访存优化
首先，可以将矩阵A进行转置以达到访存合并，其次单线程可以计算更多的输出点，以减少数据的重复加载，向量化加载也可以更好的提高带宽利用率；
优化方案如下图所示：
转置后使用A的一列与B的一列乘累加，得到C的一个点；代码实现如下：
// global_work_size[] = {(N + 3)/4, (M + 3) / 4}// 单线程计算16个点；// 读者可以在不同架构的设备上尝试其他方案，虽然单线程计算点越多重复加载数据越小，但是也可能导致寄存器溢出，性能反而下降严重#pragma OPENCL EXTENSION cl_khr_fp16 : enable__kernel void gemm_opt(__global half* A, __global half* B, __global half* C, int M, int N, int K){int idx = get_global_id(0) << 2;int idy = get_global_id(1) << 2;if(idx > N || idy > M) return;half4 cval[4] = {(half4)(0), (half4)(0), (half4)(0), (half4)(0)};for(int i = 0; i < K; i++){half4 a = vload4(0, A + idy + i * M);half4 b = vload4(0, B + idx + i * N);cval[0] += a.s0 * b;cval[1] += a.s1 * b;cval[2] += a.s2 * b;cval[3] += a.s3 * b;}vstore4(cval[0], 0, C + idy * N + idx);vstore4(cval[1], 0, C + (idy + 1) * N + idx);vstore4(cval[2], 0, C + (idy + 2) * N + idx);vstore4(cval[3], 0, C + (idy + 3) * N + idx);}
该版本实现性能如下：
从数据可以看出，目前版本相对于直接实现版本提升了17倍，Adreno的实现提升8倍。其主要原因在于缺少L1 Cache的加持Adreno设备的buffer吞吐远低于Mali设备的吞吐。接下来通过使用Texture内存对两种设备做进一步的优化。
Adreno设备的Texture方案
下图是Qualcomm文档中关于纹理内存的描述，
从图中可以看出，shader在加载数据的时候，texture内存和buffer内存是通过不同的通道进行的，texture内存的加载可以使用到单独的Texture Processor/L1 Cache,而buffer内存的加载只能使用L2 Cache，因此合理的使用Texture 内存存储数据可以进一步提升上一版本性能。
Texture和buffer内存一般是通过不同的硬件单元进行加载的，所以，在使用纹理内存的时候，是选择A/B其一存储在Texture 内存，另外一个存储到Buffer内存呢？还是选择两块内存都使用Texture呢？
这里给出结论，Qualcomm上使用双Texture内存，Mali部分机型上使用两种不同的内存类型来存储数据，部分机型使用双Texture内存。感兴趣的读者可以在不同机型上测试不同的case。吐槽一下，Qualcomm的文档更新太慢，以上信息来源于5xx GPU的文档。
使用纹理内存的优化版本如下：
// global_work_size[] = {(N + 3) / 4, (M + 3) / 4}#pragma OPENCL EXTENSION cl_khr_fp16 : enable__constant sampler_t default_sampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP | CLK_FILTER_NEAREST;__kernel void gemm_opt(__read_only image2d_t A, __read_only image2d_t B, __write_only image2d_t C, int M, int N, int K){int idx = get_global_id(0);int idy = get_global_id(1);if((idx << 2) > N || (idy << 2) > M) return;half4 c[4] = {(half4)(0), (half4)(0), (half4)(0), (half4)(0)};for(int i = 0; i < K; i++){half4 a = read_imageh(A, default_sampler, (int2)(idy, i));half4 b = read_imageh(B, default_sampler, (int2)(idx, i));c[0] += a.x * b;c[1] += a.y * b;c[2] += a.z * b;c[3] += a.w * b;}idy = idy << 2;write_imageh(C, (int2)(idx, idy), c[0]);write_imageh(C, (int2)(idx, idy + 1), c[1]);write_imageh(C, (int2)(idx, idy + 2), c[2]);write_imageh(C, (int2)(idx, idy + 3), c[3]);}
该版本性能如下：
该版本相对于基础版本有3倍的提升，可以看出Texture内存的使用可以极大的提升访存性能，进而发挥GPU的计算能力。
Mali Valhall 设备优化方案
纹理内存方案
上图是Mali 各个架构下的GPU型号。Mali设备都是硬件厂商可配置的，同一GPU型号，可能存在多种配置。本文采用Valhall架构下的G77进行测试，SOC为MTK的天玑1000，设备为G77 MP9.
上文最后一个版本是针对Qualcomm架构给出的双Texture版本，那么在mali架构下是否是相同方案最优呢？Bifrost/ValHall架构相关文档中并未提及Texture内存与Buffer内存使用不同的Cache，因此这两个架构下，可以享受不同加载单元可以并行加载所带来的收益。同时， 从G76开始，Arm针对Texture内存的加载进行了加强，所以在Mali架构下，采用单Texture内存的方案进行优化。其他架构下，感兴趣的读者可以查看相应的文档或者相关测试。
实现方案如下：
// global_work_size[] = {(N + 3) / 4, (M + 3) / 4}#pragma OPENCL EXTENSION cl_khr_fp16 : enable__constant sampler_t default_sampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP | CLK_FILTER_NEAREST;__kernel void gemm_opt(__read_only image2d_t A, __global half* B, __write_only image2d_t C, int M, int N, int K){int idx = get_global_id(0);int idy = get_global_id(1);if((idx << 2) > N || (idy << 2) > M) return;half4 c[4] = {(half4)(0), (half4)(0), (half4)(0), (half4)(0)};int idx_ofs = idx << 2;for(int i = 0; i < K; i++){half4 a = read_imageh(A, default_sampler, (int2)(idy, i));half4 b = vload4(0, B + idx_ofs + i * N);c[0] += a.x * b;c[1] += a.y * b;c[2] += a.z * b;c[3] += a.w * b;}idy = idy << 2;write_imageh(C, (int2)(idx, idy), c[0]);write_imageh(C, (int2)(idx, idy + 1), c[1]);write_imageh(C, (int2)(idx, idy + 2), c[2]);write_imageh(C, (int2)(idx, idy + 3), c[3]);}
该版本性能如下：
该版本相对于基础版本有10%左右的性能提升。因为Mali设备的Image内存相对于buffer内存吞吐优势并不明显，所以从buffer版本到Texture版本，Adreno的性能提升大于Mali设备的性能提升。
以上版本仅通过调整使用的内存类型提升数据吞吐以提升GEMM的性能。在此基础上，可以进一步通过更优的tile划分，更优的LocalWorkSize的配置来进一步提升GEMM性能。这些优化手段会给当前版本带来更大的性能提升，通过更深入的优化，在当前版本基础上，两款GPU都可以有至少50%的性能提升，之后的文章中会逐步介绍。
Mali(Valhall)的FMA方案
Mali GPU的valhall架构相对于之前的biforst架构做了大幅调整，ValHall架构开始其渲染和计算使用相同的统一的计算单元进行。下图是关于ValHall架构处理单元的介绍，可以看到，一个FMA单元单个周期可以处理16个FP32的FMA和32个FP16的FMA指令。
下面为使用FMA指令的优化版本：
// global_work_size[] = {(N + 3) / 4, (M + 3) / 4}#pragma OPENCL EXTENSION cl_khr_fp16 : enable__constant sampler_t default_sampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP | CLK_FILTER_NEAREST;__kernel void gemm_opt(__read_only image2d_t A, __global half* B, __write_only image2d_t C, int M, int N, int K){int idx = get_global_id(0);int idy = get_global_id(1);if((idx << 2) > N || (idy << 2) > M) return;half4 c[4];for(int i = 0; i < 4; i++){c[i] = (half4)(0);}int idx_ofs = idx << 2;for(int i = 0; i < K; i += 2){half4 a0 = read_imageh(A, default_sampler, (int2)(idy, (i + 0)));half4 a1 = read_imageh(A, default_sampler, (int2)(idy, (i + 1)));half4 b0 = vload4(0, B + idx_ofs + (i + 0) * N);half4 b1 = vload4(0, B + idx_ofs + (i + 1) * N);c[0] = fma(a0.x, b0, c[0]);c[1] = fma(a0.y, b0, c[1]);c[2] = fma(a0.z, b0, c[2]);c[3] = fma(a0.w, b0, c[3]);c[0] = fma(a1.x, b1, c[0]);c[1] = fma(a1.y, b1, c[1]);c[2] = fma(a1.z, b1, c[2]);c[3] = fma(a1.w, b1, c[3]);}idy = idy << 2;write_imageh(C, (int2)(idx, idy), c[0]);write_imageh(C, (int2)(idx, idy + 1), c[1]);write_imageh(C, (int2)(idx, idy + 2), c[2]);write_imageh(C, (int2)(idx, idy + 3), c[3]);}
该版本对具体性能如下：
使用FMA单元后，性能提升有36%左右。而在Adreno650上使用FMA则非常的慢，这是因为在5xx的文档中高通指出其FMA内置函数是通过软件模拟的，非常慢，目前看即使到650设备为止，该指令依然是软件模拟的。
下图是本文各版本之间的性能性能对比图，可以看出不同实现之间的巨大差异，后期通过更细的优化方法，将得到更大比例的性能提升。
其他优化方案简介
前文一直使用的是单线程计算16个点，这是一种分块方案，但未必是最优的；所以在接下来的优化方案中，可以使用在各个维度上的分块策略，提升数据的复用度和cache命中率；合理的分块可以为矩阵乘法带来大幅度的性能提升。
除了分块策略之外，前文的LocalWorkSize一直是NULL，使用编译器的默认work group方案；在GPU优化中work group的划分，对资源划分以及调度都有很大影响。在adreno和mali的文档中也都有描述，默认的local work size未必是最优的。因此更好的local work group划分也将更好的提升性能。
除此之外，高通设备的local memory等其他资源也都有诸多探索空间，之后会逐步展开。下图是目前使用一些细节优化所达到的较优的优化版本性能。
本文主要根据Adreno和Mali硬件上访存策略的差异，对初始版本做了简单优化。当前最优版本无论是Adreno还是Mali上距离峰值性能还有很大差异，所以在后续的介绍中会针对具体配置，在tile划分策略，LocalMemory的使用以及Local Work Size的配置等方面进行更细致的优化，进一步提升当前版本性能。
下面是广告时间
商汤研究院-高性能异构计算研究员GPU方向地点：北京/上海/深圳
职位描述：参与深度学习计算引擎PPL的GPU架构研发，支持公司内部业务广泛落地，并参与维护运营开源社区。
职位要求
计算机基础扎实，熟悉C/C++，具备系统软件开发架构能力；
熟悉计算机体系结构以及并行计算基本技术；
具有nvidia，amd，intel，qualcomm，arm或者imagination等至少一种GPU体系结构上的性能调优经验；
熟悉CUDA，OpenCL，Vulkan，Metal或者OpenGL compute shader等至少一种开发原语；
有代数矩阵运算、信号处理、计算机视觉、图像处理或3D图形学算法在不同处理器上移植和调优经验的优先；
有GPU服务器环境调优或者RDMA高性能通信库开发经验的优先；有开发运营过开源软件或者为知名开源软件贡献过代码的优先；
对高性能计算，异构计算，体系结构和性能优化工作有强烈的兴趣。
工作地点：北京、上海、深圳
点击【阅读原文】投递履历~