第三章 基于GPU的大規模稀疏矩陣直接求解器
3.0 簡介
3.1 基于quotient graph的符号分析
3.1.1 頂點重排序
3.1.2 建構消去樹
3.1.3 尋找超結點
3.1.4 符号分解
3.2 多波前法
3.3 超節點方法
3.4 多波前+超節點方法的并行分解算法
小結
參考資料
第三章 基于GPU的稀疏直接求解器
前言
本章可能是所有章節中最難得了，如果節奏過快，對于即使有一定相關開發經驗的人來說要徹底了解其中的内容也會相當困難，是以作者會給出相比其它章節更多的說明。第一章中類似的技術将會在本章得到應用，通過閱讀本章，讀者不僅可以進一步鞏固第一章中學到的GPU計算的優化技術，還能了解到一些有趣的算法。在這個過程中讀者還會更深的體會到從算法的僞代碼到實際的可執行代碼之間的難度。
3.0 稀疏矩陣直接解法概述
很多科學計算問題中都涉及到求解線性方程組，而這些方程組不僅規模很大且通常是稀疏的。解決這類問題通常使用疊代解法或直接解法。疊代解法一般記憶體空間需求小，程式實作比較簡單；但是疊代解法的收斂速度依賴于矩陣的性态，如果矩陣過于病态，即使采用複雜的預處理技術(預處理越複雜，也意味着求解效率越低，記憶體空間需求越大，進而疊代解法優勢也會逐漸喪失)也很難收斂甚至不收斂。直接解法的缺點是記憶體空間的需求比疊代解法大很多，但是其具有解的精度高，計算時間穩定，沒有收斂性問題等優點；尤其是當矩陣不變，而具有多個右端項或右端項不斷變化時，效率比疊代法更高。但是直接法程式實作的難度非常的大，其中所用到的算法通常深度抽象，邏輯交錯龐雜，即使是作者本人，在事隔幾年後讀自己的代碼也廢了不少精力，是以希望讀者若對這一領域感興趣一定要有足夠的耐心和毅力；若工作中并不涉及這一領域則建議跳過本章以免浪費太多精力。稀疏直接求解方法繁多紛雜，但整體步驟大體一緻：
1 頂點重排
2 符号分解
3 數值分解
4 三角回代
5 疊代改善
有時也将第1步和第2步合并在一起統稱為符号分析，最後一步的疊代改善則是根據具體需求決定是否使用，不是必須的。本章以求解大規模對稱正定矩陣為例講解開發快速求解大規模稀疏矩陣的主要算法和相關優化技術。對稱正定矩陣的求解技術已經發展的非常成熟，通常使用cholesky三角分解方法，該方法數值穩定性好，且分解過程中無需選主元。在開始之前我們對一些符号進行約定：
adj(v) :頂點v的鄰接節點的集合。
deg(v) : 頂點v的度，也就是頂點為鄰接節點的個數。
reach(v) :頂點v的可達集，即通過頂點v可以到達的頂點的集合。
snode( v0, v1, v2, …) ：表示組成超節點的所有頂點編号的集合，簡寫為snode。
innz(v) :頂點v對應的波前中的非零元的行号或列号。
3.1 基于quotientgraph的符号分析
在進行真正的數值三角分解之前，我們需要事先确定存儲分解因子所需要的記憶體大小，以避免在分解的過程中動态記憶體配置設定。這一過程稱之為符号分解或符号消元，這不僅可以提高程式的運作效率，還可以節省很多記憶體并避免更多的記憶體碎片。假設矩陣A具有n個非零元，分解後的分解因子L具有m個非零元(是以有m-n個非零元填入，通常m-n>>n)，那麼按照定義進行符号分解需要O(m)的存儲空間。但是可以證明，隻需要O(n)大小的記憶體即可完成符号分解，且這一大小是可以事先确定的。實作這一方法的模型就是quotientgraph，我們後續符号分析中的頂點重排(reordering)，建構消去樹以及尋找超節點都是基于quotientgraph的思想在O(n)的複雜度下完成的。
3.1.1 頂點重排
我們知道，三角分解過程中會産生很多填入元，如果直接進行分解，其填入元數量一般非常多。當矩陣達到一定規模時，對存儲空間的需求會大幅度暴漲，求解就會非常耗時和困難甚至無法完成。是以我們需要事先對原矩陣的稀疏結構圖中的頂點進行重排序處理。重排序的過程并不改變原始圖的拓撲結構，因而具有等價關系，下面兩幅圖展示了重排序的作用，左邊是排序前的原始矩陣，右邊是排序後的矩陣：
3.png (8.78 KB, 下載下傳次數: 0)
2016-6-24 11:19 上傳
4.png (7.13 KB, 下載下傳次數: 0)
2016-6-24 11:20 上傳
fig3.1 fig3.2
可以看到，如果直接進行分解，矩陣的所有元素都程式設計了非零元素；而重排序後進行分解非零元素的數量不變，這雖然隻是最理想的情況，但實際上通過頂點的重排序仍然能大幅度減少分解過程中新産生的非零元素。稀疏矩陣頂點圖的重排序算法通常都是NP問題，進而隻能采用啟發式算法。目前最常用的算法有以縮減外形(或帶寬)為目的算法，的如RCM算法；以減少填入元為目的的類最小度算法(minimum degree)和用于分而治之的嵌套剖分(nested dissection)算法。雖然更窄的帶寬通常意味着更少的填入元，但是其效果一般不如類最小度算法。嵌套剖分目前比較成熟的技術是通過multilevel方法(類似多重網格方法中的思想)找到極小頂點分割子将整個圖剖分成多個獨立的子圖，然後再采用類最小度算法對各個子圖分别進行局部優化處理(fig3.3)。類最小度算法常用的有最小度算法(MD),多重最小度算法(MMD)和近似最小度算法(AMD)。多數情況下性能最高效果最好的是近似最小度算法，本書中我們僅采用最小度算法。、
5.png (6.9 KB, 下載下傳次數: 0)
2016-6-24 11:21 上傳
（白色部分為0元素）
最消度算法的思想是當一個頂點的度比較小時，當消去該頂點時引入新的邊的機率會更小或是引入的新邊的數量會更少，是以所有最小度算法都是從一個具有最小度的頂點開始。第一個具有最小度的頂點可以随機的選擇，也可以選擇所有具有相同最小度的頂點中具有最大僞直徑的頂點。為了簡單，我們周遊目前未消元的頂點并選擇第一個頂點度最小的頂點。
最小度算法 pesudo-code
do{
v=one vertex of graph with minimum degree
elimination:
delete all edges adjacency to ‘v’
add new edges between the adjacency verties of ‘v’(if non-existent edges betweenthem)
update graph: remove ‘v’ from graph
}while(graph.num_verties>0)
最小度算法的實作有很多種方法，這裡作者僅給出自己的實作（基于quotientgraph模型）。以下代碼都經過數百次針對不同布局和大小的矩陣的測試；但是按照嚴格的标準來說，數百次的測試遠遠不夠，是以作者不保證代碼中沒有BUG，請讀者慎用。
首先，需要選擇一個具有最小度的頂點,這裡我們稱之為主元頂點：
int search_pivot( const int* deg, const char* flag, int n )
{
int i, d, p, q;
for( i=0; i<n; ++i ){ if(flag==0 ) break; }
d=deg; p=i;
while( (++i)<n ){ if( flag==0 ){ q=deg; if( q<d){ d=q; p=i; } } }
return p;
}
然後需要得到這個頂點的可達集:
int gather_reaches( aux_param_t* params,const int64* ptr, const int* nbor, int e )
int64 p, q;
int i, id, nen, n=0;
p=ptr[e]; q=p+params->len[e];
params->flag[e]=1;
while( p<q ){
id=nbor[p++];
if( params->flag[id]==0 ){
if( params->mark[id]==0 ){
params->reach[n++]=id;params->mark[id]=1;
} else {
nen=gather_enbors(params, ptr, nbor, id );
for( i=0; i<nen; ++i ){
id=params->nbn;
if( params->mark[id]==0 ){ params->reach[n++]=id; }
params->mark[id]=2;
}
}
}
return n;
gather_reaches函數通過周遊主元頂點所有未消去的鄰接節點以及所有以消去的鄰接節點的鄰接節點來獲得主元頂點的可達集;gather_enbors用來擷取某個以消去頂點的鄰接頂點：
int gather_enbors( aux_param_t* params, const int64* ptr, const int* nbor, int e )
{
int k, id, v, n=0;
int64 p, q;
if( params->esize[e]>0 )
{
v=e;
for( k=0; k<params->esize[e]+1;++k ){
p=ptr[v];q=p+params->len[v];
while( p<q ){
id=nbor[p++];
if( params->flag[id]==0 ){ params->nbn[n++]=id; }
v=params->elink[v];
else
p=ptr[e];q=p+params->len[e];
while( p<q ){
id=nbor[p++];
if( params->flag[id]==0 ){ params->nbn[n++]=id; }
在消元過程中，可能有多個具有相同鄰接頂點的以消去頂點，為了避免重複計算，需要對頂點進行”吸收”,亦即将多個具有相同鄰接節點的以消去頂點合并到其中一個超頂點。頂點吸收可以減少備援計算，同時也是基于quotient-graph模型計算時所必須的，否則記憶體空間将無處容納頂點在消元過程中增加的邊:
void absorb_elements( int* nbor, const int64* ptr,aux_param_t* params, int nr, int e )
{
int i, id, enb, n, n_deads, n_slots, sign, v;
int64 ep, p=ptr[e], q=p+params->len[e];
enb=0;
for( enb=0; p<q; ++p ){
if( params->flag[( id=nbor[p] )]!=0 ){ params->nbn[enb++]=id; }
if( enb==0 ){
for( id=0; id<nr; ++id ){ params->mark[params->reach[id]]=0; } return;
for( p=0, i=0; i<enb; ++i ){
id=params->nbn;
params->spot[p++]=id;n=params->esize[id]; params->flag[id]=2;
for( q=0; q<n; ++q ){ params->spot[p++]=( id=params->elink[id] ); }
v=e;params->esize[e]=(int)p;
for( i=0; i<p; ++i ){ id=params->spot; params->elink[v]=id; v=id;}
i=n_slots=sign=0; v=e;
for(;;)
++n_slots;ep=ptr[v]; p=ep; q=ptr[v+1];
if( i>=nr ){ sign=1; break; }
nbor[p++]=params->reach[i++];
params->len[v]=(int)(p-ep);
if( sign!=0 ) break;
v=params->elink[v];
params->esize[e]=--n_slots;
for( params->flag[e]=2, v=0; v<nr; ++v )
{
id=params->reach[v];
if(params->mark[id]!=2 ){ params->mark[id]=0; continue; }
params->mark[id]=0;p=ptr[id]; q=p+params->len[id];
while( p<q ){ if( params->flag[nbor[p]]==2 ){nbor[p]=e; break; } ++p; }
ep=p; i=0;
while( (++p)<q ){
if( params->flag[nbor[p]]!=2 ){
params->spot[i++]=nbor[p];
n_deads=(int)(q-(++ep)-i);
if( n_deads>0 ){
for( p=0; p<i; ++p, ++ep ){ nbor[ep]=params->spot[p]; }
params->len[id]-=n_deads;
最後我們需要更新頂點的度，首先我們需要計算出目前消元頂點的可達集（是對可達集中的頂點的度進行更新），為此周遊消元頂點的鄰接節點，如果其鄰接節點是已消去的頂點，則周遊該頂點的鄰接頂點，如果其鄰接頂點未被标記，則加入可達集；如果消元頂點的鄰接頂點是未消去頂點，則将其加入可達集（如果未标記的話）。頂點v的可達集可表示為
reach(v)= adj(s)∪∑adj(e)
接着我們周遊可達集中的頂點，對于可達集中的每個頂點，周遊其鄰接頂點，如果其鄰接頂點為已消去頂點，則将該頂點的頂點度加1，否則不變。如果其鄰接頂點為非消去頂點，則将其頂點度加1,如此進行，直至可達集中的所有頂點處理完成，下面是具體的實作代碼：
voidupdate_degrees( aux_param_t* params, const int64* ptr, const int* nbor, int nr )
int k, d, s, n, i, id;
for( k=0; k<nr; ++k )
s=params->reach[k];
if( params->deg==0 ) continue;
d=0;p=ptr; q=p+params->len; params->mark=1;
do{
id=nbor[p];
if( params->flag[id]==0 ){
if( params->mark[id]==0 ){
params->spot[d++]=id;
params->mark[id]=1;
}else {
n=gather_enbors(params, ptr, nbor, id );
for( i=0; i<n; ++i ){
id=params->nbn;
if( params->mark[id]==0 ){
params->spot[d++]=id;params->mark[id]=1;
}
}while( (++p)<q );
params->deg=d;params->mark=0;
for( i=0; i<d; ++i ){ params->mark[params->spot]=0; }
現在可以組裝成完整的基于quotient-graph的最小度重排序函數了:
void mdo( int* iperm, char* aux, int64* Ap, int* Ai, int n )
aux_param_t params;
int64 start, end;
int d, i, e, m;
params.flag =(char*)aux;
params.mark =params.flag+n;
params.deg =(int*)( params.mark+n );
params.len =params.deg+n;
params.esize =params.len+n;
params.elink =params.esize+n;
params.reach =params.elink+n;
params.nbn =params.reach+n;
params.spot =params.nbn+n;
for( start=0, i=0; i<n; ++i )
params.flag=0;
params.mark=0;
params.esize=0;
end=Ap[i+1];
d=(int)(end-start);
params.deg=d;
params.len=d;
start=end;
for( i=0; i<n-1; ++i )
e=search_pivot(params.deg, params.flag, n );
iperm=e;
m=gather_reaches(&params, Ap, Ai, e );
if( m==0 ) continue;
absorb_elements(Ai, Ap, &params, m, e );
update_degrees(&params, Ap, Ai, m );
for( i=0; i<n; ++i ){ if( params.flag==0 ) break; }
iperm[n-1]=i;
aux是一個輔助記憶體數組，之是以沒有在内部配置設定，是因為在求解大規模矩陣時，通常會對多個子矩陣的頂點圖進行重排序，這時候可以配置設定與參與計算的線程數量對等數量的aux數組，每個線程中的aux數組的大小是該線程中所要計算的圖中尺寸最大的圖所對應的大小，這樣aux數組便可被各個線程中的多個頂點圖複用，而不用頻繁的開辟和釋放操作。通過重排序得到了新的頂點編号後，就可以建立消去樹了。需要說明的是，在我們的最小度算法中使用的是外度更新，而非标準的頂點度更新。由于考慮了額外更深層次頂點的影響，是以使用頂點外度的效果通常都要好于使用标準頂點度。
最小度算法有很多改進措施，第一個改進方法是将多個無差別頂點同時消去，這一方法也就是MMD算法。所謂無差別頂點就是當兩個不相鄰頂點的鄰接頂點集合與自身的集合相同的頂點，定義為
u∪adj(u) <==> v∪adj(v)
用可達集可表示為
reach(v, S)∪v <==> reach(u,S)∪u
第二個改進方法是進行不完全的度更新，也就是不需要每次都對完整的鄰接頂點集合的所有頂點進行更新，而隻更新滿足一定條件的具有最小度的頂點，這一方法稱為近似最小度(AMD)算法，AMD算法多數情況下的排序效果比MD和MMD算法好，同時算法的時間複雜度也大幅降低。
void build_etree( char* aux, int* parents, const int64* Cp, const int* Ci, int n )
int64 a, b, p, q, *Rp,*ptr;
int *Ri,*ancestor, i, k, next;
Rp=(int64*)aux;
ptr=Rp+n+1;
Ri=(int*)(ptr+n);
ancestors=Ri+Cp[n];
memset( Rp, 0, n*sizeof(int64) );
p=Cp[0];
for( k=0;k<n; ++k ){
q=Cp[k+1];
while( p<q ){ ++Rp[Ci[p++]]; }
}
a=Rp[0]; Rp[0]=0;
b=Rp[k+1]; Rp[k+1]=a;a+=b;
memcpy( ptr, Rp, n*sizeof(int64) );
for(k=0;k<n; ++k ){
while( p<q ){ Ri[ptr[Ci[p++]]++]=k; }
for( k=0;k<n; ++k )
{
parents[k]=-1;ancestors[k]=-1;
for( p=Rp[k]; p<Rp[k+1]; ++p )
i=Ri[p];
while((i>=0)&(i<k)){
next=ancestors;
ancestors=k;
if( next<0 ){ parents=k; }
i=next;
其中parents存儲了消元過程中每個頂點之間的依賴關系，這個資訊将在後面符号分解和尋找超節點的過程中被用到。
我們隊矩陣A的稀疏結構X進行重新編碼得到具有新的稀疏結構Y的矩陣B，Y與X中頂點之間的拓撲關系并未改變，但是卻具有更少的填入元。然後通過稀疏結構Y進行符号分解。符号分解的目的是擷取分解過程中分解因子中各個非零元素的索引以及整個分解因子所需要的存儲空間大小，進而避免在數值分解的過程中動态的計算索引以及動态的記憶體配置設定，這樣将具有更高的效率并避免記憶體碎片的積累。符号分解之是以有效是因為很多頂點對應的波前(後續章節會講到)可能進行多次更新，如果直接進行數值分解，就需要多次重複的計算索引以及很多低效的記憶體配置設定釋放操作，嚴重影響效率。符号分解的代碼同樣基于quotient-graph,前面的一些函數會得到複用。
首先我們需要先計算得到分解因子中每列非零元素的數量，并最終得到整個分解因子所需要的記憶體大小：
int64 symbolic_eliminate( char* aux, int* Ln, const int64* Yp, const int* Yi, int n )
int64nA=Hp[n]<<1;
int* Ai=(int*)(aux+ELIMINATION_BUFFER_BYTES(n));
int64*counter=(int64*)(Ai+nA);
int64*Ap=counter+n+1;
int64 p, q,c, nnz;
int k, m;
for( p=0,k=0; k<n; ++k ){
q=Yp[k+1];counter[k]=q-p; p=q;
for( p=0,k=1; k<=n; ++k ){
q=Yp[k];
while( p<q ){ ++counter[Yi[p++]]; }
for(Ap[0]=0, k=0; k<n; ++k ){
Ap[k+1]=Ap[k]+counter[k];counter[k]=Ap[k];
q=Yp[k+1];
while( p<q ){ Ai[counter[Yi[p++]]++]=k; }
c=counter[k];
while( p<q){ Ai[c++]=Yi[p++]; }
}
aux_param_t params;
params.flag =(char*)aux;
params.mark =params.flag+n;
params.len =(int*)( params.mark+n );
params.esize =params.len+n;
params.elink =params.esize+n;
params.reach =params.elink+n;
params.nbn =params.reach+n;
params.spot =params.nbn+n;
params.flag [k]=0;
params.mark [k]=0;
params.esize[k]=0;
q=Ap[k+1];
params.len[k]=(int)(q-p);
p=q;
nnz=0;
for( k=0;k<n; ++k ){
m=gather_reaches(&params, Ap, Ai, k );
if( m==0 ) continue;
Ln[k]=m; nnz+=m;
absorb_elements( Ai, Ap,&params, n, k );
return nnz;
symbolic_eliminate的傳回值是分解因子非零元的數量，而分解因子每列非零元的數量儲存在Ln中,然後我們可以在分解前一次性的配置設定索引數組的記憶體，并利用Ln擷取每列非零元素的行号以及存儲位置：
void symbolic_decompose( char* aux, int* Li, const int64* Lp, const int* Sp, const int* nodes, const int* superIDs, const int64* Ap, const int* Ai, int n, intn_snodes )
int64 ii, ii0, ii1, pp,qq, *pos;
char*mark;
int *spot,*a;
int i k,v, d, sid, c;
mark=aux;
pos=(int64*)(aux+n);
spot=(int*)(pos+n_snodes);
memset( mark, 0, n );
for( k=0;k<n_snodes; ++k )
v=nodes[Sp[k+1]-1];
ii0=Ap[v];
ii1=Ap[v+1];
pp=Lp[k];
while(ii0<ii1){
i=Ai[ii0++];
if(i>v){ Li[pp++]=Ai[ii0++]; }
pos[k]=pp;
ii0=Lp[k];
ii1=Lp[k+1];
d=(int)(ii1-ii0);
if((k>0)&(d>1)){
a=&Li[ii1];
sort( a,a+d );
for( ii=ii0; ii<ii1; ++ii )
v=Li[ii];
sid=superIDs[v];
if( sid==k ) continue;
pp=Lp[sid]+1;
qq=pos[sid];
c=0;
while( pp<qq ){
v=Li[pp++];mark[v]=1; spot[c++]=v;
for( pp=ii+1; pp<ii1; ++pp ){
v=Li[pp];
if( mark[v]==0 ){ Li[qq++]=v; }
for( i=0; i<c; ++i ){ mark[spot]=0; }
pos[sid]=qq;
這裡我們是基于超節點進行的符号分解，是以在執行symbolic_decompose之前還需要建構超節點，超節點方法的相關内容将在2.3節進行介紹。注意紅色部分的代碼，由于符号分解的過程中我們隻記錄了正在進行分解的列的行号，但是為了以正确的順序進行後續的分解，需要在此之前對其中的行号（或列号）從小到大排序。現在整個符号分析已經完成，将其過程總結如下：
1 對頂點進行重排序。
2 利用重排序得到的新頂點編碼建構消去樹，得到消元過程中結點依賴關系的數組parents。
3 利用新的頂點編碼和第2步中得到的parents資訊尋找并構造超節點，并根據超節點内頂點編碼的連續性把超節點标記為靜态超節點或動态超節點的。
4 利用新的頂點編碼和第3步中得到的超節點資訊進行符号分解确定整個分解過程所需要的記憶體大小和非零元素的索引資訊。
3.2 多波前法（multifrontalmethod）
多波前法發展于有限元中的波前法，其目的是提高稀疏矩陣分解過程中的并行性，是目前求解大規模及超大規模稀疏矩陣普遍使用的高效且成熟的方法。下面我們以易于了解的方式通過圖的展示來說明多波前法的思想（假設已經得到了分解後的稀疏結構）。
6.png (8.63 KB, 下載下傳次數: 0)
2016-6-24 11:25 上傳
7.png (8.41 KB, 下載下傳次數: 0)
2016-6-24 11:26 上傳
第一步：波前{0,1,2}可以同時進行分解，波前0的結果對{3，4,7}進行更新，
波前1的結果對波前{4，6}進行更新，波前2對波前{5，6}進行更新。
第二步：對波前{3}進行分解，然後更新波前4。
第三步：對波前{6，7}進行并行分解。
3.3 超結點法（supernodalmethod）
超結點法也是用來提高計算效率的一種方法，但不同于多波前法同時對多個獨立的波前進行并行計算， 而是通過提高計算密度和緩存友好來提升效率。假設一個連續的頂點序列{ vi,vi+1, vi+2, …, vi+n, … }，如果
E1 parent(vi)=vi+1, parent(vi+1)=vi+2,…； 且
E2 num_adj(vi)+1=num_adg(vi+1),num_adj(vi+1)+1=num_adj(vi+2), …
那麼将該頂點序列稱之為一個超結點，如圖所示：
8.png (415 Bytes, 下載下傳次數: 0)
2016-6-24 11:27 上傳
圖中有5個超節點（紅色cell代表分解過程中生成的填入元素，其它代表原始的稀疏結構），分别為{1}， {2}, {0，3, 6,7}，{4，5}。如果按照超節點進行計算，其順序為：
a. 首先對超節點{1},{2},{4,5}進行分解（由于沒有依賴關系，是以三個超節點可同時進行分解）
b. 使用分解得到的波前矩陣對3，6，7列進行更新
c. 對超節點{0,3,6,7}進行分解
對超節點進行分解可以利用高效的稠密矩陣計算提高計算效率。一般來說，超節點越少，則消去樹的高度越矮，計算密度越大，進而效率也越高，是以可以通過一些方法來減少超節點的數量。第一種方法是超節點融合，亦即把超節點當做普通的頂點，那麼又可以對超節點進一步融合成更大的超節點；但是超節點融合會占用更多的記憶體。第二種方法是在建立超節點時有限度的容忍，超節點中頂點的鄰接節點的數量無需嚴格滿足E2，隻要相差小于閥值即可加入超節點，仍以上圖為例，可以将超節點{1}和{0,3,6,7}合并為一個超節點：
{0，1，3，6，7}
這樣雖然浪費了一個元素的空間，但是卻可以減少超節點的數量（意味着消去樹層級更少，也是以減少了消去過程中的同步開銷，同時計算密度更大進而具有更高的效率）。
第三種方法是擴充超節點的概念，超節點中的頂點無需滿足連續性，但E1條件仍然要滿足，可将其稱之為動态超節點。現實中更有效的方法是：将超節點分為靜态超節點和動态超節點，同時配合第二種方法，這樣不僅可以避免不必要的記憶體移動(因為動态超節點在計算之前必須要将其移入一塊連續的記憶體中)，還可以最大限度的提高計算的密度，且記憶體的整體需求也一般會比較低。将超節點方法引入符号分解中還可以大幅降低存儲索引所需的記憶體并減少符号分解所需的計算量。
對于上圖來說，{1}, {2}, {4,5} 是靜态超節點，{0,3,6,7}是動态超節點。
容易看出，對于某個超節點snode{v0,v1,v2,…,vn},隻需要存儲其具有最小波前的頂點vn對應的那列的非零元的行号，而其它頂點對應的波前中非零元的行号可以表示為:
innz(v0) = { v1, v2, …,vn } + innz(vn)
innz(v1)= { v2, v3, …, vn } + innz(vn)
innz(v2)= { v3, v4, …, vn } + innz(vn)
……
進而整個超節點的非零元的行号（或列号）資訊可以通過超節點中的包含的頂點的标号以及最後一個頂點的非零元行号（或列号）完整的表示，亦即snode{ v0, v1, v2, …,vn } + innz(vn)，容易看出這樣将使索引存儲具有最小的空間需求。
建立超節點可以從正向尋找 (從第一列開始向後尋找)，也可以反向尋找(從最後一列開始向前尋找)，不過反向方法通常具有更好的效果，這裡分别給出兩種實作，在使用中可以從前向和後向方法中選擇結果最好的。
尋找超結點的反向方法：
int find_supernodals_backward( char* aux, int* Sp, int*verties, char* flags, const int64* Lp, const int* Li, int n )
int64 ii0, ii1, ii2, ii3;
int* child=(int*)aux;
int*id=child+n;
int i, k,n_snodes, snode_size, p, q, parent;
char b;
for( k=0;k<n; ++k ){ child[k]=-1; }
for( k=0;k<n-1; ++k )
parent=Li[ii0+1];
ii2=Lp[parent];
ii3=Lp[parent+1];
if(((ii1-ii0)-(ii3-ii2))==1){
child[parent]=k;
for( i=0,n_snodes=0, k=n-1; k>=0; --k )
if(child[k]!=-2)
snode_size=1;
id[i++]=k;
p=child[k];
child[k]=-2;
while(p>=0){
++snode_size;
id[i++]=p;
q=child[p];
child[p]=-2;
p=q;
Sp[n_snodes++]=snode_size;;
for( i=0;i<n; ++i ){ verties=id[n-i-1]; }
for( i=0;i<n_snodes; ++i ){ child=Sp[n_snodes-i-1]; }
for(Sp[0]=0, i=0; i<n_snodes; ++i ){ Sp[i+1]=Sp+child; }
for( k=0;k<n_snodes; ++k ){
char b=0;
for( i=Sp[k]; i<Sp[k+1]; ++i ){
if((verties+1)!=verties[i+1]){ b=1; break; }
flags[k]=b;
returnn_snodes;
尋找超結點的正向方法：
int find_supernodals_forward( char* mark, int* Sp, int*verties, char* flags, const __int64* Lp, const int* Li, int n )
int64 ii0,ii1, ii2, ii3;
int i, k,parent, n_snodes, snode_size;
const int imax=n-1;
for(Sp[0]=0, i=0, n_snodes=0, k=0; k<n; ++k )
if( mark[k]!=0) continue;
snode_size=1;
mark[k]=1;
verties[i++]=k;
ii0=Lp[k]; ii1=Lp[k+1];
parent=(k<imax)?Li[ii0+1]:n;
while(parent<n){
ii2=Lp[parent];ii3=Lp[parent+1];
if(((ii1-ii0-ii3+ii2)!=1)|(mark[parent]!=0)) break;
++snode_size;
mark[parent]=1;
verties[i++]=parent;
parent=(parent<imax)?Li[ii2+1]:n;
ii0=ii2;ii1=ii3;
Sp[n_snodes+1]=Sp[n_snodes]+snode_size;++n_snodes;
b=0;
find_supernodal_*函數的傳回值時超節點的數量；verties的大小是頂點的數量，每個超節點所包含的頂點編号在verties中是連續存儲的；Sp中包含了超節點的第一個頂點在數組中的索引位置以及超節點的大小資訊，第i個超節點中的第k個頂點為verties[Sp+k], 其中k<Sp[i+1]。可以如下周遊第i個超節點中的頂點：
p=Sp；q=Sp[i+1];
do{
v=verties[p];
…
}while((++p)<q);
flags中則标記了超節點的靜态和動态屬性，值為0表示靜态超節點，值為1表示動态超節點。
3.4 超結點的分解
對超結點的分解可以友善的利用高效的稠密矩陣算法，下圖是超結點分解的示意圖
由于大部分的計算都集中在syrk,是以我們重點對syrk的實作進行優化，第一章中的矩陣乘法所用到的技術可以得到複用，隻不過矩陣B變成了矩陣A的轉置。雖然可以通過直接使用*gemm來實作syrk操作，但是這種方法還是比較低效，因為每一步的計算量通常并不足夠大，可能無法充分利用硬體的計算資源；而且syrk比gemm具有更少的全局記憶體通路量。另一個需要自己開發這些操作而不使用cublas等通用庫的原因是它們的實作中syrk操作都是基于完整尺寸的矩陣，對于稀疏矩陣求解這種記憶體非常寶貴的計算來說過于浪費。是以我們實作自己的syrk操作通過對稱存儲且隻需要大約一半的存儲空間，其資料布局如下圖所示：
10.png (11.02 KB, 下載下傳次數: 0)
2016-6-24 11:30 上傳
3.5 基于多波前+超結點方法的并行分解
可以将多波前法和超節點方法結合使用以最大化并行度，每個超節點作為一個超波前(如無特别說明，我們将超節點波前統一簡稱為波前)，多個波前可以同時進行計算。在安排超節點的計算順序時，除了根據已經建立的消去樹外，還要考慮目前的計算資源，包括GPU的數量，每個GPU上得裝置記憶體的大小以及主機記憶體的大小。我們需要在求解性能和存儲空間之間取得一個最佳的平衡。是以每個裝置上可以允許同時進行計算的波前需要設定一個上限以避免在消去樹的某一層中的波前對存儲空間的需求超出已有的記憶體資源。同時當有多個GPU裝置時還需要綜合考慮如何将計算和所需要的裝置記憶體平均到每個裝置上。通常對消去樹的深度周遊可以最小化存儲資源使用量，而寬度周遊則可最大化并行度。另一個需要考慮的就是當矩陣規模很大時，記憶體和顯存往往隻能容下小部分資料，這時候就需要考慮核外求解，也就是将已經參與完計算的波前資料寫入磁盤。但是後續的計算可能會要求載入已經寫入磁盤的波前資料，那麼在确定波前的消去順序時又需要考慮如何最小化磁盤讀寫操作。确定各個波前消去順序的方法很多，但是波前消去的順序的多選擇性僅限于消去樹同一層的波前之間，各個層之間的波前的消去仍需按順序進行，是以可以很容易的根據自己的設想或需求設計相關算法，這裡不再叙述。下面給出基于多波前+超節點方法的分解算法：
idev= current_cpu_tid : 目前裝置的id
A: 待分解的矩陣
U[idev]: 第idev裝置上得波前更新矩陣，用來存儲每個波前的syrk的計算結果
for( level = 0; level<etree.n_levels; ++level )
start_snode =etree.snode_list[level][idev]
end_snode=etree.snode_list[level][idev+1]
get_front( front[idev],snode, A )
for( snode=start_snode+1;snode<end_snode; ++snode ){
if(next_snode.flag == dynamic ){
async{
assembly_front(front[idev].hFactor, next_snode, A )
copy(front[idev], next_snode )
} else {
front.[idev].hFactor=A[snode.Ap[next_snode.id]: snode.Ap[next_snode.id+1]];
async_copy(front, next_snode )
front[idev].dU= snode.factor( front[idev].dF )
locked_extend_add(A, front[idev], snode )
if((compute_mode==out_core)&&(snode.have_no_parents)&&(snode!=last)){
async_output_to_disk(snode )
上述算法是基于多GPU的多波前+超節點方法的實作，使用了靜态+動态超節點分類技術。根據消去樹的指導，逐層求解各個層内的波前(亦即超節點波前)。開始分解時整個待分解矩陣存儲在記憶體或磁盤中，然後在每一層計算的開始處使用異步傳輸将下一個要計算的波前資料從主機記憶體到裝置記憶體的複制過程與目前波前的計算重疊進行，進而可以有效隐藏CPU和GPU之間進行資料傳輸延遲，大部分時候這個延遲基本可以完全隐藏掉。然後通過判斷目前超節點的标志，如果是靜态超節點，隻需直接進行原地計算；如果是動态超節點，則要先組裝超節點再進行計算。計算的結果分為兩部分，一個是波前本身的分解因子F，另一個是目前波前對其它波前的更新矩陣U。由于可能有多個波前對同一個波前的計算有貢獻，是以我們需要将extend_add操作鎖定以保證計算不發生寫沖突。最後需要根據一定的條件決定是否将以完成計算的波前寫入磁盤，如果設定的計算模式是核外求解模式并且沒有後續的波前計算需要用到目前的波前，則将資料寫入磁盤。同時還要注意的是第三個判别條件，如果該波前是整個消去樹中的最後一個波前，那麼考慮到後續需要進行回代求解，是以從效率上考慮，就無需再寫回磁盤了。
本章小結
這一章我們介紹并給出了稀疏直接求解器的具體實作以及先關的優化技術，這其中的很多技術都可以直接或是稍加變通的用在很多數值算法中。特别是超節點方法特别适合用GPU進行計算，但是也需要通過各種優化盡量減少CPU和GPU之間的資料通信帶來的開銷，比如：
通過對超節點的靜态和動态分類減少資料移動的開銷，最大化原地分解。
通過動态并行對完整的三角型超節點進行分解以減少和主機通信的開銷。
通過使用記憶體映射和點對點傳輸降低CPU和GPU之間以及GPU和GPU之間的資料複制開銷。
通過對超節點波前的消去順序進行排程排序使得每個波前進出記憶體的次數最小化。
使用基于雙全局記憶體和雙主機記憶體的雙緩沖技術将待計算節點/波前的傳輸和目前節點/波前的計算重疊進行，進而充分利用CPU和GPU之間的異步并發模式。
1《Direct Methods for Sparse LinearSystem》，Timothy A. Davis
2 An Approximate Minimum Degree Ordering Algorithm，Patrick R. Amestoy Timothy A.Davisy Iain S. Du_z SIAM J, 1996
3 The Computational Complexity of theMinimum Degree Algorithm，P. Heggernesy S. C. Eisenstatz G. Kumfertx A. Pothen，2000
4 Theevolution of the minimum degree ordering algorithm, SIAM Review，1989
5 J.A. GEORGE AND J. W. H. LIU, A quotient graph model for symmetricfactorization,
in SparseMatrix Proceedings 1978, I. S. Duff and G. W. Stewart, eds., SIAM
Publications,1978
6 A NESTED DISSECTION APPROACH TO SPARSE MATRIXPARTITIONING FOR PARALLEL COMPUTATIONS, ERIK G.BOMAN_ AND MICHAEL M. WOLF
7 Minimum Degree Reordering Algorithms: A Tutorial Stephen Ingram University of British Columbia
原文釋出時間為：2016-6-24 11:36:14
原文由：nvadmin 釋出，版權歸屬于原作者
本文來自雲栖社群合作夥伴NVIDIA，了解相關資訊可以關注NVIDIA官方網站