移动端GPU矩阵乘优化

首发：NeuralTalk
作者：姚定界

移动端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，峰值数据是实际测试乘加计算的吞吐，非理论峰值.

测试数据

矩阵维度为：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的一个点；代码实现如下：

1.  // global_work_size[] = {(N + 3)/4, (M + 3) / 4}
2.  // 单线程计算16个点；
3.  // 读者可以在不同架构的设备上尝试其他方案，虽然单线程计算点越多重复加载数据越小，但是也可能导致寄存器溢出，性能反而下降严重    
4.  #pragma OPENCL EXTENSION cl_khr_fp16 : enable
    
5.  `__kernel void gemm_opt(__global half* A, __global half* B, __global half* C,  int M,  int N,  int K)`
    
6.  `{`
    
7.   `int idx = get_global_id(0)  <<  2;`
    
8.   `int idy = get_global_id(1)  <<  2;`
    
9.    
    
10.   `if(idx > N || idy > M)  return;`
    
11.    
    
12.   `half4 cval[4]  =  {(half4)(0),  (half4)(0),  (half4)(0),  (half4)(0)};`
    
13.   `for(int i =  0; i < K; i++)`
    
14.   `{`
    
15.   `half4 a = vload4(0, A + idy + i * M);`
    
16.   `half4 b = vload4(0, B + idx + i * N);`
    
17.    
    
18.   `cval[0]  += a.s0 * b;`
    
19.   `cval[1]  += a.s1 * b;`
    
20.   `cval[2]  += a.s2 * b;`
    
21.   `cval[3]  += a.s3 * b;`
    
22.   `}`
    
23.   `vstore4(cval[0],  0, C + idy * N + idx);`
    
24.   `vstore4(cval[1],  0, C +  (idy +  1)  * N + idx);`
    
25.   `vstore4(cval[2],  0, C +  (idy +  2)  * N + idx);`
    
26.   `vstore4(cval[3],  0, C +  (idy +  3)  * N + idx);`
    
27.  `}`

该版本实现性能如下：

从数据可以看出，目前版本相对于直接实现版本提升了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的文档。

使用纹理内存的优化版本如下：

1.  `// global_work_size[] = {(N + 3) / 4, (M + 3) / 4}`
    
2.  `#pragma OPENCL EXTENSION cl_khr_fp16 : enable`
    
3.  `__constant sampler_t default_sampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP | CLK_FILTER_NEAREST;`
    
4.  `__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)`
    
5.  `{`
    
6.   `int idx = get_global_id(0);`
    
7.   `int idy = get_global_id(1);`
    
8.    
    
9.   `if((idx <<  2)  > N ||  (idy <<  2)  > M)  return;`
    
10.    
    
11.   `half4 c[4]  =  {(half4)(0),  (half4)(0),  (half4)(0),  (half4)(0)};`
    
12.    
    
13.   `for(int i =  0; i < K; i++)`
    
14.   `{`
    
15.   `half4 a = read_imageh(A, default_sampler,  (int2)(idy, i));`
    
16.   `half4 b = read_imageh(B, default_sampler,  (int2)(idx, i));`
    
17.    
    
18.   `c[0]  += a.x * b;`
    
19.   `c[1]  += a.y * b;`
    
20.   `c[2]  += a.z * b;`
    
21.   `c[3]  += a.w * b;`
    
22.   `}`
    
23.    
    
24.   `idy = idy <<  2;`
    
25.   `write_imageh(C,  (int2)(idx, idy), c[0]);`
    
26.   `write_imageh(C,  (int2)(idx, idy +  1), c[1]);`
    
27.   `write_imageh(C,  (int2)(idx, idy +  2), c[2]);`
    
28.   `write_imageh(C,  (int2)(idx, idy +  3), c[3]);`
    
29.  `}`

该版本性能如下：

该版本相对于基础版本有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内存的方案进行优化。其他架构下，感兴趣的读者可以查看相应的文档或者相关测试。

实现方案如下：

1.  `// global_work_size[] = {(N + 3) / 4, (M + 3) / 4}`
    
2.  `#pragma OPENCL EXTENSION cl_khr_fp16 : enable`
    
3.  `__constant sampler_t default_sampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP | CLK_FILTER_NEAREST;`
    
4.  `__kernel void gemm_opt(__read_only image2d_t A, __global half* B, __write_only image2d_t C,  int M,  int N,  int K)`
    
5.  `{`
    
6.   `int idx = get_global_id(0);`
    
7.   `int idy = get_global_id(1);`
    
8.    
    
9.   `if((idx <<  2)  > N ||  (idy <<  2)  > M)  return;`
    
10.    
    
11.   `half4 c[4]  =  {(half4)(0),  (half4)(0),  (half4)(0),  (half4)(0)};`
    
12.    
    
13.   `int idx_ofs = idx <<  2;`
    
14.   `for(int i =  0; i < K; i++)`
    
15.   `{`
    
16.   `half4 a = read_imageh(A, default_sampler,  (int2)(idy, i));`
    
17.   `half4 b = vload4(0, B + idx_ofs + i * N);`
    
18.   `c[0]  += a.x * b;`
    
19.   `c[1]  += a.y * b;`
    
20.   `c[2]  += a.z * b;`
    
21.   `c[3]  += a.w * b;`
    
22.   `}`
    
23.    
    
24.   `idy = idy <<  2;`
    
25.   `write_imageh(C,  (int2)(idx, idy), c[0]);`
    
26.   `write_imageh(C,  (int2)(idx, idy +  1), c[1]);`
    
27.   `write_imageh(C,  (int2)(idx, idy +  2), c[2]);`
    
28.   `write_imageh(C,  (int2)(idx, idy +  3), c[3]);`
    
29.  `}`

该版本性能如下：

该版本相对于基础版本有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指令的优化版本：

1.  `// global_work_size[] = {(N + 3) / 4, (M + 3) / 4}`
    
2.  `#pragma OPENCL EXTENSION cl_khr_fp16 : enable`
    
3.  `__constant sampler_t default_sampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP | CLK_FILTER_NEAREST;`
    
4.  `__kernel void gemm_opt(__read_only image2d_t A, __global half* B, __write_only image2d_t C,  int M,  int N,  int K)`
    
5.  `{`
    
6.   `int idx = get_global_id(0);`
    
7.   `int idy = get_global_id(1);`
    
8.    
    
9.   `if((idx <<  2)  > N ||  (idy <<  2)  > M)  return;`
    
10.    
    
11.   `half4 c[4];`
    
12.   `for(int i =  0; i <  4; i++)`
    
13.   `{`
    
14.   `c[i]  =  (half4)(0);`
    
15.   `}`
    
16.    
    
17.   `int idx_ofs = idx <<  2;`
    
18.   `for(int i =  0; i < K; i +=  2)`
    
19.   `{`
    
20.   `half4 a0 = read_imageh(A, default_sampler,  (int2)(idy,  (i +  0)));`
    
21.   `half4 a1 = read_imageh(A, default_sampler,  (int2)(idy,  (i +  1)));`
    
22.    
    
23.   `half4 b0 = vload4(0, B + idx_ofs +  (i +  0)  * N);`
    
24.   `half4 b1 = vload4(0, B + idx_ofs +  (i +  1)  * N);`
    
25.    
    
26.   `c[0]  = fma(a0.x, b0, c[0]);`
    
27.   `c[1]  = fma(a0.y, b0, c[1]);`
    
28.   `c[2]  = fma(a0.z, b0, c[2]);`
    
29.   `c[3]  = fma(a0.w, b0, c[3]);`
    
30.    
    
31.   `c[0]  = fma(a1.x, b1, c[0]);`
    
32.   `c[1]  = fma(a1.y, b1, c[1]);`
    
33.   `c[2]  = fma(a1.z, b1, c[2]);`
    
34.   `c[3]  = fma(a1.w, b1, c[3]);`
    
35.   `}`
    
36.    
    
37.   `idy = idy <<  2;`
    
38.   `write_imageh(C,  (int2)(idx, idy), c[0]);`
    
39.   `write_imageh(C,  (int2)(idx, idy +  1), c[1]);`
    
40.   `write_imageh(C,  (int2)(idx, idy +  2), c[2]);`
    
41.   `write_imageh(C,  (int2)(idx, idy +  3), c[3]);`
    
42.  `}`

该版本对具体性能如下：

使用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的配置等方面进行更细致的优化，进一步提升当前版本性能。

往期回顾