【图像处理】NEON编程3 矩阵乘法

原文地址

矩阵

在这篇文章中，我们将会关注怎样高效的进行4x4矩阵相乘，这是在3D图形学中经常性被使用的。我们会假定矩阵在内存中以列的方式进行存储，这是在 OpenGL-ES中所使用的格式。

算法

这里通过展开计算过程，查看下矩阵运算的过程，并且辨别下能够用NEON指令实现的相关操作。

在上图中，我们将第一个矩阵的每一列（红色标记）乘以在第二个矩阵中对应的数值（蓝色标记），然后将结果加起来，得到结果中的一个列的数值。这个操作将被重复四次。

如果每个列都是在NEON寄存器中的一个向量，我们能够用vector-by-scalar multiplication 指令来计算每个列。即上图中显示的指令。我们可以用同样指令的累计版本来将结果都加起来组成每一列的结果。

当我们处理第一个矩阵的列，产生结果中的一列时，读取和写入数据到或者从内存中都是一个线性操作，不需要交错读取及存储指令。

代码

浮点数运算

首先，我们先关注乘以单精度浮点数矩阵的实现。

从内存中加载矩阵到NEON寄存器，该矩阵是使用列顺序进行存储，所以矩阵的列都是线性的存储在内存中。一个列能够用VLD1指令加载到NEON寄存器，用VST1写入到内存。

     vld1.32  {d16-d19}, [r1]!    @ load first eight elements of matrix 0     vld1.32  {d20-d23}, [r1]!    @ load second eight elements of matrix 0     vld1.32  {d0-d3}, [r2]!      @ load first eight elements of matrix 1     vld1.32  {d4-d7}, [r2]!      @ load second eight elements of matrix 1

NEON有32个64位寄存器，我们可以从矩阵中加载所有的数据到寄存器，并且还有剩余的寄存器可以用来作为累加用。这里d16到d23保存第一个矩阵中的16个数据，而d0到d7保存第二个矩阵中的16个数据。

回顾一下: D and Q 寄存器

大部分NEON指令能够用以下两种方式使用寄存器组：

• 32个双字寄存器，64bits大小，从d0到d31
• 16个四字寄存器，128bits大小，从q0到q15

这些寄存器用别名表示，因此，在Q寄存器，以及在相对应的两个D寄存器中，数据都是一样的。例如，q0可被分为d0及d1，数据可以在两种模式下进行访问。在C中，这种方式类似于union

对于浮点数矩阵乘法的例子，我们频繁的使用Q寄存器，因为我们处理4个32-bit浮点数的列，这将对应于一个128位的Q寄存器。

回到代码中

应用NEON乘法指令，我们能够计算矩阵相乘的一个列。

     vmul.f32    q12, q8, d0[0]      @ multiply col element 0 by matrix col 0     vmla.f32    q12, q9, d0[1]      @ multiply-acc col element 1 by matrix col 1     vmla.f32    q12, q10, d1[0]  @ multiply-acc col element 2 by matrix col 2     vmla.f32    q12, q11, d1[1]  @ multiply-acc col element 3 by matrix col 3

第一个指令实现了在矩阵乘法分解图中的操作，x0，x1，x2和x3（在q8寄存器中）分别乘以y0（d0的第0个数据），结果保存在q12.

接下来的指令用于处理第一个矩阵中的其他列，乘以第二个矩阵第一列中对应的元素。结果将会累加到q12，来构成最终结果的第一列。

需要注意的是，乘法指令中的标量是D寄存器。尽管q0[3]的数据与d1[1]相同，但是这里使用d1[1]更加说得过去，并且GNU汇编器并不支持这种格式。因此，我不得不使用D寄存器。

宏：

     .macro mul_col_f32 res_q, col0_d, col1_d     vmul.f32    \res_q, q8, \col0_d[0]      @ multiply col element 0 by matrix col 0     vmla.f32    \res_q, q9, \col0_d[1]      @ multiply-acc col element 1 by matrix col 1     vmla.f32    \res_q, q10, \col1_d[0]  @ multiply-acc col element 2 by matrix col 2     vmla.f32    \res_q, q11, \col1_d[1]  @ multiply-acc col element 3 by matrix col 3     .endm

4x4浮点矩阵相乘可以被这样实现：

     vld1.32  {d16-d19}, [r1]!            @ load first eight elements of matrix 0     vld1.32  {d20-d23}, [r1]!            @ load second eight elements of matrix 0     vld1.32  {d0-d3}, [r2]!              @ load first eight elements of matrix 1     vld1.32  {d4-d7}, [r2]!              @ load second eight elements of matrix 1     mul_col_f32 q12, d0, d1           @ matrix 0 * matrix 1 col 0     mul_col_f32 q13, d2, d3           @ matrix 0 * matrix 1 col 1     mul_col_f32 q14, d4, d5           @ matrix 0 * matrix 1 col 2     mul_col_f32 q15, d6, d7           @ matrix 0 * matrix 1 col 3     vst1.32  {d24-d27}, [r0]!            @ store first eight elements of result     vst1.32  {d28-d31}, [r0]!            @ store second eight elements of result

定点运算

使用定点数运算通常哔浮点数更快，因为它占用更小的内存带宽来读取及写入数值，并且整数的乘法通常来说更快。

然而，在使用定点数运算时，当需要保证程序要求的计算精度时，你必须仔细选择表示方法来避免溢出或者饱和。

使用定点数进行矩阵乘法与浮点数非常类似。在这个例子中个，我们使用 Q1.14 定点数格式，但是操作上与其他格式类似，可能只需要最终左移的位数。

宏定义：

     .macro mul_col_s16 res_d, col_d     vmull.s16   q12, d16, \col_d[0]   @ multiply col element 0 by matrix col 0     vmlal.s16   q12, d17, \col_d[1]   @ multiply-acc col element 1 by matrix col 1     vmlal.s16   q12, d18, \col_d[2]   @ multiply-acc col element 2 by matrix col 2     vmlal.s16   q12, d19, \col_d[3]   @ multiply-acc col element 3 by matrix col 3     vqrshrn.s32 \res_d, q12, #14            @ shift right and narrow accumulator into                                             @  Q1.14 fixed point format, with saturation     .endm

与浮点数版本的宏进行对比，你将会发现：

• 数值是16位而不是32为，我们可以使用D寄存器保存4个输入。
• 两个16位数的相乘结果是一个32位数，我们使用VMULL及VMLAL，因为它们将保存结果至Q寄存器，使用两倍于数据的大小来保存所有的位数。
• 最后的结果为16位，但是累加器是32位的。我们通过VQRSHRN获取16位结果，一个向量，饱和计算，四舍五入，右移窄运算。这将把所有正确的值加起来到各自的数据，右移并饱和计算到新的更窄的数据大小。

32位数到16位数的减少对内存访问有所影响。数据会被加载和存储使用更少的指令，代码如下：

 vld1.16  {d16-d19}, [r1]       @ load sixteen elements of matrix 0 vld1.16  {d0-d3}, [r2]         @ load sixteen elements of matrix 1 mul_col_s16 d4, d0                      @ matrix 0 * matrix 1 col 0 mul_col_s16 d5, d1                      @ matrix 0 * matrix 1 col 1 mul_col_s16 d6, d2                      @ matrix 0 * matrix 1 col 2 mul_col_s16 d7, d3                      @ matrix 0 * matrix 1 col 3 vst1.16  {d4-d7}, [r0]         @ store sixteen elements of result

如果考虑到调度的话

我们会在以后的文章中讨论调度的细节，现在的话，可以看看在代码中查看改进指令调度的影响。

在宏中，临近的乘法指令写入到相同的寄存器，使得NEON流水线在开始下个指令执行时必须等待每个乘法完成。

如果我们把这些指令移出宏然后重新安排顺序，我们能够分离出那些依赖于其他指令的部分，这部分指令能够在其他指令在后台完成时也被执行。
在这个实例中，我们重组代码来分离累加寄存器的使用。

    vmul.f32    q12, q8, d0[0]              @ rslt col0  = (mat0 col0) * (mat1 col0 elt0)     vmul.f32    q13, q8, d2[0]              @ rslt col1  = (mat0 col0) * (mat1 col1 elt0)     vmul.f32    q14, q8, d4[0]              @ rslt col2  = (mat0 col0) * (mat1 col2 elt0)     vmul.f32    q15, q8, d6[0]              @ rslt col3  = (mat0 col0) * (mat1 col3 elt0)     vmla.f32    q12, q9, d0[1]              @ rslt col0 += (mat0 col1) * (mat1 col0 elt1)     vmla.f32    q13, q9, d2[1]              @ rslt col1 += (mat0 col1) * (mat1 col1 elt1)     ...     ...

使用这个版本，在基于 Cortex-A8 的系统中，矩阵乘法的性能将加倍提升。

你可以在 Technical Reference Manual for your Cortex core找到更多的关于指令周期及延迟的细节。为了能够像上述一样提升性能，熟悉这些细节或许会很有帮助。