基于力學性能的矩陣求逆電路的研究

0基于動陣列的求解矩陣思維是科學計算和工程的基本問題。科學計算的一般部分可以概括為矩陣思維，矩陣思維的困難在于矩陣逆調。傳統的矩陣求逆算法多用處理器串行計算來實現,嚴重制約著計算速度的提高。為加快矩陣求逆的運算速度,可以采用硬件實現的方法。在這一方面,已經有很多人進行了深入研究,特別是使用基于心動陣列的方法來構造電路結構,使用十分廣泛。這些論文研究的重點是通過對相應算法的研究改進,最大限度地提高資源的利用效率及時空復雜度,而忽略電路結構本身的硬件實現難度。例如文獻中使用QR分解的方法來實現矩陣求逆,提出的電路結構非常緊湊,復用率非常高,但是基本處理單元需要完成的運算非常復雜,不但涉及比較復雜的乘除運算,還需要進行開平方操作,這不僅增加了硬件開銷,更主要的是大大增加了硬件實現的難度。事實上,很多時候工程上涉及的矩陣求逆對硬件開銷及速度并沒有非常高的要求。所以,設計一種性能上略微低一點,但是易于硬件實現的電路結構就尤為重要了。基于這種考慮,構造了基于心動陣列的矩陣求逆電路,雖然在復用率及時空復雜度上較上述文獻中的結構略微不足,但其基本處理單元只是涉及十分簡單的加減乘除,非常易于硬件實現。實現了關鍵的浮點乘除法和浮點加減法運算,通過EDA工具對模型進行仿真及綜合,驗證了其正確性。結果表明,這種心動陣列能快速實現可逆矩陣的求逆運算,且有較高的精度。1基于狀態轉換的矩陣分解這樣就把一般的可逆矩陣求逆轉化成了三角矩陣求逆。同時考慮到對于三角矩陣U、V,假設U=VT,則V-1=((VT)T)-1=(UT)-1=(U-1)T,可見上三角和下三角矩陣求逆是一回事。綜上所述,只要實現LU分解及上三角矩陣求逆,就可以實現一般可逆矩陣的求逆。下面將主要針對這兩部分的VLSI實現展開討論。文中采用Gauss消去法對矩陣進行LU分解。Gauss消去法的實質是對方程組的系數矩陣進行三角化,其一般計算過程可用下面遞推式來表示:a(k)ij(k)ij=a(k-1)ij(k?1)ij-mika(k-1)kj(k?1)kj(1)(k=1,…,n;i=k+1,…,n;j=k+1,…,n+1)mik=a(k-1)ik(k?1)ik/a(k-1)kk(k?1)kk(2)式中mik稱為乘數,a(k)ij(k)ij為經過k次Gauss消去后的第i個方程的第j個系數。上述矩陣的LU分解的算法可以映射為圖1所示電路結構。在這個陣列中,有兩種基本的處理單元PE(ProcessingElement):方形的PE用于內積步的計算;圓形PE用于實現除法運算。矩陣A的元素從陣列頂部按列輸入,其中ai(j+1)比aij晚一個時鐘周期輸入。計算后,所得到的上三角矩陣U的元素存于各個PE之中。當A的第一行元素經過到第一行PE時,它們被存入第一行相應PE之中,之后的各行元素經過第一行PE時,只進行相應的變換而不再存入PE之中。A的第二行元素經過第一行PE變換后存入第二行相應的PE當中,之后各行元素經過第二行PE時也只進行相應變換后流向下一行PE。以此類推,直至A矩陣第n行元素(經過了n-1次變換)存入第n行PE為止。第i行左端的圓形PE依次產生消去第i列下三角元素的各個乘數mij(i=j+1,j+2,…,j=n-1,j=n),然后將它們沿水平方向向右傳送,參加內積步計算,最后由陣列右端輸出。至整個計算結束時,PEij中所存數據是U的元素uij,而陣列的第i行輸出是矩陣L的第n-i列元素。在L矩陣元素輸出之后,通過控制信號的設定,緊接著輸出U的元素,陣列的第i行輸出的是矩陣U的第i行元素。1.2內積步運算上三角矩陣求逆的電路結構參考了文獻,其結構如圖2所示。整個陣列呈三角形狀,需要n(n+1)/2個PE(圖中n=4)。在這個陣列中有兩類PE:一種是六角形內積器,實現內積步運算;另一種是右上側邊上的圓形PE,它實現除法運算,除了產生輸出外,還將vij的值垂直向下傳回陣列。如圖2中所示,三角矩陣A的元素沿陣列的左上向右下方向流動,單元矩陣E的元素則沿陣列的左下至右上方向流動;反饋輸入的V的元素沿陣列垂直方向向下流動,與A、E的元素在PE中匯合,執行內積步計算。分析可知,由右上方向輸出的數據,正好是上三角矩陣A的逆矩陣V的元素。2硬件設備電路的設計和實現在硬件實現基于心動陣列的矩陣求逆設計中主要設計如下幾個關鍵技術。2.1浮點數編碼格式現在科學與工程計算中需要處理的數據多是浮點數,所以矩陣求逆是在實數范圍內進行的,矩陣元素需用實數表示,因此要設計用于浮點運算的電路結構。文中,浮點數的表示采用IEEE754標準,其編碼格式為1位符號位、8位階碼和23位尾數來表示一個浮點數。尾數的表示范圍為,1是默認存在的,小數點的位置在第22位之前,比如1.5的23位尾數用16進制表示為400000,而1的23位尾數全是0。階碼的表示范圍為-127～128,值得一提的是階碼部分沒有符號位,而將真正的階碼值加127進行編碼,如6.0的實際階碼值應該是2,但是編碼時其階碼部分的值為129。符號位為1表示該實數為負,為0則表示是正數。另外,如果階碼和尾數都為0,則表示該實數值為0。由矩陣求逆的算法可知,設計中需設計浮點乘除電路和浮點加減法電路。2.2pe、六角形pe的浮點運算處理單元是心動陣列的基本單元,也是心動陣列設計的關鍵所在,下面介紹一下陣列中PE的實現。LU分解和上三角矩陣求逆的四種PE可見圖1和圖2。其中方形PE和六角形PE實現的功能基本類似,其核心是實現result=opc-opa*opb的浮點運算;而兩種圓形PE實現的功能也類似,其核心是實現result=opa/opb的浮點運算。上述opa,opb,opc及result都是采用IEEE754標準編碼的32位浮點數。從上面的分析不難看出,設計浮點乘除電路和浮點加減法電路是實現PE功能的關鍵。2.2.1確定尾數的編碼規則圖3所示是浮點乘法的電路結構。浮點乘法的符號位和階碼求解都比較簡單,為了求解其尾數,需要從DesignWare調用了一個24位乘法器,根據其結果來判定尾數值。如果乘積最高位為1,說明在編碼規則下其值大于2,需要右移1位,相當于截取46～24位作為尾數的前22位,相應的階碼加1;否則,其值滿足尾數區間,截取乘積45～23位作為尾數的前22位。另外這里采用馮·諾依曼舍入法,截取的尾數末位總為1,這樣可以有效避免誤差的累計。2.2.2尾數的計算圖4所示為浮點除法電路的結構。和乘法類似,除法的難點也在于尾數的計算,為了提高計算精度,這里把被除數左移24位,調用48位除法器,根據原被除數和除數尾數的大小,截取所得商的22～0位或23～1位即可得到相應的尾數值。2.2.3浮點疊加運算浮點加減法電路由于比較復雜,這里暫不列出去電路結構圖,整個電路實現了對階、尾數運算、規格化及舍入等操作,最終得到結果。其中電路設計中值得一提的有:1)通過設定一個標志位mode,電路可以實現浮點加法和減法的運算,當mode為1時,電路結構實現的是減法運算,當mode為0時,其實現的是加法運算。實現的方法是把mode和opb(第二個操作數)的符號位相異或得到實際運算中opb的符號位,然后對opa(第一個操作數)和opb進行加法運算。2)對于兩個操作數中有一個為零或者都為零的情況,事先進行判斷,通過一個MUX電路把運算結果直接傳遞給輸出,提高效率。3)結構中有一個“確定左移位數的判斷邏輯”的模塊,它實現尾數規格化的功能,通過判斷操作數尾數相加所得值的第一個“1”(從高位到地位)的位置,來確定需要左規格化的次數。另外采取了四位一組的方法,有效地減小了面積。3l、u方向的轉置矩陣文中用一個4階可逆矩陣來驗證求解其逆矩陣。由于整個求逆過程需要的周期比較長,故把整個波形圖分割成四個子圖來說明,如圖5所示。輸入的原始矩陣A為:A=[4321310020101001]=[4080000040400000400000003f800000404000003f800000004000000003f80000003f800000003f800000]A=??????4321310020101001??????=??????4080000040400000400000003f800000404000003f800000004000000003f80000003f800000003f800000??????LU分解的仿真波形如圖5(a)所示。端口data-in1～4輸入A矩陣第1～4列元素,每列之間延遲一個時間單元。端口m-out14依次輸出L矩陣的第一列元素及U矩陣的第一行元素,m-out24端口輸出L矩陣的第二列元素及U矩陣的第二行元素,m-out34端口輸出L矩陣的第三列元素及U矩陣的第三行元素,m-out44端口輸出U矩陣的第四行元素。值得注意的是,L、U都是三角函數,且L矩陣的對角線元素為1,對于這些已知的0和1這里沒有輸出。由此可得L、U矩陣為:L=[10003f4000001003f0000003f999998103e8000003f1999983e638e351]?[10003/41001/26/5101/43/52/91]U=[4080000040400000400000003f8000000bfa00002bfc00001bf400001003fe666643eccccc80003f8e38e4]?[43210-5/4-3/2-3/4009/52/500010/9]LU分解后分別對U矩陣及L的轉置矩陣求逆,其仿真波形分別如圖5(b)和(c)所示。圖5(b)中,端口data-in1～4輸入U矩陣的元素,端口e-out1～4依次輸出第一行開始的對角線方向的逆矩陣元素,其中e-out1輸出對角線元素。由此可得U的逆矩陣為:圖5(c)的仿真波形和圖5(b)類似,唯一有區別的是圖5(c)得到的輸出是L的轉置矩陣的逆矩陣元素,這點在之后做矩陣乘法時需要注意。可以從圖中看出,得到矩陣如下:(LΤ)-1=[3f800000000bf4000013f800000003eccccc8bf9999993f80000003de38e3cbeaaaaaabe638e353f800000]?[1000-3/41002/5-6/5101/9-1/3-2/91]把L、U矩陣相乘即可得到A的逆矩陣C,如圖5(d)所示。端口dout1～4分別輸出C矩陣的1～4行,每行之間延遲一個時間單元輸出,即A的逆矩陣為:C=[bdccccbe3e9999953e4cccd23dcccccd3e9999953dcccce1bf19999bbe99999b3e4ccccabf19999a3f19999abe4cccc93dcccccfbe99999bbe4cccc93f666665]?[-0.09999990.29999990.20000000.10000000.29999990.1000000-0.6000000-0.30000000.1999999-0.60000000.6000000-0.19999990.1000000-0.3000000-0.19999990.8999999]整個矩陣求逆過程從輸入第一個A矩陣元素到輸出所有A的逆矩陣元素共花費了38個時鐘周期。從上面的仿真波形中可以看出,除去由浮點數精度所限帶來的截斷誤差(僅為10-7數量級),4階矩陣求逆的結果和理論所得的值完全吻合。用smic的0.18μm工藝對Verilog代碼進行綜合,可得整個電路的面積約為2.29mm×2.29mm,最長路徑約為25ns,時鐘頻率約為40MHz。由于整個求逆電路主要由LU分解、上三角矩陣求逆、矩陣相乘三部分構成,三者可以實現流水線操作,因此,考慮到工程上通常是大量數據需要處理,其實際的工作效率可以提高3倍左右。文中的最長路徑是由DesignWare中調用的除法器決定的,不考慮流水線,時鐘頻率可以達到40MHz,對于一般的工程問題完全可以滿足要求。文獻中提出的結構雖然只需要5n=20個時鐘周期就可以完成整個求逆過程,但是由于它基本處理單元中運算操作相當復雜,必然會對其時鐘頻率產生一定的影響而