FPGA 架構(gòu)上面向稀疏矩陣求解的靜態(tài)調(diào)度算法

王晞陽，陳繼林，李猛，劉首文

（1.國家超級計算無錫中心，江蘇 無錫 214072；2.中國電力科學(xué)研究院有限公司，北京 100192；3.國網(wǎng)湖北省電力有限公司，武漢 430070）

0 概述

在電力系統(tǒng)仿真中，電磁暫態(tài)仿真是主要的系統(tǒng)應(yīng)用之一，也是電力系統(tǒng)安全分析和運行的關(guān)鍵組件［1］，該應(yīng)用的核心算法是對大規(guī)模線性方程組Ax=b進(jìn)行求解。通過對實際數(shù)據(jù)的分析可知，線性方程組中的系數(shù)矩陣A通常為稀疏矩陣，具體到電力系統(tǒng)，其稠密度通常小于1%［2-3］，如果不能有效利用矩陣的稀疏性，在使用計算機(jī)處理大型稀疏矩陣時，大量的存儲和計算資源將會浪費在無效的零元上，導(dǎo)致存儲空間不足和處理效率低下［4］。因此，在實現(xiàn)大型稀疏矩陣存儲和計算時，需要利用專用的算法和數(shù)據(jù)結(jié)構(gòu)，基于特殊設(shè)計的硬件架構(gòu)，最大化地利用計算系統(tǒng)的算力和效能［5-6］。綜上，電力系統(tǒng)的電磁暫態(tài)加速求解問題最終聚焦于稀疏線性方程組的加速求解問題。

稀疏線性方程組的求解方法包括直接法和迭代法兩大類。由于迭代法存在迭代次數(shù)不可控、結(jié)果精度較低等問題，因此一般使用直接法進(jìn)行求解。直接法是指在不考慮計算舍入誤差的情況下，通過矩陣分解和三角方程［7-8］進(jìn)行求解。下三角稀疏矩陣求解是求解稀疏線性方程組的核心算法的組成部分［9］，而層次集合法［10-11］是最重要的并行性分析方法之一。近年來，越來越多的算法開始使用眾核加速硬件來實現(xiàn)并行加速。NAUMOV［12］提出基于層次集合的方法，其將多個小的層次集合進(jìn)行合并，以減少同步運算。PARK等［13］通過對同步進(jìn)行剪枝來提升效率。LIU等［14-15］提出利用圖形處理器（Graphics Processing Unit，GPU）的原子操作來實現(xiàn)無同步的算法。SU等［16］介紹大規(guī)模的線程級無同步算法。LU等［17］介紹基于循環(huán)塊優(yōu)化的算法。上述算法主要在GPU 上進(jìn)行優(yōu)化，WANG等［18］則提出了在國產(chǎn)神威眾核上實現(xiàn)的加速算法。

相比眾核架構(gòu)（GPU 或神威眾核），現(xiàn)場可編程門陣列（Field Programmable Gate Array，F(xiàn)PGA）具有片上緩沖大、傳輸帶寬可定制、調(diào)度靈活等優(yōu)勢，可以最大化地利用計算系統(tǒng)的算力［19］。吳志勇等［20］設(shè)計了一種面向FPGA 求解稀疏矩陣的硬件架構(gòu)，本文基于這一硬件結(jié)構(gòu)，提出軟件映射和調(diào)度求解算法。該算法給出數(shù)據(jù)分布和指令排布的過程，通過將下三角稀疏矩陣的求解過程靜態(tài)映射到多個FPGA 片上的處理單元（Process Elements，PEs），以軟硬件協(xié)同的方法實現(xiàn)下三角稀疏矩陣在定制化FPGA 架構(gòu)上的高速求解。

1 FPGA 稀疏矩陣求解器硬件設(shè)計

FPGA 稀疏矩陣求解器的硬件設(shè)計在文獻(xiàn)［20］中已經(jīng)詳細(xì)說明，本文只簡述該求解器的基本原理和結(jié)構(gòu)，以便進(jìn)一步介紹本文所提軟件調(diào)度算法。直接法需要遵循下三角稀疏矩陣內(nèi)部的對角元依賴關(guān)系，其求解路徑構(gòu)成一個有向無環(huán)圖（Directed Acyclic Graph，DAG）［21］，如圖1 所示。由于DAG 的條件路徑依賴關(guān)系無法直接映射為并行化算法，因此無法直接在算法層進(jìn)行并行優(yōu)化。下三角稀疏矩陣具有稀疏性，構(gòu)成的DAG 存在隱式的數(shù)據(jù)并行，意味著多個對角元和非對角元之間存在并行處理的可能。整個系統(tǒng)的基本思想是設(shè)計多個硬件處理單元，每個處理單元都能夠有效地處理對角元和非對角元。通過軟件對該稀疏矩陣進(jìn)行預(yù)處理，找到稀疏下三角的并行性，并且將劃分完成的下三角稀疏矩陣映射到這些硬件處理單元上，從而實現(xiàn)軟硬件協(xié)同的并行求解。

圖1 有向無環(huán)圖Fig.1 Directed acyclic graph

由于軟件預(yù)處理及劃分完成后的稀疏下三角各個部分（對角元及非對角元之間）存在相互依賴，因此需要為多個處理單元間提供硬件連接通路。在實際的FPGA 實現(xiàn)過程中，使用二維單向環(huán)網(wǎng)來實現(xiàn)多個處理單元間的通信連接，從而完成各單元間的數(shù)據(jù)傳輸［22］。FPGA 稀疏矩陣求解器的硬件結(jié)構(gòu)如圖2 所示。

圖2 求解器的硬件結(jié)構(gòu)Fig.2 Hardware structure of solver

從圖2 可以看出，系統(tǒng)中存在多個處理單元，每個處理單元都有一個指令控制部件M，用來存儲PE內(nèi)部控制部件的指令碼。另外，每個處理單元內(nèi)還有3 類操作部件，分別是Mul 部件、Add 部件和T/R部件，其中：Mul 部件用來處理稀疏矩陣對角元和非對角元的雙精度浮點復(fù)數(shù)乘法（包括除法）操作；Add部件用來處理非對角元對同行對角元的更新；T/R部件用來實現(xiàn)相鄰PE 間的數(shù)據(jù)傳輸。

每個PE 內(nèi)部包括5 個專用緩沖，分別是A、A'、B、B'和B"，其中：緩沖A用來存儲稀疏矩陣對角元系數(shù)；緩沖A'用來存儲稀疏矩陣非對角元系數(shù)；緩沖B用來存儲稀疏矩陣右端項和求得的未知數(shù)；緩沖B'用來存儲非對角元和同列未知數(shù)的乘積；緩沖B"用來存儲其他PE 傳輸過來的已求解的未知數(shù)x。由于每個PE對B"緩沖的寫操作一定都指向不同的地址，因此B"緩沖可以在2 個相鄰單元間共享，以減少對FPGA 資源的占用，有利于FPGA 綜合實現(xiàn)。處理單元間則通過二維單向環(huán)網(wǎng)相連來進(jìn)行通信?；谠撚布軜?gòu)，原本串行的稀疏矩陣求解過程被分解為以下操作步驟：

1）使用Mul 部件，從緩沖A中讀取對角元系數(shù)，從緩沖B中讀取右端項，求解獲得未知數(shù)x，再存放到緩沖B中。這個x不被其他PE 上的非對角元需求，無需傳輸。

2）使用Mul 部件，從緩沖A中讀取對角元系數(shù)，從緩沖B中讀取右端項，求解獲得未知數(shù)x，存放到緩沖B中。如果這個x被其他PE 的非對角元需求，則同時將其發(fā)送到R部件。

3）使用Mul 部件，從緩沖A'中讀取非對角元系數(shù)，從緩沖B中讀取求得的x，計算獲得非對角元乘積，將其存放在緩沖B'中。

4）使用Mul 部件，從緩沖A'中讀取非對角元系數(shù)，從緩沖B"中讀取求得的x，計算獲得非對角元乘積，將其存放在緩沖B'中。本次讀取到的x對應(yīng)第2步操作中傳輸?shù)絉部件中的x。

5）使用Add 部件，從緩沖B'中讀取非對角元乘積，從緩沖B中讀取對角元，計算更新對角元，并將結(jié)果存放到緩沖B中。

6）使用T部件，將左方或者下方傳輸來的x通過網(wǎng)絡(luò)傳輸?shù)较乱粋€PE，并保存到R部件和緩沖B"中。

根據(jù)上述稀疏矩陣求解過程，雙精度浮點復(fù)數(shù)乘法單元和加法單元的數(shù)據(jù)調(diào)度過程分別如圖3和圖4所示。

圖3 乘法單元的數(shù)據(jù)調(diào)度Fig.3 Data scheduling of the multiplication unit

圖4 加法單元的數(shù)據(jù)調(diào)度Fig.4 Data scheduling of the addition unit

2 軟件靜態(tài)調(diào)度算法

基于上文FPGA 稀疏矩陣求解器的硬件架構(gòu)，本文實現(xiàn)了稀疏矩陣求解算法和求解過程的映射，充分挖掘了稀疏矩陣求解過程中的并行性，這一映射過程通過軟件靜態(tài)調(diào)度算法來實現(xiàn)。中央處理器（Central Processing Unit，CPU）一般通過順序求解所有解向量來實現(xiàn)，而GPU 則會對解向量進(jìn)行分塊，在分塊后的解向量之間通過引用計數(shù)構(gòu)建依賴關(guān)系，然后逐塊觸發(fā)從而完成求解。無論是CPU 還是GPU，都可以通過動態(tài)尋址來定位解向量，但是，F(xiàn)PGA 不具備動態(tài)尋址的能力，因此，必須構(gòu)建靜態(tài)指令流，指示每個操作部件在特定步時的訪存地址和處理操作。

2.1 稀疏矩陣預(yù)處理

在正式求解之前，需要對矩陣進(jìn)行預(yù)處理。稀疏矩陣存儲通常使用壓縮稀疏列矩陣（Compressed Sparse Column matrix，CSC）或壓縮稀疏行矩陣（Compressed Sparse Row matrix，CSR）［23］。首先，將矩陣進(jìn)行重排序，以減少LU 分解增加的額外非零元，本文使用metis［24］工具進(jìn)行重排序，其能提供一組可以獨立運行的命令行程序，同時也提供應(yīng)用程序接口（Application Programming Interface，API），方便集成到C/C++或Fortran 程序中，該程序可以得到上述算法目標(biāo)的一個近似解；其次，對完成重排序的稀疏矩陣進(jìn)行LU 分解，獲得右下局部稠密的下三角稀疏矩陣，如果沒有特殊說明，下文以L代指這里的下三角稀疏矩陣。

2.2 矩陣劃分

針對上述硬件架構(gòu)，軟件調(diào)度算法需要為其每個PE 進(jìn)行任務(wù)分派和調(diào)度，將求解下三角稀疏矩陣L需要用到的數(shù)據(jù)分布到多個計算單元上。

數(shù)據(jù)分布的設(shè)計目標(biāo)是使得盡可能多的計算單元盡可能飽和地運行。第一個目標(biāo)是使得分布到多個PE 上的數(shù)據(jù)盡可能均勻，以保證每個PE 需要求解的數(shù)據(jù)量接近。由于稀疏矩陣的非零元計算過程實際上是一個有向無環(huán)圖，因此將一個稀疏矩陣映射到多個PE 上的問題可以認(rèn)為是一個圖劃分（Graph Partition）問題。第二個目標(biāo)是使得分割后的各個PE 間通信盡可能少，從而減小通信開銷，使計算單元盡可能飽和地運行。

圖的均勻劃分是一個NP 困難（NP-hard）問題，同樣使用metis 提供的接口，對L矩陣的對角元進(jìn)行劃分。劃分完成的對角元分布在不同的分塊中，這些分塊一一映射到多個不同的PE 上。原有的各個對角元之間的依賴關(guān)系被分成了PE 內(nèi)部的對角元依賴關(guān)系和跨PE 的對角元依賴關(guān)系，分別保存在pe_inside_diag_dep 和pe_outside_diag_dep 數(shù)組中。

非對角元被分布到其所在行的對角元所在的PE 上，這樣操作的好處是使得非對角元就緒后，可以直接更新其所在行的對角元，傳輸部件只需要傳輸對角元，能夠大幅簡化與傳輸部件相關(guān)的算法設(shè)計。

2.3 傳輸算法

在實際的硬件實現(xiàn)中，當(dāng)PE 數(shù)較多時，PE 間直接相連會帶來巨大的網(wǎng)絡(luò)資源消耗，因此，硬件設(shè)計通常使用二維mesh 網(wǎng)絡(luò)來替代PE 間的直接連接。但是，二維mesh 網(wǎng)絡(luò)帶來的問題是每個PE 只和周圍2 個維度上的4 個PE 直接相連，當(dāng)要和遠(yuǎn)端的PE通信時，需要經(jīng)過中間的PE 進(jìn)行連接，此時需要多個時鐘周期來完成傳輸。

當(dāng)前某個PE 計算獲得的x有可能會被最遠(yuǎn)端的PE 訪問，因此，在軟件調(diào)度算法上，一開始選擇廣播算法進(jìn)行數(shù)據(jù)共享。廣播算法可以保證當(dāng)前PE 計算求得的未知數(shù)x最終一定能夠被所有的PE 訪問到，并且對于二維單向環(huán)網(wǎng)而言，廣播算法非常容易設(shè)計和實現(xiàn)，如圖5 所示。二維單向環(huán)網(wǎng)中所有的PE 對于網(wǎng)絡(luò)是完全對稱的，從非0 號PE 發(fā)起的廣播傳輸過程也是類似的。在使用廣播算法時，同一時刻多個PE 發(fā)起的消息只要滿足PE 的行號與列號之和不同于PE 數(shù)量（mod PE），就可以在系統(tǒng)中并行傳輸。

圖5 從0 號PE 開始的16PE 廣播算法Fig.5 16PE broadcast algorithm starting from PE 0

實際上，并非所有的PE 都需要獲得當(dāng)前正在傳輸?shù)膞，廣播算法會占用不需要x的PE 的傳輸部件T/R。為了解決這一問題，可以將廣播算法改為點對點傳輸算法，點對點算法的路由與廣播算法一致，但只占用發(fā)送PE 和接收PE 之間的通路，空出的通路可以同時傳輸其他x，這使得系統(tǒng)中能夠同時傳輸更多的x，進(jìn)一步提升了傳輸效率。

2.4 指令排布和生成

在數(shù)據(jù)分布完成后，就可以對指令進(jìn)行排布，需要完整地規(guī)劃每個PE 上的部件操作，從而充分利用每個PE 的操作能力。在指令排布時需要確定以下信息：

1）需要明確硬件時鐘的頻率和延遲。

上述硬件設(shè)計中的操作部件并不是立刻得出結(jié)果的，當(dāng)頻率不同時，各個操作部件可能會產(chǎn)生一定的延遲。在一般情況下，頻率越高，延遲越大。當(dāng)頻率控制在300～400 MHz 時，Mul 部件的延遲是5 拍，而Add 部件的延遲是3 拍。在指令排布時，設(shè)定Mul_lat 和Add_lat 分別表示Mul 部件和Add 部件的延遲。

2）需要明確系統(tǒng)中各個部件間的依賴關(guān)系。

如上文所述，系統(tǒng)中主要存在2 種依賴關(guān)系，即PE 內(nèi)部對角元之間的依賴和跨PE 的對角元之間的依賴，分別被保存在pe_inside_diag_dep和pe_outside_diag_dep 數(shù)組中，這里直接使用這2 個數(shù)組即可。

3）需要明確系統(tǒng)中各個部件間的沖突關(guān)系。

通過分析可以發(fā)現(xiàn)，當(dāng)前硬件系統(tǒng)中沒有訪存沖突，所有的內(nèi)存模塊至多只有2 個寫入端口和2 個讀出端口，可以通過分頻實現(xiàn)復(fù)用。然而，從上述求解過程可以看出，系統(tǒng)中存在不同操作對部件的爭用沖突，主要包括以下3 種：

（1）對Mul 部件的爭用。當(dāng)Mul 部件用于對角元求解時，就無法用于非對角元的乘法，反之亦然。

（2）對R部件的沖突。如果T部件在某一時刻被占用，那么相應(yīng)的R部件也被占用，無法再進(jìn)行傳輸。

（3）對Add 部件的爭用。當(dāng)Add 部件正在處理某一行的右端項更新時，同行的非對角元不能同時更新右端項，否則會造成讀寫沖突。

為了解決沖突問題，需要為上述可能發(fā)生沖突的部件設(shè)置部件占用標(biāo)記。某一時刻，當(dāng)部件占用標(biāo)記被置位時，意味著某個部件已經(jīng)被占用，無法再用其排布指令。本文將上述3 個部件占用標(biāo)記分別設(shè)置為Mul_occupy、R_occupy 和Add_occupy，其中，Add_occupy 是一個長度為Add_lat 的列表，用來標(biāo)記當(dāng)前加法部件正在處理的行。

圖6 所示為指令排布的算法流程，算法的偽代碼如算法1 所示。在排布時，乘法單元會優(yōu)先處理非對角元，只有當(dāng)非對角元全部處理完成后才會處理對角元，這樣做是為了盡可能多地在PE 上生成可供處理的對角元，以便保持PE 的滿載，提高整個系統(tǒng)的并行度。

圖6 指令排布算法流程Fig.6 The procedure of instruction layout algorithm

算法1指令排布算法

指令排布完成后就可以生成指令。系統(tǒng)中主要有6 種不同的操作，其中，Mul 部件占據(jù)了4 種操作，Add 部件和T部件各占據(jù)1 種操作，每個時刻都有可能出現(xiàn)所有操作部件同時工作的情況。由于各操作的取數(shù)空間基址和偏移都不同，需要將多種不同的操作整合起來生成指令。生成的指令放在緩沖M中，每個存儲單元的位寬為128 位，指令格式中的每個字段的位、名稱和含義如表1 所示。

表1 指令及其含義Table 1 Instructions and their meanings

2.5 硬件空間限制

在硬件實現(xiàn)的過程中，通常都存在一定的空間限制。由于緩沖區(qū)A、A'、B、B'、B"都使用片上內(nèi)存實現(xiàn)，導(dǎo)致其空間相對受限，也意味著求解的矩陣大小是受限的。

緩沖區(qū)A和B分別存儲分布在某個PE 上的對角元系數(shù)和右端項（包括求解得到的x）。假設(shè)對圖的分割相對均勻，則A和B的大小約為（N/PE_num），其中，N表示矩陣階數(shù)，PE_num 表示PE 的個數(shù)。

A'存儲非對角元系數(shù)，B'暫存乘以同列未知數(shù)x后的非對角元。非對角元不作劃分，而是直接分配到所在行對角元的PE 上。假設(shè)分布均勻，則A'和B'的大小約為（nnz/PE_num)，其中，nnz 表示非零元的個數(shù)。

B"存儲所有PE 間求解得到并且需要傳輸?shù)膞。B''的大小受矩陣稀疏度和PE_num 影響，矩陣越稠密，PE 的個數(shù)越多（劃分越多），B''則越大。B"最大不會超過N，實際大小依賴于具體矩陣，一般來說會比N小得多。

硬件總的空間占用約為(N+nnz+N×PE_num）個實數(shù)或復(fù)數(shù)，在使用算法前，需要評估FPGA 上提供的空間是否能滿足實際的空間需求。

3 性能測試分析

算法求解的過程可以分為3 個階段：

1）第一階段是矩陣求解初期的稀疏部分，此時絕大部分PE（大于80%）能夠找到就緒的對角元或非對角元，用于運行乘法部件，即使有依賴元暫未求解或硬件資源存在沖突，通常也只需等待1～2 拍就能夠繼續(xù)運行。

2）第二階段是矩陣求解中期的相對稀疏部分，此時只有部分PE（大于10%而小于80%）能夠運行乘法部件，剩余的PE 由于對角元未就緒（需要等待同行非對角元處理完成）或找不到非對角元（需要等待求解出的未知數(shù)x）進(jìn)行處理，因此只能等待。

3）第三階段是矩陣求解后期的稠密部分，此時只有極少部分PE（小于10%）能夠運行乘法部件，多個未知數(shù)x求解時存在強(qiáng)相關(guān)性，求解過程串行化明顯增加。

圖7 所示為FPGA 加速器的邏輯結(jié)構(gòu)。本文對應(yīng)的項目在Xilinx XCVU37P（8 GB HBM、4 200 萬等效邏輯門）FPGA 上實現(xiàn)硬件加速器，其與主機(jī)間采用PCIE4.0 X16 接口（帶寬為64 GB/s）連接。主機(jī)采用20 核Intel Xeon Gold 5320H（4.3 GHz）處理器，128 GB 3 200 Mb/s DDR4配置，運 行Linux 操作系統(tǒng)，移植和適配了電力系統(tǒng)仿真所需的庫環(huán)境。FPGA 加速器實物如圖8 所示。本文硬件結(jié)構(gòu)設(shè)計開銷情況如表2 所示。

圖7 FPGA 加速器邏輯結(jié)構(gòu)Fig.7 Logic structure of FPGA accelerator

圖8 FPGA 加速器實物示意圖Fig.8 Physical diagram of FPGA accelerator

表2 硬件資源利用情況Table 2 Utilization of hardware resource

本文對2 個典型算例進(jìn)行測試，2 個算例均來源于實際的電網(wǎng)模型。算例1 的矩陣大小為10 188×10 188，非零元為25 720 個；算例2 的矩陣大小為21 464×21 464，非零元為121 890 個。算例測試結(jié)果如表3 所示。

表3 2 個典型算例測試結(jié)果Table 3 Test results of two typical examples

從表3 可以看出，在64 個PE 配置的FPGA 加速器中，算例1 約耗時1 251 拍，則每拍能夠處理約20 個非零元，算例2 約耗時4 672 拍，每拍能夠處理約25 個非零元。

將當(dāng)前國網(wǎng)電力系統(tǒng)電磁暫態(tài)仿真程序中計算最為密集的稀疏矩陣消元求解核心段作為實際測試對象，對基于FPGA 和基于傳統(tǒng)CPU/GPU 的2 種環(huán)境進(jìn)行對比測試和分析。當(dāng)FPGA 加速器系統(tǒng)在300 MHz 頻率、256 個PE 配置下時，針對上述核心段的處理效率能達(dá)到30.84 GFLOPs。與此同時，將該電磁暫態(tài)仿真程序移植到20 核Intel Xeon Gold 5320H（4.3 GHz）處理器平臺上，進(jìn)行多核上的多線程并行優(yōu)化，同時降低大數(shù)據(jù)量離散訪問導(dǎo)致的CACHE 緩存顛簸，減少訪存延遲。處理器中4 核用于操作系統(tǒng)和驅(qū)動運行，其余16 核用于優(yōu)化后的稀疏矩陣并行運算，其處理效率僅為5.22 GFLOPs。基于FPGA 算法的實測性能是傳統(tǒng)CPU/GPU 求解算法的5.9 倍，加速效果顯著。

雖然本文所提算法相較傳統(tǒng)CPU/GPU 算法有明顯的加速效果，但仍然存在優(yōu)化的空間。從上述3 個求解階段來看，第一階段已經(jīng)充分利用了所有的硬件計算資源，幾乎沒有可以優(yōu)化的余地；第二階段的節(jié)拍占比最大，隨著矩陣規(guī)模的增大，其增長速度也最快，存在優(yōu)化算法調(diào)度的可能，主要優(yōu)化思路是更好地進(jìn)行矩陣的劃分和映射，進(jìn)一步減小通信代價，提升算法并行度；第三階段主要是矩陣的稠密部分，可以使用逆矩陣乘法進(jìn)行優(yōu)化。在LU 分解時，通過行列調(diào)整可以使得矩陣L右下角局部稠密，對于稠密的部分，可以求其逆矩陣，通過逆矩陣乘法進(jìn)行求解。由于矩陣乘法可以在多個PE 上完全并行，因此當(dāng)矩陣有一定的稠密度時，用逆矩陣乘法進(jìn)行求解相比直接法更具性能優(yōu)勢。

4 結(jié)束語

本文提出一種基于FPGA 硬件的靜態(tài)調(diào)度優(yōu)化算法，利用該算法實現(xiàn)了下三角稀疏矩陣的求解。通過對稀疏矩陣直接法求解步驟的分解和對稀疏矩陣的解析排布，設(shè)計一種節(jié)拍級的靜態(tài)調(diào)度流程，以充分利用FPGA 的硬件資源獲取較高的求解效率和硬件利用率。性能測試結(jié)果驗證了該算法的高效性，對于在FPGA 上實現(xiàn)類似的基于圖劃分的隱式數(shù)據(jù)并行算法具有一定的參考意義。下一步擬對LU 分解中的右下角局部稠密矩陣進(jìn)行優(yōu)化，無需修改現(xiàn)有的硬件拓?fù)?，僅增加稠密懲罰操作指令，在軟件方面，對稠密部分進(jìn)行求逆并均勻劃分求得的逆矩陣，然后通過乘法來求解該矩陣。