ARM處理器NEON編程及優(yōu)化技巧——矩陣乘法的實(shí)例

矩陣

本節將介紹如何用NEON有效的處理一個(gè)4x4的矩陣乘法運算，這種類(lèi)型的運算經(jīng)常用于3D圖形，我們認為這些矩陣在內存里是按照列為主排列的，這是按照OPENGL-ES的通用格式。

矩陣乘法算法

我們首先看一下矩陣乘法的計算方式，計算的展開(kāi)，用NEON指令來(lái)進(jìn)行子操作過(guò)程。

圖1. 以列為主的矩陣乘法運算

由于數據是按照列序存儲的，因而矩陣乘法就是把第一個(gè)矩陣的每一列乘以第二個(gè)矩陣的每一行，然后把乘積結果相加。乘累加結果 作為結果矩陣的一個(gè)元素。

圖2. 矩陣乘法中的向量乘以標量的運算

假設每列元素在NEON寄存器中表示為一個(gè)向量，那么上述的矩陣乘法就是一個(gè)向量乘以標量的運算，而后續的累加也同樣可以同向量乘以標量的累加指令實(shí)現。因為我們的操作是在第一個(gè)矩陣的列，然后計算列的結果，讀列元素和寫(xiě)列元素都是線(xiàn)性的加載和存儲操作，不需要interleave的加載和存儲操作。

代碼

浮點(diǎn)運算版本

首先看一個(gè)單精度浮點(diǎn)的矩陣乘法實(shí)現。首先加載矩陣元素到NEON寄存器，然后按照列序做乘法，用VLD1做線(xiàn)性的加載數據到NEON寄存器，用VST1把計算結果保存到內存。

vld1.32 {d16-d19}, [r1]!            @ 加載矩陣0的上8個(gè)元素 

vld1.32 {d20-d23}, [r1]!            @ 加載矩陣0的下8個(gè)元素 

vld1.32 {d0-d3}, [r2]!              @ 加載矩陣1的上8個(gè)元素 

vld1.32 {d4-d7}, [r2]!              @ 加載矩陣1的下8個(gè)元素 

NEON有32個(gè)64位寄存器，因而加載所有的輸入矩陣元素到16個(gè)64-bit寄存器，我們仍然有16個(gè)64位寄存器做后續的處理。

D和Q寄存器

多數的NEON指令有兩種方法來(lái)訪(fǎng)問(wèn)寄存器組：

• 作為32個(gè)雙字寄存器，64-bit位寬，命名為d0-d31
• 作為16個(gè)四字寄存器，128-bit位寬，命名為q0-q15

圖3. NEON寄存器組

這些寄存器中一個(gè)Q寄存器是一對D寄存器的別名，如Q0是d0和d1寄存器對的別名，寄存器中的值可以用兩種方式訪(fǎng)問(wèn)。這種實(shí)現方式很類(lèi)似C語(yǔ)言里的union聯(lián)合的數據結構。對于浮點(diǎn)的矩陣乘法，我們會(huì )經(jīng)常使用Q寄存器的表達方式，因為經(jīng)常會(huì )處理4個(gè)32-bit的單精度浮點(diǎn)，這對應于128-bit的Q寄存器。

代碼部分

通過(guò)以下4條NEON乘法指令能完成一列4個(gè)結果：

vmul.f32    q12, q8, d0[0]        @ 向量乘以標量（MUL），矩陣0的第一列乘以矩陣1的每列的第一個(gè)元素0 

vmla.f32    q12, q9, d0[1]        @ 累加的向量乘以標量（MAC），矩陣0的第二列乘以矩陣1的每列的第二個(gè)元素1 

vmla.f32    q12, q10, d1[0] @ 累加的向量乘以標量（MAC），矩陣0的第二列乘以矩陣1的每列的第二個(gè)元素2 

vmla.f32    q12, q11, d1[1] @ 累加的向量乘以標量（MAC），矩陣0的第二列乘以矩陣1的每列的第二個(gè)元素3 

第一條指令是圖2中的列元素x0, x1, x2, x3 (寄存器q8)乘以y0 (d0[0])，然后結果保存到q12寄存器。接下來(lái)的指令操作類(lèi)似，就是把第一個(gè)矩陣的其他列乘以第二個(gè)矩陣的第一列的響應元素，結果累加到寄存器Q12里。需要注意的是標量元素如d1[1]也可以用q0[3]表示，但是可能編譯器如GNU匯編器會(huì )不能接受這種方式。

如果我們只需要矩陣乘以向量的運算，如很多3D圖像處理中的那樣，那么此時(shí)的計算就結束了，可以把結果向量保存 到內存了，但是為了完成矩陣相乘，還需要完成后面的迭代操作，使用寄存器Q1到Q3的y4到yF的元素。如果定義如下的宏，那么就能簡(jiǎn)化代碼結構了：

.macro mul_col_f32 res_q, col0_d, col1_d 

vmul.f32    res_q, q8, col0_d[0]      @向量乘以標量（MUL），矩陣0的第一列乘以矩陣1的每列的第一個(gè)元素0 

vmla.f32    res_q, q9, col0_d[1]      @累加的向量乘以標量（MAC），矩陣0的第二列乘以矩陣1的每列的第二個(gè)元素1 

vmla.f32    res_q, q10, col1_d[0] @累加的向量乘以標量（MAC），矩陣0的第二列乘以矩陣1的每列的第二個(gè)元素2 

vmla.f32    res_q, q11, col1_d[1] @累加的向量乘以標量（MAC），矩陣0的第二列乘以矩陣1的每列的第二個(gè)元素3 

.endm 

那么整個(gè)4x4的矩陣乘法代碼可能如下：

vld1.32 {d16-d19}, [r1]!            @ 加載矩陣0的上8個(gè)元素 

vld1.32 {d20-d23}, [r1]!            @ 加載矩陣0的下8個(gè)元素 

vld1.32 {d0-d3}, [r2]!              @ 加載矩陣1的上8個(gè)元素 

vld1.32 {d4-d7}, [r2]!              @ 加載矩陣1的下8個(gè)元素 

mul_col_f32 q12, d0, d1         @ 矩陣 0 * 矩陣1的第一列 

mul_col_f32 q13, d2, d3         @ 矩陣 0 * 矩陣1的第二列 

mul_col_f32 q14, d4, d5         @ 矩陣 0 * 矩陣1的第三列 

mul_col_f32 q15, d6, d7         @ 矩陣 0 * 矩陣1的第四列 

vst1.32 {d24-d27}, [r0]!            @ 保存結果的上8個(gè)元素 

vst1.32 {d28-d31}, [r0]!            @ 保存結果的下8個(gè)元素 

定點(diǎn)算法

定點(diǎn)算法計算往往比浮點(diǎn)計算更快，因為往往定點(diǎn)運算可能需要更少的內存帶寬，整數值的乘法也會(huì )比浮點(diǎn)算法更為簡(jiǎn)單。但是定點(diǎn)算法，你需要很仔細的選擇表示格式來(lái)避免溢出或者飽和，這些會(huì )影響你的算法最終的精度。定點(diǎn)算法實(shí)現的矩陣乘法和浮點(diǎn)算法類(lèi)似，在本例中，用Q1.14定點(diǎn)格式,但是基本的實(shí)現格式基本類(lèi)似，只是實(shí)現中可能需要對結果做一些移位調整。下面是列乘的宏：

.macro mul_col_s16 res_d, col_d 

vmull.s16   q12, d16, col_d[0] @ 向量乘以標量（MUL），矩陣0的第一列乘以矩陣1的每列的第一個(gè)元素0 

vmlal.s16   q12, d17, col_d[1] @ 累加的向量乘以標量（MAC），矩陣0的第二列乘以矩陣1的每列的第二個(gè)元素1 

vmlal.s16   q12, d18, col_d[2] @ 累加的向量乘以標量（MAC），矩陣0的第二列乘以矩陣1的每列的第二個(gè)元素2 

vmlal.s16   q12, d19, col_d[3] @ 累加的向量乘以標量（MAC），矩陣0的第二列乘以矩陣1的每列的第二個(gè)元素3 

vqrshrn.s32 res_d, q12, #14             @ 把結果右移14位，并把累加結果變成Q1.14定點(diǎn)格式，并飽和運算 

.endm 

比較定點(diǎn)和浮點(diǎn)算法的宏，你會(huì )發(fā)現如下的主要區別：

• 矩陣元素的值為16位而不是32位，因而用D寄存器來(lái)保存4個(gè)輸入元素
• 矩陣乘法的結果是把16x16=32位的數據，使用VMULL和VMLAL來(lái)吧結果保存到Q寄存器。
• 最終結果也是16位，因而需要把32位累加器結果來(lái)得到16-bit的結果，使用VQRSHRN，飽和處理把32位的結果舍入到16位的narrow右移操作。

把數據從32-bits變成16-bits也能有效的處理內存訪(fǎng)問(wèn)，加載和存儲數據都只需要更少的帶寬。

vld1.16 {d16-d19}, [r1]     @ 加載16個(gè)元素到矩陣0 

vld1.16 {d0-d3}, [r2]       @ 加載16個(gè)元素到矩陣1 

mul_col_s16 d4, d0                      @ 矩陣0乘以矩陣1的列0 

mul_col_s16 d5, d1                      @ 矩陣0乘以矩陣1的列1 

mul_col_s16 d6, d2                      @ 矩陣0乘以矩陣1的列2 

mul_col_s16 d7, d3                      @ 矩陣0乘以矩陣1的列3 

vst1.16 {d4-d7}, [r0]       @ 保存16個(gè)結果元素 

指令重排

我們先展示一下指令重排如何能提高代碼性能。在宏中，臨近的乘法指令會(huì )寫(xiě)入到相同的目標寄存器，這會(huì )讓NEON的流水線(xiàn)等待前面的乘法結果完成才能開(kāi)始下一條指令的執行。如果不使用宏定義，而合理安排指令的次序，把那些相關(guān)依賴(lài)的指令變成不依賴(lài)，這些指令就能并發(fā)而不會(huì )造成流水線(xiàn)的stall。

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處理平臺上性能提升了一倍。從文檔arm.com/help/index.jsp?topic=/com.arm.doc.set.cortexa/index.html" rel="nofollow">Technical Reference Manual for your Cortex core可以看到 各個(gè)指令的需要時(shí)間以及延遲，有這些延遲周期，能夠更為合理的安排代碼次序，提升性能。