ARM处理器NEON编程及优化技巧三—矩阵乘法的实例

原文地址：http://houh-1984.blog.163.com/blog/static/311278342011119101429988/

ARM的NEON协处理器技术是一个64/128-bit的混合SIMD架构，用于加速包括视频编码解码、音频解码编码、3D图像、语音和图像等多媒体和信号处理应用。本文主要介绍如何使用NEON的汇编程序来写SIMD的代码，包括如何开始NEON的开发，如何高效的利用NEON。首先会关注内存操作，即如何变更指令来灵活有效的加载和存储数据。接下来是由于SIMD指令的应用而导致剩下的若干个单元的处理，然后是用一个矩阵乘法的例子来说明用NEON来进行SIMD优化，最后关注如何用NEON来优化各种各样的移位操作，左移或者右移以及双向移位等。本节是一个用NEON优化矩阵乘法的实例。

矩阵

本节将介绍如何用NEON有效的处理一个4x4的矩阵乘法运算，这种类型的运算经常用于3D图形，我们认为这些矩阵在内存里是按照列为主排列的，这是按照OPENGL-ES的通用格式。

矩阵乘法算法

我们首先看一下矩阵乘法的计算方式，计算的展开，用NEON指令来进行子操作过程。



图1. 以列为主的矩阵乘法运算

由于数据是按照列序存储的，因而矩阵乘法就是把第一个矩阵的每一列乘以第二个矩阵的每一行，然后把乘积结果相加。乘累加结果 作为结果矩阵的一个元素。



图2. 矩阵乘法中的向量乘以标量的运算

假设每列元素在NEON寄存器中表示为一个向量，那么上述的矩阵乘法就是一个向量乘以标量的运算，而后续的累加也同样可以同向量乘以标量的累加指令实现。因为我们的操作是在第一个矩阵的列，然后计算列的结果，读列元素和写列元素都是线性的加载和存储操作，不需要interleave的加载和存储操作。

代码

浮点运算版本

首先看一个单精度浮点的矩阵乘法实现。首先加载矩阵元素到NEON寄存器，然后按照列序做乘法，用VLD1做线性的加载数据到NEON寄存器，用VST1把计算结果保存到内存。

[code]    vld1.32     {d16-d19}, [r1]!            @ 加载矩阵0的上8个元素

[code]    vld1.32     {d20-d23}, [r1]!            @ 加载矩阵0的下8个元素

[code]    vld1.32     {d0-d3}, [r2]!              @ 加载矩阵1的上8个元素

[code]    vld1.32     {d4-d7}, [r2]!              @ 加载矩阵1的下8个元素

NEON有32个64位寄存器，因而加载所有的输入矩阵元素到16个64-bit寄存器，我们仍然有16个64位寄存器做后续的处理。

D和Q寄存器

多数的NEON指令有两种方法来访问寄存器组：

作为32个双字寄存器，64-bit位宽，命名为d0-d31
作为16个四字寄存器，128-bit位宽，命名为q0-q15



图3. NEON寄存器组

这些寄存器中一个Q寄存器是一对D寄存器的别名，如Q0是d0和d1寄存器对的别名，寄存器中的值可以用两种方式访问。这种实现方式很类似C语言里的union联合的数据结构。对于浮点的矩阵乘法，我们会经常使用Q寄存器的表达方式，因为经常会处理4个32-bit的单精度浮点，这对应于128-bit的Q寄存器。

代码部分

通过以下4条NEON乘法指令能完成一列4个结果：

[code]    vmul.f32    q12, q8, d0[0]        @ 向量乘以标量（MUL），矩阵0的第一列乘以矩阵1的每列的第一个元素0

[code]    vmla.f32    q12, q9, d0[1]        @ 累加的向量乘以标量（MAC），矩阵0的第二列乘以矩阵1的每列的第二个元素1

[code]    vmla.f32    q12, q10, d1[0]     @ 累加的向量乘以标量（MAC），矩阵0的第二列乘以矩阵1的每列的第二个元素2

[code]    vmla.f32    q12, q11, d1[1]     @ 累加的向量乘以标量（MAC），矩阵0的第二列乘以矩阵1的每列的第二个元素3

第一条指令是图2中的列元素x0, x1, x2, x3 (寄存器q8)乘以y0 (d0[0])，然后结果保存到q12寄存器。接下来的指令操作类似，就是把第一个矩阵的其他列乘以第二个矩阵的第一列的响应元素，结果累加到寄存器Q12里。需要注意的是标量元素如d1[1]也可以用q0[3]表示，但是可能编译器如GNU汇编器会不能接受这种方式。

如果我们只需要矩阵乘以向量的运算，如很多3D图像处理中的那样，那么此时的计算就结束了，可以把结果向量保存 到内存了，但是为了完成矩阵相乘，还需要完成后面的迭代操作，使用寄存器Q1到Q3的y4到yF的元素。如果定义如下的宏，那么就能简化代码结构了：

[code]    .macro mul_col_f32 res_q, col0_d, col1_d

[code]    vmul.f32    \res_q, q8, \col0_d[0]      @向量乘以标量（MUL），矩阵0的第一列乘以矩阵1的每列的第一个元素0

[code]    vmla.f32    \res_q, q9, \col0_d[1]      @累加的向量乘以标量（MAC），矩阵0的第二列乘以矩阵1的每列的第二个元素1

[code]    vmla.f32    \res_q, q10, \col1_d[0]     @累加的向量乘以标量（MAC），矩阵0的第二列乘以矩阵1的每列的第二个元素2

[code]    vmla.f32    \res_q, q11, \col1_d[1]     @累加的向量乘以标量（MAC），矩阵0的第二列乘以矩阵1的每列的第二个元素3

[code]    .endm

那么整个4x4的矩阵乘法代码可能如下：

[code]    vld1.32     {d16-d19}, [r1]!            @ 加载矩阵0的上8个元素

[code]    vld1.32     {d20-d23}, [r1]!            @ 加载矩阵0的下8个元素

[code]    vld1.32     {d0-d3}, [r2]!              @ 加载矩阵1的上8个元素

[code]    vld1.32     {d4-d7}, [r2]!              @ 加载矩阵1的下8个元素

[code]    mul_col_f32 q12, d0, d1                 @ 矩阵 0 * 矩阵1的第一列

[code]    mul_col_f32 q13, d2, d3                 @ 矩阵 0 * 矩阵1的第二列

[code]    mul_col_f32 q14, d4, d5                 @ 矩阵 0 * 矩阵1的第三列

[code]    mul_col_f32 q15, d6, d7                 @ 矩阵 0 * 矩阵1的第四列

[code]    vst1.32     {d24-d27}, [r0]!            @ 保存结果的上8个元素

[code]    vst1.32     {d28-d31}, [r0]!            @ 保存结果的下8个元素

定点算法

定点算法计算往往比浮点计算更快，因为往往定点运算可能需要更少的内存带宽，整数值的乘法也会比浮点算法更为简单。但是定点算法，你需要很仔细的选择表示格式来避免溢出或者饱和，这些会影响你的算法最终的精度。定点算法实现的矩阵乘法和浮点算法类似，在本例中，用Q1.14定点格式,但是基本的实现格式基本类似，只是实现中可能需要对结果做一些移位调整。下面是列乘的宏：

[code]    .macro mul_col_s16 res_d, col_d

[code]    vmull.s16   q12, d16, \col_d[0]         @ 向量乘以标量（MUL），矩阵0的第一列乘以矩阵1的每列的第一个元素0

[code]    vmlal.s16   q12, d17, \col_d[1]         @ 累加的向量乘以标量（MAC），矩阵0的第二列乘以矩阵1的每列的第二个元素1

[code]    vmlal.s16   q12, d18, \col_d[2]         @ 累加的向量乘以标量（MAC），矩阵0的第二列乘以矩阵1的每列的第二个元素2

[code]    vmlal.s16   q12, d19, \col_d[3]         @ 累加的向量乘以标量（MAC），矩阵0的第二列乘以矩阵1的每列的第二个元素3

[code]    vqrshrn.s32 \res_d, q12, #14             @ 把结果右移14位，并把累加结果变成Q1.14定点格式，并饱和运算

[code]    .endm

比较定点和浮点算法的宏，你会发现如下的主要区别：

矩阵元素的值为16位而不是32位，因而用D寄存器来保存4个输入元素
矩阵乘法的结果是把16x16=32位的数据，使用VMULL和VMLAL来吧结果保存到Q寄存器。
最终结果也是16位，因而需要把32位累加器结果来得到16-bit的结果，使用VQRSHRN，饱和处理把32位的结果舍入到16位的narrow右移操作。

把数据从32-bits变成16-bits也能有效的处理内存访问，加载和存储数据都只需要更少的带宽。

[code]    vld1.16     {d16-d19}, [r1]             @ 加载16个元素到矩阵0

[code]    vld1.16     {d0-d3}, [r2]               @ 加载16个元素到矩阵1

[code]    mul_col_s16 d4, d0                      @ 矩阵0乘以矩阵1的列0

[code]    mul_col_s16 d5, d1                      @ 矩阵0乘以矩阵1的列1

[code]    mul_col_s16 d6, d2                      @ 矩阵0乘以矩阵1的列2

[code]    mul_col_s16 d7, d3                      @ 矩阵0乘以矩阵1的列3

[code]    vst1.16     {d4-d7}, [r0]               @ 保存16个结果元素

指令重排

我们先展示一下指令重排如何能提高代码性能。在宏中，临近的乘法指令会写入到相同的目标寄存器，这会让NEON的流水线等待前面的乘法结果完成才能开始下一条指令的执行。如果不使用宏定义，而合理安排指令的次序，把那些相关依赖的指令变成不依赖，这些指令就能并发而不会造成流水线的stall。

[code]    vmul.f32    q12, q8, d0[0]              @ rslt col0  = (mat0 col0) * (mat1 col0 elt0)

[code]    vmul.f32    q13, q8, d2[0]              @ rslt col1  = (mat0 col0) * (mat1 col1 elt0)

[code]    vmul.f32    q14, q8, d4[0]              @ rslt col2  = (mat0 col0) * (mat1 col2 elt0)

[code]    vmul.f32    q15, q8, d6[0]              @ rslt col3  = (mat0 col0) * (mat1 col3 elt0)

[code]    vmla.f32    q12, q9, d0[1]              @ rslt col0 += (mat0 col1) * (mat1 col0 elt1)

[code]    vmla.f32    q13, q9, d2[1]              @ rslt col1 += (mat0 col1) * (mat1 col1 elt1)

[code]    ...

[code]    ...

用以上的处理方式，矩阵乘法的性能在Cortex-A8 处理平台上性能提升了一倍。从文档Technical
Reference Manual for your Cortex core可以看到 各个指令的需要时间以及延迟，有这些延迟周期，能够更为合理的安排代码次序，提升性能。

源代码

以上两个函数的源代码可以参考如下：



matrix_asm_sched.s (4.54K)

Reference：

http://www.blog.163.com/houh-1984/

Coding for NEON - Part 3:
Matrix Multiplication

http://forums.arm.com/index.php?app=blog&module=display§ion=archive&blogid=7&st=50