一種gpu加速的電力潮流雅可比矩陣的qr分解方法
【專利摘要】本文公開了一種GPU加速的電力潮流雅可比矩陣的QR分解方法，所述方法包括：CPU中對雅可比矩陣J進行QR符號分解，得到Household變換矩陣V和上三角陣R陣的稀疏結(jié)構(gòu)；根據(jù)R陣的稀疏結(jié)構(gòu)，對矩陣J各列進行并行化分層；GPU中按層次遞增的順序計算分層QR分解內(nèi)核函數(shù)SparseQR。本發(fā)明利用CPU控制程序的流程并處理基礎(chǔ)數(shù)據(jù)和GPU處理密集的浮點運算相結(jié)合的模式提高了電力潮流雅可比矩陣QR分解的效率，解決了電力系統(tǒng)靜態(tài)安全性分析中潮流計算耗時大的問題。
【專利說明】
一種GPU加速的電力潮流雅可比矩陣的QR分解方法
技術(shù)領(lǐng)域
[0001]本發(fā)明屬于電力系統(tǒng)高性能計算應(yīng)用領(lǐng)域，尤其涉及一種GPU加速的電力潮流雅可比矩陣的QR分解的GPU線程設(shè)計方法。
【背景技術(shù)】
[0002]潮流計算是電力系統(tǒng)中應(yīng)用最廣泛、最基本和最重要的一種電氣運算。在電力系統(tǒng)運行方式和規(guī)劃方案的研究中，都需要進行潮流計算以比較運行方式或規(guī)劃供電方案的可行性、可靠性和經(jīng)濟性。同時，為了實時監(jiān)控電力系統(tǒng)的運行狀態(tài)，也需要進行大量而快速的潮流計算。因此，在系統(tǒng)規(guī)劃設(shè)計和安排系統(tǒng)的運行方式時，采用離線潮流計算;在電力系統(tǒng)運行狀態(tài)的實時監(jiān)控中，則采用在線潮流計算。
[0003]而實際生產(chǎn)過程中，無論離線潮流和在線潮流計算都對潮流的計算速度有這比較高的要求。在涉及規(guī)劃設(shè)計和安排運行方式的離線潮流中，因設(shè)備落地方案等情況復(fù)雜，需要仿真運行的種類多，潮流計算量大，單個潮流計算時間影響整體仿真時長;而在電力系統(tǒng)運行中進行的在線潮流計算對計算時間敏感度高，需要實時給出潮流計算結(jié)果，如在預(yù)想事故、設(shè)備退出運行對靜態(tài)安全的影響的潮流計算中，系統(tǒng)需要計算大量預(yù)想事故下潮流分布，并實時地做出預(yù)想的運行方式調(diào)整方案。
[0004]傳統(tǒng)的牛頓拉夫遜法潮流計算中，修正方程組求解占潮流計算時間的70%，修正方程組求解的計算速度影響程序的整體性能。而隨著CPU計算速度提升的放緩，現(xiàn)階段的單個潮流計算計算時間已經(jīng)達(dá)到一個瓶頸。目前對潮流計算的加速方法主要集中在使用集群和多核服務(wù)器對多潮流進行粗粒度加速，實際成產(chǎn)中對單個潮流內(nèi)部運算加速的研究較少。
[0005]GPU是一種眾核并行處理器，在處理單元的數(shù)量上要遠(yuǎn)遠(yuǎn)超過CPU。傳統(tǒng)上的GPU只負(fù)責(zé)圖形渲染，而大部分的處理都交給了 CPU?，F(xiàn)在的GPU已經(jīng)法陣為一種多核，多線程，具有強大計算能力和極高存儲器帶寬，可編程的處理器。在通用計算模型下，GPU作為CPU的協(xié)處理器工作，通過任務(wù)合理分配分解完成高性能計算。
[0006]稀疏線性方程組求解計算具有并行性。對方程組系數(shù)矩陣進行QR符號分解后，得到Household變換矩陣V和上三角陣R陣的稀疏結(jié)構(gòu)，根據(jù)R陣的稀疏結(jié)構(gòu)，對矩陣J各列進行并行化分層。其中每層中的列的計算相互獨立，沒有依賴關(guān)系，天然可以被并行的計算處理，適合GPU加速。因此通過CPU和GPU之間合理的調(diào)度可以快速完成方程組系數(shù)矩陣進行QR分解，并求解稀疏線性方程組，國內(nèi)外學(xué)者已經(jīng)開始對GPU進行稀疏線性方程組加速求解的方法進行了研究，但是沒有深入的優(yōu)化線程設(shè)計，單純從計算量的分配上研究計算線程設(shè)計，對線程計算方式，數(shù)據(jù)索引方式?jīng)]有進行深入研究，無法使程序充分發(fā)揮GPU的優(yōu)勢。
[0007]因此，亟待解決上述問題。

【發(fā)明內(nèi)容】

[0008]發(fā)明目的:針對現(xiàn)有技術(shù)的不足，本發(fā)明提供一種能大幅減少電力潮流雅可比矩陣QR分解計算時間并能提升潮流計算速度的一種GHJ加速的電力潮流雅可比矩陣的QR分解方法。
[0009]潮流計算:電力學(xué)名詞，指在給定電力系統(tǒng)網(wǎng)絡(luò)拓?fù)?、元件參?shù)和發(fā)電、負(fù)荷參量條件下，計算有功功率、無功功率及電壓在電力網(wǎng)中的分布。
[0010]并行計算:相對于串行運算，是一種一次可執(zhí)行多個指令的算法，目的是提高計算速度，及通過擴大問題求解規(guī)模，解決大型而復(fù)雜的計算問題。
[0011 ] GPU:圖形處理器(英語:6瓜口11；^8?1'00688；[呢1]11；[1:，縮寫:GPU)。
[0012]本發(fā)明公開了一種GPU加速的電力潮流雅可比矩陣的QR分解方法，所述方法包括:
[0013](I)CPU中對雅可比矩陣J進行QR符號分解，得到Household變換矩陣V和上三角陣R陣的稀疏結(jié)構(gòu);根據(jù)R陣的稀疏結(jié)構(gòu)，對矩陣J各列進行并行化分層；
[0014](2)GPU中按層次遞增的順序計算分層QR分解內(nèi)核函數(shù)SparseQR。
[0015]其中，所述步驟(I)中，并行化分層將矩陣J的η列歸并到MaxLevel層中，屬于同一層中的列并行進行QR分解;每層包含的列的數(shù)量為LeVelnum(k)，k表示層號;存儲第k層中所有列號至映射表Mapk。
[0016]優(yōu)選的，所述步驟(2)中，分層(^分解內(nèi)核函數(shù)定義為3口3186(^〈他1。。1?，隊111^如>，其線程塊大小Nthread固定為128，當(dāng)對k層進行計算時，線程塊數(shù)量Nb1cks = Levelnum(k)，總線程數(shù)量為:Nb1cks XNthreads;按照層次遞增的順序，調(diào)用內(nèi)核函數(shù)SparseQR〈Levelnum(k)，Nthreads〉分解屬于第k層的所有列；SparseQR〈Levelnum(k)，Nthreads>的計算流程為:
[0017]((2.1)⑶DA自動為每個線程分配線程塊索引blockID和線程塊中的線程索引threadID；
[0018](2.2)將blockID和threadID賦值給變量bid和t，之后通過bid和t來索引bid號線程塊中的t號線程；
[0019](2.3)第1^(1號線程塊負(fù)責(zé)01?分解雅可比矩陣1的第」=1^?1{(1^(1)列；
[0020](2.4)第1^(1號線程塊中，變量1從1遞增到」-1，如果1?(1，」)#0，執(zhí)行以下計算:
[0021]首先，計算中間變量0= 2￥(:[:11，；01'.J(1:n，j)，其中V(1:n，i)是Household變換矩陣V中第i列的第i?η行元素構(gòu)成的向量，J(1:n，j)是雅克比矩陣J中第j列的第i?η行元素構(gòu)成的向量;然后，采用公式氕1:11，」)=氕1:11，」)-|^(1:11，1)更新雅可比矩陣1的第」列；
[0022](2.5)計算第j列的Household變換向量；
[0023](2.6)更新了矩陣第」列:了(」，」)=&，了(」+1:11，」)=0。
[0024]優(yōu)選的，所述步驟(2.4)中:
[0025]首先，計算中間變量0= 2￥(:[:11，；01'.J(1:n，j)，具體計算步驟為，
[0026]I)判斷線程編號t是否小于n-1，否則該線程結(jié)束執(zhí)行；
[0027]2)線程中的變量Temp+ = 2V( i+t，i) X J( i+t，j);
[0028]3)t = t+128，返回 I)執(zhí)行；
[0029]4)線程中Temp的值賦給共享內(nèi)存變量Cache[threadID]；
[0030]5)對線程塊中的線程進行同步；
[0031]6)變量 M= 128/2;
[0032]7)當(dāng)M不等于O時執(zhí)行8);否則結(jié)束計算，跳轉(zhuǎn)到12);
[0033]8)判斷 thread ID 小于M執(zhí)行:Cache [threadID] + = Cache [threadID+M]，否則線程不執(zhí)行；
[0034]9)對線程塊中的線程進行同步；
[0035]10)M/ = 2；
[0036]11)返回7)執(zhí)行；
[0037]12)P = Cache(0)；
[0038]然后，采用公式J(1: η，j) = J( 1: η，j)-βν(1: η，i)更新雅可比矩陣J的第j列，具體步驟為，
[0039]I)判斷線程編號t是否小于n-1，否則線程結(jié)束執(zhí)行；
[0040]2)J(t+i,j)=J(t+i,j)-PV(t+i,i);
[0041]3)t = t+128，返回 I)。
[0042]進一步，所述步驟(2.5)中，
[0043]首先，計算變量a2= J( j:n，j)T.J( j:n，j)，具體計算步驟如下:
[0044]I)判斷線程編號t是否小于n-j，否則該線程結(jié)束執(zhí)行；
[0045]2)線程中的變量Temp，Temp+ = J( j+t, j) X J( j+t, j);;
[0046]3)t = t+128，返回 I)執(zhí)行；
[0047]4)線程中Temp的值賦給共享內(nèi)存變量Cache[threadID]；；
[0048]5)對線程塊中的線程進行同步；
[0049]6)變量 M= 128/2;
[0050]7)當(dāng)M不等于O時執(zhí)行8);否則結(jié)束計算，跳轉(zhuǎn)到12);
[0051 ]8)判斷 thread ID 小于M執(zhí)行:Cache [threadID] + = Cache [threadID+M]，否則線程不執(zhí)行；
[0052 ]9)對線程塊中的線程進行同步；
[0053]10)M/ = 2；
[0054]11)返回7)執(zhí)行；
[0055]12)a2 = Cache(0)；
[0056]其次，計算，v(j:η，j) = J( j:η，j )-aej(j:n)，其中是ej是第j個元素為I的單位向量；
[0057]具體步驟如下:
[0058]I)判斷線程編號t是否小于n-j，否則線程結(jié)束執(zhí)行；
[0059]2)V(t+j, j)=J(t+j, j)-aej(t+j)；；
[0060]3)t = t+128，返回 I)。
[0061 ]然后，計算變量b2 = V( j:η，j )τ.V( j:η，j);
[0062]具體計算步驟如下:
[0063]I)判斷線程編號t是否小于n-j，否則該線程結(jié)束執(zhí)行；
[0064]2)線程中的變量為Temp，Temp+ = V( j+t，j) XV( j+t, j);
[0065]3)t = t+128，返回 I)執(zhí)行；
[0066]4)線程中Temp的值賦給共享內(nèi)存變量Cache[threadID]；
[0067]5)對線程塊中的線程進行同步；
[0068]6)變量 M= 128/2;
[0069]7)當(dāng)M不等于O時執(zhí)行8);否則結(jié)束計算，跳轉(zhuǎn)到12);
[0070]8)判斷 thread ID 小于M執(zhí)行:Cache [threadID] + = Cache [threadID+M]，否則線程不執(zhí)行；
[0071]9)對線程塊中的線程進行同步；
[0072]10)M/ = 2；
[0073]11)返回7)執(zhí)行；
[0074]12)b2 = Cache(0)；
[0075]最后，計算V (j:n, j)=V(j:n, j)/b；
[0076]具體步驟如下:
[0077]I)判斷線程編號t是否小于n-j，否則線程結(jié)束執(zhí)行；
[0078]2)V(t+j，j)=V(t+j，j)/b;
[0079]3)t = t+128，返回 I)。
[0080]有益效果:與現(xiàn)有技術(shù)比，本發(fā)明的有益效果為:首先本發(fā)明采用CPU對電力潮流的雅可比矩陣J進行QR符號分解，根據(jù)R陣的稀疏格式，可以減少不必要的浮點計算;其次根據(jù)R陣的稀疏結(jié)構(gòu)將J陣的各列分到可以并行計算的不同層次，并將分層結(jié)果傳給GPU;再者在GPU中按層次遞增的順序計算分層QR分解內(nèi)核函數(shù)SparseQR，且GPU線程中采用了歸約算法處理點積運算，提高了計算效率;最后本發(fā)明利用CPU控制程序的流程并處理基礎(chǔ)數(shù)據(jù)和GPU處理密集的浮點運算相結(jié)合的模式提高了電力潮流雅可比矩陣QR分解的效率，解決了電力系統(tǒng)靜態(tài)安全性分析中潮流計算耗時大的問題。
【附圖說明】
:
[0081 ]圖1為本發(fā)明的實例計算時間；
[0082]圖2為本發(fā)明的實例分層結(jié)果；
[0083]圖3為本發(fā)明的流程示意圖。
【具體實施方式】
:
[0084]如圖3所示，本發(fā)明一種GPU加速的電力潮流雅可比矩陣的QR分解方法，所述方法包括:
[0085](I)CPU中對雅可比矩陣J進行QR符號分解，得到Household變換矩陣V和上三角矩陣R矩陣的稀疏結(jié)構(gòu)，符號分解之后的J的稀疏結(jié)構(gòu)等于V+R;根據(jù)R陣的稀疏結(jié)構(gòu)，對矩陣J各列進行并行化分層。
[0086](2)GPU中按層次遞增的順序計算分層QR分解內(nèi)核函數(shù)SparseQR。
[0087]一、CPU中對電力潮流雅可比矩陣J進行QR符號分解方法
[0088]首先，在CPU中對雅可比矩陣J進行QR符號分解，得到Household變換矩陣V和上三角陣R陣的稀疏結(jié)構(gòu)，符號分解之后的J的稀疏結(jié)構(gòu)等于V+R;然后，并行化分層將矩陣J的η列歸并到MaxLevel層中，屬于同一層中的列并行進行QR分解；每層包含的列的數(shù)量為LeVelnum(k)，k表示層號;存儲第k層中所有列號至映射表Mapk。最后，CPU將GPU計算所需數(shù)據(jù)傳輸給GPU，數(shù)據(jù)包括:雅可比矩陣J，其維度η，上三角陣R陣，層數(shù)MaxLeveI，每層包含的列數(shù)Levelnum以及映射表Map。
[0089]其中QR符號分解原理和并行化分層原理參見
[0090]aDirectMethodsforSparseLinearSystemsTimothyA.Davi s, SI AM,Philadelphia，2006。本專利使用的QR符號分解和并行化分層程序參見CSparse:
[0091 ] aConciseSparseMatrixpackage.VERS10N3.1.4,Copyright(c)2006-2014,T imothyA.Davis，OctlO，2014 Aouseho Id變換原理參見文獻:胡冰新，李寧，呂俊.一種采用Household變換遞歸實現(xiàn)的復(fù)矩陣QR分解算法[J].系統(tǒng)仿真學(xué)報，2004，16(11):2432-2434。
[0092]二、GPU中按層次遞增的順序啟動分層QR分解內(nèi)核函數(shù)SparseQR
[0093]分層QR分解內(nèi)核函數(shù)定義為SparseQIKNblocks，Nthreads>，其線程塊大小Nthread固定為128，當(dāng)對k層進行計算時，線程塊數(shù)量Nb1cks = Levelnum(k)，總線程數(shù)量為:Nb1cks XNthreads;按照層次遞增的順序，調(diào)用內(nèi)核函數(shù)SparseQR〈Levelnum(k)，Nthreads〉來分解屬于第k層的所有列。
[0094]SparseQR〈Levelnum(k)，Nthreads〉的計算流程為:
[0095](I)CUDA自動為每個線程分配線程塊索弓IblockID和線程塊中的線程索引threadID；
[0096](2)將blockID和threadID賦值給變量bid和t，之后通過bid和t來索引bid號線程塊中的t號線程；
[0097](3)第bid號線程塊負(fù)責(zé)QR分解雅可比矩陣J的第j=Mapk(bid)列；
[0098](4)第bid號線程塊中，變量i從I遞增到j(luò)-Ι，如果R(i，j)#0，執(zhí)行以下計算:
[0099]首先，計算中間變量0= 2￥(:[:11，；01'.J(1:n，j)，其中V(1:n，i)是Household變換矩陣V中第i列的第i?η行元素構(gòu)成的向量，J(1:n，j)是雅克比矩陣J中第j列的第i?η行元素構(gòu)成的向量;然后，采用公式氕1:11，」)=氕1:11，」)-|^(1:11，1)更新雅可比矩陣1的第」列；
[0100]首先，計算中間變量0= 2￥(:[:11，；01'.J(1:n，j)，具體計算步驟為，
[0101]I)判斷線程編號t是否小于n-1，否則該線程結(jié)束執(zhí)行；
[0102]2)線程中的變量Temp+ = 2V( i+t，i) X J( i+t，j);
[0103]3)t = t+128，返回 I)執(zhí)行；
[0104]4)線程中Temp的值賦給共享內(nèi)存變量Cache[threadID]；
[0105]5)對線程塊中的線程進行同步；
[0106]6)變量 M=128/2;
[0107]7)當(dāng)M不等于O時執(zhí)行8);否則結(jié)束計算，跳轉(zhuǎn)到12);
[0108]8)判斷 thread ID 小于M執(zhí)行:Cache [threadID] + = Cache [threadID+M]，否則線程不執(zhí)行；
[0109]9)對線程塊中的線程進行同步；
[0110]10)M/ = 2；
[0111]11)返回7)執(zhí)行;
[0112]12)P = Cache(0)；
[0113]然后，采用公式J(1: η，j) = J( 1: η，j)-βν(1: η，i)更新雅可比矩陣J的第j列，具體步驟為，
[0114]I)判斷線程編號t是否小于n-1，否則線程結(jié)束執(zhí)行；
[0115]2)J(t+i，j)=J(t+i，j)_PV(t+i，i);
[0116]3)t = t+128，返回 I)。
[0117](5)計算第j列的Household變換向量:
[0118]首先，計算變量a2= J( j:n，j)T.J( j:n，j)，具體計算步驟如下:
[0119]I)判斷線程編號t是否小于n-j，否則該線程結(jié)束執(zhí)行；
[0120]2)線程中的變量Temp，Temp+ = J( j+t, j) X J( j+t, j);;
[0121]3)t = t+128，返回 I)執(zhí)行；
[0122]4)線程中Temp的值賦給共享內(nèi)存變量Cache[threadID];;
[0123]5)對線程塊中的線程進行同步；
[0124]6)變量 M= 128/2;
[0125]7)當(dāng)M不等于O時執(zhí)行8);否則結(jié)束計算，跳轉(zhuǎn)到12);
[0126]8)判斷 thread ID 小于M執(zhí)行:Cache [threadID] + = Cache [threadID+M]，否則線程不執(zhí)行；
[0127]9)對線程塊中的線程進行同步；
[0128]10)M/ = 2；
[0129]11)返回7)執(zhí)行；
[0130]12)a2 = Cache(0)；
[0131 ]其次，計算，v( j:η，j) = J( j:η，j )-aej(j:n)，其中是ej是第j個元素為I的單位向量；
[0132]具體步驟如下:
[0133]I)判斷線程編號t是否小于n-j，否則線程結(jié)束執(zhí)行；
[0134]2)V(t+j, j)=J(t+j, j)-aej(t+j)；；
[0135]3)t = t+128，返回I)。
[0136]然后，計算變量b2= V( j:η，j )τ.V( j:η，j);
[0137]具體計算步驟如下:
[0138]I)判斷線程編號t是否小于n-j，否則該線程結(jié)束執(zhí)行；
[0139]2)線程中的變量為Temp，Temp+ = V( j+t，j) XV( j+t, j);
[0140]3)t = t+128，返回 I)執(zhí)行；
[0141 ]4)線程中Temp的值賦給共享內(nèi)存變量Cache[threadID]；
[0142]5)對線程塊中的線程進行同步；
[0143]6)變量 M=128/2;
[0144]7)當(dāng)M不等于O時執(zhí)行8);否則結(jié)束計算，跳轉(zhuǎn)到12);
[0145]8)判斷 thread ID 小于M執(zhí)行:Cache [threadID] + = Cache [threadID+M]，否則線程不執(zhí)行；
[0146]9)對線程塊中的線程進行同步；
[0147]10)M/ = 2；
[0148]11)返回7)執(zhí)行；
[0149]12)b2 = Cache(0)；
[0150]最后，計算V (j:n, j)=V(j:n, j)/b；
[0151]具體步驟如下:
[0152]I)判斷線程編號t是否小于n-j，否則線程結(jié)束執(zhí)行；
[0153]2)V(t+j, j)=V(t+j, j)/b；
[0154]3)t = t+128，返回 I)。
[0155](6)更新J 陣第 j列:J( j，j)=a，J( j+1:n，j)=0。
[0156]本發(fā)明所使用的GPU計算平臺配備一張TeslaK20CGPU卡和IntelXeonE5-2620CPU，GPU的峰值帶寬可達(dá)208GB/s，單精度浮點計算量峰值可達(dá)3.52Tf lops，CPU主頻為2GHz XPU計算平臺配備IntelCorei7-3520M2.90GHz的CPU。在CPU和GPU計算平臺上分別對16557維稀疏雅可比矩陣實例進行了測試，圖1為16557維稀疏雅可比矩陣在不同計算庫上的計算時間。在GPU上使用兩級并行策略的QR算法求解時間為67ms，落后于KLU的數(shù)值分解時間，略快與UMPACK的計算時間。GPU加速算法落后于KLU計算時間的原因主要是每層包含的列數(shù)隨著層數(shù)增加迅速減少，第I層包含1642列，前30層每層大于100列，但到了 207層之后僅剩I?2列。前面層次中可并行計算的列數(shù)量多，后面層次中的列需要串行執(zhí)行，無法充分發(fā)揮GPU的計算性能。圖2為16557維稀疏雅可比矩陣的并行度變化曲線(消去樹每層包含的列數(shù)，縱坐標(biāo)以10為底的對數(shù)表示)。
【主權(quán)項】
1.一種GRJ加速的電力潮流雅可比矩陣的QR分解方法，其特征在于:所述方法包括: (I)CPU中對雅可比矩陣J進行QR符號分解，得到Household變換矩陣V和上三角陣R陣的稀疏結(jié)構(gòu);根據(jù)R陣的稀疏結(jié)構(gòu)，對矩陣J各列進行并行化分層； (2 )GPU中按層次遞增的順序計算分層QR分解內(nèi)核函數(shù)Spar seQR。2.根據(jù)權(quán)利要求1所述的GPU加速的電力潮流雅可比矩陣的QR分解方法，其特征在于:所述步驟(I)中，并行化分層將矩陣J的η列歸并到MaxLevel層中，屬于同一層中的列并行進行QR分解;每層包含的列的數(shù)量為Levelnum(k)，k表不層號;存儲第k層中所有列號至映射表Mapk ο3.根據(jù)權(quán)利要求1所述的GPU加速的電力潮流雅可比矩陣的QR分解方法，其特征在于:所述步驟(2)中，分層QR分解內(nèi)核函數(shù)定義為SparseQtKNbIcicks，Nthreads>，其線程塊大小^1^固定為128，當(dāng)對1^層進行計算時，線程塊數(shù)量他1。。1? = 1^611111111(10，總線程數(shù)量為:Nb1cks XNthreads;按照層次遞增的順序，調(diào)用內(nèi)核函數(shù)SparseQR〈Levelnum(k)，Nthreads>分解屬于第k層的所有列；SparseQR〈Levelnum(k)，Nthreads>的計算流程為: (2.1)0]0六自動為每個線程分配線程塊索引1310(^10和線程塊中的線程索引訪^&(110; (2.2)將blockID和threadID賦值給變量bid和t，之后通過bid和t來索引bid號線程塊中的t號線程； (2.3)第1^(1號線程塊負(fù)責(zé)01?分解雅可比矩陣1的第」=1^?1{(1^(1)列； (2.4)第1^(1號線程塊中，變量1從1遞增到」-1，如果1?(1，」)#0，執(zhí)行以下計算: 首先，計算中間變量β = 2￥(:[:11，；01'*了(;[:11，」)，其中￥(;[:11，；0是Househo I d變換矩陣V中第i列的第i?η行元素構(gòu)成的向量，J(1:n，j)是雅克比矩陣J中第j列的第i?η行元素構(gòu)成的向量;然后，采用公式J(1:n，j)=J(1:n，j)-0V(1:n，i)更新雅可比矩陣J的第j列； (2.5)計算第j列的Household變換向量； (2.6)更新了矩陣第」列:了(」，」)=3，了(」+1:]1，」)=0。4.根據(jù)權(quán)利要求3所述的GPU加速的電力潮流雅可比矩陣的QR分解方法，其特征在于:所述步驟(2.4)中: 首先，計算中間變量0 = 2￥(:[:11，；01'.J(1:n，j)，具體計算步驟為， 1)判斷線程編號t是否小于n-1，否則該線程結(jié)束執(zhí)行； 2)線程中的變量1611^1+= 2￥(1+1:，；0 X J(i+t，j); 3)t= t+128，返回 I)執(zhí)行； 4)線程中Temp的值賦給共享內(nèi)存變量Cache[threadID]； 5)對線程塊中的線程進行同步； 6)變量M= 128/2; 7)當(dāng)M不等于O時執(zhí)行8);否則結(jié)束計算，跳轉(zhuǎn)到12); 8)判斷thread ID 小于 M 執(zhí)行:Cache [threadID] + = Cache [threadID+M]，否則線程不執(zhí)行； 9)對線程塊中的線程進行同步；10)M/= 2； 11)返回7)執(zhí)行；12)P= Cache(0)；然后，采用公式J(1: η，j) = J(1: η，j)-βν(1: η，i)更新雅可比矩陣J的第j列，具體步驟為， 1)判斷線程編號t是否小于n-1，否則線程結(jié)束執(zhí)行； 2)J(t+i，j)=J(t+i，j)-0V(t+i，i); 3)t= t+128，返回 I)。5.根據(jù)權(quán)利要求3所述的GPU加速的電力潮流雅可比矩陣的QR分解方法，其特征在于:所述步驟(2.5)中， 首先，計算變量a2 = J( j:n，j)T.J( j:n，j)，具體計算步驟如下: 1)判斷線程編號t是否小于n-j，否則該線程結(jié)束執(zhí)行； 2)線程中的變量Temp，Temp+= J( j+t, j) X J( j+t, j)；； 3)t= t+128，返回 I)執(zhí)行； 4)線程中Temp的值賦給共享內(nèi)存變量Cache[threadID]；； 5)對線程塊中的線程進行同步； 6)變量M= 128/2; 7)當(dāng)M不等于O時執(zhí)行8);否則結(jié)束計算，跳轉(zhuǎn)到12);8)判斷thread ID 小于 M 執(zhí)行:Cache [threadID] + = Cache [threadID+M]，否則線程不執(zhí)行； 9)對線程塊中的線程進行同步；10)M/= 2； 11)返回7)執(zhí)行； 12)a2= Cache(0)； 其次，計算，V( j:n，j)=J( j:n，j)-aej( j:n)，其中是ej是第j個元素為I的單位向量； 具體步驟如下: 1)判斷線程編號t是否小于n-j，否則線程結(jié)束執(zhí)行； 2)V(t+j，j)=J(t+j，j)_aej(t+j);； 3)t= t+128，返回 I)。 然后，計算變量b2 = V( j:η，j )τ.V( j:η，j); 具體計算步驟如下: 1)判斷線程編號t是否小于n-j，否則該線程結(jié)束執(zhí)行； 2)線程中的變量為Temp，Temp+= V( j+t，j) XV( j+t, j); 3)t= t+128，返回 I)執(zhí)行； 4)線程中Temp的值賦給共享內(nèi)存變量Cache[threadID]； 5)對線程塊中的線程進行同步； 6)變量M= 128/2; 7)當(dāng)M不等于O時執(zhí)行8);否則結(jié)束計算，跳轉(zhuǎn)到12);8)判斷thread ID 小于 M 執(zhí)行:Cache [threadID] + = Cache [threadID+M]，否則線程不執(zhí)行； 9)對線程塊中的線程進行同步；10)M/= 2；11)返回7)執(zhí)行；12)b2= Cache(0)；最后，計算 V (j: η，j) = V( j: η，j )/b;具體步驟如下:1)判斷線程編號t是否小于n-j，否則線程結(jié)束執(zhí)行；2)V(t+j，j)=V(t+j，j)/b;3)t= t+128，返回 I)。
【文檔編號】H02J3/06GK106026107SQ201610592223
【公開日】2016年10月12日
【申請日】2016年7月26日
【發(fā)明人】周贛, 孫立成, 張旭, 柏瑞, 馮燕鈞, 秦成明, 傅萌
【申請人】東南大學(xué)