楊淑丹 董方敏

（三峽大學(xué)計(jì)算機(jī)與信息學(xué)院 宜昌 443002）

1 引言

電力系統(tǒng)潮流計(jì)算在電力系統(tǒng)規(guī)劃預(yù)測(cè)、生產(chǎn)運(yùn)行、調(diào)度管理、科學(xué)計(jì)算中有著廣泛的應(yīng)用。隨著電力網(wǎng)絡(luò)大規(guī)模、強(qiáng)耦合的特點(diǎn)日益突出，傳統(tǒng)的串行潮流計(jì)算已經(jīng)無(wú)法滿足大電網(wǎng)仿真和實(shí)時(shí)模擬的要求，高效的并行潮流算法正成為研究的熱點(diǎn)。從數(shù)學(xué)角度來(lái)看，潮流計(jì)算是求解一組由有功、無(wú)功功率、電壓和功角構(gòu)成的非線性方程組，求解時(shí)間占達(dá)整個(gè)潮流計(jì)算時(shí)間的70%以上，因此對(duì)潮流并行計(jì)算中的方程組求解方法的研究是實(shí)現(xiàn)實(shí)時(shí)仿真計(jì)算的關(guān)鍵。

本文針對(duì)電力系統(tǒng)潮流并行計(jì)算中的方程組求解這一核心問(wèn)題，分析了潮流并行計(jì)算大型矩陣常用算法，通過(guò)對(duì)稀疏矩陣壓縮為稠密矩陣的算法以及矩陣分塊算法進(jìn)行性能比較，給出了不同算法的適用范圍；同時(shí)，對(duì)常用的潮流計(jì)算線性方程組求解并行算法進(jìn)行探討，對(duì)比分析了共軛梯度法、廣義極小殘差法等不同算法的優(yōu)點(diǎn)和局限性；并對(duì)進(jìn)一步的研究方向進(jìn)行了展望。

2 電力系統(tǒng)潮流計(jì)算模型

電網(wǎng)穩(wěn)定狀態(tài)特性通常采用電力系統(tǒng)潮流計(jì)算進(jìn)行描述。設(shè)網(wǎng)絡(luò)節(jié)點(diǎn)數(shù)為n，PQ節(jié)點(diǎn)數(shù)為m-1，PV節(jié)點(diǎn)數(shù)為n-m，則可以將潮流方程表示成為以下直角坐標(biāo)形式：

式（1）是第i個(gè)電力系統(tǒng)有功功率方程式，共有n-1個(gè)，式（2）是第i個(gè)無(wú)功功率方程式，共有n-m-1個(gè)。其中，Pi為PQ和PV節(jié)點(diǎn)注入的有功功率，Qi為PV節(jié)點(diǎn)注入的無(wú)功功率，ei和 fi分別為迭代過(guò)程中節(jié)點(diǎn)電壓的實(shí)部、虛部，ej和 fj分別為平衡點(diǎn)電壓的實(shí)部、虛部，Gij為節(jié)點(diǎn)的自導(dǎo)，Bij為節(jié)點(diǎn)間的互導(dǎo)。

求解潮流計(jì)算的方法有很多，最常用的算法為牛頓拉夫遜法。其修正方程式如下

式中，ΔPi、ΔQi為節(jié)點(diǎn)的注入有功功率和無(wú)功功率，Δei、Δfi為節(jié)點(diǎn)電壓不平衡量的實(shí)部和虛部，i=1，…，n。修正方程也可簡(jiǎn)寫(xiě)為

其中，Δf為不平衡量的列向量，J為雅可比矩陣，Δx為節(jié)點(diǎn)電壓列向量。根據(jù)上述公式，可以將電力系統(tǒng)潮流計(jì)算轉(zhuǎn)化為對(duì)線性方程組的求解問(wèn)題，即

在大規(guī)模電網(wǎng)中，由于節(jié)點(diǎn)之間的連接線遠(yuǎn)少于電力系統(tǒng)，大型系數(shù)矩陣 A一般為高度稀疏陣。因此，通常采用不同的大型稀疏陣壓縮算法來(lái)節(jié)省存儲(chǔ)空間，以及通過(guò)矩陣分塊方法提高潮流并行計(jì)算效率；同時(shí)，可對(duì)線性方程組求解算法進(jìn)行優(yōu)化，針對(duì)不同情況選用合適的并行算法進(jìn)行求解，提高運(yùn)行效率。本文將從大型稀疏矩陣處理算法以及方程組并行求解算法兩個(gè)方面，來(lái)探討分析電力系統(tǒng)潮流并行計(jì)算中方程組求解的問(wèn)題。

3 潮流并行計(jì)算矩陣處理算法

3.1 稀疏矩陣壓縮為稠密矩陣算法

為了減少存儲(chǔ)空間，便于方程組并行求解，通常采用稀疏陣壓縮算法改進(jìn)方程組的求解，下面將依次介紹COO，ELL，HYB三種壓縮方法。

1）COO（Coordinate form）存儲(chǔ)

COO存儲(chǔ)［1～2］是采用三元組對(duì)稀疏矩陣的非零元進(jìn)行存儲(chǔ)，其中row，column，value分別代表元素所在的行、列和元素的值，三元組分別用三個(gè)數(shù)組進(jìn)行存儲(chǔ)。COO存儲(chǔ)不存在數(shù)據(jù)填充的問(wèn)題，壓縮比較徹底，但是缺乏行、列的規(guī)整性，因此并行時(shí)的性能一般。

2）ELL（ELL-pack form）存儲(chǔ)

ELL存儲(chǔ)［3］是使用一個(gè)n×k的矩陣進(jìn)行存儲(chǔ)，k指非零元素最多行的非零元素?cái)?shù)目，ELL存儲(chǔ)不適用于行之間非零元相差很大和非零元很少的稀疏矩陣。

3）HYB（Hybrid form）存儲(chǔ)

HYB存儲(chǔ)［4］是上述兩種存儲(chǔ)方式的綜合，可實(shí)現(xiàn)GPU+CPU同時(shí)進(jìn)行運(yùn)算，來(lái)減少運(yùn)算時(shí)間。其實(shí)現(xiàn)算法可描述如下：

輸入：稀疏矩陣A，稀疏矩陣行數(shù)N，列數(shù)M，行非零元素的下限值min；

若每行中非零元數(shù)目＜min

則非零元存儲(chǔ)在ELL格式中，在GPU中運(yùn)算；

否則

將其存儲(chǔ)在COO格式中，在CPU中運(yùn)算。

但該算法存在任務(wù)分配不均、GPU和CPU的訪問(wèn)存儲(chǔ)器和調(diào)度模式不同、而容易造成時(shí)間浪費(fèi)等問(wèn)題，通?？蓮臄?shù)據(jù)分割、數(shù)據(jù)共享、存儲(chǔ)算法等方面進(jìn)行改進(jìn)。

表1對(duì)上述三種矩陣壓縮存儲(chǔ)方法進(jìn)行了對(duì)比。

表1 稀疏矩陣壓縮為稠密矩陣的算法對(duì)比

3.2 矩陣分塊算法

采用矩陣分塊算法，能夠?qū)⒋笮途仃嚪指畛啥鄠€(gè)小型矩陣，有利于進(jìn)行并行運(yùn)算，減輕運(yùn)算壓力。下面分別介紹對(duì)角加邊法，單、雙對(duì)角加邊法這三種分塊算法：

1）對(duì)角加邊法矩陣A進(jìn)行對(duì)角加邊變換（block bordered diag?onal form，BBDF）［5～7］的方式大致可以分為三類(lèi)：節(jié)點(diǎn)撕裂法，支路切割法，統(tǒng)一網(wǎng)絡(luò)分塊法。下面分別介紹這三種分塊方法：

（1）節(jié)點(diǎn)撕裂法

節(jié)點(diǎn)撕裂法［8］是將兩個(gè)系統(tǒng)中交集部分取一合適節(jié)點(diǎn)，進(jìn)行撕裂，使得兩個(gè)系統(tǒng)都多一個(gè)節(jié)點(diǎn)，使得研究的系統(tǒng)簡(jiǎn)單方便，如圖1所示。

圖1 節(jié)點(diǎn)撕裂法

其中S1和S2為兩個(gè)區(qū)域電網(wǎng)，SA為由邊界節(jié)點(diǎn)構(gòu)成的上級(jí)電網(wǎng)，S1in和S2in為兩個(gè)區(qū)域電網(wǎng)的內(nèi)部節(jié)點(diǎn)，A為邊界節(jié)點(diǎn)的集合，A1，A2，A3為同一邊界節(jié)點(diǎn)在不同區(qū)域的分區(qū)，節(jié)點(diǎn)撕裂后平衡狀態(tài)方程為

式中，uA1，uA2為兩個(gè)分裂區(qū)域的電壓，θ1，θ2為分裂區(qū)的相角，PA1，QA1，PA2，QA2分別為兩個(gè)分裂區(qū)的有功和無(wú)功功率。

（2）支路切割法

支路切割法［9～10］即在子網(wǎng)絡(luò)之間的聯(lián)絡(luò)線處加入電流源來(lái)代替各子網(wǎng)絡(luò)間的耦合狀態(tài)，電流源矢量作為該方法的協(xié)調(diào)變量，如圖2～3所示。

根據(jù)推導(dǎo)，最終得到第i個(gè)子網(wǎng)絡(luò)的注入電流的不平衡方程：

圖2 互聯(lián)電力系統(tǒng)

圖3 互聯(lián)電力系統(tǒng)的分解

（3）統(tǒng)一網(wǎng)絡(luò)分塊法

統(tǒng)一網(wǎng)絡(luò)分塊法是對(duì)大規(guī)模電力系統(tǒng)進(jìn)行網(wǎng)絡(luò)分塊，利用節(jié)點(diǎn)編號(hào)優(yōu)先順序的理論來(lái)對(duì)網(wǎng)絡(luò)進(jìn)行劃分。通常有靜態(tài)優(yōu)化法，半動(dòng)態(tài)優(yōu)化法，動(dòng)態(tài)優(yōu)化法三種編號(hào)方法［11］，半動(dòng)態(tài)優(yōu)化法使用最為廣泛。下面以五節(jié)點(diǎn)網(wǎng)絡(luò)為例依次介紹這三種編號(hào)方法：

圖4 有向A圖

①靜態(tài)優(yōu)化法：按照每個(gè)節(jié)點(diǎn)出線度由小到大的順序進(jìn)行編號(hào)，則有向A圖編號(hào)順序即為靜態(tài)優(yōu)化編號(hào)方法。

②半動(dòng)態(tài)優(yōu)化法：先編號(hào)出線度最小的節(jié)點(diǎn)，進(jìn)行消去，其次為消去節(jié)點(diǎn)后剩余節(jié)點(diǎn)出線度最少的節(jié)點(diǎn)進(jìn)行編號(hào)，依次類(lèi)推…，圖中1，2節(jié)點(diǎn)編號(hào)不用更改，消去這兩個(gè)節(jié)點(diǎn)節(jié)點(diǎn)以后，剩余節(jié)點(diǎn)的出線度相同，所以可任意編號(hào)。

③動(dòng)態(tài)優(yōu)化法：對(duì)節(jié)點(diǎn)先不進(jìn)行編號(hào)，對(duì)任意節(jié)點(diǎn)進(jìn)行消去，為消去后出現(xiàn)支路最少的節(jié)點(diǎn)先編號(hào)，依次類(lèi)推…，圖中1，2，3，4節(jié)點(diǎn)，由于消去不產(chǎn)生注入元，可任意順序消去。

對(duì)角加邊法中三種方法各有特點(diǎn)，表2對(duì)這三種方法特點(diǎn)進(jìn)行比較分析。

2）單、雙邊塊對(duì)角法

采用對(duì)角加邊分塊法，能夠?qū)⒋笮途仃嚪指畛啥鄠€(gè)小塊矩陣，減輕計(jì)算壓力，便于并行求解。但由于電力網(wǎng)絡(luò)的不確定性，該方法無(wú)法滿足求解要求，因此在該方法基礎(chǔ)上提出單、雙邊塊對(duì)角變換法（singly/doubly bordered block diagonal form，SBBD/DBBD）［12］，下面依次介紹這兩種分塊方法。

表2 BBDF中矩陣分塊方法的對(duì)比

（1）對(duì)A進(jìn)行雙邊塊對(duì)角DDBD變換

對(duì)矩陣A進(jìn)行變換后轉(zhuǎn)置，則DBBD的表現(xiàn)形

式為

其中塊 A11，A22，…，ANN是nN×nN，邊界矩陣C1是n1×p矩陣，R1是 p×n1矩陣，且 p≤n。邊界塊中的每一列叫做耦合列，每一行叫做耦合行，對(duì)A進(jìn)行DDBD變換時(shí)，只有當(dāng)耦合行和列中的 p很小時(shí)，才能滿足協(xié)調(diào)矩陣All之間的協(xié)調(diào)性。其中A矩陣中的元素aij需要滿足的邊界條件是， 門(mén) 檻 值u∈(0，1)。采用該方法，容易忽略列邊界條件，導(dǎo)致u無(wú)法滿足條件，矩陣變換失敗。因此，只有在不消除邊界條件的前提下，通過(guò)移動(dòng)行、列的位置，DBBD的變換方式才有效。

（2）對(duì)A進(jìn)行SBBD變換

采用DBBD變換時(shí)，移動(dòng)矩陣行列會(huì)增加矩陣的維數(shù)，增加方程組求解的迭代次數(shù)。為了簡(jiǎn)化求解方法，對(duì)矩陣進(jìn)行更進(jìn)一步變換，即SBBD變換，其表現(xiàn)形式如下：

其中塊 All是 ml×nl矩陣，ml＞nl，且邊界矩陣為ml×p(p≤nl)，矩陣 A中的元素aij需要滿足的邊界條件是，在進(jìn)行變換時(shí)，應(yīng)盡可能減少耦合列的數(shù)量以保證分解塊的穩(wěn)定性。

4 潮流計(jì)算線性方程組求解并行算法

線性方程組并行求解最常用的處理方法有共軛梯度法，不完全Cholesky預(yù)處理共軛梯度法，重開(kāi)始廣義極小殘差法等。下面分別對(duì)這些并行算法進(jìn)行介紹。

4.1 共軛梯度法

共軛梯度法（Conjugate Gradient，CG）［13］也稱(chēng)共軛斜量法，是求取正定線性方程組首選方法。該方法最大的特點(diǎn)為當(dāng)求解的線性方程組階數(shù)很高時(shí)，通常只需經(jīng)過(guò)比階數(shù)小很多的迭代次數(shù)，就可以取得方程組的解。其算法描述如下：

對(duì)于 Ax=b線性方程，輸入 A，b，其中 A為n×n矩陣，b為n階向量，輸出x

1）給定初值 x(0)=0，d(0)=0，r(0)=-b，迭代精度 ε＞0；

2）在第k次迭代時(shí)，求取x(k)分五步：

（1）計(jì)算殘向量：r(k)=Ax(k-1)-b；

（2）計(jì)算修正方向向量：

（3）計(jì)算迭代步長(zhǎng)：

（4）求出新近似量：x(k)=x(k-1)+ξ(k)d(k)；

（5）如不滿足迭代精度要求，則重復(fù)以上步驟。

該算法適用于A為正定矩陣的情況，若A為非對(duì)稱(chēng)非正定矩陣，LDLT分解矩陣不一定存在，使用該算法可能會(huì)導(dǎo)致程序中斷。

為了滿足實(shí)際需求，求解更大規(guī)模的線性方程組，可在共軛梯度算法上加入不同種類(lèi)預(yù)處理。下面介紹一種常見(jiàn)的預(yù)處理方法，即不完全Cholesky預(yù)處理共軛梯度法（Modified Incomplete Cholesky Preconditioner Conjugate Gradient，MIC PCG）［14］，算法和PCG基本一致，這里不再贅述。

MIC PCG較PCG相比，迭代次數(shù)顯著減少，在進(jìn)行向量?jī)?nèi)積運(yùn)算時(shí)，共享內(nèi)存和全局內(nèi)存達(dá)到最優(yōu)，能夠提高2倍以上的加速比。

4.2 廣義極小殘差法

重開(kāi)始廣義極小殘差法（Generalized Minimal Residual Method，GMRES）［15～16］是 krylov子空間法的一種算法，是求解大規(guī)模線性方程組最受歡迎的一種，具有收斂速度快且比較穩(wěn)定的特征。該算法能很好地在GPU上并行求解方程組，通過(guò)對(duì)顯存中線程塊、線程數(shù)的合理分配，使得每個(gè)線程能夠同時(shí)并行求解方程，達(dá)到性能最優(yōu)。

其算法描述如下：

對(duì)于 Ax=b線性方程，輸入 A，b，其中 A為n×n矩陣，b為n階向量，輸出x

1）給定初值 x(0)=0，r(0)=-b，迭代精度ε＞0；

2）選取合適的m，使用Arnoldi迭代方法求?得Vn和，其中m是標(biāo)準(zhǔn)正交基向量的個(gè)數(shù)，Vn是由標(biāo)準(zhǔn)正交基組成的m×n矩陣，H~n是(n +1)×n階的上Hessenberg矩陣；

3）極小化殘余量范數(shù)，即 ‖H~nyn-βe1‖，求取使殘余量范數(shù)最小的 yn。其中 β=‖b- Ax0‖，e1=(1 ，0，…0 )；

4）根據(jù) xn=Vnyn，求得 xn；

5）如不滿足迭代要求，則重復(fù)以上步驟。

根據(jù)以上算法，可以得出每次迭代所需要的計(jì)算次數(shù)以及它的復(fù)雜度，設(shè)矩陣A為n×n階稀疏陣，m個(gè)標(biāo)準(zhǔn)正交基，每行有k個(gè)非零元，則可以得出表3。

從表中可以看出，稀疏矩陣與向量乘積、向量?jī)?nèi)積與范數(shù)在整體運(yùn)算中是最耗時(shí)間的。當(dāng)k比m 2+1大時(shí)，稀疏矩陣與向量乘積運(yùn)算主要影響整個(gè)程序的性能，當(dāng)k比m 2+1小時(shí)，向量乘積與范數(shù)運(yùn)算主要影響整個(gè)程序的性能?？紤]到在GPU運(yùn)算過(guò)程中，向量?jī)?nèi)積與范數(shù)運(yùn)算的瓶頸比稀疏矩陣與向量乘積運(yùn)算少很多，因此加速效果主要受m大小所影響，m越大加速效果越好。但受浮點(diǎn)精度的影響，m不宜取值過(guò)大，取值過(guò)小又會(huì)導(dǎo)致收斂停滯。通常根據(jù)GPU的性能、系數(shù)矩陣A的特征、浮點(diǎn)精度的大小等因素折中選取合適的m值。表4對(duì)上述并行算法進(jìn)行了比較。

表3 每次迭代中不同計(jì)算所需的計(jì)算次數(shù)和復(fù)雜度

表4 并行算法性能比較

5 結(jié)語(yǔ)

本文討論對(duì)比了電力系統(tǒng)潮流并行計(jì)算中方程組的常用求解方法。其中，稀疏矩陣壓縮算法適用于節(jié)點(diǎn)間物理連線較少的電力系統(tǒng)，可通過(guò)壓縮矩陣來(lái)節(jié)省顯存存儲(chǔ)空間，提高并行速度；矩陣分塊算法能進(jìn)一步降低問(wèn)題粒度，便于在GPU上實(shí)現(xiàn)多線程塊并行求解線性方程組，提高并行化效率；線性方程組并行算法，通過(guò)預(yù)處理降低條件數(shù)來(lái)提升收斂速度，并在迭代求解過(guò)程中實(shí)現(xiàn)矩陣與向量?jī)?nèi)積等的并行計(jì)算，實(shí)現(xiàn)求解效率的有效提升。

目前已有很多較為有效的電力系統(tǒng)潮流并行計(jì)算方法，但由于潮流并行計(jì)算方程組的求解是一個(gè)比較復(fù)雜的問(wèn)題，也還存在許多尚待進(jìn)一步完善和解決的問(wèn)題。如求解線性方程組時(shí)，在迭代次數(shù)減少的同時(shí)，如何解決收斂性能變差的問(wèn)題；如何通過(guò)線程塊的合理劃分，保證二者之間的平衡性，達(dá)到性能最優(yōu)的問(wèn)題；如何在降低時(shí)間復(fù)雜度的情況下，進(jìn)行網(wǎng)絡(luò)中大型矩陣壓縮算法的合理混合使用的問(wèn)題；以及對(duì)方程組并行求解添加預(yù)處理算法來(lái)提高并行效率等等，這些將是未來(lái)需要進(jìn)一步研究的方向。

［1］Hoang-Vu Dang，Bertil Schmidt.The Sliced COO Format for Sparse Matrix Vector Multiplication on CUDA enabled GPUs［J］.Procedia Computer Science，2012，9：57-61.

［2］李佳佳，張秀霞，譚光明，等.選擇稀疏矩陣乘法最優(yōu)存儲(chǔ)格式的研究［J］.計(jì)算機(jī)研究與發(fā)展，2014，51（4）：882-888.LI Jiajia，ZHANG Xiuxia，TAN Guangming，et al.Study of Choosing the Optimal Storage Format of Sparse Matrix Vector Multiplication［J］.Journal of Computer Research and Development，2014，51（4）：882-888.

［3］商磊.大規(guī)模線性方程組求解的并行算法及應(yīng)用［D］.西安：西北工業(yè)大學(xué)，2007，20-23.SHANG Lei.Parallel Algorithms and Applications for Solv?ing Large scale Linear Equations［D］.Xi'an：Northwest Polytechnic University，2007，20-23.

［4］陽(yáng)王東，李肯立.基于HYB格式稀疏矩陣與向量乘在CPU+GPU異構(gòu)系統(tǒng)中的實(shí)現(xiàn)與優(yōu)化［J］.計(jì)算機(jī)工程與科學(xué)，2016，38（2）：205-206.YANG Wangdong，LI Kenli.Implementation and Optimiza?tion of HYB based SPMV on CPU+GPU Heterogeneous Computing Systems［J］.Computer Engineering&Science，2016，38（2）：205-206.

［5］薛巍，舒繼武，王心豐，等.電力系統(tǒng)潮流并行算法的研究進(jìn)展［J］.清華大學(xué)學(xué)報(bào)（自然科學(xué)版），2002，42（9）：1194-1199.XUE Wei，SHU Jiwu，WANG Xinfeng，et al.Advance of parallel algorithms for power flow simulation［J］.J Tsingh ua Univ（Sci&Tech），2001，42（9）：1193-1199.

［6］楊挺，向文平，王洪濤，等.電力系統(tǒng)對(duì)角加邊模型的數(shù)據(jù)中心求解方法［J］.中國(guó)電機(jī)工程學(xué)報(bào)，2015，35（3）：512-514.YANG Ting，XIANG Wenping，WANG Hongtao，et al.An Algorithm for Solving the Block Bordered Diagonal Form of Electrical Power System in Data Center［J］.Proceedings of the CSEE，2015，35（3）：512-514.

［7］朱永興，張步涵.電力系統(tǒng)潮流分解協(xié)調(diào)并行計(jì)算［J］.電力系統(tǒng)及其自動(dòng)化學(xué)報(bào)，2010：22（5）.ZHU Yongxing，ZHANG Buhan.Decomposition and Coor?dination Power Flow Parallel Calculation［J］.Proceedings of the CSUEPSA，2010：22（5）.

［8］張伯明，陳壽孫.高等電力系統(tǒng)網(wǎng)絡(luò)分析［M］.北京：清華大學(xué)出版社，1996：144-165.ZHANG Boming，CHEN Shousun.Network Analysis of Ad?vanced Power System［M］.Beijing：Tsinghua University Press，1996：144-165.

［9］黃彥全，肖建，劉蘭，等.基于支路切割方法的電力系統(tǒng)潮流并行協(xié)調(diào)算法［J］.電網(wǎng)技術(shù)，2006，30（4）：21-22.HUANG Yanquan，XIAO Jian，LIU Lan，et al.A Coordina?tional Parallel Algorithm for Power Flow Calculation Based on Branch Cutting［J］.Power System Technology，2006，30（4）：21-22.

［10］龔仁敏，張永浩，翟學(xué)峰，等.面向整定計(jì)算的大型電力網(wǎng)絡(luò)分塊計(jì)算的實(shí)用化通用計(jì)算方法［J］.電力信息化，2007年增刊：42-44.GONG Renmin，ZHANG Yonghao，ZHAI Xuefeng，et al.A Practical General Method of Large Scale Power Net?work Partition Computing for Setting Calculation［J］.Electric power Information，2007 Supplement：42-44.

［11］陳珩.電力系統(tǒng)穩(wěn)態(tài)分析［M］.3版.北京：中國(guó)電力出版社，2007：160-163.CHEN Heng.Steady-State Analysis of Power System［M］.The Third Release.Beijing：China Electric Power Press，2007：160-163.

［12］Iain S.Duff，Jennifer A.Scott.Stabilized block diagonal forms for parallel sparse solvers［R］.Council for the Cen?tral Laboratory of the Research Councils，2004.

［13］吳長(zhǎng)江.基于CUDA的大規(guī)模線性稀疏方程組求解器的設(shè)計(jì)［D］.成都：中國(guó)電子科技大學(xué)，2013：27-30.WU Changjiang.Design of the of the solver for large scale linear sparse equations based on CUDA［D］.Chengdou：University of Electronic Science and Technol?ogy of China，2013：27-30.

［14］Jiaquan Gao，Ronghua Liang，Jun Wang.Research on the conjugate gradient algorithm with a modified incom?plete Cholesky preconditioner on GPU［J］.J.Parallel Dis?trib.Comput.74，2014.

［15］Kai He，Sheldon X.-D.Tan，Hengyang Zhao，et al.Par?allel GMRES solver for fast analysis of large linear dy?namic systems on GPU platforms［J］.INTEGRATION，the VLSI journal，2016（52）.

［16］張慧，于春肖，白雪婷，等.Krylov子空間E-變換GMRES（m）算法［J］.遼寧工程技術(shù)大學(xué)學(xué)報(bào)（自然科學(xué)版），2014，33（9）：1290-1291.ZHANG Hui，YU Chunxiao，BAI Xueting，et al.E-transform GMRES（m）algorithm based on krylov sub?space［J］.Journal of Liaoning Technical University（Natu?ral Science），2014，33（9）：1290-1291.