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

矩陣

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

矩陣乘法算法

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

圖1. 以列為主的矩陣乘法運(yùn)算

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

圖2. 矩陣乘法中的向量乘以標(biāo)量的運(yùn)算

假設(shè)每列元素在NEON寄存器中表示為一個(gè)向量，那么上述的矩陣乘法就是一個(gè)向量乘以標(biāo)量的運(yùn)算，而后續(xù)的累加也同樣可以同向量乘以標(biāo)量的累加指令實(shí)現(xiàn)。因?yàn)槲覀兊牟僮魇窃诘谝粋€(gè)矩陣的列，然后計(jì)算列的結(jié)果，讀列元素和寫列元素都是線性的加載和存儲(chǔ)操作，不需要interleave的加載和存儲(chǔ)操作。

代碼

浮點(diǎn)運(yùn)算版本

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

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寄存器，我們?nèi)匀挥?6個(gè)64位寄存器做后續(xù)的處理。

D和Q寄存器

多數(shù)的NEON指令有兩種方法來訪問寄存器組：

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

圖3. NEON寄存器組

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

代碼部分

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

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

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

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

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

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

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

.macro mul_col_f32 res_q, col0_d, col1_d 

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

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

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

vmla.f32    res_q, q11, col1_d[1] @累加的向量乘以標(biāo)量（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]!            @ 保存結(jié)果的上8個(gè)元素 

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

定點(diǎn)算法

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

.macro mul_col_s16 res_d, col_d 

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

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

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

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

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

.endm 

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

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

把數(shù)據(jù)從32-bits變成16-bits也能有效的處理內(nèi)存訪問，加載和存儲(chǔ)數(shù)據(jù)都只需要更少的帶寬。

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è)結(jié)果元素 

指令重排

我們先展示一下指令重排如何能提高代碼性能。在宏中，臨近的乘法指令會(huì)寫入到相同的目標(biāo)寄存器，這會(huì)讓NEON的流水線等待前面的乘法結(jié)果完成才能開始下一條指令的執(zhí)行。如果不使用宏定義，而合理安排指令的次序，把那些相關(guān)依賴的指令變成不依賴，這些指令就能并發(fā)而不會(huì)造成流水線的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í)間以及延遲，有這些延遲周期，能夠更為合理的安排代碼次序，提升性能。