本發(fā)明屬于電力系統(tǒng)高性能計(jì)算應(yīng)用領(lǐng)域，尤其涉及一種GPU加速的批處理同構(gòu)稀疏矩陣乘滿(mǎn)向量的處理方法。

背景技術(shù)：

潮流計(jì)算是電力系統(tǒng)中應(yīng)用最廣泛、最基本和最重要的一種電氣運(yùn)算。在電力系統(tǒng)運(yùn)行方式和規(guī)劃方案的研究中，都需要進(jìn)行潮流計(jì)算以比較運(yùn)行方式或規(guī)劃供電方案的可行性、可靠性和經(jīng)濟(jì)性。同時(shí)，為了實(shí)時(shí)監(jiān)控電力系統(tǒng)的運(yùn)行狀態(tài)，也需要進(jìn)行大量而快速的潮流計(jì)算。因此，在系統(tǒng)規(guī)劃設(shè)計(jì)和安排系統(tǒng)的運(yùn)行方式時(shí)，采用離線(xiàn)潮流計(jì)算；在電力系統(tǒng)運(yùn)行狀態(tài)的實(shí)時(shí)監(jiān)控中，則采用在線(xiàn)潮流計(jì)算。

而實(shí)際生產(chǎn)過(guò)程中，無(wú)論離線(xiàn)潮流和在線(xiàn)潮流計(jì)算都對(duì)潮流的計(jì)算速度有這比較高的要求。在涉及規(guī)劃設(shè)計(jì)和安排運(yùn)行方式的離線(xiàn)潮流中，因設(shè)備落地方案等情況復(fù)雜，需要仿真運(yùn)行的種類(lèi)多，潮流計(jì)算量大，單個(gè)潮流計(jì)算時(shí)間影響整體仿真時(shí)長(zhǎng)；而在電力系統(tǒng)運(yùn)行中進(jìn)行的在線(xiàn)潮流計(jì)算對(duì)計(jì)算時(shí)間敏感度高，需要實(shí)時(shí)給出潮流計(jì)算結(jié)果，如在預(yù)想事故、設(shè)備退出運(yùn)行對(duì)靜態(tài)安全的影響的潮流計(jì)算中，系統(tǒng)需要計(jì)算大量預(yù)想事故下潮流分布，并實(shí)時(shí)地做出預(yù)想的運(yùn)行方式調(diào)整方案。

傳統(tǒng)的牛頓拉夫遜法潮流計(jì)算中，修正方程組求解占潮流計(jì)算時(shí)間的70％，修正方程組求解的計(jì)算速度影響程序的整體性能。而隨著CPU計(jì)算速度提升的放緩，現(xiàn)階段的單個(gè)潮流計(jì)算計(jì)算時(shí)間已經(jīng)達(dá)到一個(gè)瓶頸。目前對(duì)潮流計(jì)算的加速方法大多集中在單個(gè)潮流內(nèi)部運(yùn)算加速，實(shí)際生成中對(duì)批處理加速潮流計(jì)算的研究較少。

GPU是一種眾核并行處理器，在處理單元的數(shù)量上要遠(yuǎn)遠(yuǎn)超過(guò)CPU。傳統(tǒng)上的GPU只負(fù)責(zé)圖形渲染，而大部分的處理都交給了CPU。現(xiàn)在的GPU已經(jīng)法陣為一種多核，多線(xiàn)程，具有強(qiáng)大計(jì)算能力和極高存儲(chǔ)器帶寬，可編程的處理器。在通用計(jì)算模型下，GPU作為CPU的協(xié)處理器工作，通過(guò)任務(wù)合理分配分解完成高性能計(jì)算。

稀疏線(xiàn)性方程組求解計(jì)算具有并行性。使用迭代法求解大型稀疏線(xiàn)性方程已經(jīng)得到廣泛的關(guān)注與應(yīng)用。在迭代法求解稀疏線(xiàn)性方程過(guò)程中有大量的稀疏矩陣乘滿(mǎn)向量的計(jì)算需求。因?yàn)榫仃囍忻啃信c滿(mǎn)向量相乘的計(jì)算相互獨(dú)立，沒(méi)有依賴(lài)關(guān)系，天然可以被并行的計(jì)算處理，適合GPU加速。

N-1靜態(tài)安全分析需要計(jì)算大量的潮流，在使用迭代法求解潮流修正方程時(shí)稀疏矩陣乘滿(mǎn)向量是其中重要的模塊，耗時(shí)很多。但由于斷線(xiàn)最多只影響稀疏矩陣四個(gè)位置的元素，因此通過(guò)冗余設(shè)計(jì)，可將所有斷線(xiàn)故障的節(jié)點(diǎn)導(dǎo)納陣冗余設(shè)計(jì)成和基態(tài)電網(wǎng)稀疏結(jié)構(gòu)相同的一系列矩陣，然后再按相同的稀疏存儲(chǔ)格式合并存儲(chǔ)，以實(shí)現(xiàn)GPU讀取數(shù)據(jù)的合并訪問(wèn)。再通過(guò)CPU和GPU之間合理的調(diào)度可以批處理同構(gòu)稀疏矩陣乘滿(mǎn)向量。國(guó)內(nèi)外學(xué)者已經(jīng)開(kāi)始對(duì)GPU進(jìn)行稀疏線(xiàn)性方程組加速迭代求解的方法進(jìn)行了研究，但是沒(méi)有對(duì)其中重要的模塊稀疏矩陣乘滿(mǎn)向量的加速求解做過(guò)專(zhuān)門(mén)研究，也沒(méi)做過(guò)批處理同構(gòu)稀疏矩陣乘滿(mǎn)向量的工作，無(wú)法使程序充分發(fā)揮GPU的優(yōu)勢(shì)。

因此，亟待解決上述問(wèn)題。

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

發(fā)明目的：針對(duì)現(xiàn)有技術(shù)的不足，本發(fā)明提供了一種能提高了算法并行度和訪存效率，大幅減少N-1靜態(tài)安全分析中迭代法求解潮流修正方程時(shí)龐大稀疏矩陣乘滿(mǎn)向量計(jì)算需求耗時(shí)的GPU加速的批處理同構(gòu)稀疏矩陣乘滿(mǎn)向量的處理方法。

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

N-1靜態(tài)安全分析：對(duì)運(yùn)行中的網(wǎng)絡(luò)或某一研究態(tài)下的網(wǎng)絡(luò),按N-1原則,研究一個(gè)個(gè)運(yùn)行元件因故障退出運(yùn)行后,元件有無(wú)過(guò)負(fù)荷及母線(xiàn)電壓有無(wú)越限。

GPU：圖形處理器(英語(yǔ)：GraphicsProcessingUnit，縮寫(xiě)：GPU)。

技術(shù)方案：本發(fā)明所述的一種GPU加速的批處理同構(gòu)稀疏矩陣乘滿(mǎn)向量的處理方法，大量同構(gòu)稀疏矩陣A₁～A_bs的乘滿(mǎn)向量操作：A₁x₁＝b₁，…，A_bsx_bs＝b_bs，其中x₁～x_bs為被乘的滿(mǎn)向量，b₁～b_bs為結(jié)果滿(mǎn)向量，bs為批處理的矩陣數(shù)量，所述方法包括：

(1)在CPU中將所有矩陣A₁～A_bs存儲(chǔ)為行壓縮存儲(chǔ)格式，矩陣A₁～A_bs共享相同的行偏移數(shù)組CSR_Row和列號(hào)數(shù)組CSR_Col，行偏移數(shù)組元素CSR_Row[k]中存儲(chǔ)的是矩陣第k行之前的非零元總數(shù)，k取值范圍從1到n+1；每個(gè)矩陣具體的數(shù)值存儲(chǔ)于各自的數(shù)值數(shù)組CSR_Val₁～CSR_Val_bs中，被乘滿(mǎn)向量存儲(chǔ)于數(shù)組x₁～x_bs中，結(jié)果滿(mǎn)向量存儲(chǔ)于數(shù)組b₁～b_bs中，數(shù)組CSR_Val₁～CSR_Val_bs、數(shù)組x₁～x_bs和數(shù)組b₁～b_bs都按照矩陣編號(hào)遞增順序連續(xù)存儲(chǔ)；

(2)CPU將GPU內(nèi)核函數(shù)所需數(shù)據(jù)傳輸給GPU；

(3)將矩陣A₁～A_bs的乘滿(mǎn)向量任務(wù)分配給GPU線(xiàn)程，一個(gè)線(xiàn)程塊負(fù)責(zé)矩陣A₁～A_bs中特定一行的乘滿(mǎn)向量計(jì)算，并優(yōu)化內(nèi)存訪問(wèn)模式；

(4)GPU中執(zhí)行批處理同構(gòu)稀疏矩陣乘滿(mǎn)向量?jī)?nèi)核函數(shù)spmv_batch，批處理同構(gòu)稀疏矩陣乘滿(mǎn)向量?jī)?nèi)核函數(shù)定義為spmv_batch<N_blocks，N_threads>，其線(xiàn)程塊大小N_thread為bs，線(xiàn)程塊數(shù)量N_blocks＝n，總線(xiàn)程數(shù)量為：N_blocks×N_threads；調(diào)用內(nèi)核函數(shù)spmv_batch<N_blocks，N_threads>來(lái)批處理并行計(jì)算同構(gòu)稀疏矩陣乘滿(mǎn)向量。

其中，所述步驟(2)中，其所需數(shù)據(jù)包括：稀疏矩陣的維數(shù)n，批處理的矩陣數(shù)量bs，稀疏結(jié)構(gòu)數(shù)組CSR_Row和CSR_Col，稀疏矩陣數(shù)值數(shù)組CSR_Val₁～CSR_Val_bs，被乘滿(mǎn)向量x₁～x_bs，結(jié)果滿(mǎn)向量b₁～b_bs。

優(yōu)選的，所述步驟(3)中，所述數(shù)值數(shù)組CSR_Val₁～CSR_Val_bs、被乘滿(mǎn)向量x₁～x_bs、結(jié)果滿(mǎn)向量b₁～b_bs的存儲(chǔ)區(qū)域均為bs行的矩陣，并對(duì)三個(gè)矩陣進(jìn)行轉(zhuǎn)置操作，以實(shí)現(xiàn)合并訪問(wèn)。

再者，所述步驟(4)中，所述內(nèi)核函數(shù)spmv_batch<N_blocks，N_threads>的計(jì)算流程為：

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

(4.2)將blockID和threadID賦值給變量bid和t，之后通過(guò)bid和t來(lái)索引bid號(hào)線(xiàn)程塊中的t號(hào)線(xiàn)程；

(4.3)第bid號(hào)線(xiàn)程塊負(fù)責(zé)所有矩陣第bid行的乘滿(mǎn)向量操作；

(4.4)第bid號(hào)線(xiàn)程塊中，每個(gè)線(xiàn)程負(fù)責(zé)一個(gè)矩陣的第bid行乘滿(mǎn)向量，具體步驟如下：

1)讀取當(dāng)前行以及下一行的第一個(gè)非零元的索引CSR_Row[bid]和CSR_Row[bid+1]，定義變量j＝CSR_Row[bid]；

2)判斷j是否小于CSR_Row[bid+1]，否則線(xiàn)程結(jié)束執(zhí)行；

3)讀取當(dāng)前列CSR_Col[j]，并把它賦值給變量cur_col；

4)更新結(jié)果滿(mǎn)向量b_t的bid號(hào)元素值b_t[bid]+＝CSR_Val_t[j]×x_t[cur_col]；

5)j＝j(luò)+1，返回2)。

有益效果：與現(xiàn)有技術(shù)比，本發(fā)明具有以下顯著優(yōu)點(diǎn)：本發(fā)明利用CPU負(fù)責(zé)控制程序的整體流程和準(zhǔn)備數(shù)據(jù)，GPU負(fù)責(zé)計(jì)算密集的向量乘法，利用批處理模式提高了算法并行度和訪存效率，大幅降低了批量稀疏矩陣乘滿(mǎn)向量的計(jì)算時(shí)間，可以解決N-1靜態(tài)安全性分析中迭代法求解潮流修正方程時(shí)龐大的稀疏矩陣乘滿(mǎn)向量計(jì)算需求耗時(shí)多的問(wèn)題。

附圖說(shuō)明：

圖1為本發(fā)明的實(shí)例電網(wǎng)的數(shù)據(jù)表；

圖2為本發(fā)明的實(shí)例計(jì)算時(shí)間；

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

具體實(shí)施方式：

如圖3所示，本發(fā)明一種GPU加速的批處理同構(gòu)稀疏矩陣乘滿(mǎn)向量的處理方法，大量同構(gòu)稀疏矩陣A₁～A_bs的乘滿(mǎn)向量操作：A₁x₁＝b₁，…，A_bsx_bs＝b_bs，其中x₁～x_bs為被乘的滿(mǎn)向量，b₁～b_bs為結(jié)果滿(mǎn)向量，bs為批處理的矩陣數(shù)量，所述方法包括：

(1)在CPU中將所有矩陣A₁～A_bs存儲(chǔ)為行壓縮存儲(chǔ)格式，矩陣A₁～A_bs共享相同的行偏移數(shù)組CSR_Row和列號(hào)數(shù)組CSR_Col，行偏移數(shù)組元素CSR_Row[k]中存儲(chǔ)的是矩陣第k行之前的非零元總數(shù)，k取值范圍從1到n+1；每個(gè)矩陣具體的數(shù)值存儲(chǔ)于各自的數(shù)值數(shù)組CSR_Val₁～CSR_Val_bs中，被乘滿(mǎn)向量存儲(chǔ)于數(shù)組x₁～x_bs中，結(jié)果滿(mǎn)向量存儲(chǔ)于數(shù)組b₁～b_bs中，數(shù)組CSR_Val₁～CSR_Val_bs、數(shù)組x₁～x_bs和數(shù)組b₁～b_bs都按照矩陣編號(hào)遞增順序連續(xù)存儲(chǔ)；

(2)CPU將GPU內(nèi)核函數(shù)所需數(shù)據(jù)傳輸給GPU，所需數(shù)據(jù)包括：稀疏矩陣的維數(shù)n，批處理的矩陣數(shù)量bs，稀疏結(jié)構(gòu)數(shù)組CSR_Row和CSR_Col，稀疏矩陣數(shù)值數(shù)組CSR_Val₁～CSR_Val_bs，被乘滿(mǎn)向量x₁～x_bs，結(jié)果滿(mǎn)向量b₁～b_bs；

(3)將矩陣A₁～A_bs的乘滿(mǎn)向量任務(wù)分配給GPU線(xiàn)程，一個(gè)線(xiàn)程塊負(fù)責(zé)矩陣A₁～A_bs中特定一行的乘滿(mǎn)向量計(jì)算，并優(yōu)化內(nèi)存訪問(wèn)模式；其中，所述數(shù)值數(shù)組CSR_Val₁～CSR_Val_bs、被乘滿(mǎn)向量x₁～x_bs、結(jié)果滿(mǎn)向量b₁～b_bs的存儲(chǔ)區(qū)域均為bs行的矩陣，并對(duì)三個(gè)矩陣進(jìn)行轉(zhuǎn)置操作，以實(shí)現(xiàn)合并訪問(wèn)；

(4)GPU中執(zhí)行批處理同構(gòu)稀疏矩陣乘滿(mǎn)向量?jī)?nèi)核函數(shù)spmv_batch，批處理同構(gòu)稀疏矩陣乘滿(mǎn)向量?jī)?nèi)核函數(shù)定義為spmv_batch<N_blocks，N_threads>，其線(xiàn)程塊大小N_thread為bs，線(xiàn)程塊數(shù)量N_blocks＝n，總線(xiàn)程數(shù)量為：N_blocks×N_threads；調(diào)用內(nèi)核函數(shù)spmv_batch<N_blocks，N_threads>來(lái)批處理并行計(jì)算同構(gòu)稀疏矩陣乘滿(mǎn)向量。

本發(fā)明內(nèi)核函數(shù)spmv_batch<N_blocks，N_threads>的計(jì)算流程為：

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

(4.2)將blockID和threadID賦值給變量bid和t，之后通過(guò)bid和t來(lái)索引bid號(hào)線(xiàn)程塊中的t號(hào)線(xiàn)程；

(4.3)第bid號(hào)線(xiàn)程塊負(fù)責(zé)所有矩陣第bid行的乘滿(mǎn)向量操作；

(4.4)第bid號(hào)線(xiàn)程塊中，每個(gè)線(xiàn)程負(fù)責(zé)一個(gè)矩陣的第bid行乘滿(mǎn)向量，具體步驟如下：

1)讀取當(dāng)前行以及下一行的第一個(gè)非零元的索引CSR_Row[bid]和CSR_Row[bid+1]，定義變量j＝CSR_Row[bid]；

2)判斷j是否小于CSR_Row[bid+1]，否則線(xiàn)程結(jié)束執(zhí)行；

3)讀取當(dāng)前列CSR_Col[j]，并把它賦值給變量cur_col；

4)更新結(jié)果滿(mǎn)向量b_t的bid號(hào)元素值b_t[bid]+＝CSR_Val_t[j]×x_t[cur_col]；

5)j＝j(luò)+1，返回2)。

本發(fā)明所使用的GPU計(jì)算平臺(tái)配備一張TeslaK20CGPU卡和IntelXeonE5-2620CPU，GPU的峰值帶寬可達(dá)208GB/s，單精度浮點(diǎn)計(jì)算量峰值可達(dá)3.52Tflops，CPU主頻為2GHz。CPU計(jì)算平臺(tái)配備IntelCorei7-3520M2.90GHz的CPU。GPU計(jì)算平臺(tái)上對(duì)圖1中的四個(gè)實(shí)例電網(wǎng)的電納矩陣進(jìn)行了測(cè)試，圖2為對(duì)應(yīng)四個(gè)電網(wǎng)的不同bs下稀疏矩陣乘滿(mǎn)向量的測(cè)試時(shí)間。