高階矩陣并行求逆的FPGA設(shè)計(jì)與實(shí)現(xiàn)

伍松，毛宇河

（1.廣西科技大學(xué) 機(jī)械與交通學(xué)院，廣西 柳州 545006；2.廣西汽車零部件與整車技術(shù)重點(diǎn)實(shí)驗(yàn)室（廣西科技大學(xué)），廣西 柳州 545006）

1 引言

在現(xiàn)代數(shù)字信號(hào)處理中，許多先進(jìn)的技術(shù)方法都涉及高階矩陣運(yùn)算，例如在聲全息領(lǐng)域里的波疊加法近場(chǎng)聲全息技術(shù)[1]中或在對(duì)碳纖維復(fù)合材料的轉(zhuǎn)動(dòng)軸進(jìn)行模態(tài)分析時(shí)都會(huì)用到傳遞矩陣或傳遞矩陣的逆變換過(guò)程[2]，這些均與高階矩陣運(yùn)算息息相關(guān)。矩陣求逆過(guò)程亦是如此，其在數(shù)據(jù)分析、信號(hào)處理、系統(tǒng)理論、多元分析、現(xiàn)代控制等領(lǐng)域都有重要的使用價(jià)值，可見矩陣求逆的應(yīng)用范圍十分廣泛。

但由于求逆過(guò)程自身的繁瑣性和計(jì)算復(fù)雜性，因此大多數(shù)求逆運(yùn)算都過(guò)于依賴通用型計(jì)算機(jī)中的軟件開發(fā)平臺(tái)（如MAT?LAB、LabVEIW等），并且在需要進(jìn)行矩陣求逆等相關(guān)運(yùn)算時(shí)，還需先掌握開發(fā)平臺(tái)的配套語(yǔ)言和使用方式，增加了研究人員的設(shè)計(jì)難度。

另外矩陣求逆的傳統(tǒng)方法大多都是通過(guò)最小均方算法或最小二乘法求解近似值，這樣的求解過(guò)程往往精度低，收斂慢，難以達(dá)到設(shè)計(jì)要求[3]。同時(shí)在硬件實(shí)現(xiàn)方面，矩陣求逆也會(huì)受到資源及傳輸復(fù)雜程度的限制，相關(guān)的研究大都是針對(duì)特殊矩陣或者小規(guī)模矩陣進(jìn)行實(shí)現(xiàn)的，例如，文獻(xiàn)[4]是針對(duì)三角矩陣求逆進(jìn)行硬件實(shí)現(xiàn)的過(guò)程，而針對(duì)高階任意矩陣或大規(guī)模矩陣的相關(guān)文獻(xiàn)卻較少，因此對(duì)于眾多計(jì)算密集型問題的處理中，研究高階矩陣及大規(guī)模求逆的硬件實(shí)現(xiàn)過(guò)程十分具有挑戰(zhàn)性[5]。

即針對(duì)高階矩陣求逆問題，采用并行算法，基于現(xiàn)場(chǎng)可編程門陣列（Field Programmmable Gate Array，F(xiàn)PGA）的高速并行處理能力，采用流水線法以及LU分解求逆方法來(lái)完成基于FPGA[6]的8階滿秩方陣快速求逆處理器的設(shè)計(jì)與實(shí)現(xiàn)。這種快速實(shí)現(xiàn)高階矩陣求逆過(guò)程，不僅解決了傳統(tǒng)方法收斂性差、精度低的問題，而且還解決了矩陣求逆自身的計(jì)算復(fù)雜性以及高階矩陣難以求逆的問題，同時(shí)也避免了求逆運(yùn)算對(duì)于通用性計(jì)算機(jī)的過(guò)度依賴，能為很多使用到矩陣求逆過(guò)程的研究領(lǐng)域提供便利，提高運(yùn)算效率，具有廣泛的應(yīng)用價(jià)值，對(duì)工程應(yīng)用也有一定的幫助。

2 理論分析

LU分解是矩陣求逆過(guò)程中使用較為廣泛的一種分解求解辦法。其基本原理是先將n階方陣A分解成一下三角矩陣L同一上三角矩陣U的乘積形式后，利用三角矩陣易求逆的原則，分別求出下三角矩陣L和上三角矩陣U的逆形式，再通過(guò)矩陣相乘得到最終的逆矩陣A?1。即：

該方法可較大程度地減少求逆過(guò)程的運(yùn)算量和復(fù)雜程度，同時(shí)能解決最小均方、最小二乘等方法求近似解精度低收斂性差的問題。

2.1 矩陣的LU分解

設(shè)n階方陣的各階順序主子式均不為零，則A的LU分解則能僅包含一種結(jié)果。若以矩陣乘法原理為支撐，可推導(dǎo)得出LU分解的迭代結(jié)果算法[7]，即：

由上面的結(jié)果式可知，矩陣的LU分解是一種循環(huán)迭代過(guò)程，U矩陣和L矩陣中各結(jié)果元素的求得過(guò)程均是使用遞推迭代，過(guò)程類似。但由式（4）、式（5）分析可知，L矩陣實(shí)際上要比U矩陣少一步，因此在此處的設(shè)計(jì)上要注意小心防止出錯(cuò)。

2.2 上三角及下三角矩陣求逆

假設(shè)下三角矩陣L的逆矩陣為R，上三角矩陣U的逆矩陣為S，雖然主要計(jì)算8階矩陣的逆結(jié)果，但對(duì)于n階上下三角逆矩陣而言，其結(jié)果形式可規(guī)律呈現(xiàn)。例如，n階下三角逆矩陣R的一般形式為：

上三角逆矩陣S的結(jié)果式，其表達(dá)式同式（6）。

由式（6）結(jié)果式可以看出，計(jì)算上三角或下三角矩陣求逆過(guò)程基本可歸納為以下幾步：

（1）首先由分解出的L矩陣和U矩陣首先可以求得逆矩陣R和逆矩陣S的對(duì)角元素；

（2）由對(duì)角元素，在j=1的情況下，以式（6）中遞推式計(jì)算出R矩陣第1列或S矩陣第1行的結(jié)果元素；

（3）后再在j取2到n的順序條件下，用通過(guò)（1）～（2）及此步驟計(jì)算得到的結(jié)果元素，按i=j?1至i=1的先后順序倒序遞推，便可得出R矩陣或S矩陣的其他元素，完成該過(guò)程。

由上面的計(jì)算過(guò)程得知，LU分解法的迭代過(guò)程多，且迭代次數(shù)會(huì)隨著矩陣階數(shù)的提升而成倍增加，如果想要利用該方法實(shí)現(xiàn)高階矩陣的求逆過(guò)程，其迭代計(jì)算的復(fù)雜程度直接影響硬件設(shè)計(jì)過(guò)程的難易度。因此準(zhǔn)確控制求逆中的迭代過(guò)程，降低其在設(shè)計(jì)上的復(fù)雜性，是設(shè)計(jì)的關(guān)鍵。此外，本設(shè)計(jì)在該過(guò)程中將采用并行計(jì)算設(shè)計(jì)，因此保證中間過(guò)程計(jì)算的準(zhǔn)確性是設(shè)計(jì)的難點(diǎn)。

2.3 矩陣相乘

該過(guò)程是將上三角逆矩陣S與下三角逆矩陣R進(jìn)行矩陣相乘，依靠公式A?1=U?1L?1=SR，來(lái)求得原始方陣A最終逆矩陣，進(jìn)而完成求逆的計(jì)算過(guò)程。

那么由矩陣乘法原公式容易得到逆矩陣A的結(jié)果式為：

因此由式（7）得知，該過(guò)程所需的輸入數(shù)據(jù)數(shù)量較多，因此控制好輸入數(shù)的準(zhǔn)確性和輸入及儲(chǔ)存的時(shí)效性將是本過(guò)程的關(guān)鍵。

2.4 并行運(yùn)算

因本處理器中矩陣LU分解的輸入輸出過(guò)程及L矩陣與U矩陣的求逆過(guò)程是相對(duì)獨(dú)立且互不干涉，若將每個(gè)過(guò)程都串行連接，在時(shí)鐘控制下依次執(zhí)行，雖能得到計(jì)算結(jié)果，卻運(yùn)算時(shí)間長(zhǎng)，計(jì)算過(guò)程繁瑣。

因此為了提高運(yùn)算速度，考慮將獨(dú)立互不干涉的過(guò)程并行執(zhí)行，即對(duì)上述兩個(gè)過(guò)程采用并行結(jié)構(gòu)（Parallel Structure）[8]，形成并行形的數(shù)據(jù)輸入過(guò)程及矩陣求逆過(guò)程，該過(guò)程將原串行過(guò)程改為兩次并行過(guò)程，即可在原串行基礎(chǔ)上提速四倍；同時(shí)并行結(jié)構(gòu)還具有并發(fā)性和同時(shí)性，是各領(lǐng)域尤其是大數(shù)據(jù)領(lǐng)域達(dá)到高性能及高實(shí)時(shí)性計(jì)算需求的理想選擇。并行結(jié)構(gòu)的示意圖，如圖1所示。

圖1 LU分解輸入輸出并行結(jié)構(gòu)示意圖Fig.1 Parallel Structural Diagram of Input and Output Process of LU Decomposition

2.5 數(shù)的定標(biāo)

我們知道，在Quartus II開發(fā)平臺(tái)上，浮點(diǎn)數(shù)及小數(shù)是無(wú)法進(jìn)行綜合的，因此在使用vhdl語(yǔ)言描述硬件程序時(shí)，是不能在程序中使用和出現(xiàn)帶小數(shù)點(diǎn)的數(shù)據(jù)的，但如將浮點(diǎn)型近似為整型后，其數(shù)據(jù)的精度將大大降低。因此為了能夠在設(shè)計(jì)中使用浮點(diǎn)數(shù)和小數(shù)，采用定標(biāo)法給予解決。

所謂定標(biāo)法（Calibration Method），就是確定一個(gè)數(shù)小數(shù)點(diǎn)位置的方法[9]，而由定標(biāo)數(shù)確定出來(lái)的數(shù)，稱為定點(diǎn)數(shù)。若設(shè)定標(biāo)數(shù)為ξ，定點(diǎn)數(shù)為Φ，原浮點(diǎn)數(shù)設(shè)為X，則在二進(jìn)制中，浮點(diǎn)數(shù)與定點(diǎn)數(shù)的關(guān)系為：

由關(guān)系式看出，浮點(diǎn)數(shù)變?yōu)槎c(diǎn)數(shù)后，其精度將隨ξ的增大而增大，但隨著ξ的增加，在位數(shù)一定的情況下，其原數(shù)據(jù)的涵蓋范圍也將減小，即數(shù)據(jù)的精度和范圍是對(duì)立的。由式（8），也能推斷出，該方法的適用范圍并不局限于輸入數(shù)為浮點(diǎn)數(shù)或小數(shù)，整數(shù)同樣可以，若輸入數(shù)為整數(shù)，則式（8）也可變?yōu)椋?/p>

那么當(dāng)整型數(shù)需要進(jìn)行除法或開方運(yùn)算想要保住小數(shù)位數(shù)不丟失時(shí)，可采用定標(biāo)法，這便是思路，輸入數(shù)據(jù)為整數(shù)，但因過(guò)程數(shù)據(jù)會(huì)產(chǎn)生小數(shù)，為保留住過(guò)程數(shù)據(jù)的精度不流失，保留小數(shù)位數(shù)，因此采用定標(biāo)法。

2.6 任意高階滿秩方陣求逆法

雖然在設(shè)計(jì)上僅完成了8階滿秩矩陣的求逆過(guò)程，但最終目的是要依靠本設(shè)計(jì)去完成解決任意階滿秩矩陣的求逆工作。經(jīng)研究，將本設(shè)計(jì)與高階矩陣分解求逆方法結(jié)合起來(lái)，可以解決任意階矩陣的求逆。下面將對(duì)此方法進(jìn)行闡述。

首先需要簡(jiǎn)要介紹在解決高階滿秩矩陣求逆時(shí)所要用到的分解求逆方法。假設(shè)矩陣為N階滿秩方陣，且其主對(duì)角線上的K（K

式中：A、D—σ 階和（N?σ）階矩陣；B、C—σ×(N?σ)階和(N?σ)×σ階矩陣，矩陣H?1=(D?CA?1B)?1—δ階矩陣[10]，若令σ=8，則矩陣A的逆可由設(shè)計(jì)計(jì)算得到，因此只需求得H?1即可；再將H繼續(xù)按式(10）分解，又可得到下一個(gè)需求逆的8階矩陣和一個(gè)（N?16）階矩陣，8階矩陣求解方式同上，同樣又可按照式(10）進(jìn)行分解；以此類推，最后必然分解得到m個(gè)8階矩陣和1個(gè)階數(shù)小于8且大于0的矩陣R，然后再將矩陣R繼續(xù)按照σ=4，σ=2及σ=1的次序繼續(xù)進(jìn)行分解求逆過(guò)程后便可得到最終結(jié)果。

其中，4階矩陣、2階矩陣及1階矩陣的逆，可以簡(jiǎn)單通過(guò)修改8階矩陣程序參數(shù)來(lái)獲取，十分方便，因此綜上可知對(duì)于N階方陣的分解公式可為：

若令4α+2β+γ=t，則t可視為N除8后的余值，m視為商，便可快速得到分解次數(shù)，再按照上述的分解求逆方法便可得到任意階滿秩矩陣的求逆結(jié)果了。

3 FPGA 實(shí)現(xiàn)

本并行矩陣求逆的實(shí)現(xiàn)主要采用Altera公司提供的Quartus II 13.0軟件開發(fā)平臺(tái)來(lái)編寫求逆IP核的設(shè)計(jì)文件，再通過(guò)設(shè)計(jì)好的IP核搭建出最終的處理器，該過(guò)程便是設(shè)計(jì)思路。

與理論算法對(duì)應(yīng)，硬件實(shí)現(xiàn)電路包括4個(gè)主要模塊：總控制模塊（Total FSM，TFSM）、LU 分解分模塊（LU Decomposition，LU Decp）、上下三角矩陣求逆分模塊（LU Inversion，LU Inv）和矩陣乘法分模塊（Matrix Multiplication，Matrix Mult）。其中，總控制模塊包含總體狀態(tài)控制器（Total State Controller，TSC）和總體地址控制及狀態(tài)轉(zhuǎn)換計(jì)數(shù)器（Total Address&State Counter，TASC）；LU分解分模塊包含LU分解器（LU Resolver，LU RS）及其RAM存儲(chǔ)器，其中RAM存儲(chǔ)器包含L矩陣存儲(chǔ)器L_RAM和U矩陣存儲(chǔ)器U_RAM；上下三角矩陣求逆分模塊包含其自身求逆控制器（In?version FSM，IFSM）、求逆模塊（Inversion Model，Inv Model）及其RAM存儲(chǔ)器，其中求逆模塊和存儲(chǔ)器模塊分別包含L矩陣求逆部分和U 矩陣求逆部分；矩陣乘法分模塊包含（8×8）矩陣乘法器（Matrix Multiplication，Mx Mult）及其RAM存儲(chǔ)器；硬件系統(tǒng)的總體結(jié)構(gòu)，如圖2所示。

圖2 矩陣求逆系統(tǒng)總體結(jié)構(gòu)示意圖Fig.2 Overall Structure Diagram of Matrix Inversion System

其中，三個(gè)分模塊中的主要模塊：LU分解器、矩陣求逆模塊及矩陣乘法器，在設(shè)計(jì)上對(duì)其程序都進(jìn)行了參數(shù)化設(shè)計(jì)，因此容易根據(jù)實(shí)際情況，通過(guò)簡(jiǎn)單程序更改，實(shí)現(xiàn)不同階數(shù)的矩陣求逆過(guò)程，提升了程序自身的可重構(gòu)能力和實(shí)用性。

另外，對(duì)于輸入矩陣A 進(jìn)入本系統(tǒng)的輸入辦法，采用的是TXT文本輸入方式，即將寫著輸入數(shù)據(jù)的TXT文本文件由程序進(jìn)行識(shí)別，通過(guò)時(shí)鐘信號(hào)的具體控制，使文件中的數(shù)據(jù)準(zhǔn)確穩(wěn)定的傳輸?shù)匠绦虍?dāng)中；且由于其獨(dú)特的輸入方式，輸入數(shù)據(jù)可根據(jù)實(shí)際情況隨時(shí)進(jìn)行更改，也可將其他軟件得出來(lái)需要求逆的數(shù)據(jù)文件，通過(guò)格式轉(zhuǎn)換，轉(zhuǎn)換為TXT文件格式，之后便可通過(guò)本程序快速得出求逆結(jié)果。因此本設(shè)計(jì)的實(shí)用性、穩(wěn)定性和便利性都是十分突出的。根據(jù)由IP核的總體結(jié)構(gòu)示意圖以及要實(shí)現(xiàn)的功能，此設(shè)計(jì)的總體工作流程，如圖3所示。

（1）首先，在復(fù)位信號(hào)的控制下，總系統(tǒng)歸于初始狀態(tài)，且TSC處于S0的穩(wěn)定復(fù)位狀態(tài)。

（2）在時(shí)鐘信號(hào)的總體控制下，TASC中計(jì)數(shù)信號(hào)加一，并發(fā)送轉(zhuǎn)換信號(hào)q，TSC接收信號(hào)后使總體狀態(tài)轉(zhuǎn)換為S1，待狀態(tài)平穩(wěn)后，控制L_RAM從TXT文件中穩(wěn)定接收輸入矩陣A的各元素?cái)?shù)據(jù)并加以存儲(chǔ)，為接下來(lái)的計(jì)算做準(zhǔn)備。

（3）在TASC計(jì)數(shù)信號(hào)的控制下，TSC繼續(xù)按S2至S4的順序發(fā)生循環(huán)狀態(tài)轉(zhuǎn)換，LU分解過(guò)程便依序穩(wěn)定進(jìn)行，其中LU分解的輸入輸出過(guò)程采用并行結(jié)構(gòu)，以增加數(shù)據(jù)傳輸?shù)男市?；且由于需要?jì)算的是8階矩陣的逆，因此該過(guò)程需循環(huán)8次，當(dāng)然也可根據(jù)實(shí)際情況，通過(guò)對(duì)本程序進(jìn)行簡(jiǎn)要修改來(lái)更改矩陣階數(shù)，以實(shí)現(xiàn)不同階數(shù)矩陣的求逆，從而增加本設(shè)計(jì)的可重構(gòu)性能。

此步驟便是LU分解分模塊的工作過(guò)程，同時(shí)當(dāng)LU分解結(jié)束時(shí)，TSC的狀態(tài)信號(hào)轉(zhuǎn)換成為S4；

（4）之后，系統(tǒng)進(jìn)入上下三角求逆過(guò)程，由TSC發(fā)出的使能信號(hào)，控制啟動(dòng)LU 求逆分模塊的內(nèi)部控制器（IFSM），使其以內(nèi)部狀態(tài)轉(zhuǎn)換穩(wěn)定控制內(nèi)部RAM接收并存儲(chǔ)LU分解分模塊求得的各數(shù)據(jù)結(jié)果；之后，在IFSM狀態(tài)循環(huán)轉(zhuǎn)換的控制下，內(nèi)部的求逆運(yùn)算過(guò)程循環(huán)執(zhí)行，且其循環(huán)次數(shù)同步驟（3），同樣與矩陣階數(shù)一致，且求逆過(guò)程同樣采用并行結(jié)構(gòu)，使L矩陣求逆和U矩陣求逆雙過(guò)程同時(shí)進(jìn)行，提高了計(jì)算過(guò)程的效率性。當(dāng)執(zhí)行循環(huán)過(guò)程后，求逆控制器將給總控制器發(fā)送信號(hào)，一方面使總控制器繼續(xù)執(zhí)行下面步驟，一方面使矩陣乘法分模塊再次復(fù)位清空數(shù)據(jù)，使之后矩陣乘法過(guò)程能穩(wěn)定進(jìn)行；

（5）當(dāng)LU求逆分模塊工作完成后，TSC的內(nèi)部狀態(tài)繼續(xù)由S4轉(zhuǎn)換為S5，待狀態(tài)穩(wěn)定后，矩陣乘法分模塊在其控制下開始接收LU求逆的結(jié)果數(shù)據(jù)；之后，狀態(tài)變量再次自S6至S12循環(huán)轉(zhuǎn)換，控制矩陣乘法分模穩(wěn)定進(jìn)行計(jì)算過(guò)程，由式（7），因此乘法器計(jì)算采用邊相乘累加邊儲(chǔ)存的方法，按行方向依次依序得到結(jié)果數(shù)據(jù)，并存進(jìn)結(jié)果矩陣，循環(huán)往復(fù)；待計(jì)算全部完成后，通過(guò)輸出端口輸出最終結(jié)果，進(jìn)而完成整個(gè)求逆過(guò)程。其流程步驟，如圖3所示。

圖3 矩陣求逆系統(tǒng)總體工作流程圖Fig.3 Overall Work Flow Chart of Matrix Inversion System

其中，由于LU分解過(guò)程及LU矩陣求逆過(guò)程都包含了除法過(guò)程，因此為保證數(shù)據(jù)的精度不流失，采用了定標(biāo)法來(lái)進(jìn)行計(jì)算，過(guò)程數(shù)據(jù)利用定標(biāo)數(shù)采用先放大后縮小的方式來(lái)盡量減少中間過(guò)程因小數(shù)而產(chǎn)生的誤差損失，以提升結(jié)果的有效性。且為了保證快速性，在LU分解過(guò)程和LU求逆過(guò)程都采用了并行結(jié)構(gòu)，最大程度的保證求逆過(guò)程的時(shí)效性和本設(shè)計(jì)的實(shí)用性。

4 IP核的RTL級(jí)描述

設(shè)計(jì)的求逆處理器，是由8×8并行滿秩方陣求逆的IP核搭建而成的。其IP核的RTL級(jí)描述原理圖，可由四個(gè)分模塊共同構(gòu)成，其RTL級(jí)描述簡(jiǎn)圖，如圖4所示。

圖4 矩陣求逆系統(tǒng)RTL級(jí)描述簡(jiǎn)圖Fig.4 RTL Level Description Diagram of Matrix Inversion System

①為總體控制單元，是本系統(tǒng)的主控制單元，主要通過(guò)狀態(tài)轉(zhuǎn)換，為其他單元提供使能信號(hào)和地址信號(hào)。

①的上方TSC模塊表示總狀態(tài)轉(zhuǎn)換控制單元，用于執(zhí)行狀態(tài)轉(zhuǎn)換，提供使能信號(hào)，下方TASC模塊表示總地址控制及狀態(tài)轉(zhuǎn)換計(jì)數(shù)單元，用于通過(guò)計(jì)數(shù)信號(hào)控制狀態(tài)轉(zhuǎn)換控制單元執(zhí)行狀態(tài)轉(zhuǎn)換，并為其他模塊提供地址信號(hào)；

②為L(zhǎng)U分解總單元，用于執(zhí)行LU分解求逆中的LU分解過(guò)程工作，其共分為三個(gè)模塊，L_ram和U_ram兩模塊分別用于接收分解得到的L矩陣與U矩陣的結(jié)果，其中模塊L_ram同時(shí)還包含從文件接收輸入矩陣A及輸出輸入矩陣A兩項(xiàng)工作，模塊LU_Decp為L(zhǎng)U分解模塊，是LU分解單元的主要功能模塊，同時(shí)對(duì)其采用了并行結(jié)構(gòu)與定標(biāo)法，以提升運(yùn)算速度及運(yùn)算精確度；

③共同形成LU求逆總單元，其中ISC和IASC為該單元中的內(nèi)部控制部分，其原理與①相同，③中其他模塊為該單元中的主要執(zhí)行部分，其中Inv_ModelU與ram_inv_u兩模塊和Inv_ModelL與ram_inv_l是并行執(zhí)行關(guān)系，他們分別作為U 矩陣和L 矩陣的求逆執(zhí)行部分，以L 矩陣求逆部分為例，ram_inv_l 為RAM 存儲(chǔ)模塊，用于對(duì)運(yùn)算數(shù)據(jù)和結(jié)果數(shù)據(jù)進(jìn)行存儲(chǔ)、輸入和輸出，Inv_ModelL 為L(zhǎng) 求逆部分的主要功能模塊—求逆模塊，用于完成其求逆工作；

④中三模塊共同形成為乘法器總單元，其中mult_ram 為RAM存儲(chǔ)模塊，r_sct為行選擇模塊，matrix8×8為（8×8）矩陣乘法模塊，其中，RAM存儲(chǔ)器的功能除了包含數(shù)據(jù)的傳輸存儲(chǔ)，還增添了為乘法模塊提供地址的第二用途，行選擇模塊是為了保證矩陣乘法過(guò)程中在RAM循環(huán)提供右乘矩陣列方向元素時(shí)，其提供的左乘矩陣的行方向元素能準(zhǔn)確按時(shí)供應(yīng)，矩陣乘法模塊為該總單元的主要功能模塊。

5 仿真測(cè)試與分析

為了驗(yàn)證該8×8滿秩矩陣求逆處理器的功能性要求，將利用Modelsim仿真平臺(tái)對(duì)設(shè)計(jì)出來(lái)的IP核進(jìn)行仿真測(cè)試。主要測(cè)試過(guò)程、測(cè)試結(jié)果及對(duì)結(jié)果的分析結(jié)論簡(jiǎn)述如下：

設(shè)計(jì)選用Stratix高端系列芯片為目標(biāo)開發(fā)板，以方便在設(shè)計(jì)上降低對(duì)編譯速度、I/O口數(shù)、記憶位數(shù)（Memory Bits）及LEs邏輯元使用數(shù)等各參數(shù)的限制，增加在功能實(shí)現(xiàn)方面的有效性。芯片選用Altera Stratix III 系列的EP3SL340F1417C2 為目標(biāo)芯片，通過(guò)在Quartus II 13.0開發(fā)平臺(tái)上編譯、例化、綜合及在Modelsim上以10MHz為時(shí)鐘頻率進(jìn)行時(shí)序仿真分析后，得到的最終結(jié)果，如圖5所示。

圖5 矩陣求逆系統(tǒng)仿真計(jì)算結(jié)果Fig.5 Simulation Results of Matrix Inversion System

圖中的結(jié)果的顯示方式為，右斜列方向8個(gè)相鄰數(shù)據(jù)為結(jié)果矩陣的一行，例如圖中右斜第一列數(shù)：?131731，138855，?30116，88090，16053，?45763，?3422，?74240即為實(shí)際逆矩陣的第一行結(jié)果，其他行結(jié)果類似可得。圖中的所示結(jié)果是將原結(jié)果擴(kuò)大了5122倍且截?cái)嗔诵?shù)位數(shù)后的顯示結(jié)果，即相當(dāng)于結(jié)果作為了以18為定標(biāo)的定點(diǎn)數(shù)。這樣顯示的目的是為了解決在硬件系統(tǒng)無(wú)法處理浮點(diǎn)數(shù)的情況下，浮點(diǎn)數(shù)能以更高的精度在硬件中顯示和運(yùn)算的問題。其中，輸入矩陣：

若將輸入矩陣擴(kuò)大5122后放入MATLAB中，以LU分解的求逆方法求解結(jié)果，得到的結(jié)果，如圖6所示。

圖6 MATLAB中LU分解求逆計(jì)算結(jié)果Fig.6 Calculation Results of Inverse Method of LU Decomposition in MATLAB

若將本設(shè)計(jì)的結(jié)果矩陣與MATLAB求解矩陣的各元素按行方向順序展開成行向量后，在MATLAB中進(jìn)行誤差比較，可得到結(jié)果對(duì)比圖和誤差百分比圖，如圖7所示。計(jì)算出的百分誤差的平均值為?0.43%。可見用本文設(shè)計(jì)求得的結(jié)果值與MATLAB按此方法求得結(jié)果值的誤差十分小。若用MATLAB的INV函數(shù)去獲得結(jié)果，得到的結(jié)果值，如圖8所示。

圖7 MATLAB中LU分解求逆結(jié)果與設(shè)計(jì)求逆結(jié)果誤差對(duì)比Fig.7 The Error Comparison Between the Result of LU Decomposition in MATLAB and that of the Design

圖8 MATLAB中INV函數(shù)求逆計(jì)算結(jié)果Fig.8 Inverse Calculation Results of INV Function in MATLAB

將此結(jié)果按行方向順序展開成行向量后，與本設(shè)計(jì)結(jié)果進(jìn)行誤差比較，同樣可得到結(jié)果對(duì)比圖和誤差百分比圖，如圖9所示。

圖9 MATLAB中INV函數(shù)求逆結(jié)果與設(shè)計(jì)求逆結(jié)果誤差對(duì)比Fig.9 The Error Comparison Between the Result of INV Function in MATLAB and that of the Design

計(jì)算出的百分誤差的平均值為5.94%?？梢娪么嗽O(shè)計(jì)方法求得的結(jié)果值與實(shí)際結(jié)果值的誤差較小。此外，使用本設(shè)計(jì)求得結(jié)果耗費(fèi)的時(shí)間與MATLAB 按兩種方法求解花費(fèi)的時(shí)間對(duì)比，如表1所示。

表1 本設(shè)計(jì)仿真計(jì)算時(shí)間與MATLAB求解時(shí)間比較Tab.1 The Comparison Between the Simulation Calcula?tion Time of the Design and the Time of MATLAB Solution

由上表可知，本設(shè)計(jì)的計(jì)算效率已遠(yuǎn)超于MATLAB。通過(guò)圖5的數(shù)據(jù)結(jié)果可知，本設(shè)計(jì)實(shí)現(xiàn)了處理器能夠計(jì)算高階矩陣的逆，證明了其具有高速計(jì)算能力；同時(shí)通過(guò)圖5、圖6、圖8中結(jié)果數(shù)據(jù)的對(duì)比分析證明了使用定標(biāo)法計(jì)算出來(lái)的數(shù)據(jù)在擴(kuò)大了512倍后與MATLAB以兩種方式計(jì)算得出的數(shù)據(jù)十分接近；且由圖7、圖9中各結(jié)果元素的結(jié)果對(duì)比圖和誤差百分比圖，可以證明該設(shè)計(jì)的計(jì)算結(jié)果誤差較小精度較高；且再通過(guò)對(duì)表1中仿真計(jì)算時(shí)間與MATLAB 用兩種方法計(jì)算所花費(fèi)的時(shí)間進(jìn)行對(duì)比，可以得知，使用FPGA及并行結(jié)構(gòu)設(shè)計(jì)出的矩陣求逆，其速度已遠(yuǎn)超MATLAB求解的速率，從0.598s縮短到0.645μs，達(dá)到了高速求解的目的。通過(guò)以上分析可知，本設(shè)計(jì)能夠達(dá)到了預(yù)期的高速性和高效性。該設(shè)計(jì)雖然滿足求解的高精度要求，但從圖7與圖9也看出，所求得數(shù)據(jù)還存在偏差，圖9中有極個(gè)別數(shù)據(jù)的偏差較大。存在偏差的原因在于：

在LU分解和LU矩陣求逆過(guò)程中，雖然使用了定標(biāo)法盡量保證了小數(shù)位數(shù)的保留，但仍有個(gè)別較低位的小數(shù)部分丟失，遇到較小浮點(diǎn)數(shù)與較大整型數(shù)的乘積運(yùn)算時(shí)，其結(jié)果誤差就偏大。因此在后續(xù)的設(shè)計(jì)中，會(huì)對(duì)其做進(jìn)一步優(yōu)化，減小因定標(biāo)法截?cái)嗾`差產(chǎn)生的計(jì)算誤差，使本設(shè)計(jì)更加具有實(shí)用價(jià)值。

6 結(jié)論

采用LU分解、并行結(jié)構(gòu)、定標(biāo)等設(shè)計(jì)方法實(shí)現(xiàn)了一個(gè)高精度高效率的基于FPGA 支持的8 階滿秩矩陣求逆處理器的研究設(shè)計(jì)。該處理器利用了FPGA的并行能力，使矩陣求逆過(guò)程得以快速實(shí)現(xiàn)；同時(shí)該設(shè)計(jì)還具有計(jì)算穩(wěn)定性高，設(shè)計(jì)結(jié)構(gòu)清晰，可重構(gòu)性能強(qiáng)等特點(diǎn)。該設(shè)計(jì)不僅充分利用了FPGA芯片的高速并行性能，解決了高階矩陣難以求逆的問題，同時(shí)不再依賴繁瑣的通用計(jì)算機(jī)運(yùn)算，對(duì)實(shí)際工程具有一定的意義和價(jià)值。