ARMv8入门且实用的C语言汇编混合编程案例：NEON SIMD优化的通用矩阵乘

最新推荐文章于 2023-12-24 09:22:33 发布

withinmiaov

最新推荐文章于 2023-12-24 09:22:33 发布

阅读量3.8k

点赞数 2

文章标签： c语言矩阵 asm arm simd

本文链接：https://blog.csdn.net/qq_41933225/article/details/122402064

版权

1.简介

在ARMv8-A架构下，通过矩阵乘这样一个在图像处理、神经网络等各种领域都十分常见的算法，来帮助大家了解如何实现汇编与C语言混合编程。同时也会使用到NEON SIMD向量化技术，实现矩阵乘的性能优化。

实现汇编与C语言混合编程的方式有两种，一种是C语言中嵌入汇编的方式，适用于汇编代码量较小的情况，另一种是C语言程序调用汇编程序函数的方式。本案例使用C语言程序调用汇编程序函数的方式进行混合编程。

在目前的计算机领域中，通用矩阵乘（GEMM）可以说是随处可见，优化矩阵乘的性能对于很多应用程序都十分有价值，所以一直都有许多关于矩阵乘性能优化的研究。目前已经有非常成熟的高性能通用矩阵乘算法，而稀疏矩阵乘成为了新的研究热门。

相信大多数人都对矩阵乘有一定的了解，这里用4*4矩阵相乘来简要说明矩阵乘算法原理：

矩阵乘最简单的C语言实现方式（最基本的三层循环方式）：

for (int m = 0; m < M_len; m++)
	for (int n = 0; n < N_len; n++)
		for (int k = 0; k < K_len; k++)
			maxtrix_result[m * N_len + n] += maxtrix_A[m * K_len + k] * maxtrix_B[k * N_len + n];

这样一个最为原始的三重循环矩阵乘实现性能很低，我们想通过C语言与汇编语言（汇编调用NEON SIMD）混合编程的方式来优化矩阵乘的性能，实现一个稍简易版的性能更好的通用矩阵乘，即实现result[m,n] = matrix1[m,k] * matrix2[k,n]。

2.C语言汇编混合编程案例：双层矩阵乘

矩阵由单精度浮点数（float）组成，即矩阵中的每个数都是32位的。matrix1是行数为m、列数为k的矩阵，matrix2是行数为k、列数为n的矩阵，它们相乘的结果为result矩阵，行数为m，列数为n。

为了简化案例难度，我们将m、n、k规定为4的倍数。因为NEON向量寄存器长度为128位，可以放得下4个32位float浮点数。

本案例目标实现的矩阵乘算法分为两层，首先通过汇编语言实现4*4的矩阵相乘函数，再通过C语言调用汇编函数的方式实现第二层矩阵乘，通过这样两层的设计实现任意规模大小的矩阵乘，同时也通过汇编及NEON SIMD向量化优化矩阵乘程序的性能。

2.1外层C语言矩阵乘（调用汇编函数）

以16*16规模的矩阵乘作为示例说明，如图所示，如果将4*4的小矩阵看作一个单元，那么这个大矩阵也就是一个由小矩阵构成的4*4矩阵。可以看出最终结果的一个4*4单元需要4次4*4小矩阵相乘得到，即需要在C语言中调用4次汇编函数。

本实验程序包含两个文件，project.c、project.S，程序入口main()在C语言文件中，S汇编文件中只包含了4*4矩阵乘函数。

.text
.global asm_neon
asm_neon:
    ……
    RET

在project.c中调用该汇编函数可以通过以下方式，也就是将汇编函数声明为extern外部函数。

extern int asm_neon(float *matrix0, float *matrix1, float *out, int offset0, int offset1);

需要注意的是函数参数的传递方式，函数参数依次存放在从X0开始的寄存器中，即

float *matrix0 -> X0
float *matrix1 -> X1
float *out -> X2
int offset0 -> X3
int offset1 -> X4
五个参数中，三个矩阵都以地址形式传参，而两个跨步偏移量参数则是传值。

C语言实现的外层矩阵乘只需要用最基础的三层循环式矩阵乘实现方式即可，只不过每次计算都需要调用4*4矩阵乘汇编函数来实现。这里将外层矩阵乘包装为一个函数。

void process(float *maxtrix_A, float *maxtrix_B, float *maxtrix_out, int m, int n, int k)
{
	int M_len = m / 4;
	int N_len = n / 4;
	int K_len = k / 4;
	for (int x = 0; x < M_len; x++)
		for (int y = 0; y < N_len; y++)
			for (int z = 0; z < K_len; z++)
				asm_neon(&maxtrix_A[x * 4 * k + z * 4], &maxtrix_B[y * 4 + z * 4 * n], &maxtrix_out[x * 4 * n + y * 4], k * 4, n * 4); 
                //偏移量是按字节计算的，所以float的偏移需要乘4
}

2.2内层汇编4*4矩阵乘

因为每次汇编计算的4*4矩阵只是外层大矩阵中的一部分，所以需要考虑到跨步问题；并且因为外层大矩阵乘要通过多次4*4矩阵乘才能得出同一块4*4局部矩阵的最终结果，所以需要将每次计算累加，也就是需要先读入计算位置的4*4矩阵，将本次计算结果与之叠加，再存出。跨步和累加这两个因素都会在4*4矩阵乘汇编函数代码中体现出来。

4*4矩阵乘需要做什么？NEON SIMD向量化用在哪？如下所示，将数据读入128位NEON寄存器，如V4 {Y0，Y1，Y2，Y3} * V0[0] X0 -> V8 {Y0 * X0，Y1 * X0，Y2 * X0，Y3 * X0}这样的操作就是用向量操作实现的。

1.按行将matrix1、matrix2 两个4*4矩阵读入8个128位的寄存器中。
V0 {X0，X1，X2，X3} 
V1 {X4，X5，X6，X7} 
V2 {X8，X9，X10，X11} 
V3 {X12，X13，X14，X15} 
V4 {Y0，Y1，Y2，Y3} 
V5 {Y4，Y5，Y6，Y7} 
V6 {Y8，Y9，Y10，Y11} 
V7 {Y12，Y13，Y14，Y15} 

2.通过下面一系列计算得出答案
V4 {Y0，Y1，Y2，Y3} * V0[0] X0 -> V8 {Y0 * X0，Y1 * X0，Y2 * X0，Y3 * X0}
V4 {Y0，Y1，Y2，Y3} * V1[0] X4 -> V9 {Y0 * X4，Y1 * X4，Y2 * X4，Y3 * X4}
V4 {Y0，Y1，Y2，Y3} * V2[0] X8 -> V10 {Y0 * X8，Y1 * X8，Y2 * X8，Y3 * X8}
V4 {Y0，Y1，Y2，Y3} * V3[0] X12 -> V11 {Y0 * X12，Y1 * X12，Y2 * X12，Y3 * X12}

以与上面同样的相乘方式，下面的部分简写
V8 = V8 + V5 * V0[1]
V9 = V9 + V5 * V1[1]
V10 = V10 + V5 * V2[1]
V11 = V11 + V5 * V3[1]

V8 = V8 + V6 * V0[2]
V9 = V9 + V6 * V1[2]
V10 = V10 + V6 * V2[2]
V11 = V11 + V6 * V3[2]

V8 = V8 + V7 * V0[3]
V9 = V9 + V7 * V1[3]
V10 = V10 + V7 * V2[3]
V11 = V11 + V7 * V3[3]

此时V8-V11中存放的就是结果(Z0-Z15的值参考上面的原理图)
V8 {Z0，Z1，Z2，Z3} 
V9 {Z4，Z5，Z6，Z7} 
V10 {Z8，Z9，Z10，Z11} 
V11 {Z12，Z13，Z14，Z15} 

3.把存在寄存器V8-V11中的计算结果存回内存

V8 = V8 + V5 * V0[1]浮点乘加操作使用FMLA指令，与FMUL浮点乘(V8 = V5 * V0[1])不同。

下面是实际写出的汇编代码，考虑到跨步和需要与之前的计算结果累加两个因素。

    LD1 {V0.4S}, [X0], x3
    LD1 {V1.4S}, [X0], x3
    LD1 {V2.4S}, [X0], x3
    LD1 {V3.4S}, [X0], x3
    //读取矩阵2，每次偏移保存在x4寄存器中的n
    LD1 {V4.4S}, [X1], x4
    LD1 {V5.4S}, [X1], x4
    LD1 {V6.4S}, [X1], x4
    LD1 {V7.4S}, [X1], x4
    //保存结果矩阵的地址用于计算完后的存回
    mov x5, x2
	//读取结果矩阵，每次偏移保存在x4寄存器中的n
    LD1 {V8.4S}, [X2], x4
    LD1 {V9.4S}, [X2], x4
    LD1 {V10.4S}, [X2], x4
    LD1 {V11.4S}, [X2], x4
	//计算4*4的两个矩阵相乘，计算的时候通过[]选择哪个通道参与计算
    FMLA V8.4S, V4.4S, V0.S[0] 
    FMLA V9.4S, V4.4S, V1.S[0] 
    FMLA V10.4S, V4.4S, V2.S[0]
    FMLA V11.4S, V4.4S, V3.S[0] 
	
    FMLA V8.4S, V5.4S, V0.S[1]  
    FMLA V9.4S, V5.4S, V1.S[1]  
    FMLA V10.4S, V5.4S, V2.S[1] 
    FMLA V11.4S, V5.4S, V3.S[1] 

    FMLA V8.4S, V6.4S, V0.S[2]
    FMLA V9.4S, V6.4S, V1.S[2] 
    FMLA V10.4S, V6.4S, V2.S[2]
    FMLA V11.4S, V6.4S, V3.S[2]

    FMLA V8.4S, V7.4S, V0.S[3]
    FMLA V9.4S, V7.4S, V1.S[3]
    FMLA V10.4S, V7.4S, V2.S[3]
    FMLA V11.4S, V7.4S, V3.S[3]
	//将计算结果存回，偏移和读入的时候一样
    ST1 {V8.4S}, [X5], x4
    ST1 {V9.4S}, [X5], x4
    ST1 {V10.4S}, [X5], x4
    ST1 {V11.4S}, [X5], x4

3.编译运行

此处的编译运行均在ARMv8架构的鲲鹏处理器环境下执行。

通过gcc对汇编源码project.s、project.c进行编译，编译完成后运行。

gcc -o project project.c project.S
./project

编译运行结果：

程序通过正确的C语言实现矩阵乘结果来与混合汇编矩阵乘产生的结果相比较来判断结果是否正确，可以看到10个测例的正确数量为10/10，说明程序运行正确。

测试代码中所给的10个测例中的m、k、n分别为
{4, 4, 4},
{8, 12, 4},
{20, 40, 16},
{128, 36, 36},
{44, 4, 12},
{4, 48, 48},
{16, 8, 200},
{64, 64, 64},
{100, 8, 100},
{128, 256, 128}
矩阵中的数值为随机数

验证了本实验实现的矩阵乘程序适用于不同规模大小的矩阵相乘，当然前提是需要满足m、k、n都为4的倍数。

简单的性能评估与分析

另外可以看到在不加-O2、-O3编译器自动优化的情况下，我们实现的混合编程矩阵乘只用0.001s（在图中输出为ASM Total Time），而最为基础的纯C语言实现矩阵乘用了0.0398s（在图中输出为C Total Time），实现了接近40倍的性能优化。

通过这个矩阵乘案例的运行结果，可以帮助大家更好地理解到使用包括SIMD向量化在内的汇编混合编程具有优化程序性能的实际意义。除了调用SIMD向量化、减少编译器产生的额外操作外，汇编级的优化在某些处理器架构下还有更为细节的利用双指令通道实现的指令级并行优化。

从实际来讲，目前的编译器优化已经做得很好了，它也包含了自动向量化等操作，人工实现的汇编代码都很难做到在性能上比编译器生成的汇编代码好。并且向量化操作也可以通过Intrinsic这样的方式直接在C语言中调用，实际用到汇编优化的地方确实很少，除非你的程序需要追求极致的高性能。比如在OpenSSL中，为了追求极致的性能，很多加密算法都是通过汇编实现的。

4.完整的程序代码

该程序只含有两个文件：project.c、project.S

首先是project.c

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
//汇编实现的4*4矩阵乘
extern int asm_neon(float *matrix0, float *matrix1, float *out, int offset0, int offset1);
//学员需要自己实现混合编程矩阵乘函数，通过调用汇编实现的4*4矩阵乘来实现任意规模的矩阵乘
void process(float *maxtrix_A, float *maxtrix_B, float *maxtrix_out, int m, int n, int k)
{
	int M_len = m / 4;
	int N_len = n / 4;
	int K_len = k / 4;
	for (int x = 0; x < M_len; x++)
		for (int y = 0; y < N_len; y++)
			for (int z = 0; z < K_len; z++)
				asm_neon(&maxtrix_A[x * 4 * k + z * 4], &maxtrix_B[y * 4 + z * 4 * n], &maxtrix_out[x * 4 * n + y * 4], k * 4, n * 4); //偏移量是按字节计算的，所以float的偏移需要乘4
}

int main(int argc, char *argv[])
{
	//正确的测例个数
	int right_num = 0;
	//用于统计计算时间
	clock_t start_time, end_time;
	double total_time = 0, c_total_time = 0;
	//相乘的两个矩阵的边长
	int M_len, N_len, K_len;
	int tests[10][3] = {{4, 4, 4},
						{8, 12, 4},
						{20, 40, 16},
						{128, 36, 36},
						{44, 4, 12},
						{4, 48, 48},
						{16, 8, 200},
						{64, 64, 64},
						{100, 8, 100},
						{128, 256, 128}};
	//10个测例
	for (int i = 0; i < 10; i++)
	{
		M_len = tests[i][0];
		N_len = tests[i][1];
		K_len = tests[i][2];
		//相乘的两个矩阵，汇编输出结果，c语言实现的正确结果
		float *maxtrix_A = (float *)malloc(sizeof(float) * M_len * K_len);
		float *maxtrix_B = (float *)malloc(sizeof(float) * K_len * N_len);
		float *maxtrix_out = (float *)malloc(sizeof(float) * M_len * N_len);
		float *maxtrix_result = (float *)malloc(sizeof(float) * M_len * N_len);
		//用随机矩阵作为测例
		for (int m = 0; m < M_len * K_len; m++)
			maxtrix_A[m] = rand() % 10;
		for (int m = 0; m < K_len * N_len; m++)
			maxtrix_B[m] = rand() % 10;
		for (int m = 0; m < M_len * N_len; m++)
		{
			maxtrix_out[m] = 0;
			maxtrix_result[m] = 0;
		}
		//正确结果，因为矩阵为随机数值，此处用c计算得出正确结果
		start_time = clock();
		for (int m = 0; m < M_len; m++)
			for (int n = 0; n < N_len; n++)
				for (int k = 0; k < K_len; k++)
					maxtrix_result[m * N_len + n] += maxtrix_A[m * K_len + k] * maxtrix_B[k * N_len + n];
		end_time = clock();
		c_total_time += (double)(end_time - start_time);
		//汇编计算
		start_time = clock();
		process(maxtrix_A, maxtrix_B, maxtrix_out, M_len, N_len, K_len);
		end_time = clock();
		total_time += (double)(end_time - start_time);
		//比较汇编计算结果是否正确，并统计正确的测例数量
		int j = 0;
		for (; j < M_len * N_len; j++)
			if (maxtrix_result[j] != maxtrix_out[j])
				break;
		if (j == M_len * N_len)
			right_num++;

		free(maxtrix_A);
		free(maxtrix_B);
		free(maxtrix_out);
		free(maxtrix_result);
		maxtrix_A = NULL;
		maxtrix_B = NULL;
		maxtrix_out = NULL;
		maxtrix_result = NULL;
	}
	//输出评价结果
	printf("Right Num: %d/10 |", right_num);
	printf("C Total Time: %lfs | ", c_total_time / CLOCKS_PER_SEC);
	printf("ASM Total Time: %lfs\n", total_time / CLOCKS_PER_SEC);
	return 0;
}

然后是project.S

.text
.global asm_neon 
asm_neon:
    //读取矩阵1，每次偏移保存在x3寄存器中的k
	  LD1 {V0.4S}, [X0], x3
    LD1 {V1.4S}, [X0], x3
    LD1 {V2.4S}, [X0], x3
    LD1 {V3.4S}, [X0], x3
    //读取矩阵2，每次偏移保存在x4寄存器中的n
    LD1 {V4.4S}, [X1], x4
    LD1 {V5.4S}, [X1], x4
    LD1 {V6.4S}, [X1], x4
    LD1 {V7.4S}, [X1], x4
    //保存结果矩阵的地址用于计算完后的存回
    mov x5, x2
	//读取结果矩阵，每次偏移保存在x4寄存器中的n
    LD1 {V8.4S}, [X2], x4
    LD1 {V9.4S}, [X2], x4
    LD1 {V10.4S}, [X2], x4
    LD1 {V11.4S}, [X2], x4
	//计算4*4的两个矩阵相乘，计算的时候通过[]选择哪个通道参与计算
    FMLA V8.4S, V4.4S, V0.S[0] 
    FMLA V9.4S, V4.4S, V1.S[0] 
    FMLA V10.4S, V4.4S, V2.S[0]
    FMLA V11.4S, V4.4S, V3.S[0] 
	
    FMLA V8.4S, V5.4S, V0.S[1]  
    FMLA V9.4S, V5.4S, V1.S[1]  
    FMLA V10.4S, V5.4S, V2.S[1] 
    FMLA V11.4S, V5.4S, V3.S[1] 

    FMLA V8.4S, V6.4S, V0.S[2]
    FMLA V9.4S, V6.4S, V1.S[2] 
    FMLA V10.4S, V6.4S, V2.S[2]
    FMLA V11.4S, V6.4S, V3.S[2]

    FMLA V8.4S, V7.4S, V0.S[3]
    FMLA V9.4S, V7.4S, V1.S[3]
    FMLA V10.4S, V7.4S, V2.S[3]
    FMLA V11.4S, V7.4S, V3.S[3]
	//将计算结果存回，偏移和读入的时候一样
    ST1 {V8.4S}, [X5], x4
    ST1 {V9.4S}, [X5], x4
    ST1 {V10.4S}, [X5], x4
    ST1 {V11.4S}, [X5], x4
    RET

5.NEON寄存器与SIMD指令介绍

完整的A64 SIMD指令集功能与用法参考官方网站，链接为：

Documentation – Arm Developer

5.1.128位NEON寄存器

在ArmV8的Aarch64状态下，NEON提供了专门的32个128bit的向量寄存器v0-v31，用于支持浮点运算和SIMD运算。这些寄存器是独立于通用的Xn寄存器的。

如果需要使用v0-v31这些寄存器的全128位、低64位、低32位、低16位和低8位用于标量操作，那么也可以记名为Q0-Q31（全128位）、D0-D31（低64位）、S0-S31（低32位）、H0-H31（低16位）和B0-B31（低8位）。

在实际的使用过程中，可以通过在v0-v31之后添加后缀的方式来达到，既使用寄存器的全部128位又使用这128位中的一部分，这一目的。后缀的典型格式为：Vn.T（n为0-31之间的一个整数）。例如一条SIMD加法指令：ADD Vd.T, Vn.T, Vm.T。在该指令中，Vd、Vn和Vm都是Neon寄存器的名字，T是Neon寄存器的数据位宽表示方式。

标识	位宽	示例
b	8 bit	v1.b，v1.8b, v1.16b：v1寄存器的低8位，低64位，全部的128位；
h	16 bit	v1.4h, v1.8h：v1寄存器的低64位，全部的128位
s	32 bit	v1.2s, v1.4s：v1寄存器的低64位，全部的128位
d	64 bit	v1.1d, v1.2d：v1寄存器的低64位，全部的128位

这些128bit寄存器的低64位也可以当做64位向量寄存器来使用。如果将这些128位向量寄存器的低64位或更低位作为小位宽的向量寄存器的话，那么其高位部分将被清零（除了lane插入操作之外）。

5.2.SIMD指令

完整的A64 SIMD指令集功能与用法参考官方网站，链接为：

Documentation – Arm Developer

5.2.1.SIMD指令前缀修饰符

在Aarch64状态下，SIMD指令中不再有“V”这种助记符前缀，但是使用S/U/F/P等四种指令数据类型的前缀修饰符，可以被添加用来说明SIMD指令的操作数是有符号、无符号、浮点数、多项式中的某一种数据类型，例如，SMAX，UMAX，FMAX指令等。如果没有这些修饰符，则表明是与数据类型无关的SIMD指令。

指令前缀Q可以指定饱和指令。饱和指令限定了各种类型的数据范围。比如加法 SQADD 和 UQADD分别表示有符号饱和加以及无符号饱和加，如果结构数据超过了最大最小界限，饱和运算会使得结果不过超过最大或最小。比如，SQADD V0.16B, V0.16B, V1.16B实现有符号饱和加功能，其三个操作数都是16个字节。

5.2.2.SIMD指令后缀修饰符

Normal指令，对相同类型的数据进行操作，返回结果的数据类型与源类型相同。但是A64的SIMD还提供了两个源操作数和一个目标操作数之间长度不一致情况的指令后缀：

（1）长指令，使用L作为后缀，结果数据的位数是源数据位数的两倍，比如：SADDL V0.4S, V1.4H, V2.4H。

（2）Wide宽指令，对一个双字数据和一个单字数据进行操作，结果将都是双字数据，使用W作为后缀，比如：SADDW V0.4S, V1.4H, V2.4S。

（3）Narrow指令，操作两个四字向量，得到双字向量，结果数据是源数据的一半长，使用N作为后缀，比如：SUBHN V0.4H, V1.4S, V2.4S。

后缀修饰符除了定义操作数的宽度情况之外，还将指明其它一些情况：

（1）后缀P，表示矢量寄存器内部组对操作(pairwise)，即矢量寄存器中将相邻两个lane作为一组进行算法运算，并产生一个矢量结果。比如：addp v1.8B, v2.8B, v3.8B

（2）后缀V，表示Across，对寄存器中跨所有lane的操作，比如：addv B0, v1.8B // 将v1寄存器中的低64位中8个8位数据（跨8个B类型的lane）相加求和后，赋给v0的最低8位。

（3）指令后缀”2”，表示对源寄存器的高64位进行操作，比如：smull2 v0.4S, v1.4H, v2.4H //将v1和v2高64位中每一个16位相乘，并将结果放在v0的每个32位中。

5.2.3.SIMD操作数下标

NEON寄存器用于矢量运算时, 支持将16字节的NEON寄存器平均划分成若干个lane（通道），每个lane有相同类型/大小的位宽. 如v0.8h, 表示分成8个lane, 每个lane是两字节大小。Neon中含有向量功能单元，每个单元都完全实现流水化，每个时钟周期都可以开启一个新的操作。每个功能单元都可能只有单条流水线（单通道），也可以有多个并行的流水线（多通道）来并行执行相关操作。

在A64 SIMD指令中，操作数的下标(index)：表明操作数的哪一个lane被使用。比如：插入指令：INS Vd.Ts[index1], Vn.Ts[index2]。

其中，Vd为目标寄存器，Vn为源寄存器，Ts指明从最低位起需要使用的lane的宽度和个数，index1和index2表明在多个lane中仅使用哪一个lane，其上限为T-1。index为0则表示最低的lane，index的最大值为Ts中表示lane的最大数目。例如Ts为b类型，则128bit的寄存器最多一共可以有16个b类型的lane；Ts为d类型，则128bit的寄存器最多一共可以有2个d类型的lane。如果Ts为4b，则index值不能超过3。

Ts	lane宽度	下标index的总范围	index的实际范围
b	8bit	0-15	Ts为8b，则index<=7；Ts为2b，则index为0或1；
h	16bit	0-7	Ts为6h，则index<=5；Ts为2b，则index为0或1；
s	32bit	0-3	Ts为4s，则index<=3；Ts为2b，则index为0或1；
d	64bit	0或1	Ts为2d，则index<=1；

例如：ld1 {v18.S}[0], [x0], x1 // 将x0地址里面的数据取32位加载到v18的最低32位，然后x0=x0+x1。

withinmiaov

关注

2
点赞
踩
48

收藏

觉得还不错? 一键收藏
2
评论
ARMv8入门且实用的C语言汇编混合编程案例：NEON SIMD优化的通用矩阵乘

在ARMv8-A架构下，通过矩阵乘这样一个在图像处理、神经网络等各种领域都十分常见的算法，来帮助大家了解如何实现汇编与C语言混合编程。同时也会使用到NEONSIMD向量化技术，实现矩阵乘的性能优化。
复制链接

扫一扫