NEON_16
矩阵乘法示例
矩阵乘法是在许多数据密集型应用程序中执行的操作。 它由以简单方式重复的算术运算组组成:
矩阵乘法过程如下:
A-在第一个矩阵中进行一行
B-执行该行的点积与第二个矩阵中的一列
C-将结果存储在新矩阵的相应行和列中
对于32位浮点矩阵,乘法可以写为:
void matrix_multiply_c(float32_t *A, float32_t *B, float32_t *C, uint32_t n, uint32_t m, uint32_t k) {
for (int i_idx=0; i_idx < n; i_idx++) {
for (int j_idx=0; j_idx < m; j_idx++) {
C[n*j_idx + i_idx] = 0;
for (int k_idx=0; k_idx < k; k_idx++) {
C[n*j_idx + i_idx] += A[n*k_idx + i_idx]*B[k*j_idx + k_idx];
}
}
}
}
我们假设内存中矩阵的主要列布局。 即,将n×m矩阵M表示为阵列M_array,其中Mij = M_array [n * j + i]。
该代码不是最佳代码,因为它没有充分利用Neon。 我们可以开始使用内在函数对其进行改进,但是让我们先解决一个简单的问题,即先查看固定大小的小型矩阵,然后再处理较大的矩阵。
以下代码使用内部函数将两个4x4矩阵相乘。 由于我们要处理的值很小且固定,所有这些值都可以一次放入处理器的Neon寄存器中,因此我们可以完全展开循环。
void matrix_multiply_4x4_neon(float32_t *A, float32_t *B, float32_t *C) {
// these are the columns A
float32x4_t A0;
float32x4_t A1;
float32x4_t A2;
float32x4_t A3;
// these are the columns B
float32x4_t B0;
float32x4_t B1;
float32x4_t B2;
float32x4_t B3;
// these are the columns C
float32x4_t C0;
float32x4_t C1;
float32x4_t C2;
float32x4_t C3;
A0 = vld1q_f32(A);
A1 = vld1q_f32(A+4);
A2 = vld1q_f32(A+8);
A3 = vld1q_f32(A+12);
// Zero accumulators for C values
C0 = vmovq_n_f32(0);
C1 = vmovq_n_f32(0);
C2 = vmovq_n_f32(0);
C3 = vmovq_n_f32(0);
// Multiply accumulate in 4x1 blocks, i.e. each column in C
B0 = vld1q_f32(B);
C0 = vfmaq_laneq_f32(C0, A0, B0, 0);
C0 = vfmaq_laneq_f32(C0, A1, B0, 1);
C0 = vfmaq_laneq_f32(C0, A2, B0, 2);
C0 = vfmaq_laneq_f32(C0, A3, B0, 3);
vst1q_f32(C, C0);
B1 = vld1q_f32(B+4);
C1 = vfmaq_laneq_f32(C1, A0, B1, 0);
C1 = vfmaq_laneq_f32(C1, A1, B1, 1);
C1 = vfmaq_laneq_f32(C1, A2, B1, 2);
C1 = vfmaq_laneq_f32(C1, A3, B1, 3);
vst1q_f32(C+4, C1);
B2 = vld1q_f32(B+8);
C2 = vfmaq_laneq_f32(C2, A0, B2, 0);
C2 = vfmaq_laneq_f32(C2, A1, B2, 1);
C2 = vfmaq_laneq_f32(C2, A2, B2, 2);
C2 = vfmaq_laneq_f32(C2, A3, B2, 3);
vst1q_f32(C+8, C2);
B3 = vld1q_f32(B+12);
C3 = vfmaq_laneq_f32(C3, A0, B3, 0);
C3 = vfmaq_laneq_f32(C3, A1, B3, 1);
C3 = vfmaq_laneq_f32(C3, A2, B3, 2);
C3 = vfmaq_laneq_f32(C3, A3, B3, 3);
vst1q_f32(C+12, C3);
}
由于某些原因,我们选择将固定大小的4x4矩阵相乘:
一些应用程序特别需要4x4矩阵,例如图形或相对论物理学。
Neon向量寄存器具有四个32位值,因此将程序与体系结构进行匹配将使其更易于优化。
我们可以采用这种4x4内核,并在更通用的内核中使用它。
让我们总结一下这里使用的内在函数:
现在我们可以将一个4x4矩阵相乘,可以将较大的矩阵视为4x4矩阵的块来相乘。 这种方法的一个缺点是,它仅适用于二维尺寸均为四的倍数的矩阵大小,但是通过将任何矩阵都填充为零,您可以使用此方法而无需对其进行更改。
下面列出了用于更通用的矩阵乘法的代码。 内核的结构变化很小,主要的变化是增加了循环和地址计算。 与在4x4内核中一样,即使我们可以使用一个变量并重新加载,我们也为B列使用了唯一的变量名。 这提示编译器将不同的寄存器分配给这些变量,这将使处理器能够在等待另一列加载的同时完成一列的算术指令。
void matrix_multiply_neon(float32_t *A, float32_t *B, float32_t *C, uint32_t n, uint32_t m, uint32_t k) {
/*
* Multiply matrices A and B, store the result in C.
* It is the user's responsibility to make sure the matrices are compatible.
*/
int A_idx;
int B_idx;
int C_idx;
// these are the columns of a 4x4 sub matrix of A
float32x4_t A0;
float32x4_t A1;
float32x4_t A2;
float32x4_t A3;
// these are the columns of a 4x4 sub matrix of B
float32x4_t B0;
float32x4_t B1;
float32x4_t B2;
float32x4_t B3;
// these are the columns of a 4x4 sub matrix of C
float32x4_t C0;
float32x4_t C1;
float32x4_t C2;
float32x4_t C3;
for (int i_idx=0; i_idx<n; i_idx+=4 {
for (int j_idx=0; j_idx<m; j_idx+=4){
// zero accumulators before matrix op
c0=vmovq_n_f32(0);
c1=vmovq_n_f32(0);
c2=vmovq_n_f32(0);
c3=vmovq_n_f32(0);
for (int k_idx=0; k_idx<k; k_idx+=4){
// compute base index to 4x4 block
a_idx = i_idx + n*k_idx;
b_idx = k*j_idx k_idx;
// load most current a values in row
A0=vld1q_f32(A+A_idx);
A1=vld1q_f32(A+A_idx+n);
A2=vld1q_f32(A+A_idx+2*n);
A3=vld1q_f32(A+A_idx+3*n);
// multiply accumulate 4x1 blocks, i.e. each column C
B0=vld1q_f32(B+B_idx);
C0=vfmaq_laneq_f32(C0,A0,B0,0);
C0=vfmaq_laneq_f32(C0,A1,B0,1);
C0=vfmaq_laneq_f32(C0,A2,B0,2);
C0=vfmaq_laneq_f32(C0,A3,B0,3);
B1=v1d1q_f32(B+B_idx+k);
C1=vfmaq_laneq_f32(C1,A0,B1,0);
C1=vfmaq_laneq_f32(C1,A1,B1,1);
C1=vfmaq_laneq_f32(C1,A2,B1,2);
C1=vfmaq_laneq_f32(C1,A3,B1,3);
B2=vld1q_f32(B+B_idx+2*k);
C2=vfmaq_laneq_f32(C2,A0,B2,0);
C2=vfmaq_laneq_f32(C2,A1,B2,1);
C2=vfmaq_laneq_f32(C2,A2,B2,2);
C2=vfmaq_laneq_f32(C2,A3,B3,3);
B3=vld1q_f32(B+B_idx+3*k);
C3=vfmaq_laneq_f32(C3,A0,B3,0);
C3=vfmaq_laneq_f32(C3,A1,B3,1);
C3=vfmaq_laneq_f32(C3,A2,B3,2);
C3=vfmaq_laneq_f32(C3,A3,B3,3);
}
//Compute base index for stores
C_idx = n*j_idx + i_idx;
vstlq_f32(C+C_idx, C0);
vstlq_f32(C+C_idx+n,Cl);
vstlq_f32(C+C_idx+2*n,C2);
vstlq_f32(C+C_idx+3*n,C3);
}
}
}
编译和反汇编此函数,并将其与我们的C函数进行比较,结果显示:
给定矩阵乘法的算术指令较少,因为我们利用具有完整寄存器打包功能的Advanced SIMD技术。 纯C代码通常不这样做。
FMLA代替FMUL指令。 如内在函数所指定。
更少的循环迭代。 如果正确使用内在函数,则可以轻松展开循环。
但是,由于内存分配和数据类型(例如,float32x4_t)的初始化而导致不必要的加载和存储,而纯C代码中未使用这些数据类型。
可以使用以下命令在Arm机器上编译和反汇编上面的完整源代码:
gcc -g -o3 matrix.c -o exe_matrix_o3
objdump -d exe_ matrix _o3 > disasm_matrix_o3
如果您无权使用基于Arm的硬件,则可以使用Arm DS-5 Community Edition和Armv8-A Foundation Platform。
上一篇: usr/bin/ld: i386 architecture of input file XXX.a( xxx.o) is incompatible with i386:x86-64 output
推荐阅读