移动端GPU矩阵乘优化
转自如下链接，并结合自身进行丰富
https://zhuanlan.zhihu.com/p/353975826
移动端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次，取平均值作为执行时间。
刘文志《OpenCL异构并行计算》
OpenCL 2.0 异构计算 [第三版] （Heterogeneous Computing with OpenCL 2.0）
cl_event event; errNum = clEnqueueNDRangeKernel(commandQueue, kernel, 2, NULL, globalWorkSize, NULL, 0, NULL, &event); clWaitForEvents(1, &event); clFinish(commandQueue); cl_ulong time_start; cl_ulong time_end; clGetEventProfilingInfo(event, CL_PROFILING_COMMAND_START, sizeof(time_start), &time_start, NULL); clGetEventProfilingInfo(event, CL_PROFILING_COMMAND_END, sizeof(time_end), &time_end, NULL); double times = time_end - time_start; printf("OpenCL Kernel Execution time is: %0.4f\n", times / 1000000.0);
基础优化版本
直接实现版本
首先按照矩阵乘法的计算公式，实现最简单版本作为base，如下如图所示，A矩阵的第一行乘以B矩阵的第一列得到C矩阵对应行列的一个元素。
代码实现如下：
// global_work_size = {N, M}
//A的维度为M x K, B的维度为K x N, C维度为M x N
#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;
}half数据类型用16位来表示浮点数。这比更为常用的float型的示数范围要小，但他是浮点数家族的新成员。有效位数（0-9）10bits,指数位+(9-15)5bits,符号位
https://zhuanlan.zhihu.com/p/112564372
该实现版本性能如下：
该版本可以看出，对于矩阵A的访问步长为 K * sizeof(float), 显然不满足GPU访存合并的原则。其次，计算过程中存在大量的数据重复加载，例如A矩阵的第一行数据，会在计算第一行每一列数据的时候被反复加载。
合并访存优化
首先，可以将矩阵A进行转置以达到访存合并，其次单线程可以计算更多的输出点，以减少数据的重复加载，向量化加载也可以更好的提高带宽利用率；
优化方案如下图所示：
转置后使用A的一列与B的一列乘累加，得到C的一个点；代码实现如下：
// global_work_size[] = {(N + 3)/4, (M + 3) / 4}
// A的维度为M x K, B的维度为K x N, C维度为M x N
// 单线程计算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;//0到N-1
int idy = get_global_id(1) << 2;//0到M-1
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);
}vloadn
Read vectors from a pointer to memory.
gentypen vloadn (size_t offset , const __global gentype *p)gentypen vloadn (size_t offset ,const __local gentype *p)gentypen vloadn (size_t offset ,const __constant gentype *p)gentypen vloadn (size_t offset ,const __private gentype *p)Description
Return sizeof (gentype
n) bytes of data read from location (p+ (offset*n)).The read address computed as (
p+ (offset*n)) must be 8-bit aligned if gentype is charn, ucharn; 16-bit aligned if gentype is shortn, ushortn; 32-bit aligned if gentype is intn, uintn, floatn; 64-bit aligned if gentype is longn, ulongn.If the double extension is enabled, then in addition to the above the read address must be 64-bit aligned if gentype is long
n, ulongn, or doublen.If the half extension is enabled, the read address computed as (
p+ (offset*n) must be 16-bit aligned.vstoren
Write sizeof bytes (gentype
n) given by data to address.
void vstoren (gentypen data ,size_t offset ,const __global gentype *p)gentype vstoren (gentypen data,size_t offset,const __local gentype *p)gentype vstoren (gentypen data,size_t offset,const __private gentype *p)Description
Write
sizeof(gentypebytes given byn)datato address (p+ (offset*n)). The write address computed as (p+ (offset*n)) must be 8-bit aligned if gentype is charn, ucharn; 16-bit aligned if gentype is shortn, ushortn; 32-bit aligned if gentype is intn, uintn, floatn; 64-bit aligned if gentype is longn, ulongn.If the double extension is enabled, then in addition to the above the address must be 64-bit aligned if gentype is long
n, ulongn, or doublen.If the half extension is enabled, the address computed as (
p+ (offset*n) must be 16-bit aligned.1.加载和保存同类型的数据:如果发送端和接收端的数据类型相同，例如将一个区域的int4型向量加载到另一区域中，只知道“=”即可。除了对变量赋值外，等号还可以将一个区域的数据赋值传输到另一个区域中。例如，内核想将全局内存中的数据加载到局部内存中，并将处理完的结果保存回全局内存。如果要将局部内存改为私有内存，只需要将数据声明中的_local标识符去掉即可，默认的标识符便是_private,因此，如果没有标识符，编译器会将变量自动保存在私有内存中。
2.将标量数组加载保存到向量中：大多数情况下，可能需要用向量来处理数据，但数据输入却一般都保存在标量数组中。这时，“=”就不管用了，你不可以将一个地址空间中的数据类型转换到另外一个地址空间中。所以OpenCL提供了函数vloadn，用来将标量数据加载保存到向量中：
vector vloadn(size_t offset, const __(g|c|l|p) scalar *mem)；
n表示的是向量所包含的向量分数量，它必须被设定为2，3，4，8或16.vector可以是任意一种向量类型，而scalar则必须和vector中的分量类型一致。vector数据所在的内存地址可以和scalar数据所在的地址不同。例如，假设float_vec是一个float4型向量，而float_array是一个float型数组。下面代码的功能就是将float_array数组中的四个float型数据保存到float_vec向量中：
float_vec = vload4(0,float_array);
参数offset决定了哪些数组元素将被放到向量之中。这个参数和向量的大小（而不是它的标量分量的大小）有关。
3.将向量保存到标量数组中：和函数vloadn将标量数组保存到向量类似，函数vstoren做的是反向操作，将向量中的数据保存到标量数组中：
void vstoren(vector vec, size_t offset, __(g|l|p) scalar *mem)；
这个函数将向量vec中的数据保存到mem指向的标量数组保存到offset确定的地址空间中。和函数vloadn一样，n必须等于向量vec中的分数量，其值必须设定为2，3，4，8或16.但是和函数vload不同，标量数组不能保存到一个常数地址空间之中，下面代码将int4型向量int_vec保存到指针int_array所指向的整形数组中;
vstore4(int_vec, 0, int_array);
通过修改offset参数，可以将向量保存到数组中的不同位置，位置差为（向量大小*offset）。通过指针运算，可以将向量保存在数组中非对齐的内存位置。
OPenCL标准还分别提供了函数vload_halfn将半精度数据保存到浮点型向量中，以及vstore_halfn将浮点型向量中的数据保存到半精度数组中。https://man.opencl.org/read_imageh3d.html
https://blog.csdn.net/INFINALGEORGE/article/details/105198959
该版本实现性能如下：
从数据可以看出，目前版本相对于直接实现版本提升了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]);
}read_imageh
Use the coordinate (coord.x, coord.y, coord.z) to do an element lookup in the 3D image object specified by
image.coord.wis ignored.
- 如果读取的时候使用的是
float2坐标，假设为坐标为(w, h)，那么，其返回的值为(w - 0.5, h - 0.5)处的插值结果，插值的方式为我们常规意义，或者在CPU代码中对该图像进行双线性插值。当然这也和采样器sampler_t对象设置为CLK_FILTER_LINEAR有关。如果其设置为CLK_FILTER_NEAREST，那么肯定就是为最近邻插值了。举例来说，对于(float2)(1.0, 1.0)坐标，其插值目标为(1.0 - 0.5, 1.0 - 0.5)，位于(0,0), (0, 1), (1, 0), (1,1)四个像素点中间，根据双线性插值计算。其结果即为0.75- 如果读取的时候使用的是
int2坐标，那么其坐标与值的关系就和CPU中处理该image一样。
该版本性能如下：
该版本相对于基础版本有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指令。
Valhall架构的本质在新的执行核心。前代Bifrost架构是4-wide和8-wide设计，G72核心部分的执行模块就包含4-wide标量SIMD单元，warp size为4；G76则增加到两个4-wide单元，warp size为8。warp是GPU的最基本可调度单元，SIMD过程中数据处理的最小单位；在所有线程中，同时执行同一指令。
下面为使用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的配置等方面进行更细致的优化，进一步提升当前版本性能。
参考文章：
https://blog.csdn.net/qq_31049419/article/details/111565316
https://zhuanlan.zhihu.com/p/112564372
https://news.mydrivers.com/1/643/643900.htm
https://blog.csdn.net/qq_36533552/article/details/105885714
640
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
