ARMv7 Neon指令集

半块方便面

已于 2023-09-14 12:52:31 修改

阅读量294

点赞数

分类专栏：性能优化文章标签： arm开发系统架构 c++ 算法

于 2023-08-28 12:26:06 首次发布

本文链接：https://blog.csdn.net/qq_35227495/article/details/132536278

版权

性能优化专栏收录该内容

12 篇文章 0 订阅

订阅专栏

1、概念

neon 单元不能直接访问 S 寄存器，需要经过 VFP 单元转换

S 寄存器代表 single register
D 寄存器代表 double-precision register
Q 寄存器代表 ARM NEON 指令集中的一种寄存器类型，用于存储多个数据元素，通常用于并行执行向量操作。每个 Q 寄存器可以包含多个 D 寄存器。

寄存器映射关系：

Neon指令集的数据类型

在这里插入图片描述

VFP指令的数据类型

利用编译器优化代码时的一些规则（使用向量化）

简短
不要使用 breadk 退出循环
循环次数为 2 的幂
让编译器明确知道循环的次数

例如，一个累加程序，如果循环的次数是 4 的倍数，则可以在循环中明确的指出，二进制最后三位为 0，从而使得编译器安全的向量化。
```
int accumulate(int * c, int len)
{
    int i, retval;
    for(i = 0, retval = 0; i < (len & ~3) ; i++)
    {
        retval = retval + c[i]
    }
    return retval;
}
```
本文档最后展示了一个循环次数未知时的解决方案。
循环中调用的函数应该为内联函数
使用带索引的数组比带指针的数组更好
使用 restrict 关键字告诉编译器指针不引用重叠区域
间接寻址会导致无法向量化
避免当前循环依赖上一轮的结果
避免在循环中使用条件语句
使用合适的数据结构，neon不支持 64 位浮点数，不要使用 long long 类型
如果结构体的数据类型不同，neon不会对其执行向量化操作

struct pixel
{
char r;
short g; /* Green channel contains more information */
char b;
}screen[10];

编译器仅尝试对使用结构中所有数据项的加载操作进行矢量化,由于存在未使用的字节，编译器不会尝试对该结构的加载进行矢量化。NEON加载指令可以加载不对齐的结构。因此，在这种情况下，最好去掉填充，以便编译器可以对加载进行矢量化。

struct aligned_pixel
{
char r;
char g;
char b;
char not_used; /* Padding used to keep r aligned to a 32-bit word */
}screen[10];

Neon汇编器和 ABI 限制

在ARM架构中，子程序（函数）在使用寄存器S16-S31（D8-D15，Q4-Q7）之前必须保护（保存）这些寄存器的内容。这意味着如果在子程序中要使用这些寄存器，需要在进入子程序之前将其内容保存在内存中，以便在子程序执行完毕后可以恢复原来的值。

相反，寄存器S0-S15（D0-D7，Q0-Q3）不需要被保护，可以在标准过程调用的变体中用于传递参数或返回结果。这意味着这些寄存器的内容可以在子程序内部修改，而不需要在子程序之前将其保存。

同样，寄存器D16-D31（Q8-Q15）也不需要被保护，这意味着在使用这些寄存器之前无需将其内容保存，可以在子程序中自由使用它们。

在过程调用标准中，浮点参数可以通过两种方式传递：

对于软件浮点，浮点参数将使用ARM寄存器R0-R3传递，并在需要时存储在堆栈中。
如果处理器中存在浮点硬件，可以将参数传递到NEON寄存器中。这意味着可以直接在寄存器中传递浮点参数，而不必通过堆栈或其他方式传递。

编译选项

可以通过下面的命令查看内核是否支出 NEON

zcat /proc/config.gz | grep NEON

可以通过下面指令查看处理器是否支出NEON扩展

cat /proc/cpuinfo | grep neon

在 arm 工具链中的编译选项

–vectorize 运行向量化
–cpu=Cortex-A8
-O3
-Otime 时间优化为主
-Ospace 代码空间减小为主
–fpmode=fast 优化浮点运算以向量化

例如：

armcc --cpu=Cortex-A8 -O3 -Ospace

GCC 编译器中的编译选项

-ftree-vectorize
-mcpu
-mfpu=neon
-ffast-math(或者用-Ofast 在GCC 4.6或更高的版本) 优化浮点运算以向量化
-fdump-tree-vect 和 -ftree-vectorizer-verbose=level(1-9)可以输出GCC编译信息
例如：

arm-gcc -O3 -mcpu=cortex-a15 -mfpu=neon-vfpv4 -mfloat-abi=hard -ffast-math -o 
myprog.exe myprog.c

-mcpu

-mfpu

-mfloat-abi

-mfloat-abi=soft：不使用任何FPU和NEON指令，仅使用核心寄存器集。通过库调用模拟所有浮点操作。(例如，输出.o文件)

-mfloat-abi=softfp：与-mfloat-abi=soft相同的调用约定，但根据需要使用浮点和NEON指令。使用这个选项编译的应用程序可以链接到一个软浮点库。如果相关的硬件指令可用，那么你可以使用这个选项来提高代码性能，同时仍然让代码符合软浮点环境。（例如，没有neon单元）

-mfloat-abi=hard：根据需要使用浮点和NEON指令，并且还更改了ABI调用约定，以便生成更有效的函数调用。浮点和向量类型可以在NEON寄存器之间传递给函数，从而大大减少了复制的数量。这也意味着在堆栈上传递参数时需要更少的调用。（例如，输出可执行程序，且具备neon单元）

2、Neon指令集

支持的操作

加法操作
逐对相加，将相邻的向量元素相加
带有倍增和饱和选项的乘法操作
乘法和累加操作
乘法和减法操作
左移、右移和插入操作
常见的逻辑操作（与、或、异或、与非、或非）
选择最小值和最大值的操作
计算前导零、计算有符号位和计算设置位的操作

不支持的操作

除法操作：可以使用VRECPE和VRECPS指令来执行牛顿-拉弗森迭代，从而近似实现除法操作。
平方根操作：可以使用VRSQRTE和VRSQRTS指令来执行平方根的倒数估算，然后通过乘法操作来近似实现平方根操作。

标量

neon指令集支持标量乘法，标量可以访问寄存器文件，但是只支持16位和32位的乘法运算，乘法指令只能访问前 32 个寄存器文件中的标量。

Neon指令的格式

V{mod} <op> {<shape>}{<cond>}{.<dt>}<dest1>{,{dest2}}, <src1>{,<src2>}

其中，

<mod> 是一个修饰符

Q(Saturating) 饱和修饰符，用于执行饱和运算。在饱和运算中，结果会被限制在一个特定的范围内，通常是最小值和最大值之间。这可以防止数值溢出。
H(Halving) 减半修饰符，用于将运算结果减半。这在某些情况下可以用来实现效率更高的计算。
D(Doubling) 加倍修饰符，用于将运算结果加倍。类似于减半，加倍也可以用来优化特定的计算。
R(Rounding) 舍入修饰符，用于执行舍入运算。舍入可以将浮点数结果舍入到最接近的整数值。

shape 是一个修饰符，操作数和结果的类型可能不一样

L(Long) 结果比操作数的位数更大
W(Wide) 其中一个操作数位数较小，结果和其中一个操作数类型一致
N(Narrow) 结果比操作数的位数更小

<cond> 条件，与 IT 指令一起使用
<op> 操作（ADD，SUB，MUL）
<src1>,<src2> 输入寄存器
<dest1>,<dest2> 存储寄存器
<.dt> 数据类型

cond代表条件
在这里插入图片描述

条件标志对应于不同的flag位（NZCV,参考FPSCR），如果该标志成立则执行相应的操作。

ADD r0, r1, r2; // r0 = r1 + r2, don't update flags

ADDS r0, r1, r2; // r0 = r1 + r2, and update flags

ADDSCS r0, r1, r2; // If C flag set then r0 = r1 + r2, and update flags

CMP  r0, r1; // update flags based on r0-r1.


ADD r0, r1, r2; // r0 = r1 + r2, don't update flags
 
ADDS r0, r1, r2; // r0 = r1 + r2, and update flags
 
ADDSCS r0, r1, r2; // If C flag set then r0 = r1 + r2, and update flags
 
CMP  r0, r1; // update flags based on r0-r1.

在neon intrinsics中提供了compare运算符，并将比较结果输出至目标寄存器。

Pack

NEON指令集中包含了一些打包（pack）和解包（unpack）指令，用于在打包和解包数据之间进行转换。这样可以实现从内存中加载连续的打包数据，使用字（word）或双字（doubleword）加载进行高效加载，然后解包成单独的寄存器值，进行操作，最后再将数据打包回寄存器，以便高效写回到内存中。这种方式能够提高数据在寄存器和内存之间的传输效率。

Alignment

NEON架构为NEON数据访问提供了完全的非对齐支持。然而，指令的操作码包含一个对齐提示（alignment hint），允许在地址对齐且指定了对齐提示的情况下实现更快的访问速度。然而，如果指定了对齐但地址未正确对齐，将会导致数据异常（Data Abort）。

[<Rn>:<align>]

# 例子
VLD1.8 {D0}, [R1:64]
# !表示地址偏移，常用于写回模式
VLD1.8 {D0,D1}, [R4:128]! 
VLD1.8 {D0,D1,D2,D3}, [R7:256], R2

neon将非常小的数视为 0，即flush-to-zero model

shift operations

shifting vectors

可以通过指令中编码的文字量或附加的移位向量来指定移位的数量。当使用移位向量时，应用于输入向量的每个元素的移位取决于移位向量中相应元素的值。移位向量中的元素被视为有符号值，因此可以在每个元素的基础上进行左移、右移和零移位。

shifting and inserting

NEON 的移位与插入操作提供了一种将两个向量的位组合在一起的方式。例如，左移并插入（VSLI）操作将源向量的每个元素向左移动。插入在每个元素右侧的新位来自目标向量的相应位。

Imm 表示使用立即值进行移位。Reg 表示使用寄存器指定移位量的情况。

Polynomial arithmetic over {0,1}

在 {0,1} 上对两个多项式进行乘法，首先会像常规乘法一样确定部分乘积，然后将这些部分乘积进行异或操作，而不是常规的加法。多项式乘法的结果与常规乘法不同，因为它需要对部分乘积进行多项式加法。在这里，{0,1} 上的乘法实际上是异或操作，{0,1} 上的加法实际上是异或操作(VEOR)，这是在二进制域上进行多项式运算的特点之一。这种操作在数字电路、编码和密码学等领域中常常用于处理二进制数据。如下所示：

3、Neon内联函数（intrinsics）

neon向量类型的命令语法：

<type><size>x<number_of_lanes>_t

ps:ARMv8架构已经支持float64类型

我们可以使用以下指令创建向量类型的数组：

<type><size>x<number_of_lanes>x<length_of_array>_t

将上述基本的数据组合成一个结构体构成结构化数据，通常被映射到一组向量寄存器中，例如：

struct int16x4x2_t{
    int16x4_t val[2];
} <var_name>;

这些类型仅被加载、存储、转置、交错和解交错指令使用；要对实际数据执行操作，可以选择从各个寄存器中选择元素，例如，<var_name>.val[0] 和 <var_name>.val[1]。

neon intrinsics原型

<opname><flags>_<type>

neon内联函数定义在头文件 arm_neon.h 中，该文件也定义了上述的向量数据类型

部分 intrinsics 函数

如果函数带 q 后缀代表使用 Q 寄存器，否则通常代表使用 D 寄存器

## 输入和输出向量使用 D 寄存器
uint8x8_t vadd_u8(uint8x8_t a, uint8x8_t b);

## 输入和输出向量使用 Q 寄存器
uint8x16_t vaddq_u8(uint8x16_t a, uint8x16_t b);

## 输入使用 64 位向量，输出使用 128 位向量（long 代表 shape中的Long）
uint16x8_t vaddl_u8(uint8x8_t a, uint8x8_t b);

创建内存：

vcreate_u8 # 从字面值创建一个向量
vdup_n_u8  # 创建 n 个相同的值
vget_lane_u8 # 提取单个值从向量中
vset_lan_u8 # 设置单个值为一个向量

声明变量

uint32x2_t vec64a, vec64b; // create two D-register variables

声明常量

uint8x8 start_value = vdup_n_u8(0);

uint8x8 start_value = vreinterpret_u8_u64(vcreate_u64(0x123456789ABCDEFULL));

从Q寄存器访问D寄存器

使用128位寄存器的高、低64位

vec64a = vget_low_u32(vec128); // split 128-bit vector 
vec64b = vget_high_u32(vec128); // into 2x 64-bit vectors

类型转换

uint8x8_t byteval;
uint32x2_t wordval;
byteval = vreinterpret_u8_u32(wordval);

uint8x16_t byteval2;
uint32x4_t wordval2;
byteval2 = vreinterpretq_u8_u32(wordval2);

加载交错内存

例如声道数据、RGB数据

#include <arm_neon.h>
int main (void)
{
    uint8x8x3_t v; // This represents 3 vectors. 
    // Each vector has eight lanes of 8-bit data.
    unsigned char A[24]; // This array represents a 24-bit RGB image. 
    v = vld3_u8(A); // This de-interleaves the 24-bit image from array A 
    // and stores them in 3 separate vectors
    // v.val[0] is the first vector in V. It is for the red channel
    // v.val[1] is the second vector in V. It is for the green channel
    // v.val[2] is the third vector in V. It is for the blue channel.
    //Double the red channel
    v.val[0] = vadd_u8(v.val[0],v.val[0]);
    vst3_u8(A, v); // store the vector back into the array, with the red channel doubled.
    return 0;
}

一个优化例子

#include <arm_neon.h>
uint32_t vector_add_of_n(uint32_t* ptr, uint32_t items)
{
uint32_t result,* i;
uint32x2_t vec64a, vec64b;
uint32x4_t vec128 = vdupq_n_u32(0); // clear accumulators
for (i=ptr; i<(ptr+(items/4));i+=4)
{
    uint32x4_t temp128 = vld1q_u32(i); // load four 32-bit values
    vec128=vaddq_u32(vec128, temp128); // add 128-bit vectors
}
vec64a = vget_low_u32(vec128); // split 128-bit vector
vec64b = vget_high_u32(vec128); // into two 64-bit vectors
vec64a = vadd_u32 (vec64a, vec64b); // add 64-bit vectors together
result = vget_lane_u32(vec64a, 0); // extract lanes and
result += vget_lane_u32(vec64a, 1); // add together scalars
return result;
}

以下Neon指令没有对应的内敛函数

VSWP
VLDM
VLDR
VMRS
VMSR
VPOP
VPUSH
VSTM
VSTR
VBIT
VBIF.

VSWP 指令没有本征，因为编译器可以在必要时生成 VSWP 指令。

VLDM、VLDR、VSTM 和 VSTR 主要用于上下文切换，这些指令具有对齐限制。在编写内敛函数时，使用 vldx 内敛函数更为简单。不要求对齐，除非明确指定。

VMRS 和 VMSR 访问 NEON 的条件标志。使用 NEON 内联函数进行数据处理时不需要这些标志

VPOP 和 VPUSH 用于函数的参数传递。减少变量重复使用或使用更多的
更多的 NEON 固有变量，使寄存器分配器能够跟踪活动寄存器。

5、优化

note:
vmov 指令可以在 arm 寄存器和 neon 寄存器之间移动数据，但传输的效率较慢

Result-use需要注意下个周期数据是否能计算好，例如vmla需要5个周期，大多数NEON指令具有3个周期的结果延迟

重用变量

下面函数是一个 4x4 矩阵的乘法，重用了变量 r，导致后面的计算依赖于前面的变量，则后面的指令不得不等待前面指令计算完成

void altneonmult(const float *matrixA, const float *matrixB, float *matrixR)
{
float32x4_t a, b0, b1, b2, b3, r;
a0 = vld1q_f32(matrixA); /* col 0 of matrixA */
a1 = vld1q_f32(matrixA + 4); /* col 1 of matrixA */
a2 = vld1q_f32(matrixA + 8); /* col 2 of matrixA */
a3 = vld1q_f32(matrixA + 12); /* col 3 of matrixA */
b = vld1q_f32(matrixB); /* load col 0 of matrixB */
r = vmulq_lane_f32(a0, vget_low_f32(b), 0);
r = vmlaq_lane_f32(r, a1, vget_low_f32(b), 1);
r = vmlaq_lane_f32(r, a2, vget_high_f32(b), 0);
r = vmlaq_lane_f32(r, a3, vget_high_f32(b), 1);
vst1q_f32(matrixR, r); /* store col 0 of result */
b = vld1q_f32(matrixB + 4); /* load col 1 of matrixB */
r = vmulq_lane_f32(a0, vget_low_f32(b), 0);
r = vmlaq_lane_f32(r, a1, vget_low_f32(b), 1);
r = vmlaq_lane_f32(r, a2, vget_high_f32(b), 0);
r = vmlaq_lane_f32(r, a3, vget_high_f32(b), 1);
vst1q_f32(matrixR + 4, r); /* store col 1 of result */
b = vld1q_f32(matrixB + 8); /* load col 2 of matrixB */
r = vmulq_lane_f32(a0, vget_low_f32(b), 0);
r = vmlaq_lane_f32(r, a1, vget_low_f32(b), 1);
r = vmlaq_lane_f32(r, a2, vget_high_f32(b), 0);
r = vmlaq_lane_f32(r, a3, vget_high_f32(b), 1);
vst1q_f32(matrixR + 8, r); /* store col 2 of result */
b = vld1q_f32(matrixB + 12); /* load col 3 of matrixB */
r = vmulq_lane_f32(a0, vget_low_f32(b), 0);
r = vmlaq_lane_f32(r, a1, vget_low_f32(b), 1);
r = vmlaq_lane_f32(r, a2, vget_high_f32(b), 0);
r = vmlaq_lane_f32(r, a3, vget_high_f32(b), 1);
vst1q_f32(matrixR + 12, r); /* store col 3 of result */
}

多定义几个变量，消除了数据依赖，在计算之前，就可以完成数据加载。

void neonmult(const float *matrixA, const float *matrixB, float *matrixR)
{
float32x4_t a0, a1, a2, a3, b0, b1, b2, b3, r0, r1, r2, r3;
a0 = vld1q_f32(matrixA); /* col 0 of matrixA */
a1 = vld1q_f32(matrixA + 4); /* col 1 of matrixA */
a2 = vld1q_f32(matrixA + 8); /* col 2 of matrixA */
a3 = vld1q_f32(matrixA + 12); /* col 3 of matrixA */

b0 = vld1q_f32(matrixB); /* col 0 of matrixB */
b1 = vld1q_f32(matrixB + 4); /* col 1 of matrixB */
b2 = vld1q_f32(matrixB + 8); /* col 2 of matrixB */
b3 = vld1q_f32(matrixB + 12); /* col 3 of matrixB */

/* compute all the cols in the order specified by compiler */
r0 = vmulq_lane_f32(a0, vget_low_f32(b0), 0);
r0 = vmlaq_lane_f32(r0, a1, vget_low_f32(b0), 1);
r0 = vmlaq_lane_f32(r0, a2, vget_high_f32(b0), 0);
r0 = vmlaq_lane_f32(r0, a3, vget_high_f32(b0), 1);

r1 = vmulq_lane_f32(a0, vget_low_f32(b1), 0);
r1 = vmlaq_lane_f32(r1, a1, vget_low_f32(b1), 1);
r1 = vmlaq_lane_f32(r1, a2, vget_high_f32(b1), 0);
r1 = vmlaq_lane_f32(r1, a3, vget_high_f32(b1), 1);

r2 = vmulq_lane_f32(a0, vget_low_f32(b2), 0);
r2 = vmlaq_lane_f32(r2, a1, vget_low_f32(b2), 1);
r2 = vmlaq_lane_f32(r2, a2, vget_high_f32(b2), 0);
r2 = vmlaq_lane_f32(r2, a3, vget_high_f32(b2), 1);

r3 = vmulq_lane_f32(a0, vget_low_f32(b3), 0);
r3 = vmlaq_lane_f32(r3, a1, vget_low_f32(b3), 1);
r3 = vmlaq_lane_f32(r3, a2, vget_high_f32(b3), 0);
r3 = vmlaq_lane_f32(r3, a3, vget_high_f32(b3), 1);

vst1q_f32(matrixR, r0);
vst1q_f32(matrixR + 4, r1);
vst1q_f32(matrixR + 8, r2);
vst1q_f32(matrixR + 12, r3);
}

vld and vst语法格式

元素的大小会影响指针的对齐(内存地址能够整除元素大小)。通常情况下，按照元素的大小进行对齐会提供更好的性能，并且这可能是目标操作系统的要求。例如，当加载32位元素时，应将第一个元素的地址对齐到32位。这种对齐方式有助于提高内存访问的效率，特别是在处理涉及硬件向量化指令集（如NEON）的情况下。在许多体系结构中，对齐的数据访问能够提高数据传输速度和处理器操作的效率。

除加载RGB数据外，如果arm寄存器r0中的值一样，会有另外两种用法：

在这里插入图片描述

arm寄存器的一些结构化操作

[ $R_{n}$ :{,:align}] 按指定对齐格式加载数据
[ $R_{n}$ :{,:align}]！寄存器地址在被load后向后移动（移动位数等于加载的数据大小）
[ $R_{n}$ :{,:align}], $R_{m}$
在内存访问之后，指令会将指针增加由寄存器 Rm 中的值表示的偏移量。这在读取或写入一组按固定宽度分隔的元素时非常有用。一个示例是从图像中读取垂直线的数据。通过使用寄存器 Rm 来指定偏移量，可以有效地读取或写入一系列数据元素，从而避免手动计算和操作内存指针。

其他的加载和存储指令

VLDR and VSTR : 加载或存储单个寄存器为64位值。
VLDM and VSTM : 将多个寄存器加载为64位值。这对于从堆栈中存储和检索寄存器很有用。

数组元素不是2的幂次倍

解决方案有：

padding，补白的位置可以初始为对计算结果没有影响的值。缺点：需要过多的内存；补白的值可能会造成结果错误
overlapping，该方法可以用于查找数组最大最小元素等操作
单独处理，对于未对齐的向量可以采用该方法，缺点：速度慢；代码多

矩阵乘操作

floating-point

矩阵乘操作，看下图，可以分解为以下操作

以列为主存，y0被第一列的所有元素乘了一遍:

上图的代码实现

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

vmul.f32 q12, q8, d0[0] @ multiply element 0 (y0) by matrix col 0 (x0-x3)
vmla.f32 q12, q9, d0[1] @ multiply-acc element 1 (y1) by matrix col 1 (x4-x7)
vmla.f32 q12, q10, d1[0] @ multiply-acc element 2 (y2) by matrix col 2 (x8-xb)
vmla.f32 q12, q11, d1[1] @ multiply-acc element 3 (y3) by matrix col 3 (xc-xf)

note:
标量的存储一般使用 D 寄存器，因为 GNU 不接受使用 Q 寄存器去读取一个标量。

完整的实现如下：

首先创建一个宏定义，用于行列相乘（即封装了上面的代码）

.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 col2
vmla.f32 \res_q, q11, \col1_d[1] @ multiply-acc col element 3 by matrix col 3
.endm

调用宏定义函数，执行计算流程

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

fixed point

定点运算一般更快，因为 CPU 可以使用较少的内存带宽，访问较小位数的数据，但是需要注意：

认真选择表示方法避免饱和或溢出
保持应用程序所要求的结果的精确度

宏定义函数，注意两个 16 位相乘，数值可能会扩大到 32 位，所有要使用 vmull （Long）：

.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

Q1.14 定点数的格式如下：

其中，VQRSHRN的意思为：

Q 饱和度转换
R 四舍五入
N narrow，计算结果的数据类型位数减少
SHR 右移

调用宏定义函数的代码：

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)

cross product

$i, j, k$ 是相互正交的单位向量，其中： $a=[a_{i},a_{j},a_{k}],b=[b_{i},b_{j},b_{k}]$ ,交叉积表示为：

$\times b = i(a_{j}b_{k} - a_{k}b_{j}) + j(a_{k}b_{i} - a_{i}b_{k}) + k(a_{i}b{j} - a_{j}b_{i})$

结构向量为：

$a_{j}b_{k} - a_{k}b_{j}, a_{k}b_{i} - a_{i}b_{k}, a_{i}b{j} - a_{j}b_{i}]$

代码实现如下：（注释有问题，向量最后一个元素为起始位置)

void cross_product_s(float32_t *r, float32_t* a, float32_t* b) 
{
// Vector a is stored in memory such that ai is at the lower address and
// ak is at the higher address. Vector b is also stored in the same way.
// Registers are shown as vectors containing elements, for example:
// [element3, element2, element1, element0] where element3 uses the higher bits
// and element0 uses the lower bits of the register.
float32x2_t vec_a_1 = vld1_f32(a + 1); //D register = [ak, aj]
float32x2_t vec_a_2 = vld1_f32(a); //D register = [aj, ai]
float32x2_t vec_b_1 = vld1_f32(b + 1); //D register = [bk, bj]
float32x2_t vec_b_2 = vld1_f32(b); //D register = [bj, bi]
float32x4_t vec_a = vcombine_f32(vec_a_1, vec_a_2); //Q register = [aj, ai, ak, aj]
float32x4_t vec_b = vcombine_f32(vec_b_1, vec_b_2); //Q register = [bj, bi, bk, bj]
float32x4_t vec_a_rot = vextq_f32(vec_a, vec_a, 1);
float32x4_t vec_b_rot = vextq_f32(vec_b, vec_b, 1);
// vec_a = [ aj, ai, ak, aj ]
// vec_b_rot = [ bj, bj, bi, bk ]
// vec_a_rot = [ aj, aj, ai, ak ]
// vec_b = [ bj, bi, bk, bj ]
float32x4_t prod = vmulq_f32(vec_a, vec_b_rot);
// prod = [ ajbj, aibj, akbi, ajbk ]
prod = vmlsq_f32(prod, vec_a_rot, vec_b);
// prod = [ ajbj-ajbj, aibj-ajbi, akbi-aibk, ajbk-akbj ]
vst1_f32(r, vget_low_f32(prod)); // Store the lower two elements to address r
vst1_lane_f32(r + 2, vget_high_f32(prod), 0); // Store the 3rd element
}

使用 vld3q_f32() 和 vst3q_f32() 可以更高效的计算多维叉乘：

void cross_product_q(float32_t* r, float32_t* a, float32_t* b) 
{
float32x4x3_t vec_a = vld3q_f32(a);
float32x4x3_t vec_b = vld3q_f32(b);
float32x4x3_t result;
result.val[0] = vmulq_f32(vec_a.val[1], vec_b.val[2]);
result.val[0] = vmlsq_f32(result.val[0], vec_a.val[2], vec_b.val[1]);
result.val[1] = vmulq_f32(vec_a.val[2], vec_b.val[0]);
result.val[1] = vmlsq_f32(result.val[1], vec_a.val[0], vec_b.val[2]);
result.val[2] = vmulq_f32(vec_a.val[0], vec_b.val[1]);
result.val[2] = vmlsq_f32(result.val[2], vec_a.val[1], vec_b.val[0]);
vst3q_f32(r, result);
}

以下代码可以处理任意长度向量的叉乘运算：

void cross_product(float32_t* r, float32_t* a, float32_t* b, int count) 
{
int count4 = count / 4;
count &= 3;
while(count4 > 0) 
{
cross_product_q(r, a, b);
r += 12; a += 12; b += 12; count4--;
}
if(count >= 2) 
{
cross_product_d(r, a, b);
r += 6; a += 6; b += 6; count -= 2;
}
if(count == 1) 
{
cross_product_s(r, a, b);
}
}

Neon Code优化实例

RGB565 转换到 RGB888

图像位深转换是一种常见的操作，对于 RGB565 格式的图像数据，neon单元可以通过偏移操作高效的处理

下面的代码是 RGB565 格式类型通过 neon 偏移转换为 RGB888 格式的数据，其中该格式的数据存在 d2®,d3(G),d4(B)寄存器中

vshr.u8 q1, q0, #3  @ shift red elements right by three bits,
                    @ discarding the green bits at the bottom of
                    @ the red 8-bit elements.
vshrn.i16 d2, q1, #5    @ shift red elements right and narrow,
                        @ discarding the blue and green bits.
vshrn.i16 d3, q0, #5    @ shift green elements right and narrow,
                        @ discarding the blue bits and some red bits
                        @ due to narrowing.
vshl.i8 d3, d3, #2  @ shift green elements left, discarding the
                    @ remaining red bits, and placing green bits
                    @ in the correct place.
vshl.i16 q0, q0, #3 @ shift blue elements left to most-significant
                    @ bits of 8-bit color channel.
vmovn.i16 d4, q0    @ remove remaining red and green bits by
                    @ narrowing to 8 bits

但是该转换有一个缺点：白色的表示会存在问题，原因在于在 RGB565 格式中白色的表示是[0x1F,0x3F,0x1F],上述转换会将这段数字转换为[0xF8,0xFC,0xF8],这可以使用移位操作修复（将一些重要的数据放置在低位）

使用内联函数优化上述代码：

uint16_t *src = image_src;
uint8_t *dst = image_dst;
int count = PIXEL_NUMBER;
while (count >= 8) {
uint16x8_t vsrc;
uint8x8x3_t vdst;
vsrc = vld1q_u16(src);
vdst.val[0] = 
vshrn_n_u16(vreinterpretq_u16_u8(vshrq_n_u8(vreinterpretq_u8_u16(vsrc), 3)), 5);
vdst.val[1] = vshl_n_u8(vshrn_n_u16(vsrc, 5) ,2);
vdst.val[2] = vmovn_u16(vshlq_n_u16(vsrc, 3));
vst3_u8(dst, vdst);
dst += 8*3;
src += 8;
count -= 8;
}

编译后的指令为：

vshr.u8 q1, q0, #3
vshrn.i16 d3, q0, #5
vshl.i16 q0, q0, #3
vshrn.i16 d2, q1, #5
vshl.i8 d3, d3, #2
vmovn.i16 d4, q0

RGB88 转换到 RGB565

vshll.u8 q2, d0, #8 @ shift red elements left to most-significant
                    @ bits of wider 16-bit elements.
vshll.u8 q3, d1, #8 @ shift green elements left to most-significant
                    @ bits of wider 16-bit elements.
vsri.16 q2, q3, #5 @ shift green elements right and insert into
                    @ red elements.
vshll.u8 q3, d2, #8 @ shift blue elements left to most-significant
                    @ bits of wider 16-bit elements.
vsri.16 q2, q3, #11 @ shift blue elements right and insert into
                    @ red and green elements.

使用内联函数优化上述代码：

uint8_t *src = image_src;
uint16_t *dst = image_dst;
int count = PIXEL_NUMBER;
while (count >= 8) {
uint8x8x3_t vsrc;
uint16x8_t vdst;
vsrc = vld3_u8(src);
vdst = vshll_n_u8(vsrc.val[0], 8);
vdst = vsriq_n_u16(vdst, vshll_n_u8(vsrc.val[1], 8), 5);
vdst = vsriq_n_u16(vdst, vshll_n_u8(vsrc.val[2], 8), 11);
vst1q_u16(dst, vdst);
dst += 8;
src += 8*3;
count -= 8;
}

经过编译后的代码：

vshll.u8 q2, d0, #8
vshll.u8 q3, d1, #8
vshll.u8 q8, d2, #8
vsri.16 q2, q3, #5
vsri.16 q2, q8, #11

中值滤波

Bitonic排序

Bitonic序列是一个非严格增非严格减的序列，Bitnoic是一个基于比较的排序。将任意一个长为2n的双调序列A分为等长的两半X和Y，将X中的元素与Y中的元素一一按原序比较，即a[i]与a[i+n] (i < n)比较，将较大者放入MAX序列，较小者放入MIN序列。则得到的MAX和MIN序列仍然是双调序列，并且MAX序列中的任意一个元素不小于MIN序列中的任意一个元素。如图所示：

在neon中的实现流程如下图所示：

static VFP inline uint16x8x2_t bitonic_merge_16(uint16x8_t a, uint16x8_t b)
{
b = vrev128q_u16(b);
vminmaxq(a, b);
return bitonic_resort_16(a, b);
}
static VFP inline uint16x8x2_t bitonic_resort_16(uint16x8_t a, uint16x8_t b)
{
/* Start with two vectors:
* +---+---+---+---+---+---+---+---+
* | a | b | c | d | e | f | g | h |
* +---+---+---+---+---+---+---+---+
* +---+---+---+---+---+---+---+---+
* | i | j | k | l | m | n | o | p |
* +---+---+---+---+---+---+---+---+
* All the elements of the first are guaranteed to be less than or equal to
* all of the elements in the second, and both vectors are bitonic.
* We need to perform these operations to completely sort both lists:
* vminmax([abcd],[efgh]) vminmax([ijkl],[mnop])
* vminmax([abef],[cdgh]) vminmax([ijmn],[klop])
* vminmax([acef],[bdfh]) vminmax([ikmo],[jlmp])
* We can align the necessary pairs for mixing with transpose operations,
* like so:
*/
vtrn64q(a ,b);
/* We now have:
* +---+---+---+---+---+---+---+---+
* | a | b | c | d | i | j | k | l |
* +---+---+---+---+---+---+---+---+
* +---+---+---+---+---+---+---+---+
* | e | f | g | h | m | n | o | p |
* +---+---+---+---+---+---+---+---+
*/
vminmaxq(a, b);
vtrn32q(a, b);
/* We now have:
* +---+---+---+---+---+---+---+---+
* | a | b | e | f | i | j | m | n |
* +---+---+---+---+---+---+---+---+
* +---+---+---+---+---+---+---+---+
* | c | d | g | h | k | l | o | p |
* +---+---+---+---+---+---+---+---+
*/
vminmaxq(a, b);
vtrn16q(a, b);
/* We now have:
* +---+---+---+---+---+---+---+---+
* | a | c | e | g | i | k | m | o |
* +---+---+---+---+---+---+---+---+
* +---+---+---+---+---+---+---+---+
* | b | d | f | h | j | l | n | p |
* +---+---+---+---+---+---+---+---+
*/
vminmaxq(a, b);
/* Since we now have separate vectors of odd and even lanes, we can simply
* use vzip to stitch them back together producing our original
* positioning. Then we're done.
*/
return vzipq_u16(a, b);
}
/* Macros */
/* In-place swap top 64 bits of “a” with bottom 64 bits of “b” -- one operation (vswp)
*/
#define vtrn64q(a, b) \
do { \
uint16x4_t vtrn64_tmp_a0 = vget_low_u16(a), vtrn64_tmp_a1 = vget_high_u16(a); \
uint16x4_t vtrn64_tmp_b0 = vget_low_u16(b), vtrn64_tmp_b1 = vget_high_u16(b); \
{ \
uint16x4_t vtrn64_tmp = vtrn64_tmp_a1; \
vtrn64_tmp_a1 = vtrn64_tmp_b0; \
vtrn64_tmp_b0 = vtrn64_tmp;\
} \
(a) = vcombine_u16(vtrn64_tmp_a0, vtrn64_tmp_a1); \
(b) = vcombine_u16(vtrn64_tmp_b0, vtrn64_tmp_b1); \
} while (0)
/* In-place exchange odd 32-bit words of “a” with even 32-bit words of “b” -- one 
operation
*/
#define vtrn32q(a, b) \
do { \
uint32x4x2_t vtrn32_tmp = \
vtrnq_u32(vreinterpretq_u32_u16(a), vreinterpretq_u32_u16(b)); \
(a) = vreinterpretq_u16_u32(vtrn32_tmp.val[0]); \
(b) = vreinterpretq_u16_u32(vtrn32_tmp.val[1]); \
} while (0)
#define vtrn32(a, b) 
do { \
uint32x2x2_t vtrn32_tmp = vtrn_u32(vreinterpret_u32_u16(a),vreinterpret_u32_u16(b));\
(a) = vreinterpret_u16_u32(vtrn32_tmp.val[0]); \
(b) = vreinterpret_u16_u32(vtrn32_tmp.val[1]); \
} while (0)
/* In-place exchange odd 16-bit words of “a” with even 16-bit words of “b” -- one 
operation */
#define vtrn16q(a, b) \
do { \
uint16x8x2_t vtrn16_tmp = vtrnq_u16((a), (b)); \
(a) = vtrn16_tmp.val[0]; \
(b) = vtrn16_tmp.val[1]; \
} while (0)
#define vzipq(a, b) \
do { \
uint16x8x2_t vzip_tmp = vzipq_u16((a), (b)); \
(a) = vzip_tmp.val[0]; \
(b) = vzip_tmp.val[1]; \
} while (0)
#define vminmaxq(a, b) \
do { \
uint16x8_t minmax_tmp = (a); \
(a) = vminq_u16((a), (b)); \
(b) = vmaxq_u16(minmax_tmp, (b)); \
} while (0)
#define vrev128q_u16(a) \
vcombine_u16(vrev64_u16(vget_high_u16(a)), vrev64_u16(vget_low_u16(a)))

GRB Bitonic Sort

sort

奇数 bitonic 排序的流程图（不是标志的 bitonic 排序），前 4 个元素已经排序好，后 3 个元素已经排序好，然后两部分进行比较交换（可以确定最小元素），此时 2，3，4，5 号元素的大小未知，经过比较可以获得这四个元素中的最大值和最小值，然后 1和2， 3和4， 5和6 大小未知，比较交换就行了。
在这里插入图片描述

void loadblock(uint16_t dst[8][8], uint16_t const *src, int spitch)
{
uint16x8_t q0, q1, q2, q3, q4, q5, q6, q7;
spitch >>= 1;
q0 = vld1q_u16(src); src += spitch;
q1 = vld1q_u16(src); src += spitch;
q2 = vld1q_u16(src); src += spitch;
q3 = vld1q_u16(src); src += spitch;
q4 = vld1q_u16(src); src += spitch;
q5 = vld1q_u16(src); src += spitch;
q6 = vld1q_u16(src);
/* vminmaxq() performs a pair of vminq and vmaxq operations to put the two arguments in 
order */
vminmaxq(q0, q1);
vminmaxq(q2, q3);
vminmaxq(q4, q5);

vminmaxq(q0, q2);
vminmaxq(q1, q3);
vminmaxq(q4, q6);

vminmaxq(q1, q2);
vminmaxq(q5, q6);

vminmaxq(q0, q4);
vminmaxq(q1, q5);
vminmaxq(q2, q6);

vminmaxq(q2, q4);
vminmaxq(q3, q5);

vminmaxq(q1, q2);
vminmaxq(q3, q4);
vminmaxq(q5, q6);
}

tranpose

为了方便处理，padding一行，其值填充为 99

q7 = vdupq_n_u16(UINT16_MAX);
vzipq(q0, q1);
vzipq(q2, q3);
vzipq(q4, q5);
vzipq(q6, q7);
uint32x4x4_t tmp;
tmp.val[0] = vreinterpretq_u32_u16(q0); tmp.val[1] = vreinterpretq_u32_u16(q2);
tmp.val[2] = vreinterpretq_u32_u16(q4); tmp.val[3] = vreinterpretq_u32_u16(q6);
vst4q_u32((uint32_t *)&dst[0], tmp);
tmp.val[0] = vreinterpretq_u32_u16(q1); tmp.val[1] = vreinterpretq_u32_u16(q3);
tmp.val[2] = vreinterpretq_u32_u16(q5); tmp.val[3] = vreinterpretq_u32_u16(q7);
vst4q_u32((uint32_t *)&dst[4], tmp);

find median

/* Take a set of sorted lists (32, 16, and 8 elements) and return the median.
*
* We arrange for the second and third lists to be combined in reverse order so
* that they merge directly with the first list without reversal.
*
* We return the first list unmodified because it means the caller can write it
* back to its origin and not worry about saving that data for subsequent
* calls.
*/
static VFP uint16x8x4_t bitonic_median_56(uint16_t *dst, uint16x8x4_t a, 
uint16x8x2_t b0, uint16x8_t b1)
{
uint16x8x3_t b = bitonic_merge_24r(b0, b1); /* MUST inline */
b.val[0] = vminq_u16(a.val[1], b.val[0]);
b.val[1] = vminq_u16(a.val[2], b.val[1]);
b.val[2] = vminq_u16(a.val[3], b.val[2]);
b.val[1] = vmaxq_u16(a.val[0], b.val[1]);
b.val[2] = vmaxq_u16(b.val[0], b.val[2]);
b.val[2] = vmaxq_u16(b.val[1], b.val[2]);
{
uint16x4_t tmp;
tmp = vmin_u16(vget_low_u16(b.val[2]), vget_high_u16(b.val[2]));
tmp = vpmin_u16(tmp, tmp);
tmp = vpmin_u16(tmp, tmp);
vst1_lane_u16(dst, tmp, 0);
}
return a;
}
static VFP inline uint16x8x3_t bitonic_merge_24r(uint16x8x2_t a, uint16x8_t b)
{
uint16x8x3_t result;
b = vrev128q_u16(b);
vmaxminq(a.val[0], b);
vmaxminq(a.val[0], a.val[1]);
a = bitonic_resort_16r(a.val[0], a.val[1]);
vmaxmin_half(b);
b = bitonic_resort_8r(b);
result.val[0] = a.val[0];
result.val[1] = a.val[1];
result.val[2] = b;
return result;
}
static VFP inline uint16x8_t bitonic_resort_8r(uint16x8_t a)
{
uint16x4_t a0 = vget_low_u16(a), a1 = vget_high_u16(a);
vtrn32(a0, a1);
vmaxmin(a0, a1);
vtrn16(a0, a1);
vmaxmin(a0, a1);
vzip(a0, a1);
return vcombine_u16(a0, a1);
}
/* Macros */
#define vmaxminq(a, b) \
do { \
uint16x8_t maxmin_tmp = (a); \
(a) = vmaxq_u16((a), (b)); \
(b) = vminq_u16(maxmin_tmp, (b)); \
} while (0)
#define vmaxmin_half(a) \
do { \
uint16x4_t minmax_tmp_lo = vget_low_u16(a), minmax_tmp_hi = vget_high_u16(a); \
vmaxmin(minmax_tmp_lo, minmax_tmp_hi); \
(a) = vcombine_u16(minmax_tmp_lo, minmax_tmp_hi); \
} while (0)
#define vmaxmin(a, b) \
do { \
uint16x4_t maxmin_tmp = (a); \
(a) = vmax_u16((a), (b)); \
(b) = vmin_u16(maxmin_tmp, (b)); \
} while (0)

re-use

void filter_row_bs(uint16_t *dst, uint16_t const *src, int spitch, int count)
{
    uint16_t scratch[24][8];
    struct
    {
        uint16x8x2_t ab, ef;
        uint16x8_t b, d, f, h;
    } state[3];
    int i;
    loadblock(scratch + 0, src + 0, spitch);
    loadblock(scratch + 8, src + 8, spitch);
    loadblock(scratch + 16, src + 16, spitch);
    for (i = 0; i < 3; i++)
    {
    uint16x8_t c, d, e, f;
    c = vld1q_u16(scratch[i + 3]);
    d = vld1q_u16(scratch[i + 6]);
    e = vld1q_u16(scratch[i + 9]);
    f = vld1q_u16(scratch[i + 12]);
    state[i].ab = bitonic_merge_16(e, f);
    state[i].ef = bitonic_merge_16(c, d);
    state[i].d = d;
    state[i].f = f;
    state[i].b = vld1q_u16(scratch[i + 0]);
    state[i].h = vld1q_u16(scratch[i + 15]);
    }
    src += 18;
    /* We maintain three sets (components) of partially combined lists and we recombine
    * these for a full bitonic sort in various permutations in an eight-phase rota.
    *
    * Register allocation collapses in a heap, here, and we spend a lot of our
    * time marshalling things through the stack. We could probably do much
    * better with hand-written assembly.
    */
    while (count > 0)
    {
        int i;
        uint16_t const *pft, *pft2;
        for (i = 0; i < 3; i++)
        {
        loadblock(scratch + i * 8, src, spitch);
        src += 8;
        }
        pft = src;
        if (count < 8) pft -= 24;
        pft2 = pft + 6 * spitch;
        for (i = 0; i < 3; i++, dst++)
        {
            int cnt = count;
            uint16x8x2_t ab, ef;
            uint16x8_t h;
            uint16x8_t x, y;
            uint16x8x2_t xx;
            uint16x8x4_t xxxx;
            ef = state[i].ef;
            ab = state[i].ab;

            /* stage 0 */
            x = vld1q_u16(scratch[i + 0]);
            __builtin_prefetch(pft + 0);
            xx = bitonic_merge_16(state[i].h, x);
            xxxx = bitonic_merge_32(ab, xx);
            h = x;
            ab = xx;
            xxxx = bitonic_median_56(dst + 0, xxxx, ef, state[i].b);
            if (--cnt <= 0)
            continue;

            /* stage 1 */
            y = vld1q_u16(scratch[i + 3]);
            __builtin_prefetch(pft + 16);
            xxxx = bitonic_median_56(dst + 3, xxxx, ef, y);
            if (--cnt <= 0)
            continue;

            /* stage 2 */
            x = vld1q_u16(scratch[i + 6]);
            __builtin_prefetch(pft + 23); pft += spitch;
            xx = bitonic_merge_16(y, x);
            state[i].b = x;
            ef = xx;
            xxxx = bitonic_median_56(dst + 6, xxxx, ef, state[i].d);
            if (--cnt <= 0)
            continue;

            /* stage 3 */
            y = vld1q_u16(scratch[i + 9]);
            __builtin_prefetch(pft + 0);
            (void)bitonic_median_56(dst + 9, xxxx, ef, y);
            if (--cnt <= 0)
            continue;

            /* stage 4 */
            x = vld1q_u16(scratch[i + 12]);
            __builtin_prefetch(pft + 16);
            xx = bitonic_merge_16(y, x);
            xxxx = bitonic_merge_32(ef, xx);
            state[i].d = x;
            state[i].ef = xx;
            xxxx = bitonic_median_56(dst + 12, xxxx, ab, state[i].f);
            if (--cnt <= 0)
            continue;

            /* stage 5 */
            y = vld1q_u16(scratch[i + 15]);
            __builtin_prefetch(pft + 23); pft += spitch;
            xxxx = bitonic_median_56(dst + 15, xxxx, ab, y);
            if (--cnt <= 0)
            continue;

            /* stage 6 */
            x = vld1q_u16(scratch[i + 18]);
            __builtin_prefetch(pft2); pft2 += 8;
            xx = bitonic_merge_16(y, x);
            state[i].f = x;
            ab = xx;
            xxxx = bitonic_median_56(dst + 18, xxxx, ab, h);
            if (--cnt <= 0)
            continue;

            /* stage 7 */
            y = vld1q_u16(scratch[i + 21]);
            state[i].h = y;
            state[i].ab = ab;
            (void)bitonic_median_56(dst + 21, xxxx, ab, y);
        }
        dst += 24 - 3;
        count -= 8;
    }
}

FIR filter

有限脉冲响应（FIR）滤波器是一个常见的可矢量化循环示例。代码如下：

/*fir.c*/ 
void fir(short * y, const short *x, const short *h, int n_out, int n_coefs) 
{ 
    int n; 
    for (n = 0; n < n_out; n++) 
    { 
        int k, sum = 0; 
        for(k = 0; k < n_coefs; k++) 
        { 
            sum += h[k] * x[n - n_coefs + 1 + k]; 
        } 
        y[n] = ((sum>>15) + 1) >> 1; 
    } 
}

上述内循环代码中是一个矩阵乘加运算，可以被矢量化，但是循环次数 $n\_coefs$ 未知，但是我们可以对其除 4：

for(k = 0; k < n_coefs/4; k++) 
{ 
    sum0 += h[k*4] * x[n - n_coefs + 1 + k*4]; 
    sum1 += h[k*4+1] * x[n - n_coefs + 1 + k*4 + 1]; 
    sum2 += h[k*4+2] * x[n - n_coefs + 1 + k*4 + 2]; 
    sum3 += h[k*4+3] * x[n - n_coefs + 1 + k*4 + 3]; 
} 
sum = sum0 + sum1 + sum2 + sum3;

if (n_coefs % 4) 
{ 
    for(k = n_coefs - (n_coefs % 4); k < n_coefs; k++)
    { 
        sum += h[k] * x[n - n_coefs + 1 + k]; 
    } 
}

使用 neon 指令进行优化：

#include <arm_neon.h> 
void fir(short * y, const short *x, const short *h, int n_out, int n_coefs) 
{
    int n, k; 
    int sum; 
    int16x4_t h_vec; 
    int16x4_t x_vec; 
    int32x4_t result_vec; 
    for (n = 0; n < n_out; n++) 
    { /* Clear the scalar and vector sums */ 
        sum = 0; 
        result_vec = vdupq_n_s32(0); 
        for(k = 0; k < n_coefs / 4; k++) 
        { /* Four vector multiply-accumulate operations in parallel */ 
            h_vec = vld1_s16(&h[k*4]); 
            x_vec = vld1_s16(&x[n - n_coefs + 1 + k*4]); 
            result_vec = vmlal_s16(result_vec, h_vec, x_vec); 
        } 
        /* Reduction operation - add each vector lane result to the sum */ 
        sum += vgetq_lane_s32(result_vec, 0); 
        sum += vgetq_lane_s32(result_vec, 1); 
        sum += vgetq_lane_s32(result_vec, 2); 
        sum += vgetq_lane_s32(result_vec, 3); 
        /* consume the last few data using scalar operations */ 
        if(n_coefs % 4) 
        { 
            for(k = n_coefs - (n_coefs % 4); k < n_coefs; k++) 
                sum += h[k] * x[n - n_coefs + 1 + k]; 
        }
        /* Store the adjusted result */ 
        y[n] = ((sum>>15) + 1) >> 1; 
    }
}

如果我们明确知道 $n\_coefs$ 是4的倍数，那么可以使用以下方法矢量化：

for (n = 0; n < ((n_coefs/4)*4); n++)

or

for (n = 0; n < n_coefs_by_4*4; n++)

其中， $n\_coefs\_by\_4 = n\_coefs / 4$ ，矢量化代码如下：

void fir(short * y, const short *x, const short *h, int n_out, int n_coefs_by_4) 
{
    int n; 
    for (n = 0; n < n_out; n++) 
    { 
        int k, sum = 0; 
        for(k = 0; k < n_coefs_by_4*4; k++) 
        { 
            sum += h[k] * x[n - n_coefs_by_4*4 + 1 + k]; 
        } 
        y[n] = ((sum>>15) + 1) >> 1; 
    } 
}