本發(fā)明屬于電力系統(tǒng)高性能計(jì)算應(yīng)用領(lǐng)域，尤其涉及一種大量同構(gòu)稀疏矩陣的GPU加速Q(mào)R分解方法。

背景技術(shù)：

潮流計(jì)算是電力系統(tǒng)中應(yīng)用最廣泛、最基本和最重要的一種電氣運(yùn)算。在電力系統(tǒng)運(yùn)行方式和規(guī)劃方案的研究中，都需要進(jìn)行潮流計(jì)算以比較運(yùn)行方式或規(guī)劃供電方案的可行性、可靠性和經(jīng)濟(jì)性，在電力系統(tǒng)運(yùn)行狀態(tài)的實(shí)時(shí)監(jiān)控中，需要采用在線潮流計(jì)算。傳統(tǒng)的牛頓拉夫遜法潮流計(jì)算中，修正方程組求解時(shí)間占潮流計(jì)算時(shí)間的70％，修正方程組求解的計(jì)算速度影響程序的整體性能。

靜態(tài)安全分析中的故障潮流是在基態(tài)電網(wǎng)上發(fā)生若干斷線產(chǎn)生的，任何一個(gè)N-1斷線故障都對(duì)應(yīng)一個(gè)潮流，這些故障潮流的稀疏格式都是基態(tài)潮流稀疏格式的一個(gè)子集，其稀疏線性方程組具有統(tǒng)一的稀疏結(jié)構(gòu)，求解計(jì)算具有并行性。對(duì)批量方程組系數(shù)矩陣的相同稀疏結(jié)構(gòu)J矩陣進(jìn)行QR符號(hào)分解后，得到Household變換矩陣V和上三角矩陣R矩陣的稀疏結(jié)構(gòu)，根據(jù)R矩陣的稀疏結(jié)構(gòu)，對(duì)矩陣J各列進(jìn)行并行化分層。其中每層中的列的計(jì)算相互獨(dú)立，沒有依賴關(guān)系，天然可以被并行的計(jì)算處理，適合GPU加速。因此通過CPU和GPU之間合理的調(diào)度可以快速完成方程組系數(shù)矩陣進(jìn)行QR分解，并求解稀疏線性方程組，國(guó)內(nèi)外學(xué)者已經(jīng)開始對(duì)GPU進(jìn)行稀疏線性方程組加速求解的方法進(jìn)行了研究，但是沒有深入的優(yōu)化線程設(shè)計(jì)，單純從計(jì)算量的分配上研究計(jì)算線程設(shè)計(jì)，對(duì)線程計(jì)算方式，數(shù)據(jù)索引方式?jīng)]有進(jìn)行深入研究，無法使程序充分發(fā)揮GPU的優(yōu)勢(shì)。

因此，亟待解決上述問題。

技術(shù)實(shí)現(xiàn)要素：

發(fā)明目的：針對(duì)現(xiàn)有技術(shù)的不足，本發(fā)明提供了一種適用于靜態(tài)安全分析中批量潮流修正方程組的雅可比矩陣QR分解，可提高潮流計(jì)算速度，為在線分析提供基礎(chǔ)的大量同構(gòu)稀疏矩陣的GPU加速Q(mào)R分解方法。

潮流計(jì)算：電力學(xué)名詞，指在給定電力系統(tǒng)網(wǎng)絡(luò)拓?fù)洹⒃?shù)和發(fā)電、負(fù)荷參量條件下，計(jì)算有功功率、無功功率及電壓在電力網(wǎng)中的分布。

GPU：圖形處理器(英語：Graphics Processing Unit，縮寫：GPU)。

本發(fā)明公開了一種大量同構(gòu)稀疏矩陣的GPU加速Q(mào)R分解方法，所述方法包括如下步驟：

(1)稀疏結(jié)構(gòu)相同的一系列n階矩陣A₁～A_N構(gòu)成大量同構(gòu)稀疏矩陣，在CPU上對(duì)稀疏矩陣A₁進(jìn)行QR符號(hào)分解，得到Household變換矩陣V₁和上三角矩陣R₁的稀疏結(jié)構(gòu)；根據(jù)R₁矩陣的稀疏結(jié)構(gòu)，對(duì)矩陣A₁各列進(jìn)行并行化分層，且A₁～A_N具有相同的Household變換矩陣稀疏結(jié)構(gòu)V₁、上三角矩陣稀疏結(jié)構(gòu)R₁和并行化分層結(jié)果；

(2)CPU將QR分解所需數(shù)據(jù)傳輸給GPU；

(3)任務(wù)分配和設(shè)備內(nèi)存優(yōu)化：將對(duì)矩陣A₁～A_N的QR分解任務(wù)分配到GPU上的大量線程中執(zhí)行，并根據(jù)合并訪問原則優(yōu)化內(nèi)存使用；

(4)GPU中按層次遞增的順序計(jì)算分層QR分解內(nèi)核函數(shù)Batch_QR。

其中，所述步驟(1)中，并行化分層將矩陣A₁的n列分配到M層中，屬于同一層中的列可并行QR分解；每層包含的列的數(shù)量為L(zhǎng)(k)，k表示層號(hào)；存儲(chǔ)第k層中所有列號(hào)至映射表Map_k。

優(yōu)選的，所述步驟(2)中，所述的QR分解所需數(shù)據(jù)包括矩陣A₁～A_N、矩陣維度n、矩陣V₁和矩陣R₁的稀疏結(jié)構(gòu)、矩陣A₁的并行化分層結(jié)果。

再者，所述步驟(3)中，將N個(gè)同構(gòu)稀疏矩陣A₁～A_N的同一列的QR分解工作分配給同一個(gè)線程塊的不同線程處理；為保證合并內(nèi)存訪問，將矩陣A₁～A_N在內(nèi)存中連續(xù)存放組成一個(gè)邏輯上為N行的大矩陣，然后進(jìn)行轉(zhuǎn)置操作；將矩陣V₁～V_N在內(nèi)存中連續(xù)存放組成一個(gè)邏輯上為N行的大矩陣，然后進(jìn)行轉(zhuǎn)置操作。

進(jìn)一步，所述步驟(4)中，GPU中的內(nèi)核函數(shù)定義為Batch_QR<N_blocks，N_threads>，其中線程塊大小N_threads固定為128；線程塊數(shù)量N_blocks設(shè)計(jì)成變量，當(dāng)對(duì)第k層進(jìn)行計(jì)算時(shí)，線程塊數(shù)量N_blocks＝L(k)，總線程數(shù)量為：N_blocks×N_threads；按照層次遞增的順序，啟動(dòng)內(nèi)核函數(shù)Batch_QR<L(k)，N_threads>來分解屬于第k層的所有列；Batch_QR<L(k)，N_threads>的具體計(jì)算流程為：

(4.1)CUDA自動(dòng)為每個(gè)線程分配線程塊索引blockID和線程塊中的線程索引threadID；

(4.2)將blockID賦值給變量bid，第blockID號(hào)線程塊中的128個(gè)線程負(fù)責(zé)分解矩陣A₁～A_N的第j＝Map_k(bid)列，其中：第threadID號(hào)線程負(fù)責(zé)計(jì)算矩陣A_t的第j列QR分解，t＝threadID+m×128，(m＝0，1，…，N/128)；

(4.3)第bid號(hào)線程塊的threadID號(hào)線程中，判斷t是否小于N，小于繼續(xù)執(zhí)行，否則該線程退出運(yùn)行；

(4.4)變量i從1遞增到j(luò)-1，如果R₁(i，j)≠0，執(zhí)行以下計(jì)算：

1)計(jì)算變量β，計(jì)算公式為β＝2V_t(i：n，i)^T·A_t(i：n，j)；其中V_t(i：n，i)是Household變換矩陣V_t中第i列的第i～n行元素構(gòu)成的向量，A_t(i：n，j)是矩陣A_t中第j列的第i～n行元素構(gòu)成的向量；具體計(jì)算步驟：變量c從i遞增到n計(jì)算：β+＝2V_t(c，i)×A_t(c，j)；

2)采用公式A_t(i：n，j)＝A_t(i：n，j)–β×V_t(i：n，i)，更新矩陣A_t的第j列，具體計(jì)算步驟如下：變量c從i遞增到n計(jì)算：A_t(c，j)＝A_t(c，j)–β×V_t(c，i)；

(4.5)計(jì)算Household變換矩陣V_t的第j列：首先，采用公式a²＝A_t(j：n，j)^T·A_t(j：n，j)計(jì)算中間變量a，具體計(jì)算步驟：變量c從j遞增到n計(jì)算：a²+＝A_t(c，j)×A_t(c，j)；接著，計(jì)算，V_t(j：n，j)＝A_t(j：n，j)–ae_j(j：n)，其中是e_j是第j個(gè)元素為1的單位向量，具體計(jì)算步驟：變量c從j遞增到n計(jì)算：V_t(c，j)＝A_t(c，j)–ae_j(c)；然后，采用公式b²＝V_t(j：n，j)^T·V_t(j：n，j)計(jì)算中間變量b，具體計(jì)算步驟：變量c從j遞增到n計(jì)算：b²+＝V_t(c，j)×V_t(c，j)；最后，計(jì)算，V_t(j：n，j)＝V_t(j：n，j)/b，具體計(jì)算步驟：變量c從j遞增到n計(jì)算：V_t(c，j)＝V_t(c，j)/b；

(4.6)更新矩陣A_t的第j列：A_t(j，j)＝a，A_t(j+1：n，j)＝0；

(4.7)t＝t+128，返回(4.3)。

有益效果：與現(xiàn)有技術(shù)比，本發(fā)明的有益效果為：首先，本發(fā)明利用CPU控制程序的流程并處理基礎(chǔ)數(shù)據(jù)和GPU處理密集的浮點(diǎn)運(yùn)算相結(jié)合的模式提高了大量同構(gòu)稀疏矩陣的QR分解速度，解決了電力系統(tǒng)靜態(tài)安全性分析中潮流計(jì)算耗時(shí)大的問題，且采用CPU對(duì)大量同構(gòu)稀疏矩陣的相同稀疏格式A₁進(jìn)行QR符號(hào)分解，根據(jù)R₁矩陣的稀疏格式，可以減少不必要的浮點(diǎn)計(jì)算；其次，在CPU中將矩陣A₁進(jìn)行并行化分層，并將結(jié)果傳給GPU，減少GPU對(duì)邏輯操作的運(yùn)算；再者，將批量矩陣的QR分解工作分配到大量的線程中執(zhí)行，并根據(jù)GPU的訪存模式優(yōu)化設(shè)備內(nèi)存使用，使GPU實(shí)現(xiàn)合并訪存，內(nèi)存操作速度提高了接近16倍；最后GPU中按層次遞增的順序啟動(dòng)內(nèi)核函數(shù)Batch_QR，取得了單個(gè)稀疏矩陣QR分解的平均計(jì)算時(shí)間為1.67ms的效果。

附圖說明：

圖1為本發(fā)明內(nèi)核函數(shù)任務(wù)分配示意圖；

圖2為本發(fā)明的實(shí)例測(cè)試結(jié)果；

圖3為本發(fā)明的實(shí)例性能分析；

圖4為本發(fā)明的流程示意圖。

具體實(shí)施方式：

下面結(jié)合附圖對(duì)本發(fā)明的技術(shù)方案作進(jìn)一步說明。

如圖4所示，本發(fā)明一種大量同構(gòu)稀疏矩陣的GPU加速Q(mào)R分解方法，所述方法包括如下步驟：

(1)大量同構(gòu)稀疏矩陣指稀疏結(jié)構(gòu)相同的一系列n階矩陣A₁～A_N，在CPU上對(duì)其中稀疏矩陣A₁進(jìn)行QR符號(hào)分解，得到Household變換矩陣V₁和上三角矩陣R₁的稀疏結(jié)構(gòu)，符號(hào)分解之后的A₁矩陣的稀疏結(jié)構(gòu)等于V₁+R₁；根據(jù)R₁矩陣的稀疏結(jié)構(gòu)，對(duì)矩陣A₁各列進(jìn)行并行化分層；因?yàn)锳₁～A_N的稀疏結(jié)構(gòu)相同，所以A₁～A_N具有相同的Household變換矩陣稀疏結(jié)構(gòu)V₁、上三角矩陣稀疏結(jié)構(gòu)R₁和并行化分層結(jié)果；

(2)CPU將QR分解所需數(shù)據(jù)傳輸給GPU；

(3)任務(wù)分配和設(shè)備內(nèi)存優(yōu)化：將對(duì)矩陣A₁～A_N的QR分解任務(wù)分配到GPU上的大量線程中執(zhí)行，并根據(jù)合并訪問原則優(yōu)化內(nèi)存使用；

(4)GPU中按層次遞增的順序啟動(dòng)分層QR分解內(nèi)核函數(shù)Batch_QR。

一、CPU中對(duì)稀疏矩陣A₁進(jìn)行QR符號(hào)分解和并行化分層

大量同構(gòu)的一系列n階稀疏矩陣A₁～A_N具有相同的稀疏結(jié)構(gòu)，在CPU上對(duì)其中稀疏矩陣A₁進(jìn)行QR符號(hào)分解，得到Household變換矩陣V₁和上三角矩陣R₁的稀疏結(jié)構(gòu)，符號(hào)分解之后的A₁矩陣的稀疏結(jié)構(gòu)等于V₁+R₁；根據(jù)R₁矩陣的稀疏結(jié)構(gòu)，對(duì)矩陣A₁各列進(jìn)行并行化分層；因?yàn)锳₁～A_N的稀疏結(jié)構(gòu)相同，所以A₁～A_N具有相同的Household變換矩陣稀疏結(jié)構(gòu)V₁、上三角矩陣稀疏結(jié)構(gòu)R₁以及并行化分層結(jié)構(gòu)。并行化分層將矩陣A₁的n列分配到M層中，屬于同一層中的列可并行QR分解；每層包含的列的數(shù)量為L(zhǎng)(k)，k表示層號(hào)；存儲(chǔ)第k層中所有列號(hào)至映射表Map_k。

其中QR符號(hào)分解原理和并行化分層原理參見“Direct Methods for Sparse Linear Systems”Timothy A.Davis，SIAM，Philadelphia，2006。本專利使用的QR符號(hào)分解和并行化分層程序參見CSparse：a Concise Sparse Matrix package.VERSION 3.1.4，Copyright(c)2006-2014，Timothy A.Davis，Oct 10，2014。

二、CPU將QR分解所需數(shù)據(jù)傳輸給GPU

CPU讀取電網(wǎng)基礎(chǔ)數(shù)據(jù)，并將矩陣A₁的分層結(jié)果和電網(wǎng)基礎(chǔ)數(shù)據(jù)在內(nèi)核函數(shù)開始執(zhí)行之前一次性傳輸給GPU，減少CPU與GPU之間的數(shù)據(jù)交互。所需數(shù)據(jù)包括：矩陣A₁～A_N、矩陣維度n、矩陣V₁和R₁的稀疏結(jié)構(gòu)、分層數(shù)量M層，每層包含的列的數(shù)量為L(zhǎng)(k)，映射表Map_k。

三、任務(wù)分配和設(shè)備內(nèi)存優(yōu)化

將N個(gè)同構(gòu)稀疏矩陣A₁～A_N的同一列的QR分解工作分配給同一個(gè)線程塊的不同線程處理，具體的任務(wù)分配模式，如圖1所示；為保證合并內(nèi)存訪問，將矩陣A₁～A_N在內(nèi)存中連續(xù)存放組成一個(gè)邏輯上為N行的大矩陣，然后進(jìn)行轉(zhuǎn)置操作；將矩陣V₁～V_N在內(nèi)存中連續(xù)存放組成一個(gè)邏輯上為N行的大矩陣，然后進(jìn)行轉(zhuǎn)置操作。

四、GPU中按層次遞增的順序啟動(dòng)分層QR批處理分解內(nèi)核函數(shù)

GPU中的內(nèi)核函數(shù)定義為Batch_QR<N_blocks，N_threads>，其中線程塊大小N_threads固定為128；線程塊數(shù)量N_blocks設(shè)計(jì)成變量，當(dāng)對(duì)第k層進(jìn)行計(jì)算時(shí)，線程塊數(shù)量N_blocks＝L(k)，總線程數(shù)量為：N_blocks×N_threads；按照層次遞增的順序，調(diào)用內(nèi)核函數(shù)Batch_QR<Ln(k)，N_threads>來分解屬于第k層的所有列。

Batch_QR<Ln(k)，N_threads>的計(jì)算流程為：

(1)CUDA自動(dòng)為每個(gè)線程分配線程塊索引blockID和線程塊中的線程索引threadID；

(2)將blockID賦值給變量bid，第blockID號(hào)線程塊中的128個(gè)線程負(fù)責(zé)分解矩陣A₁～A_N的第j＝Map_k(bid)列，其中：第threadID號(hào)線程負(fù)責(zé)計(jì)算矩陣A_t的第j列QR分解，t＝threadID+m×128，(m＝0，1，…，N/128)；

(3)第bid號(hào)線程塊的threadID號(hào)線程中，判斷t是否小于N，小于繼續(xù)執(zhí)行，否則該線程退出運(yùn)行；

(4)變量i從1遞增到j(luò)-1，如果R₁(i，j)≠0，執(zhí)行以下計(jì)算：

1)計(jì)算變量β，計(jì)算公式為β＝2V_t(i：n，i)^T·A_t(i：n，j)；其中V_t(i：n，i)是Household變換矩陣V_t中第i列的第i～n行元素構(gòu)成的向量；A_t(i：n，j)是矩陣A_t中第j列的第i～n行元素構(gòu)成的向量；具體計(jì)算步驟如下：變量c從i遞增到n計(jì)算：β+＝2V_t(c，i)×A_t(c，j)；

2)采用公式A_t(i：n，j)＝A_t(i：n，j)–β×V_t(i：n，i)，更新矩陣A_t的第j列，具體步驟如下：變量c從i遞增到n計(jì)算：A_t(c，j)＝A_t(c，j)–β×V_t(c，i)；

(5)計(jì)算Household變換矩陣V_t的第j列：

首先，采用公式a²＝A_t(j：n，j)^T·A_t(j：n，j)計(jì)算中間變量a，具體計(jì)算步驟：變量c從j遞增到n計(jì)算：a²+＝A_t(c，j)×A_t(c，j)；

接著，計(jì)算，V_t(j：n，j)＝A_t(j：n，j)–ae_j(j：n)，其中是e_j是第j個(gè)元素為1的單位向量，具體計(jì)算步驟：變量c從j遞增到n計(jì)算：V_t(c，j)＝A_t(c，j)–ae_j(c)；

然后，采用公式b²＝V_t(j：n，j)^T·V_t(j：n，j)計(jì)算中間變量b，具體計(jì)算步驟：變量c從j遞增到n計(jì)算：b²+＝V_t(c，j)×V_t(c，j)；

最后，計(jì)算，V_t(j：n，j)＝V_t(j：n，j)/b，具體計(jì)算步驟：變量c從j遞增到n計(jì)算：V_t(c，j)＝V_t(c，j)/b；

(6)更新矩陣A_t的第j列：A_t(j，j)＝a，A_t(j+1：n，j)＝0‘’

(7)t＝t+128，返回(3)。

在CPU和GPU混合計(jì)算平臺(tái)上分別對(duì)四個(gè)不同電網(wǎng)中的稀疏線性方程組集合的雅可比矩陣進(jìn)行了批量QR分解，具體計(jì)算時(shí)間，如圖2所示。算例4中，批處理數(shù)量N從1到變化到100時(shí)，計(jì)算時(shí)間只增加了30％(從383.5ms到500.1ms)，而當(dāng)N＝400時(shí)的平均單個(gè)雅可比矩陣的計(jì)算時(shí)間僅為1.47ms，遠(yuǎn)快于KLU和UMPACK的計(jì)算速度。該算例的性能分析，如圖3所示，全局內(nèi)存請(qǐng)求由于QR求解設(shè)計(jì)中較好的合并訪問模式僅僅隨批處理數(shù)量增長(zhǎng)而緩慢增長(zhǎng)，當(dāng)N從1變化到25時(shí)，全局內(nèi)存請(qǐng)求僅僅增長(zhǎng)了10％；當(dāng)N變化到400時(shí)，全局內(nèi)存請(qǐng)求也只增長(zhǎng)了67％，計(jì)算時(shí)間也只相應(yīng)增長(zhǎng)了53％，設(shè)備內(nèi)存帶寬達(dá)到了49.6GB/s，相較與批處理數(shù)量N為1時(shí)發(fā)生了質(zhì)變，批處理性能的獲取即源自于內(nèi)存帶寬提升。而此時(shí)的內(nèi)存帶寬和計(jì)算帶寬都只達(dá)到了峰值的20％，這也說明了K20C有足夠的能力完成此規(guī)模下的批量潮流計(jì)算。